Multiresolution Modeling of High-Latitude Ionospheric Electric Field Variability and Impact on Joule Heating Using SuperDARN Data.

The most dynamic electromagnetic coupling between the magnetosphere and ionosphere occurs in the polar upper atmosphere. It is critical to quantify the electromagnetic energy and momentum input associated with this coupling as its impacts on the ionosphere and thermosphere system are global and major, often leading to considerable disturbances in near-Earth space environments. The current general circulation models of the upper atmosphere exhibit systematic biases that can be attributed to an inadequate representation of the Joule heating rate resulting from unaccounted stochastic fluctuations of electric fields associated with the magnetosphere-ionosphere coupling. These biases exist regardless of geomagnetic activity levels. To overcome this limitation, a new multiresolution random field modeling approach is developed, and the efficacy of the approach is demonstrated using Super Dual Auroral Radar Network (SuperDARN) data carefully curated for the study during a largely quiet 4-hour period on February 29, 2012. Regional small-scale electrostatic fields sampled at different resolutions from a probabilistic distribution of electric field variability conditioned on actual SuperDARN LOS observations exhibit considerably more localized fine-scale features in comparison to global large-scale fields modeled using the SuperDARN Assimilative Mapping procedure. The overall hemispherically integrated Joule heating rate is increased by a factor of about 1.5 due to the effect of random regional small-scale electric fields, which is close to the lower end of arbitrarily adjusted Joule heating multiplicative factor of 1.5 and 2.5 typically used in upper atmosphere general circulation models. The study represents an important step toward a data-driven ensemble modeling of magnetosphere-ionosphere-atmosphere coupling processes.

Introduction

The most dynamic electromagnetic coupling between the magnetosphere and ionosphere occurs in the Earth's polar upper atmosphere. In particular, collisions between neutrals and ions drifting under the effect of the elevated high‐latitude ionospheric electric field are a major source of heating and momentum transfer, making a global impact on the upper atmosphere. The resulting energy and momentum deposition leads to the acceleration of neutral winds and Joule dissipation, triggering dramatic global upper atmosphere responses, for example, global temperature and neutral mass density enhancements, pole‐to‐equator general circulation, and atmospheric traveling disturbances (e.g., Burns et al., 2014; Fuller‐Rowell, 2014; Schunk, 2014). Practical effects include altered drag force on low‐Earth‐orbit satellites and debris by sudden changes in neutral mass density, aggravating our ability to track these objects to mitigate potential collisions; radio wave propagating disruption affected by ionospheric density changes, deteriorating reliability of communication, navigation and positioning systems; and geomagnetically induced currents in the ground resulting from intensified ionospheric currents, affecting power transmission systems, oil and gas pipelines, railway systems, and any other extended ground‐based conductor systems (e.g., Groves & Carrano, 2016; Marcos et al., 2010; Pulkkinen et al., 2017). Accurate knowledge of this energy and momentum source in the polar ionosphere is therefore of great scientific interest and has important economic and societal benefits.

The current general circulation models of the upper atmosphere exhibit systematic biases that can be attributed to the underestimation of the high‐latitude energy sources, likely resulting from an inadequate representation of the Joule heating rate (e.g., Codrescu et al., 1995, 2000; Deng et al., 2009; Matsuo & Richmond, 2008; Zhu et al., 2018). These biases are known to exist regardless of geomagnetic activity levels. The Joule heating rate is proportional to the square of the electric field magnitude and scales linearly with the Pedersen conductivity. Both of these ionospheric electrodynamics state variables are highly variable and heavily influenced by magnetosphere‐ionosphere coupling processes that are not usually self‐consistently solved in the upper atmosphere general circulation models and thus have to be empirically parameterized and/or specified as boundary conditions. The empirical models of high‐latitude ionospheric plasma convection designed to characterize the climatological behavior of the global large‐scale electric fields are not suited to representing highly variable localized multi‐scale electric fields and result in residual fields with a magnitude as large as the modeled global fields themselves (e.g., Codrescu et al., 2000; Cousins & Shepherd, 2012; Matsuo et al., 2002, 2003). Even with data assimilative procedures, the instantaneous states of the localized electric fields on scales smaller than 500 km and 5 min associated with highly transient and regional magnetosphere‐ionosphere coupling processes are difficult to capture (Matsuo & Richmond, 2008; Matsuo et al., 2005). As pointed out originally by Codrescu et al. (1995) and elaborated in later work (Codrescu et al., 2000; Deng et al., 2009; Matsuo & Richmond, 2008; Zhu et al., 2018), the underrepresented electric field variability in the upper atmosphere general circulation models is considered as one of the primary causes of the underestimation of Joule heating rate.

