Geomagnetic secular variation forecast using the NASA GEMS ensemble Kalman filter: A candidate SV model for IGRF-13

Abstract We have produced a 5-year mean secular variation (SV) of the geomagnetic field for the period 2020–2025. We use the NASA Geomagnetic Ensemble Modeling System (GEMS), which consists of the NASA Goddard geodynamo model and ensemble Kalman filter (EnKF) with 400 ensemble members. Geomagnetic field models are used as observations for the assimilation, including gufm1 (1590–1960), CM4 (1961–2000) and CM6 (2001–2019). The forecast involves a bias correction scheme that assumes that the model bias changes on timescales much longer than the forecast period, so that they can be removed by successive forecast series. The algorithm was validated on the time period 2010-2015 by comparing with CM6 before being applied to the 2020–2025 time period. This forecast has been submitted as a candidate predictive model of IGRF-13 for the period 2020–2025. Graphical abstract


Introduction
It has long been observed that the Earth's main magnetic field changes slowly in time (Bullard et al. 1950;Hide 1967;Bloxham and Gubbins 1985;Jackson et al. 2000;Lund 2018). This change, called the secular variation (SV), originates dominantly from the Earth's iron-rich fluid outer core that is vigorously convecting, driven by thermo-chemical buoyancy released from Earth's secular cooling and differentiation (geodynamo) (e.g. Lund and Olson 1987;Braginsky and Roberts, 1996). Over the past 25 years, great progress has been made in numerically modeling the geodynamo (Glatzmaier and Roberts 1995;Kagyama and Sato, 1997;Kuang and Bloxham 1997;Christensen et al. 2001;Jones et al. 2011;Matsui et al. 2016). These models can reproduce much of the qualitative aspects of the Earth's magnetic field (Christensen et al. 2010), including dipole dominance, westward drift (Aubert et al. 2013) and occasional reversals (Olson 2007). Although it remains computationally prohibitive to simulate the geodynamo in the parameter regimes appropriate for the Earth's core, the asymptotic properties emerged from numerical solutions with wide range of dynamo parameter values open the door for quantitative applications of numerical models to the geomagnetic field (e.g. Christensen and Aubert 2006;Christensen et al. 2010;Aubert and Fournier 2017;Kuang et al. 2017).
One such application is to use a numerical model to predict geomagnetic SV on time scales of several years and longer, if proper initialization is made using data assimilation techniques, as is done in numerical weather prediction (NWP). Forecast results from this geomagnetic data assimilation approach have already been used as predictive field candidates for previous IGRF models (Kuang et al. 2010;Fournier et al. 2015).
Obviously, this approach is different from other SV forecasts made without utilization of dynamo models. For example, Olsen, et al. (2010) modeled the time dependence of the Gauss coefficients with order 6 B-splines so that the predictive SV for 2010-2015 can be extrapolated from the coefficients derived from observations prior to 2010. While utilization of geodynamo models can take advantage of adding dynamically consistent motional induction in the outer core to forecast the SV, it also carries additional complications, such as model bias that results from the mismatch between the parameter values used in the models and those appropriate for the Earth's core, and propagation of observational errors in both time and space. These difficulties have been analyzed through assimilation experiments with simplified systems (Sun et al. 2007;Fournier et al 2007) and observation system simulation experiments (OSSEs) (Liu et al 2007;Aubert and Fournier, 2011;. It is expected that improvements to both modeling and assimilation techniques will make geomagnetic forecast a growing contributor of the IGRF candidate models. In this paper, we report our mean SV forecast for the period from 2020 to 2025. This effort differs from our earlier work for IGRF candidate predictive field model (Kuang et al. 2010) in two major areas: we assimilate Gauss coefficients from the CM6 geomagnetic field model (Sabaka et al in this special issue) for the period from 2000 to 2019; and the error covariance in our assimilation is determined through a large ensemble of assimilation solutions, and is updated in each analysis cycle. The two improvements can help produce more accurate SV forecasts and uncertainty estimates, since, for example, the utilization of a single field model can avoid any potential disagreement in different field models that may lead to inconsistent observation error estimates, and because the time varying error covariance can provide more consistent analysis for each cycle.
This paper is organized as follows: the algorithm of our forecasts is given in Section 2; testing and validation are presented in Section 3, followed by the forecasts and discussion.

