# Data Assimilation for Ionospheric Space-Weather Forecasting in the Presence of Model Bias

- School of Mathematical and Statistical Sciences, Arizona State University, Tempe, AZ, United States

The dynamics of many models of physical systems depend on the choices of key parameters. This paper describes the results of some observing system simulation experiments using a first-principles model of the Earth’s ionosphere, the Thermosphere Ionosphere Electrodynamics Global Circulation Model (TIEGCM), which is driven by parameters that describe solar activity, geomagnetic conditions, and the state of the thermosphere. Of particular interest is the response of the ionosphere (and predictions of space weather generally) during geomagnetic storms. Errors in the overall specification of driving parameters for the TIEGCM (and similar dynamical models) may be especially large during geomagnetic storms, because they represent significant perturbations away from more typical interactions of the earth-sun system. Such errors can induce systematic biases in model predictions of the ionospheric state and pose difficulties for data assimilation methods, which attempt to infer the model state vector from a collection of sparse and/or noisy measurements. Typical data assimilation schemes assume that the model produces an unbiased estimate of the truth. This paper tests one potential approach to handle the case where there is some systematic bias in the model outputs. Our focus is on the TIEGCM when it is driven with solar and magnetospheric inputs that are systematically misspecified. We report results from observing system experiments in which synthetic electron density vertical profiles are generated at locations representative of the operational FormoSat-3/COSMIC satellite observing platforms during a moderate (G2, *Kp* = 6) geomagnetic storm event on September 26–27, 2011. The synthetic data are assimilated into the TIEGCM using the Local Ensemble Transform Kalman Filter with a state-augmentation approach to estimate a small set of bias-correction factors. Two representative processes for the time evolution of the bias in the TIEGCM are tested: one in which the bias is constant and another in which the bias has an exponential growth and decay phase in response to strong geomagnetic forcing. We show that even simple approximations of the TIEGCM bias can reduce root-mean-square errors in 1-h forecasts of total electron content (a key ionospheric variable) by 20–45%, compared to no bias correction. These results suggest that our approach is computationally efficient and can be further refined to improve short-term predictions (∼1-h) of ionospheric dynamics during geomagnetic storms.

## 1 Introduction

The ionosphere is a portion of the Earth’s mesosphere, thermosphere, and exosphere, corresponding to altitudes from approximately 60–1,000 km, in which interactions with solar radiation create a plasma consisting of neutral and ionized gases with free electrons. The plasma tends to be concentrated in layers whose composition and thickness vary diurnally. At night, the so-called E and F layers extend, respectively, from altitudes from about 90 to 150 km and from 150 to 500 km. During the daytime, the D layer appears (from about 60 to 90 km) and the F layer bifurcates into two thinner layers,

The development of a space-weather forecasting capability, comparable to that of operational tropospheric models, is a long-term research goal [1]. No systematic overview of ionospheric models or the associated open problems is attempted here; interested readers are referred to [3, 4] and references therein. Present-day models include both empirical descriptions (e.g., the International Reference Ionosphere, or IRI) and physics-based dynamical models (e.g., the Thermosphere Ionosphere Electrodynamics Global Circulation Model, or TIEGCM). The IRI has been under development since the late 1960’s to provide, among other things, monthly means of electron density profiles and plasma temperatures during quiet-time conditions [5], i.e., in the absence of geomagnetic storms that arise from shock waves in the solar wind and/or coronal mass ejections. The TIEGEM, developed by the High-Altitude Observatory at the National Center for Atmospheric Research [6], attempts to provide self-consistent solutions for the magnetohydrodynamics of the upper atmosphere. Dynamical models will be important for space-weather forecasting efforts during storm- as well as quiet-time conditions.

Our paper builds upon recent achievements in the understanding of transport and mixing using theoretical analysis, large-scale numerical simulations, and data analysis, and focuses on multi-scale ionospheric dynamics [1, 2, 7]. The ionosphere is a dynamic mixture of ions, electrons and neutral gases surrounding the earth. Understanding and modeling multi-scale ionospheric dynamics is a daunting but important challenge because of its critical role in space sciences and space weather. Significant progress has been made in developing computational models over the past several decades. However, current modeling capabilities need further development, especially with regard to neutral-plasma interactions driven by neutral fluid dynamics and electrodynamics physical processes over the large range of temporal and spatial scales that characterize the ionosphere. Advances are needed in physics-based predictive modeling and data assimilation to make accurate high-resolution forecasts. In particular, novel computational and data assimilation capabilities are required for operational communication, navigation, remote sensing, imaging and space weather applications.

Non-equilibrium processes, transport, and mixing are exceedingly challenging to study. Their dynamics often involve sharp changes of vector and scalar fields and may also include radiation transport and chemical reactions, and diffusion of species and electric charges, among other effects. Non-equilibrium dynamics of transport and mixing are inhomogeneous, anisotropic, non-local, and statistically unsteady. At macroscopic scales, non-equilibrium transport may lead to self-organization and order, thus offering new opportunities for diagnostics and control. Capturing properties of emerging coherent structures, interfaces and mixing, enabling their accurate description can aid better understanding of the Eulerian and Lagrangian dynamics, and developing new methods of control of non-equilibrium transport in nature and technology. Significant success was recently achieved in understanding of transport and mixing on the sides of theoretical analysis, large-scale numerical simulations, and data analysis. This success opens new opportunities for studies of fundamentals of non-equilibrium dynamics across the scales including Rayleigh-Taylor (RT) instabilities and turbulence [8–10], hydrodynamic RT turbulence mixing [11, 12], and plasma dynamics in the ionosphere [13–18].

Turbulent mixing induced by RT instabilities occurs in settings as varied as exploding stars (supernovae), inertial confinement fusion, and macroscopic flows in fluid dynamics such as ionospheric plasmas [19]. Since the discovery of the plasma instability phenomenon that occurs in the nighttime equatorial F-region ionosphere, and which is revealed by rising plumes identified as large-scale depletions or bubbles, considerable efforts have been made in the development of computer models that simulate the generation and evolution of the convective equatorial ionospheric storms or Equatorial Spread F (ESF) Dynamics [2, 3, 6, 11, 12, 14, 20–32]. Analyses and simulations of primary and secondary RT instabilities in the ESF triggered by the response of plasma density to neutral turbulent dynamics and wave breaking in the lower region of ionosphere were investigated for coupled systems (ions, electrons, neutral winds), thus enabling studies of mesoscale/microscale dynamics for a range of altitudes encompassing ionospheric D, E and F layers. At smaller scales, turbulent mixing induced by hydrodynamic instabilities plays an important role and considerable efforts has focused on discovering universal properties of turbulence in various settings. Neutral-plasma interactions and physical processes at ionospheric meso/microscales is currently a grand challenge for data-driven computational studies of ionospheric dynamics [2].

To provide useful forecasts, dynamical models require an accurate and computationally efficient data assimilation capability. Traditional approaches assume that the associated model produces unbiased forecasts. Several approaches have been suggested to handle model bias in data assimilation systems. Dee and da Silva [33, 34] suggested a two-stage approach that estimates separate gain matrices for the forecast state vector and for the bias correction. Biases in observations, such as satellite radiances, can be handled through Kalman filtering (e.g., [35]). Eyre [20] shows that forecast and analysis biases are functions of the relative weights given to bias-corrected and non-bias-corrected observations and the rate at which the forecast model relaxes to its climatology. Canter *et al.* [36] have explored the possibility of applying stochastic forcing to a forecast model, with an augmented state vector to estimate the model bias, with an application to a global ocean circulation model. Machine-learning approaches also have been tried to observation bias correction [26]. The approach applied in this paper is based on one suggested by Baek et al. [37], which attempts to compensate for the model bias in the handling of the observations, as described below.

