© M.V. Codrescu et al., Published by EDP Sciences 2022

1 Introduction

Tools for ensemble modeling and data assimilation in terrestrial weather and ocean science have been developed and are in operational use. The use of an ensemble Kalman Filter (enKF) in space weather is also not new. Codrescu et al. (2004) published a paper on neutral composition enKF assimilation 2004. Although the neutral composition was recognized to be one of the most important factors in ionospheric simulations during storms (Chartier et al., 2013), the lack of neutral composition measurements has prevented the operational implementation of enKFs in space weather products and services. However, the importance of enKF for space weather research has been recognized, enKFs have been used in research, and papers have been published (Solomentsev et al., 2012; Morozov et al., 2013; Hsu et al., 2014; Chartier et al., 2016).

Today, other kinds of assimilation models based on Gauss Markov (GM) Kalman Filter (KF) processes are more popular in operational settings (Schunk et al., 2004; Spencer et al., 2004; Fuller-Rowell et al., 2006; Jee et al., 2010; Jakowski et al., 2011; Borries et al., 2015). GM KF assimilation schemes are based on stationary predefined covariance matrixes that work best if large amounts of data are available to overwhelm the empirical background model. The sudden availability of large amounts of total electron content (TEC) measurements from Global Navigation Satellite Systems (GNSS) signals made the ionosphere GM KF assimilation schemes feasible to implement. GM KF based assimilation schemes can be very good at ionosphere specification for past events, especially during quiet or moderately disturbed geomagnetic conditions when large amounts of data are available. However, in real-time environments, they can suffer from data starvation and do not have forecasting capabilities beyond persistence with a predefined evolution toward climatology.

During disturbed conditions, GM KF schemes have difficulty because their predefined quiet-time covariance matrixes do not keep pace with the changing system. To obtain the appropriate covariance matrix during disturbed conditions, it would be necessary to perform variational analysis (Rockafellar & Wets, 1998) during every assimilation time step. However, for assimilation schemes with hundreds of thousands to millions of state elements, performing variational analysis every assimilation time step (15–30 min) is not practical, and the covariance matrix needs to be estimated in some other way. An estimation of the covariance matrix using Monte Carlo methods was first proposed by Evensen (1994) as the enKF.

In a previous paper (Codrescu et al., 2018), we discussed assimilation results for total mass density and showed that assimilating measurements from one satellite improves the model results globally during quiet conditions. Our goal in this study is two-part. The first goal is to show that the Thermosphere Ionosphere Data Assimilation (TIDA) system responds to measurements and uncertainty during a significant geomagnetic disturbance. The second goal is to highlight that the proposed assimilation scheme can use neutral density measurements to estimate the forcing even in the absence or degradation of L1 measurements.

In this paper, in Section 2, we discuss the dominant processes that make the TI covariance matrix non-stationary during disturbed geomagnetic conditions. The paper continues with an experiment using TIDA to assimilate GRACE-A, GRACE-B, and CHAMP neutral density measurements during the extreme geomagnetic 2003 Halloween storm and two additional storms: September 27 – October 2, 2002, and July 26–30, 2004. Section 3 overviews the experiment setup. The measurement sources are described in Section 4, results are presented in Section 5, and finally, we conclude in Section 6.

2 The thermosphere ionosphere system

The global neutral density and composition of the thermosphere depend on system forcing and the interaction with the ionosphere (Fuller-Rowell et al., 1994). The global electron and ion density structure, roughly from 50 to 1000 km altitude, are at any given time the result of a dynamic equilibrium between plasma production, loss, and transport, processes controlled to a large extent by neutral composition and neutral winds (Fuller-Rowell et al., 1997). The processes that affect the neutral composition, density, and winds and the production, loss, and transport of plasma are highly variable on timescales of minutes to years, and their relative importance can change as a function of location on the globe, universal time, storm commencement time, season, solar cycle, waves propagating from below, and the previous state of the ionosphere–thermosphere–magnetosphere system (Sarris, 2019). On short timescales, the variations are controlled by a set of external energy inputs that include solar radiation absorption at a variety of wavelengths, solar energetic proton deposition, solar wind energy transfer through the magnetosphere that depends on the density and speed of the solar wind and the magnitude and orientation of the interplanetary magnetic field (IMF) (Codrescu et al., 2012), and waves propagating from below (Heelis & Maute, 2020). The influence of waves propagating from below will not be discussed further in this paper as their amplitudes, and phases change slowly relative to the duration of a geomagnetic storm, and their influence can be taken into account by an appropriate lower boundary condition in the assimilation background model.

