WU ET AL.
10.1029/2022SW003146
2 of 19 field-aligned currents (FACs) in response to the different phases of day and night-side reconnection (Nishimura et al., 2021). These ionospheric processes further affect the neutral atmosphere through the exchange and transport of momentum, energy, and composition in the coupled magnetosphere-ionosphere-thermosphere (MIT) system (Thayer & Semeter, 2004). Eventually, the change of the I-T system will lead to the change of conductance, current systems, and ion outflows, which in turn poses a feedback effect to the magnetosphere (Merkin & Lyon, 2010;Merkin et al., 2003;Tanaka, 2007;Yau & André, 1997).
Due to the dynamic and turbulent nature of magnetospheric processes, the high-latitude forcing of the I-T system such as aurora, electric fields, and FACs is highly multiscale and mesoscale structures that play an important role in the MIT coupling (Nishimura et al., 2021). For example, empirical models usually give the large-scale (>1,000 km) morphology of the auroral oval, which can differ greatly from the ground-based all-sky imager (ASI) observations such as those from the Time History of Events and Macroscale Interactions during Substorms (THEMIS) ASIs (Donovan et al., 2006). THEMIS ASIs depict rich mesoscale (10-100s km) structures, while a narrow field-of-view imaging can even resolve small-scale (<10 km) structures. During an expansion phase of a substorm, the auroral structure with scales smaller than 500 km contributes to 50% of the total energy flux, and mesoscale auroral processes such as poleward moving auroral forms, polar cap patches, auroral arcs, and streamers can feedback to the large-scale dynamics and impose net effects on the global distribution of electron densities (Gabrielse et al., 2021). Similar to aurora, electric fields also show multiscale features. Using Super Dual Auroral Radar Network (SuperDARN) measurements, Cousins and Shepherd (2012) found a large ratio (75%) of mesoscale to large-scale electric fields in terms of magnitude under a southward interplanetary magnetic field (IMF) condition. The scale analyses by Cousins et al. (2015) and Shi et al. (2020) show that mesoscale FACs contribute to nearly 60% of the spatial variability of FACs. These magnetosphere-originated processes are highly correlated and mesoscale auroral structures such as auroral arcs are often associated with enhanced FACs and electric fields (Nishimura et al., 2021), which have profound effects on the I-T system.
For I-T models, aurora and electric fields are the two most important drivers at high latitudes, thus it is critical to capture these two drivers realistically. In this study, we focus on the assimilation of auroral particle precipitation, specifically, energy flux and mean energy. Even though the empirical auroral models derived from historical data can capture large scales reasonably well (Newell et al., 2009;Wu et al., 2021;Zhang & Paxton, 2008;Zhu et al., 2021), they still miss the important mesoscale features. Wu et al. (2020) showed that only when the empirical auroral model is replaced by auroral observations from Special Sensor Ultraviolet Spectrographic Imagers (SSUSI) (Paxton & Meng, 1999;Paxton et al., 2002) onboard the Defense Meteorological Satellite Program (DMSP) satellites to drive the I-T model, the Thermospheric Temperature Enhancement and Inversion Layer (TTEIL) observed by the Fe-Boltzmann lidar at McMurdo, Antarctica, can be reproduced, and neutral densities in the F region match the Gravity Recovery and Climate Experiment (GRACE) observations. Similarly, Sheng et al. (2020) implemented THEMIS ASI auroral observations into the Global Ionosphere Thermosphere Model (GITM) and compared them with the simulations driven by the empirical model. The authors found that the magnitude of TIDs in GITM is almost doubled when driven by realistic THEMIS ASI observations and more consistent with observations. These previous studies indicate the necessity of developing data-driven auroral maps for the high-latitude drivers, especially when we focus on specific storms. Such efforts have been rarely made in the past, and the current work aims to address this challenge.
The existing techniques for auroral measurements include satellite and ground-based imagers, which have distinct spatial coverage and temporal samplings. SSUSI/DMSP measures global auroral emissions with a high spatial resolution and a revisit time of ∼30 min (three satellites) to the same magnetic latitude (MLAT) and magnetic local time (MLT). Ground-based instruments such as THEMIS ASIs provide both high temporal (3 s) and spatial resolution observations in North America. Empirical auroral models are built upon the statistics of a large number of historical observations and provide global auroral maps with highly smoothed patterns (Hardy et al., 1985;Roble & Ridley, 1987). They usually deviate from real-time observations, especially for mesoscale features. These deviations can often lead to systematic biases for the estimation of the general auroral activity level.
Even so, the empirical model can still provide the sensible information for large-scale features such as auroral boundaries. These data sources provide complementary information on auroral activities but are rarely used synergistically. One way to combine all data sources is by simply padding different types of auroral observations, but this method usually leads to discontinuous boundaries among different data sources and introduces unphysical gradients, which could lead to artificial perturbations in I-T models. Another approach to synthesize various data sources is the Assimilative Mapping of Ionospheric Electrodynamics (AMIE, Lu, 2017;Richmond, 1992;Richmond & Kamide, 1988), but its resolution is limited by the order of spherical cap harmonics (Matsuo, 2020).
In this paper, we apply a novel multiresolution spatial Gaussian process model (Lattice Kriging, Nychka et al., 2015) to incorporate auroral observations from satellite and ground-based data, as well as an empirical model where observations are unavailable. It uses a range limited basis function that better serves localized auroral assimilation. The mesoscale features in the satellite and ground-based observations are mostly kept in the assimilation results. In addition, the multiresolution modeling capability is fulfilled by locating multiple layers of basis functions with different resolutions. This method, therefore, provides a useful tool to study the multiscale processes and the corresponding impacts. Lattice Kriging has already been used in the lower atmospheric studies like surface temperature analysis (Heaton et al., 2019;Wiens et al., 2020). Wu and Lu (2022) have extended this model to vector fields and assimilated high-latitude electric fields using SuperDARN and Poker Flat Incoherent Scatter Radar (PFISR) data, which demonstrates its effectiveness in space weather studies. It is the first time that this model is applied to auroral assimilation.
The manuscript is organized as follows. Section 2 introduces the data sources. Section 3 describes the Lattice Kriging model including the principles and mathematical formula. Section 4 provides the detailed procedures to apply this model for auroral assimilation. Section 5 presents TIEGCM simulations driven by the empirical auroral model and the two different scales of auroral assimilation maps. Section 6 gives the conclusions and discussion.
The data sources used for the auroral assimilation include SSUSI onboard three DMSP satellites (F16, F17, F18), THEMIS ASIs, and a Kp-based empirical auroral model (Zhang & Paxton, 2008). The choice of empirical models is relatively flexible as long as the model provides full MLAT and MLT coverage. SSUSI is a remote-sensing instrument that measures ultraviolet emissions in five different wavelength bands from the Earth's upper atmosphere. The spatial resolution of SSUSI data product is ∼0.15°, which is sufficient to analyze the mesoscale structures of aurora in this study. The derived data products include the precipitating electron mean energy and energy flux. The three satellites sweep through the polar cap alternatively every 30 min, sampling through the auroral region (a swath) across the pole.
THEMIS ASIs observe the white light aurora over the North American continent from Canada to Alaska at a sampling rate of every 3 s, which provides high-resolution information about the rapid evolution of the aurora. The white light data are converted to red-green-blue colors by comparing with the nearest northern solar-terrestrial array (NORSTAR) meridian scanning photometers and multispectral ASIs, and then the color ratios are converted to energy fluxes and mean energies using the Strickland et al. (1983) formula (Mende et al., 2008). In this study, the electron mean energy and energy flux maps of spatial resolution 0.1° are used and the data are temporally down sampled to a 1 min basis.
The Zhang and Paxton (2008) Hardy et al. (1987), this model provides a more physical specification of the geo-effective energy flux and mean energy. Such information is also useful for the assessment of the statistical mean needed in the auroral assimilation (Equation 1 in Section 3.1) and for the regions where observations are not available. The empirical model can be generated at an arbitrary resolution, that is, on the satellite grids in the present study.
Owing to the noticeable auroral activity and decent data coverage on 20 February 2014, we use the auroral observations on this day as an example to demonstrate the methodology. The geomagnetic indices are shown in Figure 1. After 03:00 Universal Time (UT), a negative turning of IMF B z marks the start of geomagnetic disturbances. The Kp index reaches six and the symmetric disturbances for the magnetic H component (SYM-H) index reaches -100 nT, indicative of a moderate-intense storm. During this period, there are significant variations in the auroral electrojet (AE) indices, which reaches 1,200 nT, suggesting considerable auroral activities.
Figure 2 displays the auroral energy fluxes from the three data sources in the northern hemisphere at 11:50 UT, plotted in MLAT and MLT coordinates. This UT is chosen due to the clear auroral structures both in SSUSI and THEMIS observations. The instantaneous SSUSI observations are limited: to assimilate the auroral maps (mean energy and energy flux) for a particular time, SSUSI data falling into a 20 min time window (10 min before and 10 min after) are gathered. For example, SSUSI data from 11:40 to 12:00 UT are binned for the auroral assimilation at 11:50 UT (Figure 2a, more details in Section 4.1). SSUSI observations after the binning mainly cover the dawn and dusk sectors and show scattered auroral arc features. THEMIS ASIs provide night-time observations with several auroral enhancements spreading between 60° and 70° MLAT around midnight. The empirical model has a locally much smaller magnitude and smoother structure than the real observations but provides reasonable large-scale patterns and auroral boundaries for Kp = 6 geomagnetic condition.
In this section, we introduce the principles of the Gaussian process model adopted for Lattice Kriging (Section 3.1), and the implementation of the multiresolution data assimilation (Section 3.2).
We consider a spatial field y, whose values at all locations in a spatial domain X are assumed to follow a Gaussian process, and hence the value at every finite location {x i ,1 ≤ i ≤ n} follows a multivariate normal distribution. We will use the observations {y(x i ),1 ≤ i ≤ n}, where n is the total number of observations to predict the values for a set of new locations {x i ′,1 ≤ i ≤ n′} without observations, {y′(x i ′),1 ≤ i ≤ n′}. Krige (1951) states that the best prediction, in terms of minimizing the prediction variance, for y at any unobserved location x i ′ can be expressed as a linear superposition of the observed values, that is, ŷ(x i ′) = ∑a i y (x i )+a 0 , and these optimal coefficients, {a i ,0 ≤ i ≤ n} can be estimated from the observed data.
The spatial field of interest (the auroral map in this study) can be decomposed into a combination of a spatially varying mean μ(x), a spatially correlated field, g(x), and a spatially uncorrelated error term ϵ(x), which represents measurement uncertainties:
In our auroral application, the spatial mean function μ(x) can be retrieved from the empirical model z(x). As indicated earlier, there tends to be a systematic discrepancy between the values from the empirical model and observation. Here, we assume a scaling factor d to account for such a multiplicative bias, that is,
The majority of the spatial predictability is achieved through the spatial random field, g(x), which characterizes the detailed spatial variations of aurora. In this work, we take a spatial basis function approach by decomposing g(x) onto a series of predefined basis functions {ϕ j (x),1 ≤ j ≤ m}, that is, g(x) = ∑c j ϕ j (x), where c j is the coefficient of j th basis function and m is the total number of basis functions (see Section 3.2 for further details about the basis functions). The coefficient vector c=(c 1 ,c 2 ,…,cm) jointly follows a multivariate normal distribution with mean zero and covariance matrix Q -1 (therefore, Q represents the distribution's precision matrix). As a result, {g(x), x ∈ X} is a zero mean Gaussian process, and the covariance function takes the following form: where ρ is the spatial marginal variance of the process of interest. The detailed description for Q -1 is given in Supporting Information S1.
In terms of parameter estimation and spatial prediction, we will use matrix notation to simplify the presentation. First, we write the basis functions evaluated at the observed locations into an n × m matrix ϕ such that ϕ ij = ϕ j (x i ), the value of the j th basis function at x i . We use a vector x to denote the observed spatial locations, that is, x=(x 1 ,x 2 ,…,x n ). Then, we have g(x) = ϕc and the covariance matrix of cov (g(x), g (x′)) = ρϕQ -1 ϕ T .
Second, we stack all auroral observations {y (x i ),1 ≤ i ≤ n} and errors {ϵ(x i ),1 ≤ i ≤ n} into vectors y and ϵ, respectively. Since we assume the errors are spatially uncorrelated, the covariance matrix of ϵ is σ 2 W -1 , where W -1 is a diagonal error covariance matrix and σ 2 is a scaling factor of the error term. We also stack the empirical model at each location {z (x i ),1 ≤ i ≤ n} as vector Z. We can now write the model (Equation 1) in the following matrix form:
Here, y follows a multivariate normal (MVN) distribution with a mean of Zd and a covariance matrix
In terms of parameter estimation, we will need to estimate the fixed scaling constant d and the spatially varying effects at the observed locations c based on observations y and their spatial locations x. The best estimates of d and the conditional distribution of c can be obtained via the standard results of generalized least squares (Cressie, 1993), which are
where M λ = ϕQ -1 ϕ T +λW -1 and λ = σ 2 /ρ. Then, the estimate of c is set to the conditional mean
and the variance of ĉ is
Therefore, the predictions (conditional mean and variances) of ŷ′ at new locations are
where the primes on ŷ′, Z′, and ϕ′ indicate that the prediction can be taken at different locations from the input data. A more detailed derivation of d̂ and ĉ can be found in Supporting Information S1.
In summary, our goal is to predict the values y′ at unobserved locations x′ (with corresponding empirical model output Z′ as a predictor) and to quantify the prediction uncertainty. Equation 9gives the prediction of conditional mean ŷ′, and the associated prediction uncertainty is the square root of the diagonal terms in Equation 10. In real applications, the calculations of variances are relatively computationally expensive. Therefore, the variances at each spatial location are usually approximated by the sample variance of independent draws from the conditional distribution of ŷ′, given available observations (Monte Carlo method).
The auroral data in the magnetic latitude (ϕ) and local time (t) coordinates are mapped to the modeling coordinates, which are stretched spherical surface coordinates, using the following equations:
The setup of the basis functions is on the x-y plane. Following Nychka et al. (2015), the basis functions ϕ j (x) are chosen as compactly supported radial basis functions (RBF) φ, which are bell-shaped curves with a common width θ
where u j (1 ≤ j ≤ m) is the center of RBFs. Typically, u j is equidistant, that is, Δu = u j -u j-1 is a constant that represents the grid size (also referred to as the model resolution), and all {u j ,1 ≤ j ≤ m} of the same θ form a regular grid map covering the whole domain, which consists of one level of RBFs. The number of RBFs at this level m is related to the grid size Δu approximately by the reciprocal rule mΔu 2 = domain size. Since the latitudinal direction is directly mapped but the longitudinal direction is scaled in this coordinate, we simply refer the resolution of 1° (without differentiating between the latitudinal and longitudinal directions) to the spacing of RBFs by π/180. Take the auroral modeling for example, in the high-latitude region over 50° MLAT, the domain size is (40°×2) 2 = 6,400. In terms of 1° modeling resolution (Δu = 1°, basis functions separated by 1°), approximately 6,400 basis functions are used (m = 6,400).
For the multiresolution fitting, RBFs of different Δu and θ can be combined into a large basis set (see Figure 3 for a three-level setup of RBFs; Figure 3a shows a 1D case and Figure 3b shows a 2D case). In this sense, we relabel m with m l , Δu with Δu l , and θ with θ l with l representing the number of levels. These parameters can take different values across different levels, which lead to different resolutions. The multilevel reconstruction of the spatial variation field is then written as
where c j, l is the coefficient of the j th RBF at l th level. L is the total number of levels, which is a critical parameter in describing the multiresolution properties of the basis functions. For typical usage, θ l is set as a fixed multiple of Δu l (greater than one) to allow for an overlapping of RBFs at every point. Both L and m l (equivalently, Δu l ) can be adjusted to obtain basis function maps of different scales. Higher-resolution RBFs have more free parameters (c j ) to simulate the details of the input data, and they are expected to provide more small-scale structures of aurora than the lower resolution RBFs.
In our auroral modeling setup, the modeling domain is a 2D square over the high-latitude region. To have m l basis functions for the 2D map in the lth level, we distribute N l = √m l basis functions on each side. We choose to double N l every time as we go from a coarse to a fine level, so the overall number of basis functions (m l ) approximately increases by a factor of 4. By recalling that the grid size and the number of basis functions are related by the reciprocal rule, the grid size is approximately halved with increasing L and the fitting resolution is doubled. From Nyquist's theorem, we can simply take the smallest resolvable scale of our model approximately as the double of the grid size (2Δu l ). In this study, the number of levels (L) and the number of basis functions on each side at the coarsest level (N c ) are selected as the fundamental parameters to control the modeling resolution. In this sense, we define low-resolution modeling as L = 1, N c = 13 (N 1 = 13), which means the fitting resolution is 6.2° (modeling domain is 80°), and the resolvable scale is 12.3°. Similarly, the medium resolution is defined as L = 2, N c = 15 (N 1,2 = (15,30), 30 basis functions on each side at the finest level, the fitting resolution is 2.6°, the resolvable scale is 5.3°) and the high resolution is L = 3, N c = 25 (N 1,2,3 =(25, 50, 100), 100 basis functions on each side at the finest level, the fitting resolution is 0.8°, the resolvable scale is 1.6°).
Even though the number of RBFs roughly quadruples if we increase L, the overall computations do not grow exponentially with the number of levels. Since we formulate the precision matrix Q at each level to be sparse, the overall precision matrix is still sparse (see in Supporting Information S1 for details). Therefore, the matrix calculation does not increase cubically with the total number of the matrix elements (e.g., Gaussian elimination) but only linearly with the nonzero elements. Thus, the increase of levels of basis functions leads to a moderate increase of the whole computation.
Before feeding the SSUSI and THEMIS observations and empirical model into Lattice Kriging, two preprocessing steps are implemented. First, even after we collect 20 min of SSUSI data to form the satellite binned map (Section 2), its spatial coverage is still limited. Therefore, an interpolated satellite map using the data 1 hr before and 1 hr after the modeling time is generated to enlarge the spatial coverage and used as the fourth data source (details in Section 4.1.1). Second, we assign larger weights to higher-fidelity data (i.e., SSUSI and THEMIS observations) and smaller ones to lower-fidelity data (the interpolated satellite data and empirical model) such that the former two data sources dominate the fitting results while the latter two only play roles in the regions where observations are missing. The weighting in Lattice Kriging is realized by attributing different sampling ratios to different data sources. The low-fidelity data are downsampled to decrease their sampling rates and equivalently the weights in the fitting (details in Section 4.1.2). After the preprocessing, we use Lattice Kriging to synthesize all four data sources to generate auroral maps at all locations (Section 4.2) and produce the intermediate result. Due to the smoothing effect inherent in the fitting procedure, Lattice Kriging causes spreading and introduces nonzero values in regions with no aurora such as the polar cap and tends to smear out the auroral boundaries. This solicits a postprocessing weighting method (KNN: K nearest neighbors) for a mitigation (details in Section 4.3). Figure 4 provides a flow chart of these procedures.
Figure 5 illustrates an example of the satellite data by 20 min binning centered around 11:40 UT (Figure 5d), and the satellite map after the linear interpolation with time (Figure 5c). The data used for the interpolation are collected within a 2 hr window from 10:40-11:40 UT shown as the "1 hr before" data in Figure 5a and from 11:40-12:40 UT shown as the "1 hr after" data in Figure 5b, respectively. Compared to Figure 5d, the interpolated map (Figure 5c) shows similar results if the data being interpolated are within the 20 min window such as the region around the dawn (∼06 MLT). The similarity originates from the proximity in time for the temporal interpolation. For the cases that the satellite data are available within the 2 hr but not the 20 min window, the binning method would not show anything while the interpolation can fill up the aurora such as in a significant portion of the dusk region where a few auroral arcs are seen (∼18 MLT in Figure 5c). Even though the assumption that the aurora should change linearly during this period does not necessarily represent the truth, the interpolated results (such as their magnitude) are still closer to reality than the empirical model. Note that compared with the relatively instantaneous observations (e.g., binned satellite and ground-based data), the interpolated data are downsampled (Section 4.1.2) to ensure that they do not override the 20 min binned data when they both exist in the same regions. A similar interpolation method is used in Wu et al. (2020), which better simulates the TTEIL during the storm time than using the empirical auroral drivers. In this study, the similarity and correlation between the auroral activities separated by over 2 hr are thought to be weak, so the linear interpolation is conducted within the 2 hr window.
As discussed earlier, due to the different fidelities of the data sources (observation > satellite interpolation > empirical model), we attribute different weights to them by controlling the data sampling ratios. For simplicity, we refer to the satellite 20 min binned data as "satellite data," and the interpolated results as "satellite interpolation" or "satellite interpolated data." Sampling ratios of 1 meaning no downsampling are assigned to the satellite and ground-based data (r sat = r grd = 1). The ratios for the satellite interpolation and empirical model whose original spatial grids are the same as the satellite data are r int = 1/3 and r emp = 1/20, respectively. Considering that the spatial resolution of ground-based data is higher than that of the satellite (Section 2), the absolute sampling ratios are satellite: ground-based: satellite interpolation: empirical model = 1:2:1/3:1/20. If they overlap in the same region, their weights follow the sequence of ground-based data > satellite data > satellite interpolation > empirical model. The ratios are adjustable depending on the data quality and application purposes. 6f). Due to the lowest fidelity of the empirical model, it has the lowest data sampling density (Figure 6h), therefore its information is assimilated mainly in the regions where the other three data sources are not available.
Feeding the preprocessed data (Figures 6e-6h) as the inputs to Lattice Kriging (y in the formula of Section 3.1), we obtain the assimilation results in Figure 7 for 11:50 UT, which corresponds to the "Intermediate Results" in Figure 4 (before the KNN postprocessing method being applied). We adopt three levels (L = 3) and the numbers of basis functions from coarse to fine grids are N 1,2,3 = 25, 50, and 100, respectively. For the dawn and dusk sectors, the assimilation result mainly resembles the SSUSI observations (Figure 6e); THEMIS data (Figure 6g) contribute to the midnight sector. For the premidnight sector (21-24 MLT) where no observations are available, the assimilated aurora follows the empirical model. Figure 7a shows the prediction of the conditional mean (Equation 9 in Section 3.1). The mesoscale structures including the auroral arcs in SSUSI data and the hot spots spreading in the midnight sector measured by THEMIS ASIs are largely maintained.
The uncertainty of the assimilated energy fluxes depends on the uncertainties of the data sources. Based on the error assessment of the historical data, the uncertainty of SSUSI data can be taken as ∼15% of the measurement, the uncertainty of the THEMIS observation can be taken as ∼20% of the data itself (Gabrielse et al., 2021). Since the interpolated SSUSI data has lower fidelity, an uncertainty of 30% is assigned to the interpolated result. The uncertainty of the empirical model is chosen to be 100% of its value as a proxy since no related information is available yet. These uncertainty terms are used as inputs to ϵ in Equation 3 (Section 3.1). The standard deviation/ fitting uncertainty is then calculated following Equation 10 and shown in Figure 7b. The uncertainties are considerably smaller than the predictions of the means and smaller in the regions with observations than those without observations, reflecting the data constraints.
Despite the similarity between the assimilated auroral map and the input data (Figures 6 and7a), the former appears blurry and small energy fluxes spread into the polar cap and subauroral regions where the observations show no aurora in the input data. Even though the fine structure such as auroral arcs are retained, the peak values in the assimilated map are also lower than the real observations. In other words, the Lattice Kriging model introduces a smoothing effect, which causes leakage to the regions without aurora and reduced auroral peaks. A possible explanation is that Gaussian process models (including Lattice Kriging) use a distance weighted mean strategy to attribute contributions from input data at different locations. In addition, the specific covariance structure used in the model predicts the variances at two nearby locations with similar magnitudes. Therefore, Gaussian process models rarely predict extremely high or low values, which makes the overall spatial predictions smoother than the data. Also, the smooth spatial structures of the basis functions tend to create a smooth representation for the spatial process.
To suppress such smoothing effects, a postprocessing step relying on the KNN algorithm is applied. The same algorithm is used in Syrjäsuo andDonovan (2002, 2004) to automatically classify different types of aurorae from ASIs. We use this algorithm to identify the likelihood of having aurora for each location and eliminate the low-likelihood points. KNN is a common classification method widely used in machine learning. It relies on the assumption that the points in the same category share similar features and lie closely in the feature space. Therefore, a straight-forward way to divide points into different categories is by grouping them in the feature space justified by distance. Given a set of labeled training data, we calculate the distance of a new point to all training data points and pick up the k nearest points. We identify the category of these k points in the feature space, and the new point belongs to the category with the most members.
Based on our data set, we define the auroral activity with energy flux higher than 2 mW/m 2 as significant and set it as 1 in the feature space, otherwise it is insignificant (0). The threshold of 2 mW/m 2 is chosen based on trial and error. It is suitable to identify substantial auroral activity, avoid contamination from low fidelity data (typically on the order of 0.2 mW/m 2 ), and effectively maintain the fine structure. The preprocessed data (Figures 6e-6h) with their significant/insignificant labels (1/0) are used as the training data for KNN. For each location on the fitted map, the k-nearest points to the training data are identified. Assuming the number of data labeled as 1 is n 1 , then a ratio of n 1 /k is calculated, which represents the percentage of the k-nearest points falling into the category of significant aurora. This ratio is used as the weighting coefficient for this location. By doing so, a coefficient matrix with the same dimension as the intermediate result is formed, and their multiplication leads to a weighting process producing the final results of the auroral assimilation (Figure 4).
In Figure 8, we display the post processing results with k = 10 at 11:50 UT. The weighting coefficients from KNN are shown in Figure 8a. The intermediate results from Lattice Kriging are shown in Figure 8b (same as Figure 7a), and the final assimilation by multiplying Figures 8a and8b is given in Figure 8c, where we see the smearing of energy fluxes into the polar cap is largely suppressed. In the polar cap region where the preprocessed data clearly indicate that there are no auroral activities, the k-nearest points all fall into the feature space of 0, thus the spreading values in the polar cap region are effectively removed by multiplying a KNN weighting coefficient of 0. This can also help removing the isolated points (the ambient areas show no aurora) that may be due to measurement noise. In the auroral region around midnight, THEMIS ASI observations indicate that there is strong auroral activity, and KNN labels are mostly 1 so as the weighting coefficients, therefore, the fitting results in the auroral region are kept. In the dawn and dusk regions, the strengths of auroral activities vary so the feature space consists of both 0 and 1, and the resulting weighting coefficients are between 0 and 1. The multiplication of the weighting coefficients and the intermediate results then helps to decrease the aurora if the ambient region does not show enough significant auroral activities. This process sharpens the auroral boundary and to some extent corrects the smoothing effect caused by Lattice Kriging. The overall auroral structures become more comparable to real observations since the training data set in KNN relies on real observations. The usage of a larger k involves more points in a larger area to be weighted and introduces a smoother structure than a smaller k.
In Figure 9, we show the comparison of a simple padding of the satellite and ground-based observations with the final auroral assimilation at 11:50 UT. The padding results show an obvious discontinuity and sharp cut-off at the boundaries among the satellite data, ground-based data, and the regions without observations (Figure 9a), which largely disappear in Figure 9b. The data assimilation effectively removes the boundary discontinuity and combines different data sources more coherently than the padding. A trade-off for such coherence and continuity is the reduced peak magnitude of aurora, which cannot be corrected by the KNN postprocessing step. Nevertheless, the mesoscale aurora is largely retained including the auroral arcs, which significantly improves the reproduction of the real-time behavior of aurora compared to the empirical model. We provide a movie showing the time evolution of Figure 9 in the SI, which illustrates that the dynamic evolution of aurora with time is also captured.
As mentioned in Section 3.2, the generated auroral maps of different scales can be obtained by tuning the number of fitting levels and the number of basis functions in each level (L and N). The resolution increases and the resolvable scale becomes smaller when we increase L and N. Figure 10 shows the assimilated auroral maps at three different scales at 11:50 UT. The parameters to generate these three auroral maps are L = 1, N 1 = 13, k = 30 for large scale; L = 2, N 1,2 =(15,30), k = 20 for medium scale; and L = 3, N 1,2,3 =(25, 50, 100), k = 10 for small scale. Figures 10a-10c show the auroral energy fluxes while Figures 10d-10f show the mean energy maps from large to small scales and, equivalently, low to high resolutions. From low to high resolutions, the assimilated aurora becomes more fine-structured, and the peak values increase. The auroral arcs in the dusk sector are distinct in the high-resolution results but absent in the low-resolution ones.
To study how the data-assimilated drivers improve the simulation of I-T models, and how different scales of aurora impact the I-T system, we run TIEGCM with different auroral maps. TIEGCM is a global 3D numerical model that simulates the coupled thermosphere/ionosphere system from ∼97 to ∼600 km altitude. It self-consistently solves the fully coupled nonlinear, hydrodynamic, thermodynamic, and continuity equations of the neutral gas, the ion and electron energy equations, the O + continuity equation and ion chemistry, and the neutral wind dynamo (Qian et al., 2014;Richmond et al., 1992). In the default setup, the high-latitude drivers such as aurora and electric fields (or electric potentials) are specified as empirical models (e.g., Heelis et al., 1982;Roble & Ridley, 1987;Weimer, 2005). In our TIEGCM runs, the time-varying SuperDARN electric potential pattern, which is derived from a spherical harmonic fitting (SHF) of line-of-sight (LOS) ion velocities (Ruohoniemi & Baker, 1998) is used as a driver for electric fields. The electron precipitation pattern created in this study together with the Zhang and Paxton (2008) model are used to specify auroral particle precipitation in TIEGCM. The spatial resolution of TIEGCM is 1.25°𝐴𝐴 × 1.25°𝐴𝐴 × 1/4 scale height in latitude 𝐴𝐴 × longitude 𝐴𝐴 × altitude (Dang et al., 2018(Dang et al., , 2021)). Realistic Kp and F10.7 are used in all simulations. The time step of the TIEGCM simulation is 10 s. Diagnostic outputs are saved every 5 min.
We perform three different model runs, and the only differences among them are the auroral energy flux and mean energy inputs. These three drivers are the empirical auroral model and the assimilated aurora at low and high resolutions (Figures 11a-11c). In Run 1, the Zhang and Paxton (2008) empirical model is used as the auroral input to TIEGCM; In Runs 2 and 3, low-and high-resolution auroral patterns created in this study are used. In all runs, the high-latitude electric field input is the SuperDARN potential pattern. Since the auroral particle precipitation affects the ionization rate and therefore the electron density, we show TECs from these three runs at 11:50 UT in Figures 11d-11f and compare with the global navigation satellite system (GNSS) observations (Figure 11g). GNSS TEC is measured by the trans-ionospheric propagation time difference between two different radio frequencies from the GNSS satellite to the dual-frequency GNSS receiver. This propagation delay difference is directly proportional to the line integral of the electron density (Vierinen et al., 2016).
Compared with the TEC results driven by the empirical model (Figure 11d), the significant changes after we apply the auroral assimilation maps to drive the TIEGCM are the TEC enhancement (by a factor of ∼2) in the midnight sector where the SSUSI and THEMIS observations weight in (black rectangles in Figures 11e and11f). The changes from low to high resolutions are noticeable in TEC as more mesoscale structures are seen in the high resolution. We also compare the storm-quiet time TEC differences in Figures 11h-11j. From an observational perspective, the differential TEC is obtained by subtracting the TEC 24 hr before the targeted storm-time (Figure 11j), which corresponds to 11:50 UT on 19 February 2014. From the modeling perspective, we use the Run 1 result, which does not involve data assimilation and only shows the large-scale pattern as a proxy for the quiet-time response. Figures 11h and11i demonstrate the differential TECs from TIEGCM simulations by subtracting (d) from (e and f). The regions and magnitudes of TEC enhancements from data assimilation are in agreement with the observations, which means that the data assimilation can be used to better simulate the mesoscale ionospheric responses to auroral precipitation. The model simulations with data assimilation capture the locations of strong TEC responses more precisely than the one driven by the empirical model. The differential TECs also show comparable enhancements, which indicates the robustness of our auroral assimilation method and the resulting improvement.
It is noted that the original data resolutions of both DMSP SSUSI and THEMIS ASI data are much higher than the TIEGCM. Data assimilation can match the observation to a large extent but would be still limited by the I-T model, which incorporates it as an input. To further simulate small-scale processes and make better use of the data assimilation, the resolutions of the I-T models need to be improved. Moreover, the corresponding physics down to small scales also needs to be considered. Nevertheless, this work highlights the substantial changes from using the empirical model to data-assimilated aurora as drivers to simulate the responses of the I-T system, which is essential to better understand and predict the impacts of realistic and localized magnetospheric energy deposition.
We introduce a multiresolution Gaussian process model (Lattice Kriging) to self-consistently synthesize various data sources (satellite data, ground-based data, and empirical model) for the auroral assimilation for the first time. This model assumes that the auroral activity follows a Gaussian process. It uses the available data to estimate the fitting coefficients of the basis functions within the Kriging theory framework, and then uses these coefficients to project the estimation to the whole high-latitude region. The multilevel (or multiresolution) capability is fulfilled by distributing different levels of basis functions with different resolutions, such that different scales of aurora can be assimilated, which facilitates the study of multiscale processes like aurora.
To customize the Lattice Kriging model to auroral assimilation, we introduce two preprocessing steps and one postprocessing step. First, we interpolate the satellite data temporally to expand the spatial coverage at a particular time. The interpolated satellite data and empirical model (low-fidelity data) are then downsampled to decrease their weightings and ensure that the assimilation results are dominated by the satellite and ground-based observations (high-fidelity data) where the low-and high-fidelity data overlap. These four data sources (satellite and ground-based observations, satellite interpolation, and empirical model) are fed into Lattice Kriging to obtain the intermediate results. Due to an inherent smoothing effect of the fitting procedure, which smears out auroral boundaries and introduces nonzero values in the regions with no aurora (such as polar cap), we generate a postprocessing weighting map using KNN trained by observations to mitigate these issues. The KNN weighting coefficient indicates how likely one location has significant auroral activity. These coefficients are multiplied to the intermediate fitting results to eliminate the isolated points likely caused by measurement noises and unrealistic spreading values produced in the intermediate Lattice Kriging modeling. The reduced peak values of aurora due to the smoothing effect, however, are difficult to be compensated. Compared with the simple padding of satellite and ground-based observations, the auroral assimilation model can effectively remove the discontinuity at the boundaries of different datasets.
We use the 20 February 2014 case (a moderate geomagnetic condition) as an example to demonstrate the assimilation procedures and generate the energy flux and mean energy maps with three different scales. The large-scale maps corresponding to the low-resolution fitting miss mesoscale structures such as auroral arcs, while the small-scale maps corresponding to the high-resolution fitting show mesoscale structures that more closely resemble observations. We then apply the assimilation maps of low and high resolutions to drive TIEGCM to study the impacts of different scales on TEC. In general, the TEC in the auroral region (especially midnight sector) shows substantial enhancement that better matches observations after data assimilation due to the increased level of auroral particle precipitation and ionization. High-resolution auroral precipitation maps also produce mesoscale structures of TEC. Overall, the TIEGCM simulations highlight the importance of implementing realistic aurora as one of the magnetospheric drivers to model the mesoscale electrodynamics at high latitudes.
Despite the noticeable advantages in fusing real data to simulate the mesoscale auroral structures, the current auroral assimilation model has the following limitations, which may need further improvements. In the data preprocessing step, we combine SSUSI data over 20 min to form a snapshot, then we interpolate over a 2 hr period to expand the data coverage. One limitation from these steps is that the information of the development of aurora within that time interval is lost. This may lead to the distortion of the auroral oval if aurora changes very rapidly during the 20 min interval. For example, if a substorm onset occurred between the time when the dawnand dusk-side oval were observed, the dawn-side oval would appear expanded, while the dusk side would appear contracted. While each side of the oval might appear as narrow features, they would be coming from completely different auroral ovals. Combining observations from such a situation might lead to a double edge structure, which is purely due to the binning of SSUSI data. It is difficult to mitigate this issue by the technique itself and more data are needed to fundamentally solve it.
In the spatial modeling of aurora, when specifying the covariance structure, there are also simplified assumptions that may not represent real observations. First, the covariance matrix used here is derived from a Gaussian Markov random field (which assumes two locations are correlated only if they are adjacent, Nychka et al., 2015). In the real world, however, even distant auroral regions can be correlated if the aurora in these regions is generated from a closely connected region in the magnetotail (Nishimura, Lessard et al., 2020). An additional term indicating the medium-to-large range correlation needs to be included in the covariance matrix to describe the realistic auroral characteristics (Cousins et al., 2013;Matsuo, 2020). Second, the auroral activity may not follow the Gaussian distribution as assumed in this study. Since the different high-latitude regions connect to different regions in the magnetosphere, the auroral distributions may not be the same and they may deviate from Gaussian distribution due to the pitch angle diffusion and other wave-particle interactions (Nishimura, Lyons, et al., 2020). Therefore, the mathematical formulation may need to be modified based on a nonGaussian process model. Still, Gaussian statistics has good properties for fast computation, such as the sparse matrix calculation as aforementioned, which satisfies as a starting point. The improvements of covariance matrix and distribution type solicit statistical studies of aurora, which is beyond the scope of this study. Third, the current methodology can efficiently combine various data sources and conduct spatial fitting in a coherent way thus the boundary issue disappears, however, it is not an auroral prediction model and cannot be used to predict auroral activity for next time steps. The prediction of aurora may be achieved by the machine learning technique training a large amount of historical data. For our case, real-time observational data are still the key to drive models to produce realistic I-T responses. It is worth pointing out that there may be discrepancies between satellite and ground observations. For this event, the magnitudes from these two types of observations match to a large extent despite discrepancies in some small-scale structures. However, in case these two data sources deviate, it is necessary to examine the data quality and perform downsampling to the one with lower fidelity.
It is worth mentioning that the Lattice Kriging modeling is not limited to scalar field assimilation. Wu and Lu (2022) have extended it to assimilate vector fields such as electric fields under the curl-free condition and obtained the results with much smaller errors than the global SHF using the SuperDARN data. The fundamental principles are the same except that for the assimilation of electric fields, we need to project the basis functions of electrical potential (scalar) to electric fields (vector) and then project them onto the LOS direction, along which the observations are actually made (SuperDARN measures LOS ion drifts). Such extended capability makes the Lattice Kriging modeling appealing not only for the scalar assimilation such as GNSS TEC measurements, but also for wind measurements such as those from the Ionospheric Connection Explorer (ICON) in the future.