Our principal focus is on the TIEGCM. We consider a period corresponding to a moderately strong (G2) geomagnetic storm, which can induce significant perturbations of the ionosphere from quiet-time conditions and affect the operation of satellite communication and power distribution networks. In particular, we are interested in the case where various model driving parameters are misspecified, as may occur during periods of particularly high solar activity [38], resulting in systematic errors in model predictions of electron density and related variables. Our approach extends prior work by Baek et al. [37], who described some approaches to handle model bias in a data assimilation system.

The Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC), and its successor, COSMIC-2 are designed, built, and operated by the United States. National Oceanographic and Atmospheric Administration and other agencies in collaboration with the Taiwan National Space Organization (to whom the satellite systems are known as FormoSat-3 and -7, respectively, and abbreviated as F3/C and F7/C2). The COSMIC networks measure total electron content in the ionosphere using a method called radio occultation [29].

Variations in solar activity determine the rates of ionization and heating within the ionosphere and thermosphere. In cases of high solar activity, the ionosphere and its drivers undergo significant changes: ionization and heating rates rise sharply, geomagnetic activity undergoes significant disturbances (geomagnetic storms), and the thermosphere expands, creating enhanced pressure fields and wind circulations that modify the composition of neutral and ionized gases. (See [39] for an overview of the effects of solar flares on the ionosphere and [40] for an overview of ionospheric storms.)

The modeling of extreme ionospheric events is of practical interest due to potential adverse effects associated with space weather disturbances, which range from satellite communication interruptions [41] and Earth imaging [21] to geomagnetically induced currents, which affect power transmission systems [42]. Although the benefits of modeling the ionosphere and thermosphere and their associated electrodynamical processes with a single general circulation model have been observed [31], other ionospheric drivers are specified with auxiliary empirical models.

The ionosphere is strongly coupled to the magnetosphere in high-latitude regions. The TIEGCM does not include a model for magnetospheric dynamics; the coupling is represented by empirical models [43]. At high altitudes, the TIEGCM uses ionospheric electric fields as provided by the Heelis [23] and Weimer [43] models. The parameterization of these latter models therefore has an important effect on TIEGCM outputs. Inputs to solar and geomagnetic parameterizations usually are measured directly and are physically meaningful, but nowcasting and short-term forecasting of parameters like the *Kp* remain a topic of continuing research [44, 45]. Challenges include limited station observability coverage and limited representation of associated physical mechanisms [46].

We describe two sets of observing system simulation experiments in which synthetic observations of electron density are assimilated into the TIEGCM using various choices of driving parameters (such as for the electric fields), to assess their effect on 1-h forecasts. The dynamics of the ionosphere respond rapidly to changes in solar activity. Even if the initial state vector of an ionospheric model is accurate, its predicted evolution may not be if external, time-evolving drivers of solar activity, geomagnetic conditions [47], and the state of the dynamically coupled thermosphere are incorrectly specified [48]. We find that simple models of the time evolution of the forecast bias can provide significant improvements to 1-h predictions of electron density and related quantities.

Ensemble Kalman Filter (EnKF) approaches have been shown to be effective in estimating the state of the coupled ionosphere-thermosphere system [24, 31] and ionospheric forcing parameters [32]. Improved estimates of electron density may also help improve estimation of other systems that are strongly coupled to the ionosphere, such as the magnetosphere [19]. An overview of parameter estimation strategies in data assimilation is presented in [49], and additional parameter estimation strategies with Ensemble Kalman filters for non-global parameters are discussed in [50]. Geomagnetic storms often have significant and persistent storm-time effects that modify the state of the ionosphere relative to quiet-time conditions [51], such as by changing electron densities in ways that are difficult to predict [52]. Such effects may be another source of systematic model biases that can be difficult to remedy by adjustments to input parameters alone.

In this paper, we propose a strategy to compute spatially varying corrections to compensate for model bias introduced from suboptimal parameterized inputs. Our approach is based on methodology originally proposed in [37]. We assume that a solar and magnetospheric parameter configuration, which is specified through some parameter estimation approach, is already in place. Our strategy treats the global distribution of model bias as a collection of bias correction parameters, which we estimate using an ensemble Kalman filter. The bias corrections are not used to adjust the model state vector directly. Instead, spatially varying bias corrections are applied in the evaluation of the observation operator, with the intent of reducing model bias in its electron density predictions prior to the assimilation of observations.

The application of spatially varying bias corrections extends the degrees of freedom in the predictions made with a given parameterized configuration and may be beneficial in capturing storm-time effects. The corrections are computed independently of parameterized inputs, so this strategy may be applied together with other direct parameter estimation approaches. Furthermore, this strategy helps provide a more gentle state estimation update due to the adjusted forecast predictions, which may help to avoid the introduction of spurious dynamical artifacts from drastic adjustments to the system state that may occur due to strong model bias presence during an extreme event. To our knowledge, this bias correction strategy has not been applied previously to an ionospheric model with operational capabilities.

The proposed bias estimation strategy, which is designed for high-dimensional ionospheric systems that are sparsely observed, is applied in observing system experiments focused on the September 2011 geomagnetic storm to compensate for model bias resulting from suboptimal parameterized solar and magnetospheric inputs. We use the Local Ensemble Transform Kalman Filter (LETKF) [25] to assimilate synthetic electron density vertical profiles into the Thermosphere Ionosphere Electrodynamics Global Circulation Model (TIEGCM). The synthetic observation profile locations are the representative of the COSMIC/FORMOSAT-3 satellite constellation [53] during September 26–27, 2011. The LETKF is a type of ensemble square-root filter that computes its analysis using a low-rank estimate of the forecast covariance matrix. The analysis is computed independently grid point by grid point, by assimilating nearby observations simultaneously. The LETKF has been applied in the ionosphere with an idealized regional model [54] and for space-weather specification during an extreme event with the TIEGCM by [7]. The LETKF has also been used with the Global Ionosphere-Thermosphere Model to estimate solar parameters during periods of low and high solar activity [22, 27].

## 2 Data Assimilation Scheme

The data assimilation algorithm used in this paper is the Local Ensemble Transform Kalman Filter (LETKF). We provide a brief overview of the algorithm below and refer to [25] for a full description. Its application to the TIEGCM is also presented in [7].

The LETKF begins with an ensemble of *k* global forecast (“background”) vectors at a particular model update (“analysis”) time

Consider a model grid point indexed, say, as *L*. We associate a “local window” to *L*, consisting of the space within a prescribed horizontal and vertical distance from *L*. The subscript *L* is used to emphasize that the following quantities are associated with a specific grid point and its associated local window. Let *d* denote the number of components of the model state vector to be analyzed (Often, atmospheric models include passive scalars as well as dynamical variables, so *d* may be less than the dimension of the actual state vector. In the case of the TIEGCM, the variables that are analyzed include electron and ion density and temperature.) Let *L* be a multi-index to identify each model grid point. We denote the background global model state vector for the *j*h ensemble as *L*. The local background ensemble mean is given by *j*th column is given by

The LETKF performs the Kalman filtering step in the column space *S* of

