Predictions of the geomagnetic secular variation based on the ensemble sequential assimilation of geomagnetic field models by dynamo simulations

The IGRF offers an important incentive for testing algorithms predicting the Earth’s magnetic field changes, known as secular variation (SV), in a 5-year range. Here, we present a SV candidate model for the 13th IGRF that stems from a sequential ensemble data assimilation approach (EnKF). The ensemble consists of a number of parallel-running 3D-dynamo simulations. The assimilated data are geomagnetic field snapshots covering the years 1840 to 2000 from the COV-OBS.x1 model and for 2001 to 2020 from the Kalmag model. A spectral covariance localization method, considering the couplings between spherical harmonics of the same equatorial symmetry and same azimuthal wave number, allows decreasing the ensemble size to about a 100 while maintaining the stability of the assimilation. The quality of 5-year predictions is tested for the past two decades. These tests show that the assimilation scheme is able to reconstruct the overall SV evolution. They also suggest that a better 5-year forecast is obtained keeping the SV constant compared to the dynamically evolving SV. However, the quality of the dynamical forecast steadily improves over the full assimilation window (180 years). We therefore propose the instantaneous SV estimate for 2020 from our assimilation as a candidate model for the IGRF-13. The ensemble approach provides uncertainty estimates, which closely match the residual differences with respect to the IGRF-13. Longer term predictions for the evolution of the main magnetic field features over a 50-year range are also presented. We observe the further decrease of the axial dipole at a mean rate of 8 nT/year as well as a deepening and broadening of the South Atlantic Anomaly. The magnetic dip poles are seen to approach an eccentric dipole configuration.


Introduction
Over the past century, the geomagnetic field has undergone significant changes reflecting the vigorous dynamics of Earth's liquid iron core. The focal point of the South Atlantic Anomaly (SAA), a region of low magnetic field intensity at the Earth's surface, has moved from the shores of Southeast Brazil to its current location in Northern Argentina, while its intensity has decreased by about 26nT/year (Finlay et al. 2010). The north and south magnetic dip poles, locations where the field is vertical to Earth's surface, have both moved northwestwards, the north dip pole accelerating to a speed five times higher than that of the south dip pole during the past 20 years . The most intriguing change, however, is the decrease of the geomagnetic field dipole by about 15nT/year over the past century (Finlay 2008). These are only three particularly prominent examples of the geomagnetic secular variation (SV).
Understanding and predicting these variations is very important for modeling our space weather environment.

Open Access
*Correspondence: sanchezs@mps.mpg.de 1 Max Planck Institute for Solar System Research, Göttingen, Germany Full list of author information is available at the end of the article For instance, in the vicinity of the SAA, the radiation belts reach to much lower altitudes and can harm the electronics aboard satellites (Schaefer et al. 2016;Bourdarie et al. 2019). Predicting the evolution of the SAA is therefore important for planning satellite operations. Short-term predictions are also needed for surface exploration studies and regional magnetic anomaly surveys that depend on accurate maps of the main magnetic field. Furthermore, many commonly used electronic devices use maps of the main magnetic field as directional information.
The International Geomagnetic Reference Field (IGRF) offers a community-based model of the main magnetic field and a SV estimate. There is a new IGRF every 5 years; the newest, IGRF-13, covers the year 2020 (Alken et al. 2020a). The IGRF is a synthesis of a number of model candidates constructed using different modeling strategies. In addition to the main field (B) up to spherical harmonic degree 13, the IGRF also provides a SV ( Ḃ ) estimate that allows predicting the field 5 years later with a simple linear extrapolation: B(t 0 + 5 years) = B(t 0 ) +Ḃ × (5 years) . Though this definition calls for the mean SV over the 5-year interval, most contributions provide the SV at the year t 0 . Predicting the evolution 5 years beyond the available data requires some kind of dynamical model, which goes beyond typical field modeling studies. However, more and more candidates have started incorporating at least some key physical aspects of the core field dynamics, which could potentially improve the magnetic field forecasts. Some derive SV estimates from the interaction of the magnetic field and the inverted flow just underneath the core-mantle boundary (CMB, e.g., Hamilton et al. 2015). Others rely on full 3D-dynamo simulations within a data assimilation framework (e.g., Kuang et al. 2010;Fournier et al. 2015). Our contribution to IGRF-13 is a SV that falls into the latter category. In this paper, we explain the details of the model and discuss different approaches for estimating the SV.
Geomagnetic data assimilation and inversion schemes based on 3D-dynamo simulations have greatly improved over the past years (see the review by Fournier et al. 2010). They have evolved from early validation tests (e.g., Kuang et al. 2008;Aubert and Fournier 2011) to recent applications, mainly involving reanalyses of the core flow (Aubert 2014; Barrois et al. 2018) and long-term forecasts (Aubert 2015). A central component for the success of any sequential data assimilation scheme is the choice of a sufficiently realistic background dynamical model and the representation of its covariance. The model covariance is particularly important due to its role in propagating information from the observed to the unobserved parts of the system. The covariance can be approximated via an ensemble, such as a collection of snapshots of a given numerical model. The ensemble Kalman filter (EnKF, Evensen 1994) proposes a generalization of this concept to the time-dependent estimation problem, sequentially iterating between forecast and analysis steps. In the analysis step, the ensemble is used to infer the model covariance and to estimate the state of the model in light of the available observations. In the forecast step, the ensemble is integrated forward in time with the use of the numerical model and diverges due to its nonlinear nature. The divergence is often characterized by the e-folding time, which is related to the predictability of the system .
A first EnKF application to geomagnetism by Fournier et al. (2013) suggested that an ensemble of about N e ≈ 500 numerical dynamo simulations is required for the assimilation to converge. Since dynamo simulations at more realistic parameters are extremely costly, it would be beneficial to decrease this number. This is problematic, however, since the covariance matrix quickly degrades with decreasing ensemble size, likely causing the assimilation scheme to diverge or 'blow up' . The reason is the growing impact of spurious correlations, which can become particularly problematic when observation uncertainties are low (Gottwald and Majda 2013). Sanchez et al. (2019) demonstrated that the geomagnetic EnKF can be stabilized by keeping only the most prominent elements in the covariance matrix while setting all others to zero. Using a spectral representation, they showed that the correlations within the same spherical harmonic order clearly dominate the covariances of numerical dynamo simulations.
In this study, we show that the ensemble size can be further reduced by an additional covariance localization which exploits the dominant equatorial symmetries observed in dynamo simulations. The new scheme is used for the EnKF assimilation of geomagnetic field models covering the past two centuries, namely the COV-OBS.x1 geomagnetic field model (Gillet et al. 2015) from 1840 to 2000 and the Kalmag model (Baerenzung et al. 2020) from 2001 to 2020. We propose the resulting SV in 2020 as a candidate model for the IGRF-13 (Alken et al. 2020b) and also attempt to predict the evolution of the main field until 2070. This paper is organized as follows: we first describe the methods comprising the setup of the background model, the data assimilation scheme and the proposed covariance localization. Next, we introduce the sequence of observations used for the assimilation as well as a correction of the observational uncertainties. We further describe the results from the assimilation with additional localization, introduce the candidate SV model for the year 2020 and provide long-term predictions of the magnetic field for the next five decades. Finally, the last sections offer a discussion on the results and draw the conclusions of our study.