2.1 The system during quiet geomagnetic conditions

While the thermosphere ionosphere (TI) is never in a true steady state, it can reach quasi-steady-state conditions if the system inputs are quasi-constant over some period of time (days), as it happens during prolonged quiet geomagnetic conditions (Codrescu et al., 2008). Under steady-state conditions, the system energy input is balanced by cooling through CO₂, and NO infrared emissions and diurnally reproducible patterns can be observed in most system state variables.

While in prolonged quiet periods, the TI system reaches a quasi-steady state that can last many days (Roble, 1992). The state is said to be in a diurnally reproducible pattern. Under such conditions, statistical models of high-latitude convection electric fields (Weimer, 2005), particle precipitation (Fuller-Rowell & Evans, 1987), and solar EUV fluxes based on correlation with the F10.7 cm flux measurements (Hinteregger et al., 1981) are good enough to give acceptable model results when used in physics-based numerical models of the system. In addition, empirical, statistical models of the ionosphere (Nava et al., 2008; Bilitza, 2018) or thermosphere climatology (Picone et al., 2002) are also good during quasi-steady-state conditions. This means that during quiet conditions, the system can be modeled with a high level of confidence.

During quasi-steady-state conditions, a global equilibrium is established between heating due to solar radiation absorption on the dayside, Joule heating at high latitudes, and infrared cooling due to NO and CO₂. As a consequence, a diurnally reproducible global neutral temperature structure and circulation are established, and a relatively stable global neutral composition structure is maintained (Killeen et al., 1997). This state of the thermosphere produces a diurnally reproducible global dynamo electric field pattern (Richmond, 1989) which, in association with the stable prompt penetration electric field pattern of magnetospheric origin (Manoj & Maus, 2012), produce a diurnally reproducible ionosphere. During geomagnetically quiet conditions, the energy input from solar radiation absorption dominates the system energy input (Mlynczak et al., 2016).

In the upper atmosphere, at around 300 km altitude, where the peak electron density normally occurs, the temperature structure establishes day–night pressure gradients that drive the global circulation neutral winds (Hedin et al., 1991). The winds blow from the dayside towards the nightside, both east and west, and over the poles from the 14:00 local time sector. Mostly molecular species are present below 150 km altitude and atomic species above. The difference between the Earth’s geographic and magnetic poles contributes to the diurnal variation in the global temperature, winds, and composition structure in the sun-fixed reference frame and at any location on the globe.

2.2 The system during disturbed geomagnetic conditions

During geomagnetic storms, changes in external system forcing cause large increases in the magnitude and distribution of Joule heating, in auroral particle precipitation total energy and its distribution, and in momentum transfer to neutrals. The total energy input into the TI system at high-latitudes increases dramatically and can become larger than solar radiation heating. This has dramatic consequences for the global neutral winds and composition. Furthermore, changes in neutral winds cause changes in the dynamo-electric field pattern. Due to the tight coupling between the ionosphere and thermosphere, the changes are reflected in the ionosphere and feedback to the neutral state through ion drag, momentum transfer, heat transfer, and other mechanisms (Fuller-Rowell et al., 1994, 1996).

Empirical models are not appropriate to represent the state of the TI system during severe geomagnetic storms. Numerical models of the system also suffer during disturbed conditions because the statistical models used for forcing have very large uncertainties, and this results in unacceptable uncertainties in model simulation results. Although disturbed conditions may happen only for a small percentage of the time, it is during disturbed conditions that accurate modeling is most important.

Storm Joule heating occurs at high latitudes at about 110–115 km altitude, where molecular species (O₂ and N₂) dominate. The additional heating changes the pressure gradients and superimposes a storm circulation on top of the quiet time winds. Strong upward vertical winds are driven above the auroral zone heating area and meridional winds away from the heating area. During the storm, the effect of the divergent wind field is an upwelling of molecular species and composition change. The effect of the heating and vertical wind is the formation of a “composition bulge” where molecular constituents normally found at lower heights are moved upward. The bulge tends to spread towards lower latitudes (Proelss & von Zahn, 1978; Fuller-Rowell et al., 1994) under the influence of the meridional winds.