Denote the *j*th local forecast are denoted by the

where *j*th column is given by

The ensemble mean of the analyzed local state at the grid point *L* is an adjustment of the background ensemble mean, consisting of a linear combination of local forecast perturbations given by

where

If *S* with mean and covariance *L*. The minimizer of *J* is given by *S*,

To compensate for model nonlinearity, multiplicative covariance inflation is incorporated in the above formulation through the factor

The ensemble of local analysis perturbations is constructed from the local background perturbations as

where

## 3 Bias Estimation Methodology

We consider the data assimilation problem in the context of a forecast model for which there is some systematic difference between its predicted state

for discrete times

as a finite-dimensional projection of the ionosphere’s evolution to the grid points of the ionospheric forecast model at each time *n*, provided that the initial state of the ionosphere and its drivers is the same. In practice, however, the initial ionosphere-thermosphere state is not known precisely. Additionally, the specification of solar and magnetospheric drivers that produces each

Our approach is based on one suggested by Baek et al. [37], which they called “bias model II,” that assumes the dynamics of the ionospheric model evolve on a different attractor from the “true” dynamics of the ionosphere. In such a case, substituting

where

Generally, the correction

Here the subscript *b* refers to the “background” state vector (i.e., the forecast) and *a* to the “analyzed” state vector (i.e., the output of the data assimilation system).

The data assimilation procedure described in Section 2 proceeds as before, but it begins with an ensemble of augmented state vectors of the form *j* indexes each ensemble member. In lieu of Equation 2, the observation operator at time

In other words, we compensate for the model error by adding a correction term to the model state vector whenever a comparison to observations is made. The data assimilation procedure produces an augmented analysis ensemble of the form

The procedure described above yields a spatially varying state estimate of the model bias. The remaining component of our proposed algorithm is the propagation of the model bias estimates for which one needs a suitable model,

In principle, the bias correction procedure doubles the number of state variables that must be estimated from the same set of observations, which increases the variance in the analyzed fields. On the other hand, not all model biases contribute equally to the forecast error, and it may suffice to apply the bias correction procedure only to some of the components of the model state vector; that is, the components of

## 4 Observing System Experiment Setup

This section describes the common aspects of the observing system experiments described in the next sections. Our setup is a “perfect model” experiment, in which the results of the data assimilation system are obtained with synthetic observations generated from the TIEGCM model with known driving parameters. In this way, the sources of model bias are controlled, the performance of the data assimilation system may be compared to a known truth, and improvements due to the proposed bias correction scheme may readily be quantified. The conclusion section suggests some issues that may be considered prior to any potential operational application.

The numerical experiments simulate a scenario in which an ionospheric forecast is driven with suboptimal specification of magnetospheric and solar inputs during a period of geomagnetic disturbance. Storm-time effects can be difficult to model and predict [44, 45]. For example, the study by [51] discusses the significant enhancement in F-layer electron density peak and altitude during the geomagnetic storm of December 15, 2006. The focus of the observing system experiments in this paper is on the September 26–27, 2011 geomagnetic storm event.

### 4.1 Ionosphere Model Description

The TIEGCM is a three-dimensional nonlinear model of the coupled ionosphere-thermosphere system. It solves the momentum and energy equations for neutral and plasma species and the electrodynamical coupling between the ionosphere and thermosphere. We use version 1.94 in this paper at a geographic resolution of *p* is pressure and *e*.) The pressure levels extend from −7 to 7 with half-scale height resolution and range from about 97 to 600 km, depending on solar conditions.

The main source of ionization in the ionosphere, and thus one of its primary dynamical drivers, is the absorption of solar radiation in the thermosphere, primarily by O_{1} and neutral gases (O_{2} and N_{2}). The TIEGCM represents effects of solar irradiance and its variability through auxiliary empirical models. The default solar input model, which we use in our TIEGCM configuration, is the EUVAC irradiance model [56]. The EUVAC model specifies important ionospheric processes related to solar activity, such as ionization and dissociation rates, and heating of neutral gases, ions, and electrons. The empirical representation of solar irradiance is parameterized with the

Another key driver of ionospheric dynamics is geomagnetic activity. The ionosphere is strongly coupled to the magnetosphere in high-latitude regions and is typically approximated with auxiliary empirical models due its complexity. The magnetospheric input used in our simulations with the TIEGCM is the Heelis model [23], which specifies important high-latitude processes in the ionosphere, such as electric field patterns and auroral energy inputs, for different levels of geomagnetic disturbance.

The *Kp* index, a widely used index derived from the horizontal component of geomagnetic field disturbances, is an input to the Heelis model. It is to calculate *Hp* (Hemispheric Power), which is used in the parameterization of auroral precipitation, and the cross-tail potential (*Cp*), which specifies ion convection patterns in the polar regions. *Kp* indices are provided every 3 h and range in value from 0 to 9, to describe geomagnetic conditions ranging from quiet to extremely disturbed. Historical records of *Kp* indices are provided by the National Oceanic and Atmospheric Administration (NOAA).

### 4.2 Ensemble Generation

We generate an ensemble of 40 forecasts with normally distributed values of *Kp*, which in turn produce estimates of *Hp* and *Cp*. Each parameter is centered around the respective index values published for September 26–27, 2011. The standard deviation for the *Kp* is ±1.0. The respective standard deviations for *Cp* and *Hp* are

**FIGURE 1**. **(A)** Time series of **(B)** cross-tail potential (*Cp*) in kV, and **(C)** hemispheric power (*Hp*) in GW. The horizontal axis is in hours, starting at 00:30 UTC on 26 September and ending at 23:30 UTC on September 27, 2011. The parameters used to drive the truth are shown in green, the parameters used to drive the ensemble of forecasts are shown in pink and their respective ensemble means are shown in red. The deviation of each parameter from those used in the control simulation is shown in gray, and the ensemble mean of the deviations is shown in black.

The control simulation is driven with forcing parameters that have the same temporal evolution as the published indices but are shifted as shown by the green curves in Figures 1A–C. This shift is introduced to simulate the scenario where the solar and magnetospheric inputs used to drive the forecast are a misspecification relative to the “true” state of the ionosphere-thermosphere system. The imposed bias on *Kp* index. The corresponding bias on *Cp* and *Hp* values is shown in Figures 1B,C. Although the bias in the *Kp* index is fixed, the resulting bias in the *Cp* and *Hp* indices varies temporally. Prior to the onset of the geomagnetic storm, which is about 12:30 UTC on 26 September, the bias in the magnetospheric inputs is relatively small and constant but increases considerably over the next 6 h. The strongest bias occurs during the main phase of the geomagnetic storm, which takes place at about 16:30–19:30 UTC on 26 September. As geomagnetic conditions relax over the next 12 h, the bias in the magnetospheric inputs also decreases accordingly.

In reality, of course, the errors between the published and “optimal” model forcing parameters are more complex, but, to our knowledge, have not been fully characterized. Insofar as the model’s response varies nonlinearly with the forcing parameters, and the magnitude of the bias in the magnetospheric inputs resembles the typical uncertainty during geomagnetically disturbed periods, we believe that our choice of ensemble of forcing parameters is as good as any for this illustration of bias correction. Nevertheless, future work should address this question in more detail.

### 4.3 The TIEGCM-LETKF Data Assimilation System

Our TIEGCM-LETKF data assimilation system runs as a discrete forecast/analysis sequence, in which the LETKF is used to assimilate observations available in the time interval