Background dynamo model
The dynamo model solves for disturbances in flow U , magnetic field B and density around an adiabatic stationary ground state. The differential equations are formulated in a frame of reference that co-rotates with the outer boundary, which represents the Earth's mantle. The non-dimensional control parameters that rule the system are the Ekman number E = ν �D 2 , the modified Rayleigh number Ra = g 0 F /4πρ� 3 D 4 , the Prandtl number Pr = ν/κ , the magnetic Prandtl number Pm = ν/ , and the aspect ratio a = r i /r o . Here, g o is the outer boundary reference gravity, α the thermal expansivity, F the mass anomaly flux, ν the kinematic viscosity, κ the thermal diffusivity and the magnetic diffusivity. More details of the model are given in Christensen and Wicht (2015) or Wicht and Sanchez (2019). We used the numerical code Parody (Dormy 1997;Aubert et al. 2008) with the additional implementation of the assimilation scheme PDAF ).
The background model in our assimilation closely resembles the one used in our previous article Sanchez et al. (2019). We will refer to the article as S2019 here. All parameters are the same, except for the Rayleigh number, which we have increased by a factor two. The modified Rayleigh number is Ra = 10 −5 , the Ekman number is E = 10 −4 , the Prandtl number Pr = 1 and the magnetic Prandtl number is Pm = 10 . It seems reasonable to (1) assume that the tests and methods discussed in S2019, in particular the spectral covariance localization, still apply. The model contains an electrically conducting inner core and uses no-slip mechanical and homogeneous heat flow boundary conditions. It uses a decomposition of flow and magnetic fields in poloidal and toroidal contributions and an expansion in spherical surface harmonics (SH). The spectral representation required for our assimilation scheme is thus readily available. Recent high-end dynamo models run at Ekman numbers as low as E = 10 −7 (Schaeffer et al. 2017a). They may yield more Earth-like results, but also require extremely high spatial resolution and very short time steps. Since we run an ensemble of models in parallel, we chose a significantly larger Ekman number. We can think of the dimension-less model parameters as ratios of time scales. Particularly important for the dynamo mechanism are the magnetic diffusion time τ , the flow advection time τ U , and the rotation time τ � . The ratio of magnetic diffusion to advection time is the magnetic Reynolds number Rm = τ /τ u . The ratio of rotation to magnetic diffusion time is the magnetic Ekman number E =E/Pm= τ � /τ . Christensen et al. (2010) have shown that Earth-like dynamos occupy a wedge in the Rm/E parameter space that is accessible by numerical simulations. However, this does not mean that both parameters are Earth-like. Earth-like magnetic Reynolds number values of about 10 3 can actually be reached by the simulations, but E remains much too small even for the most extreme models where it reaches 10 −6 compared to 10 −9 for Earth. The milder parameters used here yield Rm=430 and E = 10 −6 , which is not ideal but still reasonable. The simulation lies in the Earth-like wedge defined by Christensen et al. (2010), but the simulated magnetic field is slightly less dipolar than the geomagnetic field and is too strongly concentrated in patches at the CMB.
In order to compare the non-dimensional simulations to observations, we have to choose a time scale. Being interested in the magnetic field dynamic, the secular variation time scale seems a natural choice. While the dipole varies very slowly, the typical time scale of higher order harmonics roughly decreases like τ SV /ℓ , where ℓ is the spherical harmonic degree and τ SV is known as the SV time scale. Fitting this behavior to SV observations for degrees ℓ = 2 to 14 allows estimating τ SV for the Earth and for the simulations. The secular variation time scale seems to be directly coupled to the advection time scale τ u (Christensen and Tilgner 2004). Fixing τ SV to the observed value thus means that we also fix the advective time scale. Because the parameters are not those of Earth, however, other time scales remain unrealistic. For example, the fact that the magnetic Reynolds number is about a factor two too low means that the magnetic diffusion time τ is about a factor two too short. Compromises are required here and other studies interested in different aspects could require a different time rescaling.
While the analysis of Lhuillier et al. (2011b) suggests τ SV = 415years, we find that rescaling the assimilation output with τ SV = 485year yields better hindcast results. We have therefore adopted this somewhat slower reference time scale. The sensitivity of the assimilation to the rescaling will be explored in forthcoming studies. The magnetic field was rescaled so that the mean axial dipole of the initial ensemble matches the average value of the observations over the period 1840 to 2020. Lastly, the flow is rescaled by a factor D/τ SV .