During small disturbances, the storm-induced meridional winds are overwhelmed by the quiet time circulation on the dayside but add to it on the nightside resulting in a distortion of the composition bulge relative to the shape of the auroral zone. During large storms, meridional winds can turn equatorward at mid-latitudes even on the dayside. The size, shape, and position of the auroral zone and the composition bulge are highly variable functions of solar wind density and speed, magnitude and orientation of the interplanetary magnetic field, and storm time.

During a geomagnetic storm, molecular species are transported up (upwelling) above 200 km, where they displace lighter atomic species (mostly atomic Oxygen). The lighter species are then forced equatorward by storm meridional winds. To balance the pressure gradients and close the storm circulation, the lighter species are transported down (downwelling) by storm vertical winds at some distance equatorward of the heating area, and a return poleward flow of molecular species takes place at lower altitudes (Fuller-Rowell et al., 1994, 1996). The position of the upwelling area can vary in time but depends mostly on the intensity of the storm, while the position of the downwelling area is a more complicated function of storm intensity, duration, and storm time profile.

The changes in neutral dynamics and composition cause important changes in the TI system (Proelss & von Zahn, 1978; Fuller-Rowell et al., 1994, 1996; Burns et al., 1995). In the composition bulge, plasma production decreases due to the decreased atomic Oxygen densities while the loss of plasma increases through charge exchange with molecular species followed by recombination. Poleward meridional winds and westward E-fields can also contribute to plasma loss. It is not uncommon to have less than half the quiet time plasma peak density (NmF2) in an area covered by the composition bulge following a geomagnetic storm. This is what is called the negative ionosphere storm effect. The global neutral composition can take more than 36 h to recover after a storm.

In the downwelling area, the increased atomic Oxygen causes increased plasma production and reduced loss resulting in increased plasma density. This is a positive ionospheric storm. Equatorward meridional winds and eastward E-fields can also contribute to the positive phase.

At a given location, positive storm effects are seen first, though not always, followed by negative storm effects (Codrescu et al., 1992). The ionospheric changes are most pronounced in the F2 layer but can be significant in the whole ionosphere, especially during long geomagnetic storms. Meridional winds driven by storms can cross the equator and propagate in the opposite hemisphere.

The system forcing uncertainties are greatly amplified during geomagnetic storms. Small scale electric field variability can increase dramatically, change the spatial distribution of energy input, and more than double the Joule heating that results from the convection average electric fields (Codrescu et al., 2000). The thermosphere ionosphere coupling and the dynamic changes produced by storms in each of the thermosphere and ionosphere subsystems make the modeling difficult and lead to unacceptably large simulation uncertainties.

2.3 The path forward

There are two ways to mitigate the large forcing uncertainties during storms: measure the forcing, i.e., measure the electric fields and particle precipitation at the necessary grid points every few minutes, or use any available system measurements to estimate an appropriate forcing using data assimilation scheme. As long as properly measuring the forcing is not possible, the only practical solution is a sophisticated data assimilation process, that is, a data assimilation scheme that can take advantage of all available TI measurements to reduce the external forcing uncertainty while also improving model-data comparisons.

Developing an assimilation scheme that can take advantage of a variety of TI measurements is a major challenge because the external forcing acts in multiple ways and on different timescales and because the system contains feedback loops with storm-time-dependent gains. These complications make the correlations between model variables non-stationary or, in other words, state and time-dependent. Since the covariance matrix depends on the present state of the system, it has to be calculated or estimated again during each assimilation time step.

One practical way to obtain the covariance is by Monte Carlo estimation methods (Evensen, 2003). An appropriate number of members of the background model (an ensemble) is run with representative forcing variations, and statistics of their results are used to estimate a covariance matrix. The accuracy of the estimated covariance is a function of the number of members relative to the number of degrees of freedom of the system, the forcing distribution over the ensemble members, and the error of the system estimation, at the time of the estimation.

3 Methodology