During each analysis calculation, only a subset of the TIEGCM state variables is included in the LETKF state vector and analyzed with the assimilated observations. State variables included in the LETKF state vector are electron density (Ne), neutral temperature (Tn), zonal (Un), and meridional (Vn) components of neutral winds, as well as atomic (O_{1}) and molecular (O_{2}) mass mixing ratios. The same forecast/analysis sequence and choice of analyzed state vector is used in [7]. Similar approaches are also used in other ionospheric data assimilation studies [28, 32].

The covariance inflation scheme used in this paper is the one proposed by [55], which also was used with the LETKF and TIEGCM in the observing system experiments presented in [7]. Other LETKF parameters, such as the ensemble size and localization radii, are also taken to be the same as in [7] to simplify comparisons with the results.

### 4.4 Observing Network

During each analysis step, observations are generated synthetically from the electron density component of the control simulation (“truth”),

Observations obtained through retrievals, such as COSMIC electron density profiles, have also been shown to have significant uncertainties associated with the inversion process, particularly at low altitudes [30, 58]. We include an approximation of the inversion error associated with these observations, as done in [7, 24]. In addition to standard observing errors, retrieval errors are correlated in possibly complicated ways, depending on the nature of the satellite observing platform, but here we treat them as independent. To compensate, the standard deviation of

## 5 Numerical Experiment 1: Persistent Bias Evolution

In this section, we study how well the proposed bias estimation strategy can estimate the spatiotemporal evolution of model bias in the electron density field when the bias evolution operator is persistence, that is,

At the analysis time

Figure 2A shows the global distribution of electron density bias, in units of

**FIGURE 2**. **(A)** Global maps of electron density bias in units of **(B)** Vertical structure of electron density bias at a fixed 77.5° N latitude at 11:30, 15:30, 19:30 and 23:30 UTC on September 26, 2011. The vertical partitioning of the domain is denoted with the black curves, which correspond to pressure levels −2.5, 0, 2.5, and 5. **(C)** Analogous plots of the electron density bias vertical structure at a fixed 77.5° S latitude.

Also shown in Figure 2A are four geographical regions over which we examine the model bias. The boundaries for the Northern and Southern high-latitude regions (*R*_{1} and *R*_{2}) are located at the *R*_{3}) from the nighttime (*R*_{4}).

A positive (negative) sign indicates that the background electron density overestimates (underestimates) the true electron density. Prior to the main phase of the geomagnetic storm, the most notable bias structure is in the daytime mid-latitude regions (*R*_{3}), where there is a negative bias. As geomagnetic conditions become increasingly disturbed, the misspecification in the magnetospheric input grows, resulting in the formation of a significant negative bias structure over the Southern Polar region (*R*_{1}). The negative bias structure continues to grow over regions *R*_{1} and *R*_{3} over the next few hours throughout the rest of the main storm phase. There is also a formation of negative bias in the Northern polar region (*R*_{2}) and a positive bias forming over the night-time low-to-mid latitude region (*R*_{4}). After the main phase (21:30 UTC), the negative bias in the polar regions is largely diminished but the bias over the daytime mid-latitudes remains relatively constant. This suggests that the misspecification of

Figure 2B shows the vertical structure of the electron density bias at the

### 5.1 Region-Averaged Bias Estimates

Following the procedure described in Section 3, the ensemble of global bias correction vectors at a given time *m* is the number of grid points and *d* is the number of state variables being analyzed. As discussed at the end of Section 3, we consider bias corrections only for the electron density field. For simplicity in notation, we now regard each

The forecast/analysis cycle is initialized at time

which has ensemble mean and reflects the initial spatial correlations of the electron density field. The factor of 1/4 in Eq. 14 is incorporated to reduce the variance in the electron density bias estimate, which is generally expected to be smaller than that of the electron density field. The procedure described in Section 3 is then carried out for subsequent analysis steps to produce spatially and time varying bias estimates throughout the domain that are evolved according to *L*, variance inflation is applied to the bias correction component so that its ensemble variance is at least 20% of the electron density component at *L*.

The model state vector is analyzed without the bias corrections for the first few analysis cycles after the forecast spinup so that errors associated with the unadjusted thermospheric state are reduced and the electron density bias is primarily dependent on the misspecification of the solar and magnetospheric inputs. We compute eight analysis cycles, from 16:30 to 23:30 UTC on September 25, 2011, before applying the bias correction strategy for the first time at 00:30 UTC on September 26, 2011.

To evaluate the skill of our bias estimation strategy, we examine the model bias estimates, averaged over the four geographical regions shown in Figure 2A. Consider one of the fixed bias regions *j*th forecast, *R*. Similarly, form the region-averaged electron density of the control simulation, *j*th forecast is given by

which yields an ensemble of regionally averaged model bias estimates, whose ensemble mean,

Figure 3A summarizes the temporal evolution of the region-averaged electron density bias in the Southern polar region (*R*_{1}) averaged between 1.0 and 5.5 pressure levels, which correspond to altitudes of about 260–485 km. The horizontal axis is time in hours, starting at 00:30 UTC on September 26 2011 and ending at 23:30 UTC on September 27, 2011, and the vertical axis is electron density in units of

**FIGURE 3**. **(A)** Time series of the true electron density (green), in units of **(B–D)** Analogous time series of electron density, averaged over region

The region-averaged electron density for the control simulation varies considerably within region *R*_{1}, particularly during the period of main geomagnetic disturbance, where the electron density increases sharply between 12:30 UTC and 16:30 UTC on 26 September. After the main storm phase, electron density content drops considerably and remains relatively constant throughout the rest of the simulation. The region-averaged electron density trajectories for the forecast evolve similarly. Prior to the onset of geomagnetic disturbances, the electron density bias is relatively small, but rises sharply during the main phase, peaking at around 16:30 UTC. The electron density bias is maintained until about 19:30 UTC, where it begins to drop considerably and remains relatively constant after 23:30 UTC on 26 September, until there is a small resurgence toward the end of the simulation.

Prior to the onset of geomagnetic disturbances (12:30 UTC), the bias strategy correctly detects the near-zero positive bias and continually adjusts its estimates as the model bias starts to increase during the transition into the main phase of the geomagnetic storm around 12:30 UTC. Although the temporal variations in the bias estimates evolve in a similar manner as the model bias, the temporal variations are not captured during this transitionary period in which the model bias increases and decreases sharply during the initial and relaxation phases, respectively. Most notably, the bias corrections adjust too slowly and underestimate the model bias peak during the initial storm phase. After the main phase, the model bias decreases sharply and although the bias corrections adjust, they do so too slowly and overestimate the model considerably after 20:30 UTC for the next 8 h.

The region-averaged bias over the daytime low-to-mid-latitude region (*R*_{3}), shown in Figure 3C, exhibits a relatively larger bias compared to the high-latitude regions but remains relatively constant throughout the simulation, since the misspecification of *R*_{4}), shown in Figure 3D, is considerably smaller than the bias in its surrounding regions. Its temporal evolution is relatively constant throughout the simulation and is well represented with the bias corrections.

### 5.2 Validation of Bias Correction

To validate the bias correction strategy, we compare 1-h predictions of electron density before and after the bias corrections are applied. Since the bias corrections are applied only during the evaluation of the forward operator, we evaluate the benefit of the bias corrections at the observed locations. In particular, we compare the prediction RMS error averaged over all vertical profiles within each of the geographical regions:

where *i*th observation location. The RMSE in Eq. 16 is averaged over the set of *R*.

Figure 4A shows the RMSE time series of 1 h predictions over region

**FIGURE 4**. **(A)** Time series of the root mean squared error (RMSE) of 1 h forecasted electron density predictions, in **(B–D)** Analogous time series of region-averaged RMSE values and ratios for regions

Figure 4B shows the analogous RMSE time series for region

Figure 4C shows the analogous RMSE time series for the low-to-mid-latitude daytime region (*R*_{3}). The RMSE time series for this region has a less pronounced peak during the main phase and the improvement due to the applied bias corrections is consistent throughout the simulation, although there are a few time periods where the bias corrections offer little to no improvement near the end of the simulation. The RMSE time series for the night-time region (*R*_{4}) is shown in Figure 4D. In this region, the bias is also relatively constant and the benefits of applying the bias corrections are consistent.

The numerical experiments presented in this section show that the bias correction strategy usually improves 1-h forecasts throughout the domain, particularly during the main phase of the geomagnetic storm. However, some limitations of the strategy are observed during transitionary periods of the geomagnetic storm, where model bias undergoes relatively fast temporal variations. Most notably, the high-latitude bias corrections underestimate the sharp increase in model bias during the initial storm phase and overestimate the model bias as it sharply decreases during the relaxation storm phase. The misrepresentation of model bias during these time periods is primarily due to the choice of bias evolution operator,

## 6 Numerical Experiment 2: Bias Estimation With a Time-Varying Evolution Operator

Although a simple persistence model for the bias evolution provides some improvement in 1-h forecasts of electron density, the results show that there is considerable room for improvement. The question of how to optimize the bias evolution model is a topic for future investigation, but here we outline one potential approach to the problem.

Figure 5A shows a time series (blue dots) of the RMSE between the “true” values of the electron density and those of the ensemble mean, in units of *R*_{1}) grows exponentially, and from about 19 UTC to 23 UTC decays exponentially.

**FIGURE 5**. **(A)** RMSE time series of electron density bias in units of _{∼}375 km), averaged over all grid points in region **(B)** Vertical structure of each **(left)** and relaxation **(right)** phases of the geomagnetic storm. The vertical axis is in pressure levels and the horizontal axis is the value of

By their nature, model errors tend to be of a particular sign. For example, the representative bias results in Figure 2 show regions where the model consistently underestimates or overestimates the electron density. Consequently, we propose a growth-relaxation model for the bias evolution of the form

to describe the temporal evolution of model bias estimates at each analysis time *j*th forecast is given by

### 6.1 Bias Propagation

Because the errors in the model driving parameters under geomagnetic storm conditions are not well characterized, suitable choices of the growth and decay rates for the evolution of the model bias must be determined empirically. The choices depend upon latitude, altitude, and time of year (i.e., Earth’s axial tilt). Detailed reanalyses of prior geomagnetic storms are needed, which we have not yet undertaken. For the sake of illustration, however, we consider an idealized procedure for finding time-dependent bias parameters and discuss the effect on the resulting forecast errors. Section 6.3 describes some potential future work on this topic.

We consider the RMSE of electron density, averaged over all grid points in the Southern polar region (*R*_{1}) at pressure level 2.0 as mentioned previously in Figure 5A (blue dots). Although the electron density bias a strong spatiotemporal variability (Figure 2), we perform a simple regression from 12 UTC to 15 UTC on 26 September of the RMSE (red curves in Figure 5) *vs.* time to determine the spatiotemporally averaged rate of bias increase, which we denote as *n* indexes the analysis times from 12 UTC to 15 UTC. We proceed likewise to compute *z* above 0.5, which corresponds to the lower portion of the F-layer. Below this pressure level, we set

Figure 5B shows the vertical structure of *R*_{2} and *R*_{3} and are provided in Appendix. The growth-relaxation bias evolution model is not used in the night-time region (*R*_{4}), because the bias is not observed to have a pronounced growth or decay component like the other geographical regions.

The method just described is but one of many possible parameterizations of the bias evolution model. A more sophisticated approach might further localize the

### 6.2 Results

The time series of region-averaged bias estimates, in the case where the growth-relaxation model described in Section 6.1 is used for *R*_{1} and *R*_{2} in Figures 6A,B. The quantities in these figures are computed as described in Section 5.1 and are analogous to Figures 3A,B. Regions *R*_{1} and *R*_{2} are the high-latitude regions over which rapid temporal variations occur during the transitionary periods of the geomagnetic storm. Comparison with Figures 3A,B reveals that the growth-relaxation model for *R*_{1} during the beginning of the initial phase between 11:30 and 13:30 UTC, but the application of *R*_{2} is similarly well represented throughout the simulation with the growth-relaxation model for

**FIGURE 6**. **(A,B)** Time series of bias correction parameters averaged over regions **(C,D)** Time series of the root mean squared error (RMSE) of 1 h forecasted electron density predictions, in

Figures 6C,D show the corresponding RMSE time series of 1-h electron density predictions over regions *R*_{1} and *R*_{2} in the case where the growth-relaxation model is used for *R*_{1} and *R*_{2}. The most notable benefits of the bias corrections occur during the relaxation phase, where there is improvement of about 40 and 20% in regions *R*_{1} and *R*_{2}, respectively, compared to the persistent bias evolution model.

Figure 7 summarizes the spatial distribution of the bias corrections at same representative times of the geomagnetic storm. For comparison, Figure 7A shows the same background electron density bias from Figure 2A, computed as

**FIGURE 7**. **(A)** Same global maps of electron density bias as shown in Figure 2A. **(B)** Spatial structure of the bias correction estimates, **(C)** Analogous global maps of electron density bias after the field of correction estimates is added.

Although there is considerable improvement in the overall bias distribution, there are a few locations where the bias correction scheme does not accurately represent the model bias, mainly in regions of sparse data coverage. Overall, the background prediction of electron density is improved at the observed locations, so that the data assimilation procedure may produce a more gentle update during the assimilation of the observations.

Figure 8A shows 1-h predictions of electron density averaged over all vertical profiles in region *R*_{1} at 12:30, 15:30, 18:30 and 21:30 UTC on September 26, 2011. The horizontal axis is electron density in *R*_{2} and *R*_{3} are shown in Figures 8B,C, respectively. Similar improvements in peak density and altitude are observed over these regions. The bias corrections for the most part correctly increase the electron density peak and its altitudes, but there are times where the electron density peak is correctly for improved agreement with the truth.

**FIGURE 8**. **(A)** Comparison of 1 h electron density vertical profile predictions at the indicated times on September 26, 2011. The green profile in each figure is the average of all COSMIC electron density vertical profiles located in region **(B,C)** Analogous comparison of electron density predictions over regions

### 6.3 Discussion

The main point of this exercise is to show that it may be possible, using relatively few parameters, to devise a characterization of model error that can improve short-term forecasts. Further refinement of the bias estimation procedure, particularly in the time interval following the storm peak, should be possible.

The parameterization procedure described above uses the RMSE between the predicted and “true” values of the electron density (assuming a persistent-bias model) to determine the growth and decay factors

Various authors have sought to devise useful statistical characterizations of geomagnetic storms. No comprehensive review of these efforts is attempted here, but as an illustration, we mention the “epoch analysis” of [59], who analyzed disturbance time index (Dst) data from 305 geomagnetic storms. They define 1) the start time of each storm (i.e., the commencement of the main phase) as the time at which Dst first crosses zero; 2) the peak of the storm as the time at which Dst achieves a minimum value,