Data assimilation scheme
The geomagnetic data assimilation framework employed here has been described in detail in S2019. It uses an EnKF based on 3D-dynamo simulations with the Parody-PDAF code ). An ensemble of N e dynamo simulations defines the background model. At a given epoch t i , each of the ensemble members e is characterized by its state vector where P and V are the poloidal magnetic field and flow potentials, T and U the toroidal field and flow potentials. The index k refers to radial levels in the spherical shell, ℓ and m are, respectively, the degree and order of the SH expansion and the superscript T refers to the transpose. As in S2019, given the chosen model dimensions, the state vector has a dimension of N x ∼ 2.6 × 10 6 .
From the ensemble of simulations one can evaluate the mean dynamo solution and the covariance where the † represents the complex conjugate. In our EnKF approach the ensemble of dynamo simulations defines the background (prior) model statistics, which is combined with geomagnetic field models through a sequence of forecast and analysis cycles. The forecast is the propagation of each dynamo ensemble member from time t i−1 to time t i by the model's numerical integration. When describing the dynamo state with the vector x and its time propagation with the operator M , the forecast can formally be written as (6) where indexes i and e correspond to time and ensemble member, respectively. The starting condition x a i−1,e results from the previous analysis step. Whenever observations y are available, the analysis step is performed where y i,e is a random realization of the observation from its distribution N ( y i , R i ) , where R i is the observation error covariance. It provides a correction of the forecast state based on the differences to the observations, where H i is the observation operator that screens out the observed variable from the state vector x . The Kalman Gain matrix is given by where P f i is the model forecast error covariance. The Kalman Gain matrix is essential for propagating information from observed to unobserved parts of the model state.
The diagonal elements of the model covariance matrix P f /a i provide information about the variances of the model variables, while the off-diagonal elements quantify the covariances. The variances determine the spread of the ensemble and also serve as a proxy for the model uncertainties; the covariances show the linear dependencies of the model state entries. The model covariance matrix is typically huge with N 2 x elements. Fortunately, however, in the EnKF scheme we only have to consider the elements describing the couplings of the model with the observables, that is, the low degree Gauss coefficients at the CMB. A subset of the initial model covariance P f 0 for the observable poloidal magnetic field potential at the top of the core P m ℓ (r o ) is shown in Fig. 1a. The Figure demonstrates that small ensembles (for instance with N e = 64 ) result in spurious correlations that tend to populate otherwise sparse regions of the covariance. As described in S2019, these spurious correlations can have a large impact on the performance and stability of the EnKF, but may be remedied by covariance localization.

Covariance localization
The term 'localization' originates from atmospheric data assimilation where a covariance mask is used to filter out likely spurious correlations over long distances (Houtekamer and Mitchell 2001). When C is the mask, the modified localized covariance is simply where • denotes the Schur product (element-wise). For systems like Earth's atmosphere where local effects dominate due to the thin-shell dynamics, a distance-dependent tapering makes sense. For the Earth's core, however, we expect the flows to be nearly geostrophic (i.e., to show minimal variation in the direction of the rotation axis) due to the rotation-dominated dynamics in a thick shell. This means that the flows should be dominantly equatorially symmetric, which is the property we will exploit. Earth's magnetic field, on the other hand, should be dominated by the equatorially antisymmetric axial dipole.
The preferred symmetries clearly shape the spectral covariance matrix used in geodynamo assimilation schemes. The action of the flow contribution with a given wave number m on the axial dipole causes variations in all magnetic field contributions with the same m. This leads to the block structure along the diagonal of the spectral covariance matrix shown in Fig. 1a. The respective mask accounting for this pattern is given by where k represents the index of the radial levels of the spherical shell. S2019 demonstrated that using this mask stabilized their assimilation scheme for an ensemble size of N e = 512.
Here we go a step further and also implement the checkerboard pattern evident in Fig. 1a. An equatorially symmetric flow acting on an equatorially antisymmetric field, like the axial dipole, produces again equatorially antisymmetric field. Poloidal field contributions where ℓ + m is even (odd) are equatorially symmetric (antisymmetric). It follows then that for this interaction an even (odd) order m only correlates with odd (even) degrees ℓ (Bullard and Gellman 1954). The spectral localization mask, illustrated in Fig. 1b, is therefore given by where mod is an operator defining the remainder of the division of the two arguments. The localization mask C ℓm corresponding to the subset of the poloidal field covariance is shown in Fig. 1b. The same pattern should apply to covariances between the potentials of poloidal magnetic field and the toroidal flow (S2019).
The covariance matrix between poloidal magnetic field and the poloidal flow potentials (see Fig. 12 from Appendix) obeys a shifted checkerboard pattern, consistent with The reason is that for the case above the potentials bear different dominant equatorial symmetries. The shifted pattern also applies for the covariance between the poloidal field potential and the toroidal field and codensity potentials.

Observations
We assimilate the main field component of a previous version of the Kalmag field model (Baerenzung et al. 2020) spanning the years 2001 to 2020. This version of Kalmag relies on a Kalman filter and a subsequent Kalman smoother for the assimilation of CHAMP and . Fig. 1 a Subset of the normalized model covariance P f 0 , corresponding to correlations of the poloidal magnetic field potential at the core surface P m ℓ (r o ) calculated from the ensemble of dynamo simulations. Note that for illustration purposes the matrix is a composite of two correlation matrices based on ensembles of different sizes: the upper half is calculated using N e = 512 and the lower half with N e = 64 . b The associated localization matrix C ℓm . The matrices are organized by the stacking of the spherical harmonic degrees ℓ within the corresponding order m, truncated at ℓ = 14 , and the upper insets magnify the region where m = m ′ = 3 Swarm satellite data. It models the main surface field up to degree ℓ = 20 with a classical set of Gauss coefficients g m ℓ and h m ℓ and also provides uncertainties. The Gauss coefficients are related to the state vector through eqs. (11) and (12) from Sanchez et al. (2016). The field evolution at the top of the core by Kalmag is parametrized with a second-order auto-regressive process. The model uses data until September 22 nd 2019 and provides a forecast for 2020, which is an IGRF-13 candidate.
S2019 demonstrated that the performance of the assimilation scheme improves when including data covering long time intervals. We therefore extend the assimilation interval before the Kalmag time span with the COV-OBS.x1 model for the years 1840-2000 (Gillet et al. 2015). COV-OBS.x1 uses data from magnetic observatories, surveys and satellite missions such as POGO and Magsat. It models the surface field up to degree ℓ = 14 and provides uncertainty estimates based on an ensemble approach. COV-OBS.x1 is regularized in time by temporal cross-covariances associated with a second-order auto-regressive stochastic process.
A comparison of the uncertainties σ o from COV-OBS. x1 and Kalmag models in Fig. 2 reveals a growing discrepancy over the last 60 years. COV-OBS.x1 uncertainties decrease very fast after 1990, reaching a 10 µ T RMS value at the core surface, and likely underestimate the true errors in recent epochs (as also pointed out by Barrois et al. 2018). Kalmag uncertainties, on the other hand, maintain a 30 µ T RMS level during the same period. We assume here that the Kalmag uncertainties are more realistic and use them as a lower bound. Each uncertainty coefficient from COV-OBS.x1 lower than the corresponding Kalmag value in 2001 is simply set to the corresponding Kalmag value. The resulting corrected COV-OBS.x1 RMS uncertainty estimate is shown in Fig. 2.
We choose to only assimilate Kalmag up to ℓ = 18 , since the signal-to-uncertainty ratio is small for ℓ ≥ 19 . The composite field model is assimilated every t o = 5 years for COV-OBS.x1 and every t o = 1 year for Kalmag.