In this study, we use an ensemble Kalman Filter (EnKF) implemented by the TIDA software to simulate the thermosphere–ionosphere during several geomagnetic storms: the 2003 storm: October 26–30, 2003; the 2004 storm: July 26–30, 2004; and the 2002 storm: September 27 – October 2, 2002. The 2003 storm is an extreme geomagnetic storm event known as the 2003 Halloween storm. The other two events are not as major but are provided to show that the improvement to neutral density specification is not an artifact of the lost forcing observations during the 2003 storm. TIDA is described in detail by (Codrescu et al., 2018), although the software was not called TIDA at the time. The EnKF approach to data assimilation relies on an ensemble of model realizations to derive a Monte-Carlo approximation for the state covariance matrix at each assimilation time step. A unique feature of TIDA is its ability to augment the Kalman state with the external system forcing parameters. Consequently, assimilating observations will estimate not only a correction to the typical state variables but also produces an estimated “correction” to the forcing parameters. The forcing corrections resulting from one assimilation time step are applied during the forecast phase of the following assimilation time step. This is similar to the iterative reinitialization technique (Sutton, 2018) but executed on much shorter time scales under the EnKF framework. This inference allows the system to run even in the absence of L1 system forcing measurements.

For the results presented here, we have used an ensemble with 75 members. The assimilation time step is 30 min, and the model time step is one minute. At each assimilation time step, the ensemble members are distributed by randomly sampling their forcing from Gaussian distributions centered on the previously estimated system forcing parameters. The Kalman state vector contains the following global fields: neutral temperature, constituent mixing ratios, meridional and zonal neutral winds, and mean molecular mass. The Kalman state vector is also augmented with the system forcing parameters: solar wind magnetic field magnitude in the YZ plane, the solar wind angle in the YZ plane, the solar wind speed and the solar wind number density, and the F10.7 cm radio flux value. The state vector contains over 191 thousand elements in total. In addition to the 75-member ensemble, TIDA runs two other versions of the Coupled Thermosphere Ionosphere Plasmasphere electrodynamics (CTIPe) model. The reference version, marked “Reference” on the plots, illustrates the operational model results in the absence of data assimilation. The special member is used to forecast the state using the forcing inferred during the previous assimilation time step.

TIDA uses the CTIPe general circulation model for the background physics-based numerical model. CTIPe has a long history going back to the early 1980s (Fuller-Rowell & Rees, 1980), has been running in real-time (Codrescu et al., 2012) for more than ten years, and has been tested during both quiet and disturbed conditions (Codrescu et al., 2012; Fedrizzi et al., 2012; Negrea et al., 2012; Fernandez-Gomez et al., 2019). CTIPe was transitioned into operations at the Space Weather Prediction Center (SWPC) in November 2019. Results from the SWPC real-time operational run are available at: http://ccmc-swpc.s3-website-us-east-1.amazonaws.com/plots.html.

CTIPe requires a high latitude potential pattern, a particle precipitation pattern, the average F10.7 cm radio flux values over the previous 41 days, and a seasonally dependent lower boundary condition. The high latitude potential pattern is obtained from the Weimer (Weimer, 2005) empirical model. The Weimer model requires the solar wind magnetic field magnitude in the YZ plane, the solar wind angle in the YZ plane, the solar wind speed, and the solar wind number density to produce a convection pattern. The particle precipitation is inferred from the solar wind dynamic pressure using the above variables. The intensity, characteristic energy, and spatial extent of particle precipitation are represented by the activity level and hemispheric power.

4 Data

The neutral density measurements assimilated in this experiment are derived from very sensitive accelerometers flown on the GRACE-A, GRACE-B, and CHAMP satellites (Tapley et al., 2004; Sutton, 2011) in polar orbit. For the 2002 storm, the GRACE satellites are at an altitude of ~500 km, offset by ~60° in longitude (4-h local time difference) from CHAMP at an altitude of ~400 km. For the 2003 storm, the GRACE satellites, at an altitude of ~480 km, are offset by ~50° in longitude (3.33-h local time difference) from CHAMP at an altitude of ~400 km. For the 2004 storm, the GRACE satellites were at an altitude of ~480 km, offset by ~70° in longitude (4.67-h local time difference) from CHAMP at an altitude of ~390 km. The observations are provided in 3-degree latitudinal bins averages. No bias correction or pre-processing was applied before assimilation. Furthermore, the estimated uncertainty provided with the measurements was used directly. For example, for CHAMP, the minimum and maximum uncertainty during the Halloween storm were 1.52 × 10⁻¹³ kg/m³ and 1.27 × 10⁻¹² kg/m³. For GRACE-A, the min and max uncertainties were 2.14 × 10⁻¹⁴ kg/m³ and 3.50 × 10⁻¹³ kg/m³. For GRACE-B, these values were 2.25 × 10⁻¹⁴ kg/m³ and 3.49 × 10⁻¹³ kg/m³.