The left half of the regression interval in Figure 5 corresponds approximately from the time when the Dst index first increases from its background value to when Dst first becomes negative (i.e., the start of the main storm phase). The right half corresponds approximately to the main phase of the storm, ending at the peak, when Dst reaches its relative minimum. One potential indicator of an impending geomagnetic storm is when the 1-h forecast errors begin to increase substantially; perhaps it will be possible to define a threshold value to initiate bias correction. Our initial numerical results suggest that once of the peak of the storm is reached, the forecast ensemble tracks the electron density reasonably well during the recovery phase without bias correction, although improvements may be possible.

By the criteria of [59], the September 26–27, 2011 storm would be classified as intense, because

The Dst is but one possible index to consider for this purpose. Other proposed ionospheric storm scales, such as one based on TEC and

## 7 Conclusion

A general strategy for model bias correction is presented within a data assimilation system for complex and sparsely observed dynamical systems. The proposed strategy is applied in observing system experiments corresponding to extreme conditions, in which the geomagnetic storm of September 26–27, 2011 with a global ionospheric circulation model, the TIEGEM. Synthetic electron density vertical profiles, whose locations are given by the operational COSMIC satellite observing network during the time period of the storm, are assimilated using the LETKF. Systematic model bias in electron density predictions is simulated through the misspecification of parameterized solar and magnetospheric inputs to the TIEGCM, which specify key dynamics in the ionosphere. The spatial distribution of electron density bias is estimated through a state-augmentation approach. The bias estimates are applied in the evaluation of the forward operator, with the intent to reduce the effect of model errors in the predicted electron densities prior to the assimilation of observations. This methodology permits spatially varying estimation of model bias, which may be useful during extreme events, to account for storm-time effects that are not well represented with the parameterized representation of the solar and magnetospheric drivers alone.

The results of this study suggest that relatively simple models for bias correction can provide significant improvements in 1-h electron density forecasts when the model error arises from a systematic misspecification of the ionospheric driving parameters. Two such models are considered: one in which the bias is assumed constant, and another in which an exponential growth and relaxation process is assumed in the electron density.

The constant-bias assumption underestimates the growth in the model error during the initial storm phase and overestimates it during the relaxation storm phase, particularly for the electron density over high-latitude regions. This persistence model does not adjust the bias estimates quickly enough to capture their temporal variability. Additionally, the observing network is relatively sparse over high-latitude regions during the time period under study, which further slows the adjustment of the bias estimates. Nevertheless, the bias corrections yield considerable improvements in 1-h forecasts of electron density in certain regions prior to onset of the geomagnetic storm and during the initial and main phases; however, the 1-h forecasts are worse with the bias correction during the relaxation phase.

The benefits of the second bias-correction approach, which uses a growth-relaxation model for the electron density, are seen in 1-h electron density predictions over high-latitude regions during transitionary periods of the geomagnetic storm, when the temporal variations in model bias are most prominent. Model bias corrections over the Southern polar region (*R*_{1}) yield on average about 45 and 35% improvements in 1 h electron density predictions during the main and relaxation phases of the storm, respectively. The Northern polar region (*R*_{2}) displays a similar peak in model bias evolution, and bias corrections yield an improvement of about 45 and 20% during the main and relaxation phases, respectively. Over daytime low-latitude regions, a steady improvement of about 25% is observed throughout the simulation. One-hour predictions of the global distribution of model bias are approximated reasonably well with the bias correction strategy overall. The bias corrections also yield adequate estimates of the vertical structure of the model bias and considerably improves 1-h electron density predictions of the observed vertical profiles, including the adjustment of the F-layer peak and its altitude.

The experiments presented in this paper provide a proof of concept. The results are established with synthetically generated observations, which allows us to quantify specifically the effect of the bias correction on 1-h forecasts of electron density. The choice of bias evolution operator is important when there are significant temporal variations in the model’s dynamical variables and when the observing network is sparse. Extensions to real observations will be considered in future work.

The approach described in this paper may be an effective tool for model bias estimation during geomagnetic storms, provided that there are qualitative similarities in the temporal variations of model bias during these time periods. For example, the bias model employed in our numerical experiments may be applied during other geomagnetic storms in which positive and negative storm-time effects are expected during the initial and relaxation phases, respectively. Further studies of previous geomagnetic storms may help provide improved estimates of typical temporal variations in the electron bias evolution, which may depend on solar/geomagnetic activity, latitude, pressure level, and/or season.

The extension of this approach to real observations presents some open problems. One obvious question involves model errors arising from subgrid parameterizations and uncertainties in the solar and geomagnetic driving parameters. Representativeness error in the observations must also be considered. The results presented here suggest that relatively simple dynamical characterizations of the model bias can yield a stable filter and improve 1-h forecasts of electron density. Nevertheless, better choices of the bias evolution operator

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

JD, EK, and AM contributed to conception and design of the computational study. JD performed the data analysis. JD, EK, and AM wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

The TIEGCM was obtained from the High Altitude Observatory, NCARData files containing solar and geomagnetic parameters used in this paper are published by NOAA and are packaged in the TIEGCM code. The COSMIC data information was obtained from the COSMIC Data Analysis and Archive Center (CDAAC) at UCAR (http://cdaac-www.cosmic.ucar.edu). This material is based upon work supported by the Air Force Office of Scientific Research under the award number FA9550-19-1–0064.

## Appendix

### Bias Estimation Algorithm Pseudocode

The following algorithm demonstrates the procedure used to compute the analysis at an arbitrary grid point *L*, and how the model bias corrections are applied and updated. We assume there are *d* components in the local state vector, *k*. The quantities below are related to the grid point *L* and its associated local region, where we assume there are *L*.

#### Analysis Computation

A detailed description of the analysis calculation and the associated computational complexity is provided in [25]. In the following suppose the forecast and bias estimates have been propagated to a time *L*.

a. Construct the

b. Construct the

c. Form the

d. Compute the

e. The analysis ensemble is constructed in terms of the “weight matrix”,

1) Compute the

2)Define the *k*-vector

f. Compute the analysis perturbation matrix

g. The *j*th column of the analysis ensemble is formed by adding the background ensemble mean, *j*th column of

#### Propagate Ionospheric and Model Bias State Estimates

a. Apply forecast model

b. Compute

c. Apply bias evolution operator

## References

1.N. R. C. S. S Board. Space Weather and Space Climatology: A Vision for Future Capabilities. In: *Solar and Space Physics: A Science for a Technological Society*. Washington, DC: National Academies Press, (2013). ch. 7, pp. 135–46.

2. Kelley, MC. *The Earth’s Ionosphere: Plasma Physics and Electrodynamics, Vol. 96 of International Geophysics*. 2nd ed. Cambridge, United States: Academic Press(2009).

3. Materassi, M, Forte, B, Coster, AJ, and Skone, S, eds. *The Dynamical Ionosphere: A Systems Approach to Ionospheric Irregularity*. Amsterdam, Netherlands: Elsevier (2020).

4. Schunk, R, and Nagy, A. *Ionospheres: Physics, Plasma Physics, and Chemistry, Cambridge Atmospheric and Space Science Series*. 2nd ed. Cambridge, UK: Cambridge University Press (2018).