Results
The results shown in the following explore the effects of different assimilation setups while keeping the observations (a sequence of snapshots of the Kalmag model) and dynamo model the same. The first section shows the impact of the covariance localization and the ensemble size on the forecast errors. The following section discusses the quality of 5-year forecasts based on different SV estimates using hindcast tests. We afterwards present our candidate model for the SV in 2020 and compare it to the final IGRF-13. Lastly, predictions of the Earth's main magnetic field for the next 50 years are shown.
We quantify forecast errors (sometimes called 'O-F' , observation minus forecast) by They are often compared to the forecast uncertainties, estimated by the ensemble spread Figure 3 shows the evolution of the RMS forecast errors ε f and the forecast uncertainties σ f at the top of the Earth's core for different ensemble sizes N e and localizations during the assimilation interval. The case corresponding to the scheme used in S2019 with N e = 512 and C m localization serves as a reference case and is shown by the thick black line. The forecast uncertainties σ f in Fig. 3b closely follow the evolution of the observation uncertainties σ o shown in Fig. 2, at least for the period until 2000 covered by COV-OBS.x1. The sudden decrease in σ f after 2000 is caused by the decrease in the assimilation interval from t o = 5 years to 1 year. Forecast uncertainties seem to level-off at the end of the assimilation scheme, while the errors slightly grow. The predicted uncertainties should be a representation of the discrepancies between the model and reality. Therefore discrepancies between the forecast errors and uncertainties in  Fig. 3a and b are likely related to systematic model errors, i.e., biases, not accounted for in our approach. An analysis of such bias is very important in a data assimilation approach and will be tackled in forthcoming studies. Figure 3a demonstrates that the forecast quality systematically deteriorates with decreasing ensemble size. This is accompanied by a contraction in the ensemble spread (variance) and hence a decrease of the uncertainty estimate shown in Fig. 3b. The smaller ensemble size covers a smaller fraction of the possible solution space and tends to underestimate variances and overestimate the covariances. The latter can lead to unrealistic updates of hidden (non-observed) model variables, which progressively deteriorate each analysis but further decreases the variances (Van Leeuwen 1999). The effect becomes more prominent when the assimilation interval is reduced to only 1 year after 2000. In this situation, the dynamo dynamics does not have sufficient time to spread up the ensemble in the forecast step before the next analysis cycle.

Impact of localization
When using the C m localization, the degradation becomes more noticeable for N e = 128 and unacceptable for N e = 64 , as illustrated in Fig. 3a. The more restrictive localization C ℓm helps to mitigate the problem, in particular for N e = 128 and N e = 64 . Now 50% of the remaining possibly spurious elements in the covariance matrix are simply set to zero, which alleviates artificial amplification of the unobserved state variables and mitigates the decrease in ensemble spread. The fact that the assimilation still deteriorates with decreasing N e indicates that at least some of the remaining correlations are still spurious and that the ensemble spread is too small. All in all, Fig. 3 shows that the differences in forecast error and uncertainties for the cases with N e > 100 and C ℓm remain relatively small. In the following, we use N e = 256 and the C ℓm localization to be on the safe side.