GRACE-A and B fly close together with an along track separation of approximately 220 km. It might look like including both is redundant, however, using both satellites increases the signal-to-noise ratio, and both are useful because they can, at times, occupy different voxels within the model.

The advanced composition explorer (ACE) satellite monitors the solar wind plasma parameters from the first Lagrange point (L1). These observations and observations of the F10.7 cm solar radio flux are the inputs to the model. Unfortunately, due to a solar energetic particle event during the 2003 Halloween storm, significant portions of the ACE data are bad quality and not usable. We have retrieved the available data from the NASA OMNI service (https://omniweb.gsfc.nasa.gov/form/sc_merge_min1.html), which provides solar wind values propagated to Earth’s magnetopause. During data gaps and outages, the most recent valid solar wind driver value is repeated across the gap.

5 Results

We first discuss results from the 2003 Halloween run where neutral density measurements from all three satellites, GRACE-A, GRACE-B, and CHAMP, were assimilated. During most assimilation time steps, about 30 measurements/satellites were available for assimilation. Given that fewer than 100 measurements were assimilated in each 30-min assimilation time step and that normal input parameters are not available for many hours during this period, the results are surprisingly good for such an extreme space weather event. Later in this paper, we will discuss model measurement comparisons when measurements from only one satellite are assimilated at a time to demonstrate that this is not a lucky coincidence but a consequence of the strongly forced nature of the TI system and the large scale coherence of the neutral density features in both the real system and in the CTIPe model. This conclusion is also supported by the results for the other two storms, as illustrated in Table 1. We note that all results presented in this paper are along the orbit of the moving satellites and are not orbit averaged.

Figure 1 is a scatter plot of model versus measurement results over the Halloween storms (October 27–31, 2003). The left column labeled “Reference” illustrates the CTIPe model results without data assimilation. This is approximately what would have been produced by a CTIPe real-time operational run. The middle column shows “Forecast” model results from TIDA before the present assimilation time step measurements are assimilated, while the right column shows the final “Analysis” TIDA neutral density results. The rows are for the three data sets assimilated in this run: CHAMP (top row), GRACE-A (middle row), and GRACE-B (bottom row). Forecast results can be thought of as assimilation results when only measurements older than 30 min are available for assimilation, while analysis means that measurements up to the simulation time are available.

Fig. 1 Scatter plot of model versus measurement over the 5-day Halloween 2003 storm period for the case when all three neutral density data sets are assimilated. Left column is the reference model specification with no assimilation, middle column is the forecast specification, and right column is the analysis specification.

The large overestimation of neutral density in the reference run (left column) is caused by the loss of forcing measurements on October 30. The L1 measurements needed for the convection and particle precipitation patterns were compromised by an ongoing solar energetic proton event. Since the operational run must produce results even in the absence of input measurements, the model reuses the last available forcing measurements again and again until new forcing measurements become available. The last available L1 measurements for October 30 caused a large overestimation of the Joule heating in the model when repeated, resulting in a much larger modeled neutral density for corresponding measurements.

Figure 2 illustrates the forcing parameters for the sub-period 29–30 October 2003. The reference forcing (blue) and the inferred forcing (orange) are shown. The repeated values across the L1 observation outage sustain a large magnitude for the B field in the YZ plane, and a geoeffective angle together with the large solar wind velocity results in the large overestimation of Joule heating in the reference run. In TIDA, the Kalman state is augmented with the forcing parameters, such that assimilation at each time step produces an inferred forcing. The inferred forcing from one assimilation time step is used to force the special member for the forecast step during the following assimilation time step.

Fig. 2 System forcing parameters over two of the five simulation days. “B” is the magnitude of the solar wind magnetic field in the YZ plane, “angle” is the angle of the solar wind magnetic field in the YZ plane, “SW speed” is the solar wind speed, “SW density” is the number density of the solar wind, activity level, and hemispheric power are descriptors of the intensity and spatial extent of the statistical particle precipitation pattern. The observed system forcing parameters used in the reference run is shown in blue. During large gaps of missing data, the previous valid value is repeated. The orange line shows the inferred forcing obtained from the assimilation of GRACE-A, GRACE-B, and CHAMP neutral density measurements into the special member. The estimated forcing from one assimilation time step is used over the next. The inferred drivers are not necessarily physically realistic but instead represent a set of drivers that produce the best agreement with observations.

Neutral density measurements alone do not assure a unique solution for model forcing. The estimated forcing parameters presented in Figure 2 are the best estimate for the model forcing, given the distribution of assimilated measurements and their uncertainties and the physics captured in the CTIPe model. Even when available, observations of the system forcing are not assimilated because the goal is not to reproduce realistic drivers but rather to infer forcing parameters that coerce the model into an agreement with observations.

Changes in neutral density at the height of a satellite can result from a change in temperature, a change in neutral composition, or a combination of both. Neutral density measurements alone do not contain enough information to allow TIDA to uniquely determine the cause of a model-data discrepancy and properly correct it during each assimilation time step. This and the continuous change in the position of the satellite measurements over the globe result in the ruggedness of the inferred system forcing. Additional measurements of temperature and/or neutral composition are expected to reduce the variability of the inferred forcing and further improve model-data comparisons for neutral density.

Figure 3 shows the measured (yellow), reference (red), forecast (blue), and analysis (black) neutral density values for October 29 and 30, 2003. The overestimation of density by the reference run is again obvious on October 30. On the other hand, at times, TIDA slightly underestimates the neutral density. This is most obvious for CHAMP at the end of October 30. We do not have a good explanation for this effect and plan to investigate it further.

Fig. 3 Neutral density observed (yellow), reference (red), forecast (blue), and analysis (black) for the case when all three satellite data sets were assimilated on October 29 and 30, 2003. The reference provides a benchmark for the model performance without assimilation. Forecast and analysis specifications are produced before, and after data from the current assimilation, interval have been assimilated. In this case, the reference performs poorly because of missing solar wind inputs.

Scatter plots like Figure 1 for the cases when a single satellite data set is assimilated very similar show only a little more spread compared to the 3 satellite assimilation cases, are not discussed here but can be seen in Appendix A. The difference in forcing parameters inferred by TIDA for the four assimilation cases are minor, do not bring any revelations, and again are not presented here. The difference in TIDA results when assimilating all data sets at once or one at a time can best be seen in Table 1.

Table 1 summarizes the assimilation results for three different storm periods (October 26–30, 2003, July 26–30, 2004, and September 27 – October 2, 2002), using the normalized root mean square deviation . The result of the wrong inputs forced upon the reference run by the absence of L1 measurements due to the proton event contamination in the 2003 storm is obvious and unacceptable for an operational run. Using previous neutral density measurements, i.e., measurements made before the present 30 minute assimilation time step, TIDA can infer better system forcing parameters and reduce the along the orbit NRMSD from over 150% (CHAMP reference) to below 27% (CHAMP forecast) and from over 230% (GRACE-A and B reference) to less than 36% (GRACE-A and B forecast). The assimilation of the measurements during the assimilation time step further reduces the NRMSD to 16.2% for CHAMP and to less than 16% for GRACE-A and B.

Table 1 shows that better neutral density results are obtained along the orbit of all satellites during the three storms even when measurements from only one satellite are assimilated. The fact that the TIDA results are good for any satellite when assimilating all measurements is evidence that the data sets are reasonably consistent with each other, have acceptable biases relative to each other, and that the model approximates the physics of the system well enough to be able to improve the results far away from the location of any given measurement.

The 2004 storm illustrates an extreme case of bias in the reference run, while the 2002 storm shows a more usual bias for CTIPe, in the absence of data assimilation. The L1 measurements used by the reference run for these two storms do not suffer from proton contamination, and the large NRMSDs are caused by the use of statistical particle precipitation and cross-polar cap potential patterns in the model forcing and by the missing physics in the background CTIPe model.

6 Conclusions

We have presented results over an exceptionally large “2003 Halloween” storm and two other smaller storms using TIDA, an enKF data assimilation software package adapted for strongly externally-forced systems. The non-stationary nature of the system encourages the estimation of the covariance matrix during each assimilation time step. In addition, because the forcing uncertainties are the largest source of uncertainty for model results, we estimate the system forcing parameters at each assimilation time step.

Including the external system forcing parameters in the Kalman state of TIDA allows their estimation based on all system measurements and results in considerable improvement in modeling results. Given enough system measurements, a set of forcing parameters can be inferred even in the absence of the observed forcing measurements (L1 solar wind and F10.7 values). This can assure uninterrupted modeling operations even when model inputs are not available.

TIDA, the implementation of the enKF used in this study, demonstrates the considerable improvement potential of data assimilation for TI system modeling. A small number of measurements, fewer than 100 neutral density values assimilated during each 30-min assimilation time step over 5 days, can reduce the model data NRMSD along the orbit by factors of 7–10 versus the reference with bad forcing. Furthermore, we have provided evidence that assimilation of the neutral density from a single satellite improves the specification globally during disturbed conditions.

Acknowledgments

The supporting data for this paper can be accessed at https://doi.org/10.5281/zenodo.4467407. The CTIPe model used in this paper is available from the corresponding author upon reasonable request. The TIDA system is available for license from Vector Space, LLC, contact Stefan Codrescu ssmmcc1@gmail.com. The neutral density measurements are provided by Dr. Eric Sutton and are publicly available from: https://drive.google.com/drive/folders/1oki9WOe3J4U3Pn3yebNaO1a3v-0m5Aho. The solar wind measurements used in this study were retrieved from the NASA OMNI service https://omniweb.gsfc.nasa.gov/form/sc_merge_min1.html. The authors are thankful for the constructive feedback received from the anonymous peer reviewers at the Journal of Space Weather and Space Climate. The editor thanks Liangliang Yuan and an anonymous reviewer for their assistance in evaluating this paper.

Appendix A

Individual satellite scatter plots

Fig. A.1. Scatter plot of model versus measurement when only the CHAMP satellite data set is assimilated. Left column reference state (no assimilation), middle column forecast state, right column analysis state, over the 5-day Halloween 2003 storm period.

Fig. A.2. Scatter plot of model versus measurement when only the GRACE-A satellite data set is assimilated. Left column reference state (no assimilation), middle column forecast state, right column analysis state, over the 5-day Halloween 2003 storm period.

Fig. A.3. Scatter plot of model versus measurement when only the GRACE-B satellite data set is assimilated. Left column reference state (no assimilation), middle column forecast state, right column analysis state, over the 5-day Halloween 2003 storm period.

References

All Tables

All Figures

	Fig. 1 Scatter plot of model versus measurement over the 5-day Halloween 2003 storm period for the case when all three neutral density data sets are assimilated. Left column is the reference model specification with no assimilation, middle column is the forecast specification, and right column is the analysis specification.
In the text

Fig. 2 System forcing parameters over two of the five simulation days. “B” is the magnitude of the solar wind magnetic field in the YZ plane, “angle” is the angle of the solar wind magnetic field in the YZ plane, “SW speed” is the solar wind speed, “SW density” is the number density of the solar wind, activity level, and hemispheric power are descriptors of the intensity and spatial extent of the statistical particle precipitation pattern. The observed system forcing parameters used in the reference run is shown in blue. During large gaps of missing data, the previous valid value is repeated. The orange line shows the inferred forcing obtained from the assimilation of GRACE-A, GRACE-B, and CHAMP neutral density measurements into the special member. The estimated forcing from one assimilation time step is used over the next. The inferred drivers are not necessarily physically realistic but instead represent a set of drivers that produce the best agreement with observations.

In the text

Fig. 3 Neutral density observed (yellow), reference (red), forecast (blue), and analysis (black) for the case when all three satellite data sets were assimilated on October 29 and 30, 2003. The reference provides a benchmark for the model performance without assimilation. Forecast and analysis specifications are produced before, and after data from the current assimilation, interval have been assimilated. In this case, the reference performs poorly because of missing solar wind inputs.

In the text

	Fig. A.1. Scatter plot of model versus measurement when only the CHAMP satellite data set is assimilated. Left column reference state (no assimilation), middle column forecast state, right column analysis state, over the 5-day Halloween 2003 storm period.
In the text

	Fig. A.2. Scatter plot of model versus measurement when only the GRACE-A satellite data set is assimilated. Left column reference state (no assimilation), middle column forecast state, right column analysis state, over the 5-day Halloween 2003 storm period.
In the text

	Fig. A.3. Scatter plot of model versus measurement when only the GRACE-B satellite data set is assimilated. Left column reference state (no assimilation), middle column forecast state, right column analysis state, over the 5-day Halloween 2003 storm period.
In the text