5. Bilitza, D, and Reinisch, B. International Reference Ionosphere 2007: Improvements and New Parameters. *Adv Space Res* (2007). 42:599–609. doi:10.1016/j.asr.2007.07.048

6.N. C. F. A. R. High-Altitude Observatory. Available at: https://www.hao.ucar.edu/modeling/tgcm/ (Accessed May 15, 2019).

7. Durazo, JA, Kostelich, EJ, and Mahalov, A. Local Ensemble Transform Kalman Filter for Ionospheric Data Assimilation: Observation Influence Analysis during a Geomagnetic Storm Event. *J Geophys Res Space Phys* (2017). 122:9652–69. doi:10.1002/2017JA024274

8. Abarzhi, S. Coherent Structures and Pattern Formation in Rayleigh-Taylor Turbulent Mixing. *Phys Scr* (2008). 78:015401. doi:10.1088/0031-8949/78/01/015401

9. Abarzhi, SI. Review of Theoretical Modelling Approaches of Rayleigh-Taylor Instabilities and Turbulent Mixing. *Philos Trans R Soc A* (2010). 368:1809–28. doi:10.1098/rsta.2010.0020

10. Abarzhi, S, and Rosner, R. A Comparative Study of Approaches for Modeling Rayleigh-Taylor Turbulent Mixing. *Phys Scr* (2010). T142:014012. doi:10.1088/0031-8949/2010/t142/014012

11. Mahalov, A, and Moustaoui, M. Multiscale Nested Simulations of Rayleigh-Taylor Instabilities in Ionospheric Flows. *J Fluid Eng* (2013). 136:021202. doi:10.1115/1.4025657

12. Mahalov, A. Multiscale Modeling and Nested Simulations of Three-Dimensional Ionospheric Plasmas: Rayleigh-Taylor Turbulence and Nonequilibrium Layer Dynamics at Fine Scales. *Phys Scr* (2014). 89:098001. doi:10.1088/0031-8949/89/9/098001

13. Cheng, B, and Mahalov, A. General Results on Zonation in Rotating Systems with a β-Effect and the Electromagnetic Force. In: B Galperin, and P Read, editors. *Zonal Jets: Phenomenology, Genesis, and Physics*. Cambridge, UK: Cambridge University Press (2019). p. 209–19. doi:10.1017/9781107358225.012

14. Galperin, B, and Read, P eds. *Zonal Jets: Phenomenology, Genesis, and Physics*. Cambridge, United Kingdom: Cambridge University Press (2019).

15. Tang, W, and Mahalov, A. Stochastic Lagrangian Dynamics for Charged Flows in the E-F Regions of Ionosphere. *Phys Plasmas* (2013). 20:032305. doi:10.1063/1.4794735

16. Tang, W, and Mahalov, A. The Response of Plasma Density to Breaking Inertial Gravity Waves in the Lower Regions of Ionosphere. *Phys Plasmas* (2014). 21:042901. doi:10.1063/1.4870760

17. Zhou, Y. Rayleigh-Taylor and Richtmyer-Meshkov Instability Induced Flow, Turbulence, and Mixing. I. *Phys Rep* (2017). 720–722:7201–136. doi:10.1016/j.physrep.2017.07.005 |

18. Zhou, Y. Rayleigh-Taylor and Richtmyer-Meshkov Instability Induced Flow, Turbulence, and Mixing. II. *Phys Rep* (2017). 723–725:1–160. doi:10.1016/j.physrep.2017.07.008 |

19. Merkin, VG, Kondrashov, D, Ghil, M, and Anderson, BJ. Data Assimilation of Low-Altitude Magnetic Perturbations into a Global Magnetosphere Model. *Space Weather* (2016). 14:165–84. doi:10.1002/2015SW001330

20. Eyre, JR. Observation Bias Correction Schemes in Data Assimilation Systems: A Theoretical Study of Some of Their Properties. *Q J R Meteorol Soc* (2016). 142:2284–91. doi:10.1002/qj.2819

21. Gilman, M, Smith, E, and Tsynkov, S. *Transionospheric Synthetic Aperture Imaging*. Cham: Birkäuser (2017). doi:10.1007/978-3-319-52127-5

22. Godinez, HC, Lawrence, E, Higdon, D, Ridley, A, Koller, J, and Klimenko, A. Specification of the Ionosphere-Thermosphere Using the Ensemble Kalman Filter. In: S Ravela, and A Sandu, editors. *Dynamic Data-Driven Environmental Systems Science, DyDESS*. Cham: Birkäuser (2014). p. 274–83.

23. Heelis, RA, Lowell, JK, and Spiro, RW. A Model of the High-Latitude Ionospheric Convection Pattern. *J Geophys Res* (1982). 87:6339–45. doi:10.1029/JA087iA08p06339

24. Hsu, CT, Matsuo, T, Wang, W, and Liu, JY. Effects of Inferring Unobserved Thermospheric and Ionospheric State Variables by Using an Ensemble Kalman Filter on Global Ionospheric Specification and Forecasting. *J Geophys Res Space Phys* (2014). 119:9256–67. doi:10.1002/2014ja020390

25. Hunt, BR, Kostelich, EJ, and Szunyogh, I. Efficient Data Assimilation for Spatiotemporal Chaos: A Local Ensemble Transform Kalman Filter. *Phys D Nonlinear Phenomena* (2007). 230:112–26. doi:10.1016/j.physd.2006.11.008

26. Jin, J, Lin, HX, Segers, A, Xie, Y, and Heemink, A. Machine Learning for Observation Bias Correction with Application to Dust Storm Data Assimilation. *Atmos Chem Phys* (2019). 19:10009–26. doi:10.5194/acp-2019-298

27. Koller, J, Brennan, S, Godinez, H, Higdon, D, Koller, J, and Lawrence, E. Integrated Modeling of Perturbations in Atmospheres for Conjunction Tracking. In: AMOS Conference: LA-UR-13-26931 (2013).

28. Lee, IT, Matsuo, T, Liu, JY, Wang, W, and Lin, CH. Assimilation of FORMOSAT-3/COSMIC Electron Density Profiles into a Coupled Thermosphere/ionosphere Model Using Ensemble Kalman Filtering, *J Geophys Res* (2012). 117:A10318. doi:10.1029/2012ja017700

29. Li, W, Huang, L, Zhang, S, and Chai, Y. Assessing Global Ionosphere TEC Maps with Satellite Altimetry and Ionospheric Radio Occultation Observations. *Sensors* (2019). 19:5489. doi:10.3390/s19245489 |

30. Kuo, JY, Lin, CY, Lin, CH, Tsai, HF, Solomon, SC, Sun, YY, et al. Artificial Plasma Cave in the Low-Latitude Ionosphere Results from the Radio Occultation Inversion of the FORMOSAT-3/COSMIC. *J Geophys Res* (2010). 115. doi:10.1029/2009JA015079

31. Matsuo, T, and Araujo-Pradere, EA. Role of Thermosphere-Ionosphere Coupling in a Global Ionospheric Specification. *Radio Sci* (2011). 46:RS0D23. doi:10.1029/2010rs004576

32. Matsuo, T, Lee, IT, and Anderson, J. *Upper Atmosphere Data Assimilation with an Ensemble Kalman Filter*. Washington, DC: AGU (2013). p. RS0D23.

33. Dee, DP. Bias and Data Assimilation. *Q J R Meteorol Soc* (2005). 131:3323–43. doi:10.1256/qj.05.137