The volume integrated Joule heating rate is given as where is the Pedersen conductivity which specifies the conductivity associated with ionospheric electric currents that flow perpendicular to the geomagnetic field and parallel to the electric field defined in the reference frame moving at the velocity (Jackson, 1999). Note that is the electric field in the Earth frame of reference and essentially electrostatic on time scales longer than tens of seconds and in the bottomside ionosphere where neutral species predominate over plasma, is approximately equal to the neutral wind velocity (Kelly, 2009). thus represents the dynamo fields resulting from an electromotive force induced by the neutral wind traversing the geomagnetic field. Because of very high electrical conductivity in the direction of , a geomagnetic field line is effectively equipotential where it traverses the ionosphere and therefore is nearly constant with altitude along the direction of the field line which is nearly radial at high latitudes. Note that for the sake of simplicity the radial component of the electric field is ignored from discussion in this paper. When the effect of is small, can be approximated using the height‐integrated Pedersen conductivity as where is the polar angle (i.e., magnetic co‐latitude) and is the azimuth angle (i.e., magnetic local time [MLT]). For simplicity, geomagnetic fields are here assumed strictly radial. Let us suppose that can be decomposed into global large‐scale electric fields and regional small‐scale electric fields as and represents time‐dependent mean vector fields and represents stochastic or random vector fields that belong to a certain probabilistic distribution. It is important to note that specific instances of are different as are random fields but its statistical characteristics of their randomness are assumed to be temporally stationary, thus in Equation 2 a dependence on is dropped. (Note that this assumption is made due to the necessity to aggregate data over time in the current study and should be relaxed in the future as discussed later.) It is easy to see the underestimation of the Joule heating rate could result from not accounting for effects of , which can be as large as at times, in the upper atmosphere general circulation models.

Additional sources of uncertainty in determining the Joule heating rate include neutral winds and Pedersen conductivity . The contribution of dynamo fields to can be as large as 30% especially when neutral winds are driven by elevated ionospheric plasma convection during geomagnetic storms (e.g., Lu et al., 1995; Ridley et al., 2003; Sangalli et al., 2009). Depending on the direction of the neutral wind, the dynamo field effect can increase or decrease the total Joule heating rate. For instance, Lu et al. (1995) found the neutral winds have approximately a 28% negative effect on Joule heating rate on average for the 2‐day geomagnetically disturbed period investigated. It is also important to note that is not constant with altitude, thus requiring knowledge of vertical distributions of both the neutral wind and Pedersen conductivity when computing (Thayer, 1998). In addition, ionospheric conductivity varies considerably due to ionization of the neutral species by solar extreme ultraviolet radiations and auroral energetic particle precipitations. The effects of auroral ionization can be extremely localized and transient, which are difficult to characterize with the currently available auroral models (e.g., Newell et al., 2009). Furthermore, Dimant and Oppenheim (2011) have pointed out that during geomagnetically disturbed conditions the Pedersen conductivity can be enhanced considerably due to strong anomalous electron heating and nonlinear electric currents resulting from the Farley‐Buneman instability (Buneman, 1963; Farley, 1963). Part of the instability effect was incorporated into a recent upper atmospheric general circulation modeling study by Liu et al. (2016). In spite of recent progress in modeling, considerable uncertainty still remains in representing all physical processes responsible for Joule heating in current general circulation modeling.

Spatial and temporal coherence and other properties of randomness of the electric field variability affect the estimate of Joule heating rate as they control how effectively momentum and energy are transferred from ionospheric plasma to neutral species. Matsuo and Richmond (2008) demonstrated this effect using ensemble modeling and Gaussian random fields generated with the space‐time covariance model derived from DE‐2 observations in Matsuo et al. (2002, 2005). If spatiotemporal coherence is taken into account when incorporating the effects of electric field variability into an upper atmospheric general circulation model, electric field variability becomes more effective in influencing the neutral winds and thus affecting the overall Joule heating rate. The analysis of Super Dual Auroral Radar Network (SuperDARN) plasma drift measurements by Cousins and Shepherd (2012) has furthermore revealed scale‐dependent non‐Gaussian probabilistic behaviors of the electric field variability. The observed localized transient characteristics of electric field variability are difficult to model using currently available standard statistical inferential frameworks. In response to the need for a new framework, Fan et al. (2018) have developed a multiresolution non‐Gaussian random field model by using a class of specialized needlet basis functions (Marinucci & Peccati, 2011) that has all the desired properties (including smoothness, spatial and frequency localization, frame properties), which have enabled for flexible multiresolution reconstruction of scalar electrostatic potential fields on the sphere. By using the Lyon‐Fedder‐Mobarry (LFM) magnetosphere‐ionosphere coupled model simulation results (Wiltberger et al., 2016). Fan et al. (2018) have furthermore demonstrated a measurable impact on the Joule heating rate.