Five-year prediction tests
In order to assess the applicability of our method for IGRF SV predictions, we perform a group of 5-year hindcast tests. In the EnKF context, hindcast means that the analysis steps stop at a time t i , after which observations are still available. The ensemble then freely evolves to time t f = t i + δt , where δt = 5 year is the prediction window of interest. This allows assessing the quality of the forecast by directly comparing the forecast with past or present observations, which consist of snapshots of the main field (in the form of Gauss coefficients) of the Kalmag model. Figure 4 illustrates the evolution of two observed Gauss coefficients, the axial dipole g 0 1 and h 3 5 , and their respective SVs. We have stopped the analysis at t i = 2015 and observe the behavior of the forecast until t f = 2020 . The result of the assimilation is compared with three field models: Kalmag, CHAOS-7 (Finlay et al. 2020) and IGRF-12 . While the coefficients in the field models agree very closely, the SV can differ significantly, in particular between 2010 and 2014. This period corresponds to a gap in satellite data between the CHAMP and Swarm missions, and highlights the temporal regularization in CHAOS.
The axial dipole in our assimilation scheme follows the field observations closely, as seen in Fig. 4a. The SV is also reproduced reasonably well after each analysis step, as shown in Fig. 4c, despite not being directly observed. The SV tends to deviate from the observations during the forecasts sometimes, implying that the dynamo model seems to often favor unrealistic secular accelerations g 0 1 . This is not so much of an issue during the final 5-year forecast, however, where the model dynamics successfully captures the gradual decrease of ġ 0 1 . The axial dipole in our 5-year prediction for 2020 differs by only 5 nT from the observations. The evolution of coefficient h 3 5 agrees less well with the observations (Fig. 4b and d). Particularly problematic is the epoch between 2010 and 2015 when ḣ 3 5 tends to too large values. From 2015 to 2020, the observed ḣ 3 5 first increases, but then decreases again. The 5-year assimilation forecast does not capture this variation and keeps predicting an always positive ḣ 3 5 that only gradually slows down at the end of the epoch. The predicted value for 2020 ends up deviating by about 7 nT from the observed one. The prediction quality of the other large-scale Gauss coefficients fall in between the two examples given by g 0 1 and h 3 5 , not following an obvious pattern. Figure 4 also shows the linear predictions based on the SV values from the last mean analysis step from the assimilation, coined here as LP. For g 0 1 this yields a worse result than the dynamic forecast, which captures the gradual decrease in the SV. For h 3 5 , however, the linear prediction incidentally outperforms the assimilation because of the oscillatory behavior of the SV. In order to quantify the overall quality of the forecasts, we calculate the RMS difference between the predictions and CHAOS-7 at Earth's surface truncated to ℓ = 8 between the initial and final prediction times t i and t f . Figure 5 compares this prediction error ε f for four 5-year hindcast tests from 2005 to 2020 from the dynamic forecast from our assimilation (black triangles). In Fig. 4 Comparison of the a axial dipole g 0 1 and b h 3 5 coefficients and their respective SV c ġ 0 1 and d ḣ3 5 between an assimilation hindcast with t i = 2015 and t f = 2020 (black solid curve), the Kalmag model (blue curve), CHAOS-7 (red curve) and IGRF-12 (green star). The shaded areas represent the standard deviation of the assimilation and Kalmag model. Also shown in a and b are predictions based on the linear extrapolation (LP) of the coefficients using SV in 2015 from the assimilation (dashed black line), CHAOS-7 (dashed red) and IGRF (dashed green) addition, the figure also includes the errors ε f with respect to CHAOS-7 based on the linear predictions (LP) using the SV at t i from our mean model (black squares), CHAOS-7 (red circles) and IGRF (green stars).
Since the CHAOS model serves as the reference, the LP-CHAOS prediction error represents the optimum value that can be achieved with the linear forecast in Fig 5. The error from the LP-IGRF forecast is more than twice as large for 2005 and about 20% larger than this optimum for 2015 and 2020. The fully dynamic assimilation forecasts perform slightly better for 2005 and 2020, but worse for 2010 and 2015. The error in the linear prediction based on the assimilation SV, however, comes close to the CHAOS optimum, except for 2005. Our analysis of the h 3 5 coefficient suggests that the secular acceleration (SA) is not well constrained or problematic. For the 5-year forecast, magnetic diffusion can in principle be neglected and the SA is then given by where ∇ H denotes the horizontal part of the divergence operator and U H the horizontal flow. The first term on the right-hand side depends on the SV, which is often decently realistic in our assimilation. The second term on the right-hand side of Eq. 18, however, depends on the flow changes and thus crucially on the force balance in the Navier-Stokes Eq. 1.
From analyzing a number of dynamo simulations, Christensen et al. (2012) conclude that the fast secular acceleration in geomagnetic data is mostly explained by this second term. Fournier et al. (2015) use the Coupled Earth dynamo model ) in combination with an assimilation scheme that differs from the approach chosen here. Instead of adopting a sequential approach, the authors perform independent snapshot analysis steps and assimilate not only the field, but also the SV as an additional source of information. This means that the first term in the right-hand side of Eq. 18 is less problematic. Fournier et al. (2015) show that the forecast quality increases significantly when the flow is kept constant in their dynamo simulations (what is known as a 'steady flow' approximation). They therefore argue that the second term on the right-hand side of Eq. 18 is overestimated in their dynamo-based forecasts.
We have repeated the exercise of using a steady flow in the assimilation forecast, meaning that after the last analysis step, only the induction equation (Eq. 2) is integrated forward in time by the ensemble. As seen in Fig. 5, the forecast from the assimilation only slightly improves for 2015 and 2020 by using the steady-flow approximation (white triangles), and except for 2020 it remains significantly worse than the LP-assimilation forecast, for which the SA is zero. This means that in the assimilation forecast the first term on the right-hand side of Eq. 18 is also problematic. The SV ( Ḃ r ) is often reasonably close to observations, (for instance, for the 2015-2020 hindcast case 55%, 78% and 91% of the time the error lies beneath 1 σ , 2 σ and 3 σ , respectively). Therefore, the horizontal flow U H is likely an important source for unrealistic accelerations.
However, Fig. 5 also shows that the quality of the assimilation forecasts improves over time and reaches the level close to the LP-assimilation forecast in 2020. This may indicate that the dynamo ensemble is still gradually converging to more realistic dynamics after the transition to the higher observation cadence in 2001. Since the LP forecasts perform better than dynamic forecasts for most of the hindcasts, we have decided to propose the LP forecast as a candidate model for the IGRF-13. For predictions beyond the 5-year horizon, however, we use the full dynamic ensemble forecast.