34. Dee, DP, and da Silva, AM. Data Assimilation in the Presence of Forecast Bias. *Q J R Meteorol Soc* (1998). 124:269–95. doi:10.1002/qj.49712454512

35. Aravéquia, JA, Szunyogh, I, Fertig, EJ, Kalnay, E, Kuhl, D, and Kostelich, EJ. Evaluation of a Strategy for the Assimilation of Satellite Radiance Observations with the Local Ensemble Transform Kalman Filter. *Mon Weather Rev* (2011). 139:1932–51. doi:10.1175/2010MWR3515.1

36. Canter, M, Barth, A, and Beckers, J-M. Correcting Circulation Biases in a Lower-Resolution Global General Circulation Model with Data Assimilation. *Ocean Dyn* (2017). 67:281–98. doi:10.1007/s10236-016-1022-3

37. Baek, SJ, Hunt, BR, Kalnay, E, Ott, E, and Szunyogh, I. Local Ensemble Kalman Filtering in the Presence of Model Bias. *Tellus A Dyn Meteorol Oceanogr* (2006). 58:293–306. doi:10.1111/j.1600-0870.2006.00178.x

38. Codrescu, SM, Codrescu, MV, and Fedrizzi, M. An Ensemble Kalman Filter for the Thermosphere-Ionosphere. *Space Weather* (2018). 16:57–68. doi:10.1002/2017SW001752

39. Tsurutani, BT, Verkhoglyadova, OP, Mannucci, AJ, Lakhina, GS, Li, G, and Zank, GP. A Brief Review of Solar Flare Effects on the Ionosphere. *Radio Sci* (2009). 44:RS0A17. doi:10.1029/2008rs004029

40. Buensanto, M. Ionospheric Storms - A Review. *Space Sci Rev* (1999). 88:563–601. doi:10.1029/2010RG000340

41. Dubey, S, Wahi, R, and Gwal, AK. Ionospheric Effects on Gps Positioning. *Adv Space Res* (2006). 38:2478–84. doi:10.1016/j.asr.2005.07.030

42. Stauning, P. Power Grid Disturbances and Polar Cap Index during Geomagnetic Storms. *J Space Weather Space Clim* (2013). 3:A22. doi:10.1051/swsc/2013044

43. Weimer, DR. Improved Ionospheric Electrodynamic Models and Application to Calculating Joule Heating Rates. *J Geophys Res* (2005). 110:A05306. doi:10.1029/2004JA010884

44. Tan, Y, Hu, Q, Wang, Z, and Zhong, Q. Geomagnetic Index *kp* Forecasting with LSTM. *Space Weather* (2018). 16:406–16. doi:10.1002/2017SW001764

45. Wing, S, Johnson, JR, Jen, J, Meng, C-I, Sibeck, DG, Bechtold, K, et al. Kp Forecast Models. *J Geophys Res Space Phys* (2005). 110:A04203. doi:10.1029/2004ja010500

46. Xu, W. Uncertainty in Magnetic Activity Indices. *Sci China Ser E Technol Sci* (2008). 51:1659–64. doi:10.1007/s11431-008-0261-z

47. Pi, X, Mannucci, AJ, Iijima, BA, Wilson, BD, Komjathy, A, Runge, TF, et al. Assimilative Modeling of Ionospheric Disturbances with FORMOSAT-3/COSMIC and Ground-Based GPS Measurements. *Terr Atmos Ocean Sci* (2009). 20:273–85. doi:10.3319/TAO.2008.01.04.01(F3C)

48. Chartier, AT, Jackson, DR, and Mitchell, CN. A Comparison of the Effects of Initializing Different Thermosphere-Ionosphere Model Fields on Storm Time Plasma Density Forecasts. *J Geophys Res Space Phys* (2013). 118:7329–37. doi:10.1002/2013ja019034

49. Ruiz, JJ, Pulido, M, and Miyoshi, T. Estimating Model Parameters with Ensemble-Based Data Assimilation: A Review. *J Meteorol Soc Jpn* (2013). 91:79–99. doi:10.2151/jmsj.2013-201

50. Bellsky, T, Berwald, J, and Mitchell, L. Nonglobal Parameter Estimation Using Local Ensemble Kalman Filtering. *Monthly weather Rev* (2014). 142:2150–64. doi:10.1175/mwr-d-13-00200.1

51. Pedatella, NM, Lei, J, Larson, KM, and Forbes, JM. Observations of the Ionospheric Response to the 15 December 2006 Geomagnetic Storm: Long-Duration Positive Storm Effect. *J Geophys Res Space Phys* (2009). 114. doi:10.1029/2009JA014568

52 Chartier, AT, Matsuo, T, Anderson, JL, Collins, N, Hoar, TJ, Lu, G, et al. Ionospheric Data Assimilation and Forecasting during Storms. *J Geophys Res Space Phys* (2016). 121:764–78. doi:10.1002/2014JA020799

53. Rocken, C, Kuo, YH, Schreiner, WS, Sokolovskiy, S, and McCormick, C. COSMIC System Description, *Terr Atmos Ocean Sci* (2000). 11:21–52.doi:10.3319/tao.2000.11.1.21(cosmic)

54. Durazo, J, Kostelich, E, Mahalov, A, and Tang, W. Observing System Experiments with an Ionospheric Electrodynamics Model. *Phys Scr* (2016). 91:044001. doi:10.1088/0031-8949/91/4/044001

55. Miyoshi, T. The Gaussian Approach to Adaptive Covariance Inflation and its Implementation with the Local Ensemble Transform Kalman Filter. *Mon Weather Rev* (2011). 139:1519–35. doi:10.1175/2010mwr3570.1

56. Richards, PG, Fennelly, JA, and Torr, DG. EUVAC: A Solar EUV Flux Model for Aeronomic Calculations. *J Geophys Res* (1994). 99:8981–92. doi:10.1029/94JA00518

57. Szunyogh, I, Kostelich, EJ, Gyarmati, G, Kalnay, E, Hunt, BR, Ott, E, et al. A Local Ensemble Transform Kalman Filter Data Assimilation System for the NCEP Global Model, *Tellus A Dyn Meteorol Oceanogr* (2008). 60:113–30. doi:10.1111/j.1600-0870.2007.00274.x

58. Hunt, X, Schreiner, WS, Lei, J, Sokolovskiy, SV, Rocken, C, Hunt, DC, et al. Error Analysis of Abel Retrieved Electron Density Profiles from Radio Occultation Measurements. *Ann Geophys* (2010). 28:217–22. doi:10.5194/angeo-28-217-2010

59. Yokoyama, N, and Kamide, Y. Statistical Nature of Geomagnetic Storms. *J Geophys Res* (1997). 102:14215–22. doi:10.1029/97ja00903

Keywords: non-equilibrium dynamics, data assimilation, ensemble Kalman filter, ionospheric model, model bias correction, geomagnetic storm

Citation: Durazo J, Kostelich EJ and Mahalov A (2021) Data Assimilation for Ionospheric Space-Weather Forecasting in the Presence of Model Bias. *Front. Appl. Math. Stat.* 7:679477. doi: 10.3389/fams.2021.679477

Received: 11 March 2021; Accepted: 16 April 2021;

Published: 07 May 2021.

Edited by:

Alexander Nepomnyashchy, Technion Israel Institute of Technology, IsraelCopyright © 2021 Durazo, Kostelich and Mahalov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Eric J. Kostelich, kostelich@asu.edu