By building on the statistical inferential framework developed by Fan et al. (2018), the objective of this paper is to characterize the electric field variability as multiresolution non‐Gaussian random vector fields from actual SuperDARN observations, and to evaluate its impact on the Joule heating rate. The novel elements of the data analysis method and modeling technique described in Section 3 are as follows. The work of Fan et al. (2018) is extended to vector fields in this study. This is important as existing multiresolution bases for dealing with vector fields (such as vector spherical harmonics) do not have the spatial localization property and hence are not appropriate for describing features that are spatially localized. Since needlets can be represented in terms of spherical harmonics, in particular through Legendre polynomials, the surface gradient and curl operators can be applied to them to yield vectorial needlets that inherit the spatial compactness, facilitating flexible, multiresolution representations of the curl‐free multi‐scale electrostatic fields. Furthermore, the adaptive Markov Chain Monte Carlo (MCMC) estimation approach developed in Fan et al. (2018) is used to characterize non‐Gaussian random electric fields from SuperDARN observations. The needlet‐based model is constructed on the entire spherical domain due to the spherical quadrature grid points on which needlets are placed (Fan et al., 2018; Marinucci & Peccati, 2011). Figure 1 displays the spherical quadrature grid points on the unit sphere at different resolution levels, and an illustrative example of a needlet placed at one of the points. To assure coefficients for the needlets outside the high‐latitude region are well estimated with data confined to high latitudes, the SuperDARN data points are stretched from the high‐latitude region to the entire sphere in this work by following the work by Weimer (1995) and Ruohoniemi and Baker (1998). See Section 5.1 of Fan et al. (2018) for more discussion.

Figure 1

(Top) The quadrature grid points on the unit sphere, set according to the symmetric spherical ‐designs, at different resolution levels . (Bottom) Plots of a representative spherical needlet at respective resolution levels.

An additional notable element of the study is a special pre‐processing of a standard SuperDARN FITACF data designed for the needlet‐based approach to modeling electric field variability as described in Section 2. Although similar approaches have been used in the past (Ruohoniemi & Baker, 1998), this is the first consolidated attempt to extract randomness information from the FITACF data product, and can serve as a foundation for follow‐on future studies with more data and validation with independent data. The SuperDARN data from a 4‐hour period on February 29, 2012 are selected for the study with a number of considerations including data coverage and consistency in geophysical conditions, and processed as described in Section 2. As demonstrated in Section 4, and discussed and summarized in Section 5, the study is an important cross‐disciplinary research and development effort that enables a more comprehensive data‐driven approach to modeling of magnetosphere‐ionosphere‐atmosphere coupling processes.

SuperDARN Data

The SuperDARN is an international network consisting of more than 30 low‐power HF (3–30 MHz) coherent scatter radars at middle to polar latitudes in both hemispheres that look into Earth's upper atmosphere and ionosphere (Chisham et al., 2007; Nishitani et al., 2019). The radars measure the line‐of‐sight (LOS) component of the ‐region ionospheric plasma drift velocity when decameter‐scale electron density irregularities are present and oriented favorably to produce backscatter. The irregularity motion here is due to drift. Normally, the SuperDARN radars are scheduled for 1‐min or 2‐min azimuthal sweeps in the normal mode. The step in azimuth between adjacent beams is and the range resolution is 45 km. This study uses LOS plasma drift velocity () from SuperDARN radars operating in the normal scan mode from the northern hemisphere over the 4‐hour period from 00:00 to 04:00 Universal Time (UT) on February 29, 2012. The location of these radars and their field‐of‐views (FOVs) as well as the data coverage are shown in Figure 2. This is a largely quiet period during the rising phase of Solar Cycle 24 (F107 = 100.1) with a minor geomagnetic activity of the Kp index of 3, the minimum Dst index of about −30 nT, and the Auroral Electrojet (AE) index ranging from 100 to 420 nT peaking at 02:50 UT.

Figure 2

(a) The location of Super Dual Auroral Radar Network (SuperDARN) radars in the northern hemisphere used in this study and their field‐of‐views and (b) scatter plot showing the SuperDARN measurement coverage during the 4‐hour interval on February 29, 2012.

is related linearly to the electrostatic potential , where , as described in Section 2.1. Both global large‐scale mean electric fields and regional small‐scale random electric fields , expressed in terms of the electrostatic potential, are here estimated from SuperDARN LOS velocity data. Among multiple types of SuperDARN LOS data products made available by the SuperDARN consortium for different scientific applications, the FITACF data product is used for estimating after the pre‐processing described in Section 2.2. The GRID data product that itself is a derived product of FITACF data is used for estimating using the SuperDARN Assimilative Mapping (SAM) procedure (Cousins et al., 2013b) as described in Section 2.3. The SAM analysis is conducted every two minutes using the GRID data, and global large‐scale fields' contribution to the LOS plasma drift velocity is subtracted from FITACF LOS velocity data (see Section 3.2). Note that the FITACF data is aggregated over the four hours for needlet‐based analysis as described in Section 3.

SuperDARN Line‐of‐Sight (LOS) Plasma Drift Velocity

Assuming that the geomagnetic field is strictly radial (i.e., ), the electric field is expressed as where is the radius of the ionosphere. The plasma drift velocity is thus given as where is a magnitude of the geomagnetic field that varies over the sphere. The LOS component of the velocity , which is where is a unit vector that gives the direction of the line‐of‐sight, becomes

The SuperDARN data sets include values of and for each data points.

Pre‐Processing of SuperDARN FITACF Data

