Calibrated and Enhanced NRLMSIS 2.0 Model with Uncertainty Quantification

Exospheric temperatures were obtained through a binary search, deriving the exospheric temperature required in NRLMSIS 2.0 to match satellite density estimates [29, 32]. The density estimates originate from the following sources – CHAMP 2001: Doornbos [6], CHAMP 2002-2010: Mehta et al. [18], GRACE 2002-2009: Mehta et al. [18], GRACE 2010: Sutton [24], Swarm A 2013-2018: Astafyeva et al. [1], and Swarm B 2012-2016: Astafyeva et al. [1]. Note that only GRACE-A measurements are used due to the similarity of the GRACE-A and GRACE-B orbits. The CHAMP and GRACE density estimates originate from accelerometer measurements and span an altitude range of 300–535 km while the Swarm A and B density estimates are obtained from GPS data and span an altitude range of 437–545 km. The cadence of the satellite estimates are 10 s, 5 s, 30 s, and 30 s for CHAMP, GRACE, Swarm A, and Swarm B, respectively. There are over 82 million samples total for model development and evaluation.

To run NRLMSIS 2.0, the following space weather inputs are required: F_10.7, F_81c, and Ap. F_10.7 is a measure of the 10.7 cm solar radio flux and acts as a good proxy for solar EUV emissions while F_81c is merely an 81-day centered average of the proxy. The Ap index is a measure of daily global geomagnetic activity. NRLMSIS 2.0 also has a storm-time flag that allows additional 3-hourly ap values: ap, ap₃, ap₆, ap₉, ap_12-33, and ap_33-57. The subscripts indicate the number of hours prior to epoch that ap value represents; when there are two numbers (e.g. 33-57), the 3-hourly index values are averaged over this many hours prior to epoch. NRLMSIS also has temporal inputs (e.g. time of day, day of year).

To develop a machine-learned exospheric temperature model with a high temporal cadence, we expand the input set. To account for solar activity, the model receives F_10.7, S_10.7, M_10.7, and Y_10.7, all accounting for different forms of solar emissions that affect different regions of the thermosphere [26, 2]. The ML model also uses inputs from EXTEMPLAR, particularly S_N, S_S, $\Delta T$. The two S inputs are Poynting flux totals in the Northern and Southern hemispheres [30, 31]. $\Delta T$ is a parameter derived by Weimer et al. [28, 32] and represents exospheric temperature perturbations; it is a function of Poynting flux and simulated nitric oxide emissions.

For further representation of geomagnetic activity, a time history of SYM-H is used, similar to the use of storm-time ap indices for NRLMSIS 2.0. SYM-H represents disturbances to the geomagnetic field and has similar characteristics to Dst. A benefit over Dst is its 1-min cadence [11]. The SYM-H inputs are as follows: SYM-H, SYM-H_0-3, SYM-H_3-6, SYM-H_6-9, SYM-H_9-12, SYM-H_12-33, SYM-H_33-57. Due to the model’s use of in-situ satellite estimates, it takes location as an input. The local solar time (LST) is transformed using sine and cosine functions to make it continuous about midnight,

The model also uses satellite latitude (LAT). To account for temporal dependencies, day of year (t₁ and t₂) and universal time (UT, t₃ and t₄) are transformed through similar functions,

The ML model drivers are the 21 space weather, spatial, and temporal inputs mentioned in the previous section. Each sample has an exospheric temperature estimate which will serve as the target. To reduce the variance of the output ($T_{\infty }$), we use a logarithmic transformation making the output $lo{g}_{10}{T}_{\infty }$. The final step to create ML-ready data is to perform standard normalization. This will center the data and provide a unit standard deviation for each input and output with respect to the statistics of the training data,

where x represents each input and the output, $\mu$ is the mean of that quantity over the training set, and $\sigma$ is the standard deviation of that quantity over the training set. $\mu$ and $\sigma$ must be saved for downstream use of the model.

The 81 million samples are split into training, validation, and tests sets to achieve an 80%–10%–10% distribution. Licata et al. [14] split a similar dataset to have long segments – on the order of months or years. However, the resulting model was not well-generalized across the sets. When comparing EXTEMPLAR and a previous version of MSIS over the same three time periods, the error statistics also varied. Therefore, we use a different approach to data splitting. Starting with the first sample, eight weeks are used for training, the following week is used for validation, and the subsequent week is used for testing. This pattern is repeated until the end of the dataset. This approach forces similar solar cycle and spatial coverage across the three sets. Concurrently, there is a significant number of samples within each segment providing temporally disjoint segments throughout the 17 year time-span of the dataset. Since each satellite has a different cadence and there are different numbers of satellites providing measurements at a given time, the number of samples in each segment varies. In the training, validation, and test sets, the number of samples varies from 22,454–1,450,380, 12,990–181,423, and 12,888–181,422, respectively. This means that there are between 25,878 and 362,845 samples separating training segments.