SV candidate model
Our assimilation of COV-OBS.x1 and Kalmag from the year 1840 to 2020 reveals an intricate pattern of flow and SV underneath the CMB, shown in Fig. 6. The SV is particularly strong underneath Siberia and around the west coast of Africa. The pattern mostly agrees with the SV from CHAOS-7, but the overall amplitude is slightly smaller. The flow reflects the dichotomy between the RMS errors ε f of hindcasts starting in t i and ending in t f = t i + 5 years. The predictions correspond to the full dynamic mean forecast from the assimilation (black triangles) and linear approximations (LP) based on the SV at t i from our mean model (black squares), IGRF (green stars) and CHAOS-7 (red circles). Also shown is the mean forecast error using a steady-flow approximation after the last analysis in t i (white triangles). The prediction errors are truncated at ℓ = 8 and evaluated at the Earth's surface with respect to CHAOS-7 (Finlay et al. 2020) magnetically quiet Pacific and the more active Atlantic hemispheres. Meandering westward flows at low Atlantic latitudes explain the notorious westward drift of magnetic features observed in this region. Previous studies show a clearer westward stream of larger extent over this area (Bärenzung et al. 2018;Barrois et al. 2018;Aubert 2020). It is suggested that these flows are part of an eccentric gyre that might reach deep into the core (e.g., Barrois et al. 2018;Aubert 2020). In our model, there is an overall westward flow region extending deep into the core, but comprising smaller-scale structures which reach out to the surface. They are apparent, for example as an eastward flow underneath Arabia and the Indian Ocean.
Our assimilation scheme also interprets the strong SV patches at high northern latitudes differently than other studies. For example, Livermore et al. (2017) and Aubert (2020) infer localized westward jets as the origin of this signal, while our assimilation suggests large-scale eddies. The tangent cylinder (TC) that touches the inner core Fig. 6 Instantaneous radial SV ( Ḃ r ) at the CMB and horizontal flow ( U H ) for the 2020 analysis of our assimilation. The top panel shows a Mollweide projection while the bottom panels display views on the north and south poles. The SV is truncated at ℓ = 18 , the maximum Kalmag degree assimilated in our scheme. The dashed circle represents the TC boundary is often regarded as a dynamical flow barrier. The high latitude flows by the aforementioned authors are dominated by the westward jets just outside the TC and only slower flows cross the TC. Our flows, however, cross the TC in many locations, as shown in Fig. 6. Aubert (2020) suggests that the westward jet along the TC turns into a pronounced southward-directed flow underneath China. This feature is also not present in our solution.
As previously mentioned, we have decided to contribute to the IGRF-13 SV, which estimates the variation of the magnetic field between 2020 and 2025. A table listing the coefficients from our SV candidate model is given in Table 1 of the Appendix. Figure 7a shows our candidate, truncated to the IGRF SV resolution ℓ = 8 and projected to the Earth's surface. Due to the upwards projection many of the small details in Fig. 6 have vanished and the SV underneath Siberia is now negative throughout. The pattern is reminiscent of the previous IGRF-12 since the large-scale contributions change very little over a 5-year period. In total, 14 groups have contributed with SV candidates for the IGRF-13, which is a considerable increase from only 9 candidates for IGRF-12. Figure 7b shows the difference between our candidate and the final IGRF-13 SV model, a weighted average of all 14 contributions. The largest discrepancies amount to about 10% of the peak SV and lie mainly beneath Eurasia. This encloses the region where Aubert (2020) reports the localized westward jet along the TC and the southerly flow underneath China.
The power spectra of the differences between all candidate models and the IGRF-13 are shown in Fig. 8 (gray lines). We observe that our model lies very close to the final IGRF-13 (red line). Also shown are the uncertainties in our candidate model (dashed red line), i.e., the spread of the ensemble of dynamo models, calculated using Eq. 17, and the dispersion of all SV candidates around the IGRF-13 (dashed black line). The high level of dispersion of the 14 candidates relates to the diversity of modeling approaches (Alken et al. 2020a). Our uncertainties closely follow the differences between our model and the IGRF, which points to a decent error estimation by the ensemble of dynamo models, at least on the 5-year time scale.

Long-term predictions
We have shown above that the instantaneous SV offers a better approximation to geomagnetic field changes on the 5-year time-scale than the dynamo ensemble. However, when considering decadal or even longer time scales, the dynamo dynamics starts to pay off (Aubert 2015). We have therefore used the dynamo ensemble dynamics for predicting the magnetic field evolution over the next 50 years after assimilating the observations from 1840 to 2020. Figure 9 shows the assimilation up to 2020 and the consecutive prediction for the axial dipole g 0 1 and coefficient h 3 5 . The axial dipole SV in 2020, 6.8 nT/year, is close to that suggested by CHAOS and Kalmag, but significantly larger than the IGRF-13 value of 5.7 nT/year. We predict Fig. 7 a Radial SV in 2020 of our candidate model and b the difference to IGRF-13 at the Earth's surface for ℓ ≤ 8 Fig. 8 Mauersberger-Lowes power spectra of the differences between each of the other 13 SV model candidates and the IGRF-13 (gray curves) as well as the dispersion of all the candidate models around the latter (dashed black curve). The difference between our candidate model and IGRF-13 and its corresponding uncertainty given by the spread of the ensemble (red solid and dashed curves, respectively), are also shown that the value will increase further until 2028, staying roughly constant around 9.2 nT/year until 2039, dropping and staying around 8 nT/year until 2070. Using the instantaneous SV from 2020 for a linear prediction suggests an axial dipole for 2070 that is only less than 1% smaller than for the dynamic prediction. The reason is the generally low value of the axial dipole SV. We conclude that the axial dipole strength will continue to decrease in a roughly linear fashion at a rate of 800 nT per century.
The prediction for the coefficient h 3 5 is much more dubious, as is shown in Fig. 9b and d. In 2020, the respective SV is less than half the absolute value suggested by CHAOS or Kalmag. The IGRF-13 lies somewhere in between. Our model predicts that the SV will increase significantly from −0.4nT/year in 2020 to about 2.5nT/ year in 2047 and will then remain roughly constant. The amplitude of the h 3 5 coefficient will decrease enormously from −122 nT in 2020 to only −22 nT in 2070. The linear prediction, on the other hand, suggests a drastically different value of −140 nT for 2070. We have already discussed above that h 3 5 and its SV vary on a rather short time scales from 2010 to 2020. Should this behavior continue, neither of the two predictions is trustworthy.
The evolution of the South Atlantic Anomaly is of particular practical interest. During the past century, the SAA has decreased in intensity, drifted towards the west and has grown in surface area (Finlay et al. 2010). Fig. 10 illustrates the evolution of the SAA from 1970 to 2070, using observations until 2020 and our full dynamical forecast thereafter. Contour lines show the field in 2020 from our ensemble mean, which is virtually identical to CHAOS or Kalmag. The crosses mark the positions of the two intensity minima within the SAA in a 10-year cadence. For the predictions, the crosses consist of error bars that show the location uncertainty in latitude and longitude. Currently, the deeper western minimum is located in northern Argentina. The eastern minimum, which appeared around 2010, presently sits at 42 • S and 0 • E.
We predict that the western minimum continues the westward drift it has been showing over the last decades with a rate of 0.19 • /year. After a quick move to the northeast in 2020, the eastern minimum also begins a westward drift, but at a much lower rate of 0.07 • /year. The SAA will therefore stretch in longitude. The reason for the drift is the predominantly westward flow in this region. According to our forecast, the field intensity of the SAA will decrease further. For the western minimum, we predict a rate of 16.3 nT/year, which is significantly slower than the observed value of 34.0 nT/year over the last century. Our predicted rate for the eastern minimum, on the other hand, is 42.6 nT/year. By 2070, the intensity of the eastern minimum will reach 21.2 ± 1.0 µ T and will be lower than that of its western counterpart.
Identifying the cause of the SAA structure in terms of field at the CMB and its evolution is no trivial matter. Terra-Nova et al. (2019) shows that intensity anomalies at the surface such as the SAA are due to equatorial asymmetries of both reverse as well normal magnetic flux patches at the CMB. A sequence of snapshots of flow and field at the CMB and surface intensity during our longterm forecast are available at the Appendix (Figs. 13, 14  and 15). These suggest that the relatively high-intensity decrease rate of the eastern minimum is related to the system of upwellings located beneath Southern Africa, causing the reversed flux patches and low-latitude normal flux patches in that region to diverge horizontally. Underneath the western minimum, however, the overall westward flow advects the normal and reverse flux patches alike. Figure 11 shows the observed and predicted positions of the geomagnetic poles and the magnetic dip poles. The former relates to the poles of the dipole component. The dip poles, on the other hand, are points where the total field is vertical. According to our predictions, the geomagnetic poles will continue to approach the geographic poles, which is mainly a consequence of the decreasing equatorial dipole coefficient h 1 1 . The northern magnetic dip pole has crossed to the Eastern Hemisphere, after a period of acceleration over the past 20 years (see Thébault et al. 2015) and will continue rapidly moving towards Siberia over the next 50 years. Meanwhile, the southern dip pole follows a slower path towards the northwest. Over the next decades, the dip poles will approach similar longitudes and resemble an eccentric dipole that is offset from Earth's center.