In order to estimate velocity and other parameters from autocorrelation functions calculated from the radar backscatter returns for each beam and range gate of a radar, the FITACF fitting routine is applied to estimate Doppler velocity, spectral width, and backscatter power (Ribeiro et al., 2013). Similar to the criteria used in Cousins and Shepherd (2012), several steps have been taken to ensure that only high‐quality LOS plasma velocity measurements from the F‐region ionosphere are included in the analysis. Primary selection criteria include (a) the slant range greater than 600 km, (b) the backscatter power or signal‐to‐noise ratio (SNR) greater than 8 dB, and (c) the velocity error less than 100 m/s. In addition, ground scatters are carefully excluded based on the ground scatter flag from the standard SuperDARN data processing and spectral width and velocity magnitude values. Incidental outlier data with very large velocity values can still be commonly found after the above processing, necessitating further processing on data. These outliers are excluded by applying a filtering scheme that examines the number of nearest‐neighbor points in time and space on a cell‐by‐cell basis and identifies the beam and gate cells that are isolated. The filtering is performed using a 10‐min time window and the beam‐gate cells adjacent to the cell under consideration. In general, a data point that persists for at least two minutes and has at least one spatial neighbor would be considered valid. This analysis has also produced values of median velocity and standard deviation of velocity within each valid beam‐gate cell.

Estimation of Global Large‐Scale Electric Fields by SuperDARN Assimilative Mapping (SAM)

The distribution of global large‐scale ionospheric convective electric fields is determined at 2‐min cadence from the SuperDARN GRID data over 4‐min windows using the SAM procedure (Cousins et al., 2013b). The SAM uses a set of the spherical cap harmonics functions developed by Richmond and Kamide (1988), with the spherical harmonics of the order 12 and non‐integer degrees of 72.6 for the order zonally symmetric harmonic functions that give the effective resolution of longitude and latitude in terms of the Nyquist sampling rate. The SAM solves a Bayesian spatial statistical prediction problem for ionospheric convective electric fields just as the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) (Richmond & Kamide, 1988), and computes the posterior mean given the prior mean convective electric fields specified by Cousins and Shepherd (2010). A major advantage of the SAM over the AMIE is the use of prior model error covariance developed from a large volume of SuperDARN data in Cousins et al. (2013a) for the prior model of Cousins and Shepherd (2010). The LOS plasma drift velocity due to these global large‐scale electric fields is computed according to Equation 4 and subtracted from the pre‐processed LOS velocity data explained in Section 2.2.

Due to the use of global spherical cap harmonics functions, with a limited resolution, in the SAM, it is sufficient to use the GRID data which provide the standardized LOS velocity values on an equal‐area grid over a fixed period of time of 1 or 2 min, rather than the FITACF data which contain the LOS velocity measurements recorded by individual radars as a function of beam‐azimuth range gate setting. GRID data is a highly processed data product derived from FITACF data. A median filter is first applied to the individual radar scan data to remove noise to calculate the median LOS velocities of a particular scan. The LOS vectors are then mapped within the cells of an equal‐area grid, which is defined in the geomagnetic coordinates system with each cell measuring in latitude, to eliminate biases that would derive from the much denser sampling over nearer radar range gates. The vectors contributed by a radar to a particular cell are averaged over a fixed period of time to obtain the GRID LOS data product. More details of GRID data processing can be found in Section 3 of Ruohoniemi and Baker (1998).

Needlet‐Based Approach to Modeling Electric Field Variability

The novel element of the statistical modeling approach presented here is the use of a multiresolution tight frame called needlets (Marinucci & Peccati, 2011) to represent stochastic fluctuations in the electric field vectors. In the same way wavelets facilitate analysis of transient and localized signals, needlets enable us to represent spatially localized features of the observed electric field variability in functions defined over a spherical domain. Needlets have been shown to be more efficient than spherical harmonics in representing spatially localized features on the sphere as linear combinations of spherical harmonics, through a construction involving Legendre polynomials (Scott, 2011). Furthermore, the surface gradient operators can be applied on them, thus facilitating multiresolution representations of the curl‐free multi‐scale electrostatic potential fields. In Section 3.1, spherical needlets used in this study are briefly defined, and Fan et al. (2018) should be referenced for more details.

Multiresolution Tight Frame: Spherical Needlets

Specifically, a needlet function at scale and location , , evaluated at a point on the unit sphere takes the following form: where the nonnegative function is bandlimited and enables a frequency tiling, controls the window size of the frequency tiling, are quadrature (location, weight) pairs for scale and location , ’s are the standard orthonormal, complex‐valued spherical harmonics basis functions corresponding to frequency (degree) index and phase (order) index , and is the associate Legendre polynomial of degree . The function is positive on the interval , satisfies the resolution of identity condition for , and is times continuously differentiable for some . From Equation 5, it is evident that needlets ’s are bandlimited over spherical frequencies ranging from integer index greater than or equal to to less than or equal to . In addition to being localized in the frequency domain, the needlets are localized in the spatial domain. Specifically, where is a constant independent of and . See Figures 1 and 2 of Fan (2015) for examples of a spherical needlet , and Fan (2017) for more detail. Hereafter, we choose following the prior work (e.g., Fan et al., 2018). Note that because of the linear representation of needlets in the spherical harmonic basis, needlet coefficients of a scalar function can be obtained from the spherical harmonics coefficients of the function through a linear transformation, since for any (quadratically integrable) function on ordinary sphere , where and denote the needlet and spherical harmonics coefficients, respectively.