A key consideration in ML model development is the choice for the loss (or objective) function. This will be the quantity that the algorithm will try to minimize or maximize during the training phase. For this work, we use the negative logarithm of predictive density (NLPD) loss function due to previous ML modeling results [12, 15]. NLPD is given as

where $y$, $\mu$, and $\sigma$ are the ground truth, mean prediction, and predicted standard deviation, respectively. This loss functions provides the capability for uncertainty prediction since the model can now make $\sigma$ predictions and have them incorporated in the loss function without labels. The ML model will directly predict $\mu$ and $\sigma$ for $lo{g}_{10}{T}_{\infty }$ meaning it has two outputs. The output node for $\mu$ uses a linear activation function since the normalized values can take on any value (unbounded). The output node for $\sigma$ must be positive, monotonically increasing, and have no upper-bound. Therefore, we proceed with the "softplus" activation function: $f\left(x\right)=ln(1+e^x)$.

With the inputs/output(s) selected, normalization complete, loss function chosen, and output layer determined, the only remaining task for model development is to select an architecture and train a model. To accomplish this, we defined a scheme and hyperparameter space using Keras Tuner [19] as with previous work [14, 12, 15]. Table 1 shows the parameters used and ranges provided to the tuner. Each hidden layer can have its own unique number of neurons, activation function, and dropout rate. Note: "trials" refers to number of model architectures, initial points is the number of randomly selected architectures prior to Bayesian optimization, and repeats refers to re-initializing weights and retraining.

Since there are over 65 million training samples in the dataset, we only provide the tuner with a subset of this data. The tuner uses 1 million randomly selected samples from the training set and 200,000 randomly selected samples from the validation set. Each model trained by the tuner will run for 50 training iterations (or epochs) with a batch size of 4,096. Upon completion, the 10 best models are saved based on the validation loss values. All 10 models are evaluated (see Section 2.5), and the best performing one is used as a base architecture for full training. The best architecture from the tuner is displayed in Table 2.

When comparing model to satellite densities, we use the mean absolute error (MAE) metric in percentage form. This metric is chosen due to its intuitive nature. To assess the quality of the ML uncertainty estimates, we use a calibration error score (CES) [12, 15]. This metric assesses how close the uncertainty estimates are to input prediction intervals (e.g. how close a 95% prediction interval is to containing 95% of true samples). The equations for MAE and CES are provided in Equations 5a and 5b, respectively.

In Equation 5a, ${\rho }_{satellite}$ is the estimate of density from a given satellite measurement (e.g. from CHAMP, GRACE, etc.). ${\rho }_{model}$ is the density from the model being evaluated. In Equation 5b, $m$ is the number of prediction intervals being tested, $n_{out}$ is the number of model outputs (for the ML model $n_{out}$ = 1), $p\left(z_{i,j}\right)$ is a given prediction interval, and $p\left({\hat{z}}_{i,j}\right)$ is the percentage of samples captured by the model’s uncertainty estimates with the given prediction interval. In this work, the prediction intervals span from 5% to 99% in increments of 5%.

Figure 2 contains uncertainty profiles for MSIS-UQ at CHAMP and GRACE locations on May 13, 2007. There are panels for species density, temperature, mass density, relative uncertainty, and satellite position. Please reference the figure caption and Section 2.5.2 for details.

Panels (a) and (b) show the species density profiles at CHAMP and GRACE positions, respectively. The uncertainty bounds provide valuable information on the impact of exospheric temperature uncertainty on the uncertainty of local species. For example, one can investigate the Oxygen (O) to Helium (He) transition for various locations and conditions. Panel (a) shows that at CHAMP’s position, this transition is occurring somewhere in the region of 507 to 552 km (1-$\sigma$) while at GRACE’s position, the transition may occur between 688 and 738 km (1-$\sigma$). Other insights can be gained such as which species are most impacted by exospheric temperature at a given location/altitude. Note: only 1-$\sigma$ bounds are shown here to prevent artifacts at low-values caused by the semi-logarithmic scale. The scale also causes the bounds to appear to be not-centered about the mean.

Panels (c), (d), and (e) provide information on the local temperature and mass density with uncertainty. In panel (c), MSIS-UQ severely shifts the exospheric temperature prediction and brings it closer to the estimates of CHAMP and GRACE; in both cases NRLMSIS 2.0 overpredicts temperature. The uncertainty in temperature is unobservable below 130 km and grows until it reaches the asymptotic temperature between 250 and 300 km. The uncertainty bounds and mean remain unchanged above these altitudes. Note that the CHAMP and GRACE temperature estimates are for $T_{\infty }$ but we show them at their current altitude as the temperature has converged. In panel (d), we see different trends in mass density. Again the uncertainty is minimal below approximately 200 km and begins to increase for a few hundred kilometers. The overprediction of temperature in NRLMSIS 2.0 results in higher than observed density by CHAMP and GRACE around 350 and 475 km, respectively. MSIS-UQ provides a more accurate density predictions at the satellite locations. Panel (e) shows the 1-$\sigma$ uncertainty with respect to mean density. This shows different model behavior between the two locations. At CHAMP’s location, the uncertainty increases to 24% around 460 km and decreased until around 700 km where it settles to 9%. At GRACE’s location, the uncertainty continues to increase until it reaches 18% at 600 km where it begins to decrease.