Geodynamo model and data assimilation algorithm
The geomagnetic data assimilation system used for our SV forecasts is the NASA Geomagnetic Ensemble Modeling System (GEMS). It consists of a numerical geodynamo model (Kuang and Bloxham 1999;Kuang and Chao 2003), an ensemble Kalman filter (Evensen 2009;Sun and Kuang 2015). The model solves a set of nonlinear magnetohydrodynamic equations in a spherical shell domain. Details of the model are given in Kuang and Bloxham (1999), Kuang and Chao (2003) and Jiang and Kuang (2008). The core state in this model is defined with the fluid velocity field v , the magnetic field B and the relative density anomaly �ρ that arises from temperature variations. These state variables are described by spherical harmonic expansions in the co-latitude θ and the longitude φ , with the spherical harmonic coefficients given on discrete radial grid points. If we denote by x the vector of these coefficients, then the dynamo system can be symbolically represented as where M includes all linear and nonlinear model operators. GEMS employs a sequential ensemble Kalman filter (EnKF), described below, which is the result of many years of development that started with a simple one dimensional system (Sun et al. 2007) and an optimal interpolation (OI) system (Kuang et al. 2008(Kuang et al. , 2009. These systems have shown how geomagnetic data assimilation can be used to estimate optimal geodynamo model parameters and geomagnetic field model uncertainty (Kuang et al. 2010;Kuang 2015, 2018). This was done by running assimilation experiments with various parameter values and uncertainty estimates, and comparing the resulting forecasts with geomagnetic field models. In the present work, the model uses free-slip boundary conditions for the velocity field at the inner core boundary (ICB) and the core mantle boundary (CMB), and fixed-flux boundary conditions for temperature at the ICB and the CMB. In our geodynamo model, there is also a 20-km thick electrically conducting layer, called the D ′′ -layer, at the base of the mantle, with a conductivity an order of magnitude lower that of the outer core. The D ′′ -layer is homogeneous in the present study. Since the magnetic field is no-longer a potential field in this layer, the SV in the layer differs from those of the potential field (Greff-Lefftz and Legros 1995; Kuang et al. 2017). In our model, the Rayleigh number is Ra = 1811 and the Rossby number is Ro = 1.25 × 10 −6 . A single magnetic Prandtl number Pr = 1 is used throughout.
The observations used for GEMS are the Gauss coefficients {g m l , h m l } (where l and m are spherical harmonic degrees and orders, and m ≤ l ) from geomagnetic field models. They represent only the internally generated magnetic field. In the outer core and the D ′′ -layer, the magnetic field can be decomposed into the toroidal ( B T ) and the poloidal ( B P ) components: where r is the radial unit vector, and T B , P B are the toroidal and poloidal scalars and are described in GEMS as where {Y m l } are the orthonormal spherical harmonic functions, and C.C. denotes the complex conjugate part of the expansion (similar expansions are also made for v and �ρ ). At the top of the D ′′ -layer, i.e. at , the toroidal field vanishes, i.e. j m l = 0 ; but the poloidal field coefficients b m l can be matched with the observations where r s is the mean Earth surface radius, and L obs is the highest degree of the observed geomagnetic field.
In the rest of this paper, we denote the poloidal spectral coefficients in (4) that are derived from the field models by (b m l ) o and the forecast (b m l ) f . The observations y o = (b m l ) o are assimilated into the forecast x f at analysis time t a to produce the analysis x a : where H is the observation operator matrix and K is the gain matrix, defined by: where P f is the forecast error covariance and R is the observation error covariance. The analysis x a is then used as the initial state of the dynamo system (1) for the next forecast. The forecast error covariance is computed from an ensemble of forecasts with N ens members (Sun and Kuang 2015) where µ f x is the ensemble mean forecast state. In our forecasts here, we take only the diagonal part of the observation error covariance R and is defined simply as the estimated variance of the field models, ((σ m l ) o ) 2 . The time dependent model for the observation error standard deviation is where α = 0.1 before 1900, 0.001 ≤ α ≤ 0.1 with an exponential decrease from 1900 to 2000, and where α = 0.001 afterwards. This observation error model accounts for the increased accuracy of geomagnetic models at later times, as well as the higher relative errors for the higher degrees. It has been tested through a number of numerical experiments with the goal of minimizing forecast errors, and is therefore regarded as an appropriate estimate relative to our geodynamo forecasts.
The geomagnetic field models that are used to provide the Gauss coefficients for our assimilation, and their maximum degree assimilated ( L obs ) are shown in Table 1.
Details of gufm1, CM4 and CM6 can be found in Jackson et al. (2000), Sabaka et al. (2004), and Sabaka et al. (2020) respectively. We only summarize here that the latter model uses a "comprehensive inversion" technique which co-estimates parameters for the Earth's internal field, using estimates of systematic and random errors, so as to produce an optimal separation of the internal and external fields. The values of L obs for each model here decrease further back in time, so as to avoid introducing errors from the smoothing or regularization in the earlier models. These values may in fact not be optimal, but our previous experimentation (Tangborn and Kuang 2018) has shown that reducing them for the earlier field models can lead to more accurate forecasts in the modern era.
The ensemble is initialized in 1590 using a long free model run from which uncorrelated ensemble members are generated, and subsequent analyses are computed every 20 years until 2000, when the analysis cycle is reduced to one year (the start of the forecast period using the CM6 model). P f is updated with the forecasts at each analysis time t a , but H is only updated if different L obs is chosen for making the analysis. We use an ensemble of 400 for these experiments, and we have found this size to result in forecast accuracy only slightly better than smaller ensembles of 200 or 256, and very similar to an ensemble of 512. Any further increases in the ensemble size would result in larger computational costs that outweigh small increases in accuracy.
Geomagnetic forecasts using a geodynamo model will naturally contain significant forecast errors due to mismatches in the dynamo parameter values and numerical approximations, and these can grow significantly during a multi-year forecast. But because these error sources develop on relatively long timescales, we can assume that the model error itself varies on time scales much longer than, e.g. the 5-year forecast period. This assumption allows us to implement a bias correction scheme for the SV prediction by producing a set of staggered forecasts, and then take differences between them to reduce substantially the bias (Kuang et al. 2010). The forecast y, of the observed poloidal field is defined as the projection of the state vector, x f onto the observation space:   We define the forecast error, ǫ(t) , as the difference between the forecast state, y f (t) , and the true state, y t (t): The forecast error is a combination of the model error and the growth of initial state errors, but we assume that the model errors remain the largest component. Consider a forecast that is initialized from an analysis at time t a and produces a forecast at time t a + δt . We estimate the SV from the change in the geomagnetic field after time δt: (9) The difference between these two sequences was shown to approximately remove the forecast error (Kuang et al. 2010), if the change in the model error between t a + δt and t a + δt is small compared to the SV. So we can calibrate the SV forecast with the following scheme where τ a = |t a −t a | = 1year , and t = t a + δt is the time of the SV forecast. We take δt = 5 years for these forecasts.
The details of the bias correction scheme can be found in Kuang et al. (2010).

Validation Experiments
We have validated the methodology by computing SV forecasts for the period 2009-2015. Two forecast experiments were conducted using the same assimilation run that started in 1590, and from 2000 it is carried forward with one year forecast and analysis cycles. In the first experiment (A), the assimilation is computed from 2000 to 2008; the analysis in 2008 is then used to make forecasts for the period 2009-2014. The second experiment (B) is similar, but the assimilation is computed from 2000 to 2009; and the analysis in 2009 is used to make the forecast for the period 2010-2015. The final SV forecasts are then performed with these two staggered forecasts: The SV for the Gauss coefficients ( g m l , h m l ) are then computed as follows. We show the details for g m l , while h m l follow the same procedure. In our assimilation, we utilize the scaled coefficients (12) y f (t a + δt) = y t (t a + δt) + ǫ(t a + δt) y f (t a + δt) = y t (t a + δt) + ǫ(t a + δt) where ẏ is the SV and ǫ is the SV error growth rate. We can reduce this error term by making use of the assumption of slowly varying model error and considering staggered forecasts that begin at two slightly different analysis times t a and t a , where |t a −t a | < δ . This gives us forecasts at times t a + δt and t a + δt in terms of the true state: for the SV forecasts. This is done because the parameters in the geodynamo model are far from those in the Earth's core, and assimilating g 0 1 would force the model into a regime that the current model grid cannot resolve. We first compute G m l from our forecasts, and their time variation Ġ m l via (14). The SV Gauss coefficients can then be calculated via We then approximate the axial dipole coefficient g 0 1 and its time derivative ġ 0 1 in (16) with the values from CM6 (Sabaka et al. 2020) as follows: The formulations (13)-(17) are the basis for the forecasts of the mean SV for the period 2010-2015, and are compared with CM6 for each degree l and order m up to l = 8 . In this case the years in A correspond to τ and those in B to τ . Some examples are shown in Figs. 1, which demonstrate the forecast field, ỹ f , for different l, m. The slope of each dashed line is the scaled forecast from (15) while the solid line is the scaled SV computed from CM6 coefficients. The coefficients shown here are representative of the range in accuracy obtained by the ensemble forecast and bias correction algorithm. The (l = m = 1) and (l = m = 2) coefficients in Fig. 1 are typical of coefficients that predict the CM6 field model with a higher degree of accuracy than others. The (l = 3, m = 2) in Fig. 2(a) is much less accurate accurate case (for g m l ) with a slope of opposite sign, while (l = 5, m = 3) has a typical mid-range accuracy for the dimensionless scaled SV forecasts. The final dimenional SV values determined by the validation calculations are shown in Table 2 These bias corrected SV forecasts are computed each year out to 2015, and the five year time average SV can be computed from them.
We summarize the SV forecasts by plotting mean square SV fields ( MS SV ), shown in Fig. 3 for the EnKf forecast (dashed line), CM6 (solid line) and an extrapolation of CM6 from 2009 (dotted line) as a function of degree L. The mean square SV is computed using Here it can be seen that the EnKF forecasts are more accurate than the extrapolated SV for L = 2 and shows similar accuracy for L = 3, 4 , while the accuracy for L = 5 to L = 8 is significantly worse. The individual coefficients for these three cases are given in Table 2. We discuss the possible sources of error and potential solutions in the Summary and Conclusions section.

Mean SV for 2020-2025
We repeat the same method to make the SV forecast for the period from 2020 to 2025, except that we use three staggered forecasts, again using an ensemble of 400 members each. Starting from the analysis computed for 2000, we have continued the assimilation and forecast (with the annual analysis cycle) forward in time for three sets of experiments, with the final analysis times in 2017, 2018 and 2019, and the forecasts are performed in the next six years. The experiments, labeled A, B and C; are summarized in Table 3. With these three sets of forecasts, we can use up to the following three different combinations of the differences between two staggered forecasts: where again τ a = 1 year and the AC experiments are 2 years apart. The SV for the forecast Gauss coefficients is computed as before, along with the SV for CM6.
The bias corrected SV forecast for 2025 is a combination of the different forecasts, BA, CB or CA from (19). There are two main reasons for our last analysis ending at 2019: one is that the predictive field from CM6 for 2020 will be used as the initial condition for our forecast calibration; the other is the SV from CM6 for the period 2019-2020 can be used to estimate our model bias. Some representative forecasts that demonstrate the bias correction  Tangborn et al. Earth, Planets and Space (2021) (19). The slope of each line is the time varying SV, and we confirm the consistancy of each of the experiments with the CM6 SV in each case. For example, Fig. 4(a) shows the forecasts for (l, m) = (1, 1) , and the SV for g 1 1 requires a negligble correction since the three experiments have essentially the same SV. The SV for h 1 1 shows a small change from the bias correction and the AB correction is used since its slope is closest to the CM6 estimate. The other two SV estimates are then discarded.
The final calibrated mean SV for the period 2020-2025 is computed as the following: first, the optimal annual SV forecasts ġ m(f ) l are made from our data assimilation system for 2020-2024, with the same procedure discussed in the previous section. Then, the time varying Gauss coefficients are computed via Finally, the mean annual SV Gauss coefficients are determined via The uncertainties are determined by the standard deviation of the forecast ensemble. The final results are given in Table 4 along with the associated uncertainties, which are the ensemble standard deviations from the ensemble forecast.

Summary and Conclusions
We have used an EnKF to forecast the geomagnetic secular variation (Table 4) for the time period 2020-2025, using 400 ensemble members. The development of ensemble forecasting techniques has contributed to the SV forecasts, particularly due to the ability to estimate uncertainty, which is newly available with this geomagnetic forecast. The validation forecast for 2015 shows that the SV forecasts compare favorably to the SV from CM6 (relative to an extrapolated SV from 2009-10) for degrees for k = 1, 2, ·, 5.  2-4 (degree 1 is excluded because the forecasts are scaled by (g 0 1 ) o ), but show poorer accuracy for degrees 5-8. The lower accuracy in the higher degree terms is likely due to shorter timescales for model bias for these coefficients, which could violate the assumption that the model bias is constant over a period of several years. This would make it difficult for the model to predict geomagnetic jerks (Aubert and Finlay 2019). The uncertainty estimates provided with the candidate model are seen to be too low, based on the validation experiment, in part due to limitations on the constant bias assumption. However, the 2020-25 SV used a series of 3 forecasts which helped eliminate outlier forecasts, which benefited the final SV forecast, whereas the 2010-2015 forecast only used 2 forecasts. While this is the first attempt at an SV forecast using this EnKF system, we expect that improvements can be made to the forecast in the future. In particular, we plan to gradually increase the range of geodynamo parameter values, and explore alternate timescale relationships between the dynamo model and the geomagnetic field models.
The large ensemble has placed a constraint on this forecast, because at present it is computationally challenging to run the ensemble assimilation system for 7000 years, as was done in Kuang et al. (2010) using an OI assimilation algorithm. Future work will investigate the use of localization schemes that would enable the use of a much smaller ensemble that could be run for longer time periods.