Needlet‐Based Random Electric Fields Model

Suppose that there is a total of SuperDARN LOS plasma velocity data points at locations , , in the high‐latitude region of the northern hemisphere. Note the data points shown in Figure 2b are down‐sampled as explained later in Section 3.4 before being used for the needlet model estimation. Since these data contain the observational noise, they are modeled by the following statistical model where represents the observational data, represents the underlying true velocity value, and is the observation noise or error with standard deviation . (Note that represents the normal distribution with a mean parameter and a variance parameter .) For simplicity, is henceforth assumed independent of the location (i.e., , ). As described later in Section 3.3, is one of the statistical model parameters to be estimated from SuperDARN LOS velocity data.

According to Equation 4, the velocity field can be derived from the electrostatic potential by applying the differential operators, and the electrostatic potential can be decomposed into two components: global large‐scale and spatially localized regional small‐scale components, and , which respectively correspond to and . Therefore, can also be decomposed into two components accordingly, that is,

The SAM procedure described in Section 2.3 is well suited to estimate from the SuperDARN LOS plasma velocity data . We subtract the fitted , denoted by , from , and obtain is precisely what is modeled in terms of by spherical needlets ’s.

Since SuperDARN LOS plasma velocity data points are restricted to the high‐latitude region, the data points are stretched to the entire sphere by mapping observation location points to . Since the magnitude of has a strong dependency on the latitude, a variance profile, that is, the variance of the observed LOS velocity field as a function of the latitudinal location, is introduced as a function of co‐latitude, and is modeled by the product of the variance profile and a linear combination of spherical needlets as follows: where is the variance profile function, and are needlet coefficients, which are random variables. As in Fan et al. (2018), it is assumed that ’s are distributed as scale multiples of a ‐distribution, that is, , where denotes the ‐distribution with degrees of freedom. The ‐distribution has heavier tails in comparison to the normal distribution. is used for this study following Fan et al. (2018) wherein was chosen among 2.5, 3, and 4 in their applications to the LFM model output as it yielded the best predictive performance for simulated data. Note that with infinite degrees of freedom, the ‐distribution approaches to the standard normal distribution. The assumed distribution characterizes both scale‐dependent variations and spatially localized features of the electric field variability. Moreover, ’s are assumed to be statistically independent for simplicity. Due to the non‐Gaussianity of and the spatial localization of , the resulting field is also non‐Gaussian.

The variance profile function is assumed to have the representation , where are the basis functions specified as cubic B‐splines due to their numerical stability. To avoid the non‐identifiability issue, the first B‐spline is dropped in the formula. As described later in Section 3.3, the B‐spline weights that control the variance profile and the ‐distribution population parameters , , that determines a probabilistic distribution of needlet coefficients, given a value of , are estimated from SuperDARN LOS data.

By the chain rule, we know that and where . Plugging these into Equation 4, we have

In this way, the spherical needlets are transformed into two new sets of basis functions in the domain of line‐of‐sight velocities. Recall that where is the centroid of the needlet . Then and where . Note that can be efficiently computed by using a recursive formula.

The statistical model for can be summarized as the following matrix‐vector form: where , , , , and are the design matrices constructed by the new basis functions , respectively, and . For convenience, we hereafter use to denote and to stand for so that

Adaptive Markov Chain Monte Carlo Estimation

This section describes how the parameters and that determine a probabilistic distribution of needlet coefficients and the observational noise parameter are estimated from the residual SuperDARN LOS velocity data using the adaptive MCMC method. These parameters, grouped here as the vector , are assumed a priori independent. The prior distributions of and are the non‐informative Jeffreys’ priors. The prior distribution of is assumed to be , where the hyperparameter is chosen to be sufficiently large such that the prior distribution is nearly non‐informative. Under these settings, the posterior distribution of the parameters can be computed by the following MCMC algorithm.

Since a ‐distribution can be expressed as scale mixture of Gaussians, the probability distribution of can be written in a hierarchical form by introducing an auxiliary random variable where represents the inverse gamma distribution with a shape parameter and a scale parameter . Denote by the vector consisting of the coordinates , and to be the vector comprising of . We shall employ a Gibbs sampler to obtain samples from so that the full conditional distributions of and have closed forms. The full conditional distribution of (i.e., ) is multivariate Gaussian, and hence sampling from it requires operations, where is the total number of needlets. This is computationally intractable for large . Nonetheless, numerical experiments indicate that the subblocks are weakly correlated, where , where is the number of needlets at level . Therefore, the sampling step for is achieved by successive draws from the conditional subblocks . The full conditional distribution of is not available in closed form. Therefore, we sample from using an adaptive Metropolis algorithm (Andrieu & Thoms, 2008, Algorithm 4) and incorporate it into the Gibbs sampler.

Suppose denotes the subblock of corresponding to the ‐th level of needlets, so that , and . Then the aforementioned adaptive Metropolis‐within‐Gibbs sampler can be summarized as follows:

Sample from , where and

Sample from , where are independent and distributed as