Discussion
We have used an EnKF to predict the geomagnetic field evolution based on the assimilation of the COV-OBS. x1 and Kalmag geomagnetic field models with a full 3D-numerical dynamo model. The mild parameters of the dynamo model allowed for the performing of various tests with different ensemble sizes and configurations. Since the ensembles are not large enough to sufficiently constrain the N x ∼ O(10 6 ) elements in the dynamo model covariance, the assimilation can become unstable Fig. 10 Evolution of the South Atlantic Anomaly (SAA) over the period 1970-2070 calculated from the ensemble assimilation (note that after 2020 only the forecast is shown, while before only the analyses are displayed). The background contour lines represent the field intensity based on the mean analysis in 2020 and the symbols correspond to the position and intensity (color scale) of the two minima every 10 years. The error bars correspond to the 1 σ uncertainty in latitude and longitude, while the symbol size scales with the uncertainty in the intensity of each minimum point in time when confronted with observations bearing small uncertainties. S2019 showed that correlations between model coefficients of the same wave number clearly dominate in the model covariance. Reducing the covariance matrix accordingly by setting all other elements to zero, a procedure they called spectral covariance localization, stabilized their assimilation scheme. Building on their method, we went a step further and in addition only kept covariance matrix entries that reflect the dominant equatorial symmetries (see Fig 1). This allowed a strong decrease of the ensemble size (from N e = 512 down to 128, see Fig. 3) without considerable loss of forecast quality.
The surface core flow resulting from the assimilation of COV-OBS.x1 and Kalmag (Fig. 6) bears important differences with respect to previous studies. For instance, the westward flow at low latitudes underneath the Atlantic shows considerable variations in the latitudinal component instead of a more coherent meandering flow band as in Barrois et al. (2018). However, the most prominent discrepancies of our estimation reside on the strong eastward flow underneath western Asia and the velocity at high latitudes, notably inside the TC (e.g., Aubert 2020). The former might be the cause of the peak 15 nT/ year difference between our candidate model and IGRF-13 SV underneath Asia (see Fig. 7b). The reason behind such differences in flow pattern has not been explored in depth yet, but likely resides on the characteristics of the background dynamo model. The dominant flow wave number and degree of equatorial symmetry are some characteristics of the background dynamo model which are probably introducing important errors in the flow estimation and consequently in the predictions.
Using a more advanced dynamo model is certainly desirable for improving assimilation-based forecasts. Decreasing the Ekman number (Schaeffer et al. 2017b) or imposing boundary conditions that reflect the lower mantle thermal structure ) come to mind. The higher numerical cost would require a further reduction of the ensemble size, which would in turn call for additional measures to stabilize the assimilation scheme. One possible measure is covariance inflation (Anderson and Anderson 1999). Inflation is implemented by multiplying the covariance matrix by a factor slightly larger than one, which can prevent the scheme from putting too much confidence in the model forecasts. We thus envision assimilations at more realistic dynamo setups, but somewhat smaller ensembles which retain the spectral localization discussed here and additionally employ inflation.
However, we stress that the results of our assimilation already provide realistic secular variation patterns (see Fig 6). Hindcast tests show that predictions based on instantaneous SV values slightly outperform forecasts using the dynamical model. This can be traced to problems in the SA, which depends on the flow itself and also on the flow changes. The force balance that determines the flow in the Navier-Stokes equation is therefore of crucial importance. Aubert (2020) showed that enforcing realistic force balance constraints during the analysis can improve the dynamic predictions. However, steady flow forecasts where only the induction equation is integrated forward in time while the flow field is kept constant, still provide the best forecasts (Aubert 2020).
Using steady flow is seen to deteriorate our dynamic hindcasts before 2010, but improve them after that (see Fig. 5). The non-systematic and relatively small differences between those, of around 10 nT, can mean that errors in flow estimation are comparable to the impact of errors from flow variations in our assimilation. Nonetheless, the dynamic hindcasts are observed to improve with time, suggesting that the dynamic assimilation scheme might be still converging to the forecast quality level of the linear extrapolations.
Since the linear forecasts typically provide the best 5-year hindcasts, we decided to propose the SV in 2020 as a candidate model for the IGRF-13. Our SV model is close to the final IGRF-13, with a RMS difference at the Earth's surface of about 4 nT/year (see Fig. 7). Table 1 of the Appendix lists the Gauss coefficients from our SV model candidate. Most of the other candidates also provide the SV for 2020. However, the intention of the IGRF is to linearly describe the evolution over the next 5 years. It would thus be more appropriate to use the mean SV over this period. We list the mean SV suggested by our dynamic forecast in Table 2. The values agree with the SV in 2020 within the uncertainty. The RMS difference of the average SV to the IGRF-13 is 15 nT/year and lies within the dispersion of the different candidate models (see Fig 8). Only time will tell whether this mean SV provides a better forecast for 2025.
The Earth's magnetic field is thought to bear a rather long predictability window, with an e-folding time of about 30 years (Lhuillier et al. 2011a). Prediction analyses have revealed that this window might lie somewhere between half and a full century Aubert 2015). We choose here to focus on a prediction window of 50 years using the results from our dynamic forecast from the ensemble of dynamo models. Our predictions for the axial dipole and the SAA over the next 50 years are in general agreement with the forecasts by Aubert (2015), but there are also differences. The predicted trajectories of the SAA western minimum from our ensemble point to a westward drift rate of about 0.20 • /year and a decrease of the local field intensity by about 1.0 µ T in the next 50 years. For comparison, the corresponding values from Aubert (2015) are 0.26 • /year and 1.46 µ T, respectively. Additionally, our forecast predicts that the SAA eastern minimum, which drifts towards the West at a rather slow rate of 0.07 • /year and undergoes an intensity field decrease of 2.9 µ T, will become the main minimum of the anomaly in the next 50 years (see Fig. 10).
Our predictions for the magnetic north pole correspond to a southward movement of about 472 km towards Siberia. This value is considerably lower than the frozen-flux advection-based forecast of 660 km by Livermore et al. (2020). The authors propose that most of the recent migration of the pole towards Siberia is related to the elongation of the Canadian magnetic flux lobe, which has been continuously stretched by a northwestward flow underneath that region. The differences in the pole's evolution hence likely lies on the structure of the high latitude flow between the two models.
It is important to note that these predictions do not contemplate the existence of geomagnetic jerks. Geomagnetic jerks are virtually abrupt changes in the second time derivative of the magnetic field. Recent dynamo simulations suggest that they happen where a special kind of fast magnetic waves reach the CMB (Aubert and Finlay 2019). The waves only seem to appear in small Ekman number simulations and are therefore not part of the dynamics in our dynamo model. Our scheme should fail to provide decent magnetic field prediction over epochs where the jerks play an important role in the geomagnetic field evolution.