Geomagnetic storms remain a challenge in modeling the thermosphere. Even further, post-storm characteristics vary between models. To highlight this, we show orbit-average density along CHAMP’s orbit during the 2003 Halloween storm for CHAMP, TIE-GCM, NRLMSIS 2.0, JB2008, HASDM, EXTEMPLAR, and MSIS-UQ in Figure 3. Time history SYM-H inputs for MSIS-UQ are displayed in panel (c) for reference.

In the pre-storm period (10/28–10/29), there is significant variability in model outputs. During the first peak of the storm, the models show similar trends but are mostly above the CHAMP density estimates. In the lull between the two peaks (10/30–10/31), all models show density decreases but to varying extents. NRLMSIS 2.0 shows the smallest density decay during this period, but due to the exospheric temperature corrections in MSIS-UQ, the density falls to CHAMP levels. For the second storm, all models overestimate density with the exception of TIE-GCM and JB2008. In the post-storm period (10/31–11/03), NRLMSIS 2.0 does not follow the trend observed by CHAMP – a sudden and severe decrease in density. JB2008 shows this briefly before overpredicting density to the extent of NRLMSIS 2.0. The MSIS-UQ follows the post-storm overcooling trends observed by CHAMP, showing the impact the exospheric temperature corrections to NRLMSIS 2.0 have. We now conduct the overcooling study from Licata et al. [13] to observe the characteristics in the model without effects from other drivers. These results are shown in Figure 4, and correspond to Figure 2 in Licata et al. [13].

Panels (a)–(d) in Figure 4 show the density ratios for NRLMSIS 2.0 while increasing each time-series ap value independently. We observe linear relationships between ap and density for all of the time-series inputs. The current 3-hourly ap value has the strongest impact on density and, of these four locations, never generates a density ratio greater than 2.25. These results show that in NRLMSIS 2.0, the relative importance of each ap value becomes less important the further back from epoch it is with ap_36-57 having minimal impact on density. The exception to this statement is ap_12-13. At the day equator (panel (b)), it causes the second-most positive density ratio. When ap_12-13 is very large, it causes a strong relative increase at all locations considered here. This could explain the overprediction observed in Figure 3. The second row (panels (e)–(h)) shows the results for MSIS-UQ or the change in NRLMSIS 2.0 when using SYM-H time-histories to provide $T_{\infty }$ correction. The trends shown in these panels contradict many observations when using NRLMSIS 2.0 alone. For example, at 3/4 locations, ap₃ causes the largest density ratios, even being nearly twice as large at the night equator. MSIS-UQ also shows a nonlinear relationship between geomagnetic activity and density which is not necessarily seen in NRLMSIS 2.0 alone.

MSIS-UQ also enforces the idea of negative density ratios – density decreasing while the geomagnetic drivers are increasing. This is most pronounced at the two pole locations. When the least recent geomagnetic drivers (ap_12-33 and ap_36-57) become large, the density becomes up to 25% lower than when they are set to zero. This overcooling was seen in Figure 3 and is observed in the satellite density data, therefore becoming present in MSIS-UQ. Another interesting trend is seen at low levels on geomagnetic activity particularly in panels (g) and (h). When any of the ap values increase from 0 up to 50–100, the density decreases. This seems counter-intuitive but could be caused by the approach of the study. When ap is being considered, for example, the ap and SYM-H values are set to 56 and -50 nT, respectively. When ap = 0, this would represent the time immediately after moderate geomagnetic activity while ap = 56 would represent sustained moderate geomagnetic activity. The model shows that when the conditions abruptly return to quiet values, the density increases – likely only temporarily. The bottom row (panels (i)–(l)) show the temperature ratios from MSIS-UQ corresponding to the middle row. The general trends are the same between temperature and density; however, the difference comes from the magnitude. The relative changes in temperature result in much stronger changes in density. There are negative temperature ratios, but they are much less prominent due to the consistent scaling.

In SM2, we provide a movie of global density evolution for the 2003 Halloween storm at 400 km. This contains density maps for TIE-GCM, NRLMSIS 2.0, EXTEMPLAR, and MSIS-UQ. Prior to the storm, all models have a similar diurnal structure. This is quickly disrupted at the onset of the storm when both TIE-GCM and MSIS-UQ show strong auroral density enhancements. Throughout the storm, these are the only two models that display significant global disruptions to the structure of the thermosphere at this altitude. The thermosphere expands and contracts both longitudinally and latitudinally for these models. EXTEMPLAR – a linear model based on the same dataset as MSIS-UQ – models an abnormal thermosphere, but the movement is more longitudinal or latitudinal with some wave-like oscillations about the equator. During this storm, NRLMSIS 2.0 exhibits general increases in density, but the structure is well-preserved. There are some slow-moving variations although it does not resemble the dynamics exhibited by the physics-based model.