Sample from , where are independent and distributed as where represents the gamma distribution with a shape parameter and a rate parameter .

Sample from

Sample using the adaptive Metropolis algorithm from

The proposal distribution is chosen to be where is a parameter adaptively tuned with the goal of achieving the optimal acceptance rate (Gelman et al., 1996), and is adaptively updated to approximate the covariance matrix of the full conditional distribution of .

Down‐Sampling of SuperDARN LOS Data

The SuperDARN LOS plasma velocity measurements are down‐sampled before being applied to estimation of the needlet model parameters. As shown in Figure 2b, the LOS velocity measurements are unevenly distributed over the high‐latitude region of the sphere. The aim of down‐sampling is to assure that the data are evenly distributed and the computational cost is manageable. To achieve the goal, we first partition the sphere into approximately equal‐area regions by applying Voronoi tessellation on the sphere. The number of partitioned regions is the same as the number of data points. We then calculate the area of each region. The sampling probability of each data point is proportional to the area of its corresponding region. Roughly speaking, the larger the surface area of one region, the further the data point in the region is from the neighboring data points. Therefore, a higher probability of retaining the data point is assigned. In this way, the data points after down‐sampling are approximately uniformly distributed.

Model Performance

A summary of the performance of MCMC based estimation of the needlet‐based model parameters described in Sections 3.2 and 3.3 is given here. We set to 2 since the global large‐scale components have already been subtracted from the SuperDARN LOS plasma velocity data. is set to 3 given the the computational limitation of the model with too high as well as the SuperDARN data signal‐to‐noise ratio, thus needlets at two resolution levels are used. (More discussion on the needlets resolution is given in Section 4.1 when the geophysical interpretation of SuperDARN data analysis results are presented.) The logarithm of the variance profile function is represented by a linear combination of cubic B‐splines with one interior knot . The parameter estimates are calculated as the average of 1,000 MCMC samples.

Figure 3 shows the estimated variance profile as a function of co‐latitude. The peak is around latitude, consistent with the locations of the high LOS residual velocity within auroral oval zone. As described in Section 3.4, the needlet model is fitted to a subset of the SuperDARN LOS plasma velocity data after down‐sampling. We examine the model out‐of‐sample prediction performance on the remaining data, which is shown as a scattered plot of LOS velocity magnitudes in Figure 4. The predicted values generally align with the observed values as the Pearson correlation coefficient between them is approximately 0.33. In terms of the magnitude, the predicted values are mostly smaller than the observed values. This can be explained by the following reasons: (a) The observed LOS velocity residuals show large variability and may contain some outliers. As shown in the middle plot of Figure 5, there are clearly some extreme values (indicated by sporadically appearing red arrows), not entirely captured by the model; (b) The needlet‐based model cannot represent features with resolution higher than due to the computational limitation; (c) We have assumed a simplified structure for the underlying electrostatic field, which is longitudinally (magnetic local time) symmetric with a variance profile depending on latitudes only; and (iv) Beyond these, in general, the predicted values of observations under a Bayesian paradigm (or in a random effects model) tend to shrink toward zero, even when the model represents the data perfectly. However, the observed LOS velocity residuals appear heterogeneous with respect to both longitudes and latitudes. These are the unique challenges of modeling the LOS velocity residuals, which we will discuss further in Section 5.

Figure 3

The fitted variance profile model defined in Equation 6 to residual Super Dual Auroral Radar Network line‐of‐sight plasma drift velocity data is shown as the latitudinal distribution of velocity standard deviation in m/s.

Figure 4

Out‐of‐sample prediction of residual Super Dual Auroral Radar Network line‐of‐sight plasma drift velocities versus observed values in m/s. Pearson correlation coefficient between predicted and observed values is 0.33.

Figure 5

Line‐of‐sight (LOS) plasma velocities from one 2‐min scan of Super Dual Auroral Radar Network radars from 00:04 to 00:06 UT on February 29, 2012: (Left) The FITACF LOS plasma velocity (), (middle) residual velocity (), and (right) regional small‐scale velocity modeled by the needlet model ().

Results

This section summarizes the results of the needlet‐based approach to modeling electric field variability using SuperDARN data described in Sections 2 and 3, and demonstrates how the approach can help better represent the impact of the high‐latitude ionospheric electric field variability on Joule heating rate in the upper atmosphere general circulation models. It also illustrates how uncertainty in data‐driven modeling of electromagnetic coupling between the magnetosphere and ionosphere may be represented using ensembles.

Multiresolution Non‐Gaussian Electric Fields

Figure 6 shows the electrostatic potential fields at different needlet resolution levels generated using the estimated needlet model parameters (e.g., , ) from SuperDARN LOS observations. Note that the electrostatic potential fields shown here correspond to regional small‐scale electric fields (see Equation 6). The top row displays the mean prediction conditional on the SuperDARN observations, and the middle and bottom rows show two instances of random samples conditional on the observations. These mean and two random instances are shown to illustrate that the electric field variability is in fact modeled as random fields that belong to a certain probability distribution that is, conditional on the SuperDARN observations in contrast to the past studies wherein the sample mean and standard deviations of observations have been often used. These two instances are part of a 1,000‐member ensemble set generated from 1,000 independent random draws, which are being used for the Joule heating estimation shown in the next subsection. The potential fields at two needlet resolution levels at are shown in the first two columns, and the total potential fields, which is a combination of all resolution levels, is shown in the right‐most column.