Conclusion
In this study, we used a geomagnetic sequential data assimilation scheme based on a full 3D-numerical dynamo model for forecasting the Earth's main magnetic field. The scheme adopts an EnKF approach that assimilates the main magnetic field from COV-OBS.x1 (Gillet et al. 2015) and Kalmag (Baerenzung et al. 2020) models from 1840 to 2020. A new spectral covariance localization method, extending the study by S2019, stabilized the assimilation scheme for moderate and small ensemble sizes. We chose to use an ensemble of moderate size ( N e = 256 ) to contribute to the 2020 SV prediction for the IGRF-13 as well as to predict the field evolution over the next 50 years.
Hindcasts tests focusing on the period from the year 2000 to 2020 showed that the scheme provides reasonable forecasts over a 5-year period. A linear prediction (LP) that simply uses the instantaneous secular variation yields overall better results than the dynamic forecasts from the dynamo simulations, although the latter seem to steadily improve over time. We therefore decided to propose the LP in 2020 from our assimilation scheme as an IGRF-13 SV candidate. For our forecast of the field evolution until 2070, however, we used only the prediction based on the full dynamics.
We predict that the axial dipole will keep on decreasing at a rate that will slightly increase over the next half-century. The overall SAA field intensity will drop by another 10%, continue its drift further westward and increase its longitudinal extent. The magnetic dip poles predictions show that the magnetic north pole continues rapidly drifting towards Siberia, although at a decreasing rate over the next 50 years. The dip poles hence appear to approach a configuration close to an eccentric dipole. However, these predictions should be taken with care, since a careful analysis of model errors has not yet been taken into account. In particular, the estimated flows showed important differences with respect to previous studies. Their impact on the predictability of the Earth's magnetic field will be the subject of forthcoming studies.
Abbreviations CMB: Core-mantle boundary; RMS: Root mean squared; SA: Secular acceleration; SAA: South Atlantic Anomaly; SH: Spherical harmonic; SV: Secular variation. Figure 12 shows a subset of the initial model covariance corresponding to correlations between the poloidal magnetic field potential and the poloidal flow potential, as well as the corresponding subset of the localization including the checkerboard pattern. Figures 13, 14 and 15 show snapshots of the flow and field evolution over the 50 years discussed in this study. The fields correspond to mean values from the ensemble of dynamo models during the forecast. Tables 1 and 2 list the Gauss coefficients corresponding to the SV predictions for 2020 (submitted as a candidate model for the IGRF-13) and from the dynamical trajectory of the ensemble averaged over the period 2020-2025, respectively. Note that since the latter correspond to a forecast incorporating the nonlinear dynamics of the model, the uncertainties are considerably greater than the instantaneous SV. Fig. 13 Maps of the flow just below the CMB (top), radial magnetic field at the CMB (middle) and field intensity at the Earth's surface (bottom). The horizontal flow is shown with streamlines, while the radial field is displayed by the color scale. The maps are based on the assimilation analysis in 2020 and are given in Mollweide projections