Figure 6

Electrostatic potential fields in kV at different resolutions generated from the needlet model. (Top row) The mean prediction conditional on the observations; (Middle Row) a random sample conditional on the observations and (bottom row) another random sample conditional on the observations. (Left) Field at level ; (middle) field at level ; and (right) total field summed at levels .

The smallest scales resolved at these two needlet resolution levels correspond to the spherical harmonics frequency (degree) and phase (order) of , and , , respectively, which are equivalent to the spatial scales of in latitude and in longitude for and in latitude and in longitude for , with consideration of the factor 4 latitude coordinate stretching of the analysis domain as described in Section 3.2. (Note that an approximate spatial resolution corresponding to a certain degree and order of the spherical harmonic function is obtained using the Nyquist frequency of a half wavelength.) Even though the resolutions of needlets and spherical harmonics are comparable, due to needlets’ spatial and frequency localization and frame properties that needlet‐based model can better represent localized regional features that exist in the SuperDARN observations, in comparison to spherical harmonics with a global support that are designed to capture global structures.

The SAM used to model global large‐scale electric fields from the SuperDARN observations (see Section 2.3) can resolve scales up to the spherical harmonics degree and order of about and , corresponding to the resolution of in latitude and in longitude. As shown in Figure 7 for 03:00 UT, the SAM field in fact exhibits global large‐scale features in comparison to regional small‐scale features that are present in muti‐resolution random electrostatic fields estimated from the needlet‐based modeling approach (Figure 6). Due to the variance profile, shown in Figure 3, which peaks around latitude, more distinct features appear between and in regional small‐scale fields.

Figure 7

Global large‐scale Electrostatic potential field in kV, estimated from the Super Dual Auroral Radar Network GRID data over 4‐min windows using the SuperDARN Assimilative Mapping procedure, at 3:00 UT on February 29, 2012, around the peak of Auroral Electrojet index.

Impact on Joule Heating Rate

Figure 8 shows the ensemble mean of the hemispherically integrated Joule heating rate computed with the effect of random regional small‐scale electric fields (blue) as well as the hemispherically integrated Joule heating without the effect of (black). The hemispherically integrated Joule heating rate with the effect of is integrated over the northern hemisphere high‐latitude ionosphere from to in latitude, following the definition given in Equation 1, as where is an ensemble member index, , the global large‐scale electric field is specified by the SAM as described in Section 2.3. Note that a 1,000‐member ensemble set of random regional small‐scale electric field is computed from as described in Section 4.1. The hemispherically integrated Joule heating without the effect of is given as

Figure 8

Hemispherically integrated Joule heating rate in GW from 00:00 to 04:00 UT on February 29, 2012. The ensemble mean of the Joule heating rate computed with the effect of random regional small‐scale electric fields (Equation 7) is shown in blue solid line, along with the upper and lower bounds given in terms of two standard deviations shown in blue dash lines. The Joule heating rate resulting only from global large‐scale electric field without (Equation 8) is shown in black line. As a reference, the high‐latitude geomagnetic activity, Auroral Electrojet index nT, is overlaid in red.

In Equations 7 and 8, the height‐integrated ionospheric conductivity is specified using empirical models of the solar EUV conductance and auroral conductance. The solar EUV conductance model is parameterized by solar zenith angle and the solar F10.7 index (e.g., Moen & Brekke, 1993), and the auroral conductance is based on the Ovation Prime empirical aurora model (Newell et al., 2009) and the empirical relationship of Robinson et al. (1987). Note that the Ovation Prime model is parameterized with respect to the upstream solar wind and interplanetary magnetic field conditions. Except for a minor geomagnetic activity, there is no notable geomagnetic activity during the time period of 00:00 to 04:00 UT on February 29, 2012. The overall Joule heating rate is thus small. In general, the Joule heating rate tracks temporal changes of the AE index (red) shown also in Figure 8, which is due to the changes of large‐scale electric fields. By taking the regional small‐scale electric field variability into account, the Joule heating rate increases by a factor of about 1.5 which is close to the lower end of an arbitrarily adjusted factor of 1.5 and 2.5 typically used in general circulation models. As discussed in Section 1, the biases in the upper atmosphere general circulation models attributed to an inadequate representation of the Joule heating rate exist regardless of geomagnetic activity levels. The future study needs to allow for more flexibility in the needlet model to account for temporal variation of the variance profile so that the changes of electric field variability at different geomagnetic activity levels can be better characterized.

Discussion and Conclusions

In response to the need for a new statistical inferential framework for data‐driven modeling of high‐latitude ionospheric electric field variability, Fan et al. (2018)’s spherical needlet‐based scalar random fields model is being extended for vector random fields and applied to the carefully curated SuperDARN FITACF LOS plasma velocity data set. The modeling results for the largely quiet period from 00:00 to 04:00 UT on February 29, 2012 show that the approach have the potential to rectify the underestimation of the Joule heating rate in the current upper atmosphere general circulation models due to insufficient representation of the electric field variability. The study demonstrates how data‐driven modeling of the magnetosphere‐ionosphere‐thermosphere coupling can be formulated in an ensemble modeling framework. Specific findings of the current efforts are summarized as follows.

The needlet‐based probabilistic approach to modeling regional small‐scale electric field variability can help estimate a distribution of electric field variability conditioned on actual SuperDARN LOS observations. The probabilistic modeling approach, rather than a deterministic model approach, to represent regional small‐scale electrostatic potential fields allows us to borrow information across spatial locations to model random fields globally even outside the domain of SuperDARN observation locations. The variance of random fields varies with latitudes according to the variance profile estimated from actual SuperDARN observations as shown in Figure 3. As shown in Figure 6, estimated regional small‐scale electrostatic potential fields at different resolutions exhibit considerably more localized fine‐scale features in comparison to global large‐scale potential fields modeled using the SAM procedure (Figure 7). This is enabled by the spherical needlet frames' spatial localization and overcompleteness properties and reinforces the fact that spherical harmonic basis functions, with global support, are not suited for describing features that are spatially localized.

As shown in Figure 8, the overall hemispherically integrated Joule heating rate during a largely quiet period on February 29, 2012 is increased by a factor of about 1.5 due to the effect of random regional small‐scale electric fields . This is close to the lower end of arbitrarily adjusted Joule heating multiplicative factor of 1.5 and 2.5 typically used in upper atmosphere general circulation models. It is expected that this factor will change under different geophysical conditions. The magnitude of regional small‐scale electric fields would likely become greater under southward interplanetary magnetic fields, winter conditions, and geomagnetically disturbed conditions (e.g., Matsuo et al., 2003), thus increase the Joule heating rate accordingly. As shown in Matsuo and Richmond (2008), the electric fields with the scale size smaller than vary strongly with seasons but less with interplanetary magnetic fields, suggesting a scale size dependence of regional small‐scale electric fields on the Joule heating.

The impact of the modeled electric field variability on the Joule heating rate is computed using a 1,000‐member ensemble set of regional small‐scale electric fields . This example demonstrates that how the uncertainty of the SuperDARN LOS data can be propagated to the estimate of Joule heating rate in general circulation models through the needlet‐based modeling of the ionospheric electric variability. The approach can also be applied to the output from high‐fidelity high‐resolution numerical simulations that may be computationally prohibitive to perform routinely. This study is an important step toward a data‐driven ensemble modeling of magnetosphere‐ionosphere‐atmosphere coupling processes.

Some of the methodological shortcomings identified by the current study can be addressed in future work. The needlet model can be expanded to account for non‐stationarity of the electric field variability not only with respect to magnetic latitudes but also MLT. By doing so, the electric field variability associated with specific physical processes such as convection reversal, and auroral electrojet that appear in localized locations can be better represented. As suggested by the out‐of‐sample prediction of LOS plasma drift velocities shown in Figure 4, SuperDARN LOS residual velocities with greater magnitudes ( 350 m/s) are not well predicted by the needlet model. This is evident in Figure 5. This can be addressed by increasing the needlet resolution level from to , equivalent to the spatial scales of 1.4° in latitude and 5.6° in longitude, ideally to , corresponding to the scales of 0.7° in latitude and 2.8° in longitude. Uncertainty resulting from inconsistent model assumptions associated with spatiotemporal stationarity of random fields should be better quantified using more data.

These methodological improvements will have to be accompanied with an improved uncertainty quantification in the determination of SuperDARN LOS velocity from radar backscatter. The availability of SuperDARN data with greater spatial coverage will alleviate the need to aggregate data over time allowing us to drop the assumption of spatiotemporal stationarity of random fields in the method. In order to increase the needlet resolution level to , SuperDARN data at a higher spatial resolution will be needed. Prospects for the availability of such SuperDARN LOS velocity data sets are discussed next. With more SuperDARN radars being constructed (e.g., Adak Island East and West radars in 2012, Hokkaido West radar in 2014, and Jiamusi radar in 2019 (Nishitani et al., 2019)), we could have a better spatial coverage for future work. However, there still exist a few challenges on obtaining more SuperDARN data. First, lack of ionospheric backscatter in SuperDARN data during the day would cause a data gap in MLT, particularly at mid‐latitudes (Figure 1b). Second, strong particle precipitation during geomagnetically active times could cause radar signals absorbed by the ionosphere. The Local Divergence‐Free Fitting technique from Bristow et al. (2016) is able to provide SuperDARN plasma velocity with a spatial resolution ( km) that is, comparable to the LOS velocity measurements. The LDFF technique uses all LOS velocities within a user‐defined region to produce local plasma convection which can resolve finer scale structures such as plasma flows associated with auroral arcs. This technique can be used in the future to obtain time‐dependent (mean) vector fields at much finer scales than the SAM, which is expected to improve the signal‐to‐noise ratios of residual SuperDARN LOS velocity data for the method presented in this study.