A candidate secular variation model for IGRF-13 based on MHD dynamo simulation and 4DEnVar data assimilation

We have submitted a secular variation (SV) candidate model for the thirteenth generation of International Geomagnetic Reference Field model (IGRF-13) using a data assimilation scheme and a magnetohydrodynamic (MHD) dynamo simulation code. This is the first contribution to the IGRF community from research groups in Japan. A geomagnetic field model derived from magnetic observatory hourly means, and CHAMP and Swarm-A satellite data, has been used as input data to the assimilation scheme. We adopt an ensemble-based assimilation scheme, called four-dimensional ensemble-based variational method (4DEnVar), which linearizes outputs of MHD dynamo simulation with respect to the deviation from a dynamo state vector at an initial condition. The data vector for the assimilation consists of the poloidal scalar potential of the geomagnetic field at the core surface and flow velocity field slightly below the core surface. Dimensionless time of numerical geodynamo is adjusted to the actual time by comparison of secular variation time scales. For SV prediction, we first generate an ensemble of dynamo simulation results from a free dynamo run. We then assimilate the ensemble to the data with a 10-year assimilation window through iterations, and finally forecast future SV by the weighted sum of the future extension parts of the ensemble members. Hindcast of the method for the assimilation window from 2004.50 to 2014.25 confirms that the linear approximation holds for 10-year assimilation window with our iterative ensemble renewal method. We demonstrate that the forecast performance of our data assimilation and forecast scheme is comparable with that of IGRF-12 by comparing data misfits 4.5 years after the release epoch. For estimation of our IGRF-13SV candidate model, we set assimilation window from 2009.50 to 2019.50. We generate our final SV candidate model by linear fitting for the weighted sum of the ensemble MHD dynamo simulation members from 2019.50 to 2025.00. We derive errors of our SV candidate model by one standard deviation of SV histograms based on all the ensemble members.


Introduction
The secular variation (SV) of the geomagnetic field is caused by a highly nonlinear process in the outer core of the Earth. International Association of Geomagnetism and Aeronomy (IAGA) releases a set of geomagnetic main field model and linear prediction of SV field as International Geomagnetic Reference Field (IGRF) model every 5 years, to keep the model's deviation from the actual field small. The previous generation, the 12th generation of IGRF, was released in 2015 (Thébault et al. 2015a). Each generation of IGRF has been produced by the IAGA Working Group V-MOD by integrating independent contributions from research groups all over the world, and consists of three types of models: the Definitive IGRF (DGRF) field 5 years before the release epoch, the geomagnetic main field for the release epoch, and the SV model valid for the 5 years on and after the release epoch. The truncation of spherical harmonics for IGRF and DGRF is now L MF = 13 , while for SV is L SV = 8 . Evaluations of all the submitted candidate models are unveiled altogether with the release of IGRF (e.g., Thébault et al. 2015b).
Research groups in Japan have never contributed to any IGRF models so far, although there have been important contributions of highly quality-controlled observatory data, e.g., Kakioka, operated by Japan Meteorological Agency. Since both established techniques and expertise for the fields of magnetohydrodynamic (MHD) dynamo simulation and data assimilation are present in Japan, we decided to contribute to the IGRF community by an SV candidate model first. This was actually enabled by a bilateral collaboration between Japan and the geomagnetic group in Institut de Physique du Globe de Paris (IPGP), France; this group brought to the project the required expertise to build, from raw data, a geomagnetic field model, which, in turn, has been used as input data to our assimilation scheme.
SV candidate models for IGRF have been commonly presented as linear or polynomial extrapolation of known latest SVs. The GFZ Potsdam team, Germany, for example, generated an SV candidate model for IGRF-12, namely SV from 2015.0 to 2020.0, by averaging the SV of their parent model from 2013.5 to 2014.5 (Lesur et al. 2015). Abrupt change in SV, such as geomagnetic jerks, however, may deteriorate the prediction of shorttime SVs. This suggests necessity of prediction reflecting dynamic process in the Earth's outer core, e.g., through numerical geodynamo models. There have been limited attempts to use geodynamo models for contributions to the SV models of IGRF. Kuang et al. (2010) was the first to contribute to IGRF-SV models using a geodynamobased data assimilation scheme. They adopted a Kalman filter approach and used over 7000-year data with a prediction-correction algorithm to reduce the systematic error of the geodynamo model. In addition to the full MHD prediction, the number of attempts to use controlled dynamic processes has been increasing. Fournier et al. (2015) provided an SV model for IGRF-12, using a geodynamo model (Aubert et al. 2013) with a steady fluid velocity in the outer core.
Variational data assimilation approaches also have a high possibility for providing good forecasts of the geomagnetic field, although there has been a very limited number of their applications to the geomagnetic field. The variational approaches simultaneously adjust a dynamic model to all the observations available in the assimilation window, while the sequential methods, such as the Kalman filter approach, update (or correct) the model output with given observations at each time step (see Talagrand (1997) or Fournier et al. (2010) for introduction to both). One of big advantages of the variational approaches is that the resulting trajectory of the state vector keeps the dynamic consistency of the system over the assimilation window. Fournier et al. (2007) show a possibility that the variational approach can improve historical geomagnetic data with the recent dense satellite magnetic data, from experiments using a one-dimensional MHD dynamo toy model. Li et al. (2011;2014) have developed variational assimilation schemes for more realistic threedimensional MHD systems. In particular, Li et al. (2014) performed experiments of variational data assimilation using an inertia-free MHD dynamo model with synthetic data sets, which shows a potential to handle actual magnetic observations in the variational assimilation scheme. For their application to real magnetic data, however, difficulties arise from the requirement of an adjoint model of a fully nonlinear MHD dynamo model and from computational cost for construction of the derivative of the cost function with respect to an initial state, i.e., for backward propagation of the error information (e.g., Fournier et al. 2010).
In this study, we adopt an iterative ensemble-based variational scheme based on the four-dimensional ensemble-based variational method (4DEnVar; Liu et al. 2008) with a fully nonlinear three-dimensional MHD dynamo model (Takahashi 2012(Takahashi , 2014 to predict SV in the coming 5 years from 2020.0. A noticeable feature of 4DEn-Var is that it linearizes the dynamo output variables at a given time comparable to observations with respect to the deviation of the state vector from a mean state at the initial time, assuming weak nonlinearity of the model over the assimilation window. The 4DEnVar approach, thereby, does not need an adjoint model of the forward model and numerical cost for the backward propagation of the error information, even though it is named "variational method". This strategy seems reasonable when both assimilation window and forecast period are shorter than the timescale of nonlinearity of the system. We illustrate that this approach is feasible with our iterative algorithm based on 4DEnVar (Nakano 2020) for generation of a candidate SV model for IGRF, with numerical experiments using the real geomagnetic data from 2004.50 to 2014.25 and with comparison of the forecast performances with that of IGRF-12.
The remainder of this paper is organized as follows: we first explain our assimilation method ("Method" section). Secondly, we report results of numerical experiments using the past datasets ("Numerical experiments" section), and then describe details of estimation of our candidate model for IGRF-13 SV ("The SV candidate model for IGRF-13" section). After the explanation of our error estimation ("Error estimation of the candidate SV model" section), we summarize this study ("Discussion and conclusions" section).

Method
Our intent is to forecast SV by performing the 4DEnVar data assimilation using a numerical geodynamo code and existing data with respect to the poloidal scalar potential of the geomagnetic field at the core-mantle boundary (CMB) and the toroidal and poloidal components of the core surface flow. From the variational data assimilation experiments with an inertia-free MHD dynamo model, Li et al. (2014) found that the magnetic field below the CMB is hard to reconstruct only from the magnetic data taken outside the core due to the diffusion-dominant Ekman boundary layer. In agreement with their proposal, we use the preliminary estimated core surface flow as part of observation vectors in our data assimilation as well as the geomagnetic data. Inclusion of the core surface flow in data vectors results in indirect inclusion of SV data in the data assimilation.
In this section, we first describe the assimilation theory ("Data assimilation theory" section) and details of our dynamo model and how to adjust dimensionless time to the actual time ("Geodynamo simulation: parameters and scaling of time" section). We then briefly discuss the nonlinearity of our dynamo model ("Nonlinearity of the numerical geodynamo: error growth e-folding time" section). Next, we explain details of preparation of the observational data and how to convert dimensionless simulation outputs to variables comparable to real data ("Data 1: poloidal scalar potential at the CMB obtained from the MCM model" section and "Data 2: core surface flow" section) and finally describe the way of practical implementation of the data assimilation ("Implementation of assimilation" section).

Data assimilation theory
We consider the minimization of the following cost function: where x k is the state vector of a dynamo model at time t k , y k denotes the observation vector, R k is the covariance matrix of observation noise, and h k is an observation operator which converts a state vector x k to observable variables for the comparison with y k . Given the dynamo model, x k is uniquely determined from the initial state x 0 . This allows us to represent x k as a function of x 0 , that is, x k = f k (x 0 ) . Defining a function g k as (1) , the cost function in Eq. (1) can be rewritten as follows: The minimization of this cost function is achieved by an iterative algorithm based on the 4DEnVar method (Liu et al. 2008). At the m th iteration, we approximate the cost function by using an ensemble of the simulation outputs x 0,m , which are prepared so that the ensemble mean, Σ N n=1 x (n) 0,m /N , is equal to the m th estimate x 0,m . At the m th iteration, we seek x 0 that minimizes Eq. (2), which turns out to be x 0,m+1 , with given x 0,m and x (1) 0:K ,m . As an important first step of the 4DEnVar, we express g k (x 0 ) in terms of the first-order Taylor expansion, where G k is the Jacobian of g k at x 0,m . We then approximate x 0 as a weighted sum of the ensemble members. Now we define the following matrices ⌣ X 0,m and Γ k,m for convenience: This allows us to write x 0 =x 0,m + ⌣ X 0,m w , where w consists of weight for each ensemble member. x 0,m is the mean of x (n) 0,m (n = 1, . . . , N ) . Using Eqs. (4) and (5), the function g k (x 0 ) in Eq. (3) can then be expressed: Note here that the Jacobian G k disappears in the expression of g k (x 0 ) with the aid of the relationship derived from Eq.
For the first requirement, we discuss the nonlinearity of our dynamo model using the error growth rate  in "Nonlinearity of the numerical geodynamo" section later, while we see that ensembles shrinking through iterations meet the second requirement in "Numerical experiments" section. From Eqs. (2) and (6), we introduce the following objective function; where σ m is a parameter, which is fixed to σ m = 1 in this study, while we decrease elements of R k at each step. This cost function is minimized provided that: The (m + 1) th estimate x 0,m+1 is then obtained as and we proceed to the next iteration. The first term of the right-hand side in Eq. (7) is added to ensure robustness. This iterative application of Eq. (8), which is similar to the iterative ensemble Kalman smoother algorithm (Gu and Oliver 2007;Bocquet and Sakov 2013), minimizes Eq. (2) in the subspace spanned by the ensemble members (Nakano 2020). After obtaining x 0,m+1 it is necessary to perform MHD dynamo simulations with a set of initial conditions x (1) 0,m+1 , . . . , x (N ) 0,m+1 to renew the ensemble members for the ( m + 1)th iteration. See Appendix A for how to prepare the set of initial conditions from x 0,m+1 , ⌣ X 0,m and ŵ m .
At the final (5th is chosen in this study) iteration, we also estimate the bias and trend components which correspond to model error in the dynamo model, by minimizing the following function: where b denotes the bias component, while a is the coefficient for the trend component. The bias and trend terms correspond to the offset and the linear departure in time between the observation and the model, respectively. Here, we assume that the observation is mostly explained by the dynamo model output and that the bias and trend components are minor. We thus select P a and P b as: Large norms of P −1 b and P −1 a suppress intensities of b and a while minimizing Eq. (10). w, b, a that minimize Ĵ m (w, b, a) in Eq. (10) can be obtained in a similar manner to Eq. (8) (see Appendix B for details). The minimization of Eq. (10) gives the approximate minimum of the following cost function: The final estimate and prediction are obtained by the following equation: where M indicates the final step, i.e., M = 5 , and ŵ M , b and â are solutions to Eq. (10). Then we can obtain the final estimate of Eq. (13) by the sum of a single MHD simulation starting from x 0,M , the weighted sum of the Mth ensemble members, and the trend and bias terms. Note that we can use Eq. (13) not only for the final estimate within the assimilation window, but also for the estimate in the forecast period outside the assimilation window when the future extensions of g k x 0,M and g k x (n) 0,M , (K < k) , are available, which requires only additional dynamo runs for g k x 0,M and all the ensemble members. Figure 1 shows how to prepare an SV model for IGRF-13 by our assimilation scheme and Eq. (13), where the future extensions of the ensemble members (the gray area after "Release of IGRF") generate (10) (13) g k,M = g k x 0,M + ⌣ Γ k,MŵM +b + kâ, a future prediction (the blue line after "Release of IGRF") by the weighted sum of the ensemble members via Eq. (13).
Note that the first-order approximation of the function f k allows us to approximate the state at an arbitrary time t k ′ ( 0 ≤ k ′ ≤ K ) in a similar manner to Eq. (6): where F k is the Jacobian of f k at x 0,m and ⌣ X k ′ ,m is a matrix obtained by replacing 0 by k ′ in Eq. (4) as follows: The state at t k ′ can be estimated from the following equation: without any modification to ŵ m in Eq. (8), which is the same as used in Eq. (9). To obtain x k ′ ,m+1 by Eq. (16) means that we minimize the cost function Eq. (2) (or (12)) using the first-order Taylor expansion of g �k (x k ′ ) at x k ′ ,m , instead of Eq. (3), where k = k − k ′ , and g k and G k are associated with the inverse function of f − k when �k < 0 . This enables us to choose arbitrary time, t k ′ , at which we restart MHD dynamo simulation after the data assimilation, for the purpose of better accuracy of the linear approximation at times close to t k ′ . For the future forecast, single MHD or kinematic dynamo (KD; meaning steady flow) simulation can also estimate future magnetic fields, running them from the initial condition optimized at arbitrary time t k ′ (Eq. 16) by data assimilation. We discuss differences in forecasts by MHD and kinematic dynamo simulation from different t ′ k in "Numerical experiments" section.

Geodynamo simulation: parameters and scaling of time
We adopt the geodynamo simulation code established by Takahashi (2012Takahashi ( , 2014 for our data assimilation and forecast processes. In this code, the velocity and magnetic fields are decomposed into the poloidal and toroidal parts and spatially expanded in terms of spherical harmonics in the tangential directions, while combined compact finite differencing and the Crank-Nicolson scheme are adopted for discretization in the radial and time coordinates, respectively. The temperature and composition are treated separately (Takahashi 2014). We select the following dimensionless parameters for the geodynamo simulation: the Ekman number is E = 3 × 10 −5 , the magnetic Prandtl number is Pm = 2.0 , the Prandtl number is 0.1, the compositional Prandtl (or Schmidt) number is 1.0, the modified thermal Rayleigh number is 500, and the modified compositional Rayleigh number is 2000. The magnetic Reynolds number ( R m ) of this dynamo is 181 with the standard deviation of 15, which is evaluated from the initial ensemble members for our data assimilation (see "Numerical experiments" section for details of the ensemble).
Since all the parameters and variables are dimensionless in the numerical geodynamo model, the simulation outputs should be scaled for comparison with the actual data Fig. 1 A schematic figure of the 4DEnVar and prediction of the future geomagnetic field. The weighted sum of an ensemble of dynamo run (the thick blue line) is the final estimate after the assimilation corresponding to Eq. 13. The SV candidate model for IGRF-13SV (the red dashed line) is determined by line-fitting to the MHD expectation of the next 5 years. One can use geomagnetic data up to ~ 6 months before the release of IGRF for generation of candidate models because the deadline of their submissions is ~ 3 months before the release (October 1st for  of the geomagnetic main field and the core surface flow. A time scale of the geomagnetic secular variation, τ SV , can be characterized by simple dependence of spherical harmonic degree, l (Lhuillier et al. 2011a). Christensen and Tilgner (2004) obtained a value of τ SV = 535years by fitting τ SV /l ( l ≥ 2 ) to the time scale of secular variation for each spherical harmonic degree, τ l (see Christensen and Tilgner 2004, for the definition). In the same manner, we estimate the secular variation time scale of our numerical geodynamo model to be τ * SV = 0.0411 , scaled by the viscous dissipation time; we first estimate the mean τ l of our numerical geodynamo from a longtime dynamo run (about one magnetic diffusion time, i.e., 26,000 years in actual time) and then evaluate τ * SV by fitting τ * SV /l to the estimated τ l for l ≥ 2 . Figure 2 shows the estimated τ l and the fitting result. Hence, we convert the dimensionless time of the numerical geodynamo into the actual time by using τ ν ≡ τ SV /τ * SV = 13, 049 years throughout this study.

Nonlinearity of the numerical geodynamo: error growth e-folding time
Since the 4DEnVAR approach relies heavily on the linear assumption Eq. (3) and resulting Eq. (6), we should be careful about its validity in the problem in concern, i.e., the data assimilation and forecast of the geomagnetic field for IGRF-13SV. We here investigate the timescale of nonlinearity of our numerical geodynamo by introducing the error growth rate (e.g., Hulot et al. 2010), which measures the rate at which small error in a dynamo state vector grows with time. It is known that the error growth rate is stable once a dynamo model and parameters are provided, independent of the magnitude of inserted error and types of variables into which error is contaminated (Lhuillier et al. 2011b). Its inverse, i.e., the e-folding time of error growth τ e , therefore, can be interpreted as a limit of the period over which data assimilation is feasible with the nonlinear dynamo model; when the length of assimilation window is comparable with τ e , small error in the initial state is critical to the state at the end of the window, and one cannot expect to obtain reasonable initial condition. The linear assumption for 4DEnVar (Eq. 3) clearly requires that the assimilation window should be shorter than τ e of our numerical geodynamo model.
We estimate the e-folding time τ e of our numerical geodynamo by following the method in Lhuillier et al. (2011b). We first choose ten different dynamo state vectors from a dynamo run for about one magnetic diffusion time. We then add a small perturbation to the poloidal components of the magnetic field over the first eight harmonic degree by the expression: where S m l and S m l are the poloidal scalar functions of the magnetic field with degree l and order m of spherical harmonics for the perturbed and original initial conditions at a chosen time t p , respectively. α m l is a coefficient obeying the standard normal law and ϵ is the relative magnitude of perturbation. We choose ϵ = 10 −10 in this study. MHD dynamo simulations are performed from both original and perturbed initial conditions and the discrepancy between them are measured by the magnetic error at the CMB, where t is time scaled by the magnetic diffusion time τ η , r CMB is the radial coordinate of the CMB, and �·� denotes the time average. For each pair of dynamo runs, the error growth rate and its inverse, the e-folding time τ e = −1 , are estimated by fitting the model of t + a , where a is the offset, jointly to log �y l (t) ( 1 ≤ l ≤ 8 ). The time window for the fitting is manually determined, by simultaneously considering the first eight �y l (t) . After the same procedure for ten chosen times of t p , we average ten s and τ e s for our final estimates. Figure 3 shows results of estimation of and τ e for our geodynamo model. As shown in Fig. 3a, the error growth rate (slope of the blue line) was obtained by the linear regression over the first eight spherical harmonic degrees (the gray lines). As in Fig. 3b, we finally evaluate = 187 ± 36 [ 1/τ η ] and τ e = (5.37 ± 1.02) × 10 −3 τ η , where the uncertainties are given by two standard deviations. Using the scaling law, τ η = P m τ ν = 2 · 13, 049years , τ e of our numerical dynamo is found to be τ e = 140 ± 27 years. Then we should choose the assimilation window Fig. 2 Secular variation time scale of our dynamo model. The time scale of secular variation for each spherical harmonic degree τ l is fitted by τ * SV /l for l ≥ 2 . The resulting τ * SV is 0.0411, which was utilized to scale the dimensionless numerical dynamo time shorter than ~ 110 years to avoid the effect of nonlinearity of our geodynamo model (Takahashi 2012(Takahashi , 2014. We later select a maximum assimilation window of 10 years and a prediction period of 5 years, which seem short enough to circumvent the nonlinearity effect and to hold Eq. (3). Hulot et al. (2010) claimed that τ e of the Earth's dynamo system appears to be about 30 years from the linear relationship between τ e /τ SV and R m , provided that R m > 100 , E < 10 −9 , and τ SV = 535years (Christensen and Tilgner 2004) for the Earth. Although our τ e (∼ 140 years) is much longer than τ e (∼ 30 years) by Hulot et al. (2010), it is found that τ e highly depends on numerical dynamo models even under the similar dynamo parameter settings (e.g., see Fig. 3 in Hulot et al. (2010) for the comparison between τ e by Hulot et al. (2010) and that by Olson et al. (2009) under E = 3 × 10 −4 ). We, therefore, think our long τ e is a characteristic of our geodynamo model. We here conclude that our assimilation window of maximum 10 years is short enough to adopt linear approximation of Eq. (3) and resulting Eq. (6) from the perspective of τ e of our numerical geodynamo.

Data 1: poloidal scalar potential at the CMB obtained from the MCM model
We use a version of the MCM model (Ropp et al. 2020) which spans from 2001.5 to 2019.25, referred to as the MCM3 model in this paper hereafter, to construct data vectors for our data assimilation. The MCM3 model has been built using hourly mean observatory data together with CHAMP and Swarm-A satellite magnetic data. These data have been selected depending on their local time and for magnetically quiet periods (c.f., Lesur et al. 2008). The model is made of a series of snapshot models, 3 months apart, that have been built through a Kalman filter approach combined with a correlation-based modeling technique for the analysis step (Holschneider et al. 2016). Each snapshot model is parameterized in terms of spherical harmonics and includes static core field and secular variation contributions. Lithospheric, external and induced fields are also co-estimated as well as crustal offsets at each observatory. Details of the modeling technique are described in Ropp et al. (2020). The Gauss coefficients of the core field and its secular variation can be downward continued to the CMB for use, e.g., for information on the flow in the liquid outer core.
Assuming the mantle an electrical insulator, we convert the Gauss coefficients, g m l (t) and h m l (t) , given by the MCM3 model to the poloidal scalar potential, S mc l (t) and S ms l (t) , at the CMB, respectively, as where t is the dimensional time, r e = 6371.2km and r oc = 3485km are the radii of the Earth and outer core, respectively, and S mc l (t) and S ms l (t) are poloidal scalar functions for cosine and sine terms, respectively. One of data to be assimilated is S m l (t) which denotes S mc l (t) and/ or S ms l (t) , and the truncation of spherical harmonics is L MCM = 14 . In the assimilation procedure, however, we use only coefficients with degree up to 13 assigning very small weight to coefficients for l = 14 . The magnetic data are constructed in the form of (20)  (21) to generate vectors comparable with the magnetic data vectors. We prepare the data vector by a 0.25-year interval. From the outputs of the numerical geodynamo, the poloidal scalar potentials are normalized according to Eq. (21) after the conversion of the dimensionless time to the actual time.

Data 2: core surface flow
We also include the core fluid velocity field slightly below the CMB in our data assimilation. In the present numerical simulation of geodynamo, the no-slip condition for the velocity field is imposed at the CMB (Takahashi 2012). It, therefore, is pertinent to adopt core surface flows estimated with the effect of viscous boundary layer at the CMB (Fig. 4a). Hence, we employ a method of Matsushima (2015). Before recalling the method, we define the velocity field, V , and the magnetic field, B , in the outer core. V can be expressed in terms of the poloidal and toroidal scalar functions as where (r, θ , φ) is the spherical coordinates, r is a radial unit vector, and P m l is Schmidt quasi-normalized associated Legendre function with degree l and order m . Near the core surface, the radial component of the velocity field, V r , is likely to be much smaller than the horizontal component, V h . Therefore, in core surface flow modeling, only the horizontal component, V h =θV θ +φV φ , is computed: where θ and φ are unit vectors in the θ -and φ-directions, respectively. Hence, Matsushima (2015) computed where r o is the outer radius of the rotating spherical shell in numerical MHD dynamo model, and ξ * 2 is a depth from r = r o . Therefore, taking into account the length scale r oc − r ic and the velocity scale (r oc − r ic )/τ ν (see Appendix C), we have and In the method of Matsushima (2015), (0) At the CMB (r = r oc ) , geomagnetic secular variations are caused only by magnetic diffusion due to the no-slip boundary condition; where η is the magnetic diffusivity. (1) Inside the viscous boundary layer (r = r b1 = r oc − ξ 1 ) , the viscous force plays an important role in force balance as where is the rotation rate of the mantle, ρ is the density of the core fluid, J is the electric current density, and ν is the kinematic viscosity. The magnetic diffusion, motional induction, and advection are assumed to contribute to temporal variations of the magnetic field;   Table 1 for details). Arrows denote the horizontal flow velocity field, while color contours show upwellings and downwellings given by horizontal divergence of the velocity field. The poloidal and toroidal scalar potentials of velocity are compared at r b2 in our data assimilation where core fluid is assumed to be incompressible, ∇ · V = 0.
(2) Below the boundary layer (r = r b2 = r oc − ξ 2 ) , core flow is assumed to be tangentially magnetostrophic (Matsushima 2020) as Temporal variation of the magnetic field is assumed to be caused by the motional induction and advection, and the magnetic diffusion is neglected, as in the frozen-flux approximation (e.g., Holme 2015) as The radial component of the electric current density, J r , is likely to be much smaller than the horizontal component, J h , near the core surface due to the electrically insulating mantle. Following Shimizu (2006), we assume that the electric field does not contribute to the electric current density, and we have where σ is the core electrical conductivity. Hence, the tangentially magnetostrophic constraint is given as The ratio of the Lorentz to the Coriolis forces is the Elsasser number given by We can have a local value of Λ ≈ 0.62 , using actual physical parameters for the outer core, σ = 1 × 10 6 S m −1 , Ω = 7.29 × 10 −5 rad s −1 , ρ = 1.1 × 10 4 kg m −3 , and B r = 1mT . In the numerical MHD dynamo model, we obtain Λ * = 1.778 from the initial ensemble members for our data assimilation (see "Numerical experiments" section for details of the ensemble), which is approximately three times larger than Λ for the actual outer core. Such an Elsasser number larger than 1 means relative importance of the Lorentz force in our numerical geodynamo. Hence, we consider that the tangentially magnetostrophic constraint for core surface flow is appropriate. The Elsasser number Λ * = 1.778 leads to B r ∼ 0.1 mT in our numerical geodynamo model, where we adopt σ = 1.277 × 10 5 S m −1 and Ω = 4.05 × 10 −8 rad s −1 derived from τ ν = 13049 years (see Appendix C). This shows that the magnetic field in our (32) numerical geodynamo model is rather weak compared to that of the actual Earth, possibly due to the small Ω . This is the reason why we use a relative magnetic field d S m l (t) defined in (21). To overcome this problem, we should perform more realistic numerical simulation of geodynamo for a much smaller Ekman number corresponding to much larger Ω , although this remains a problem to be solved in the future.
We compare the poloidal and toroidal scalar potentials, U m l and W m l , for the core flow field with those for the numerical flow field, up to l = 14 (the same as the truncation degree of the MCM3 model, and also practically used up to degree 13 in the assimilation by controlling the weight). The core flow data, U m l and W m l , are constructed from the main field and SV at 0.25-year intervals in the MCM3 model. To make the dimensionless numerical velocity output comparable with the core flow data, we scale the simulated velocity using the relation of τ SV /τ * SV = 13, 049years and the kinematic viscosity, which is derived from Ekman number and magnetic Prandtl number adopted in the numerical dynamo. See "Appendix C" for more details of scaling of the dimensionless velocity field. The left column of Fig. 4c shows examples of core flows at r b2 = 3428.5km , which is ξ 2 = 56.5km below the Earth's core radius of r oc = 3485.0km , corresponding to a numerical grid point in the radial direction,

Implementation of assimilation
We implement the data assimilation for the following data vectors: where d S denotes the magnetic data vector composed of the scaled poloidal component specified in Eq. (21), while d U and d W represent the velocity data vectors. Subscripts U and W denote the poloidal and toroidal scalars of the velocity field, respectively. In this study, d U and d W consist of U m l and W m l at r b2 = 3428.5 km, respectively. We consider the observation error covariance matrix, R k in Eqs.
(2) and (12), as having a simple time-independent form: where R S , R U , and R W are the diagonal covariance matrices with each degree dependence. Weights for data sets are controlled by the two scalar coefficients of α S and (38) α UW . We simply adopt a time-independent form of R k = R , since the MCM3 model is a quality-controlled smooth model. We prepare R S by the expression {R S } l = C Sl m σ 2 G l=1 (l + 1) −1 , where {R S } l is the diagonal element of R S for degree l , (l + 1) −1 is the degree dependence of the variance of Gauss coefficients based on the theory by Lowes (1975), σ 2 G l=1 is the averaged actual variance of Gauss coefficient at degree l = 1 from the MCM3 model, and C S m l [·] is a function for conversion of the covariance of Gauss coefficients to that of the poloidal scalar potentials at the CMB for degree l and order m . Similarly, we prescribe R U and R W by the expressions [·] also include the conversion from the Schmidt quasi-normalization to the normalization adopted in our MHD dynamo model (Takahashi 2012(Takahashi , 2014, since we prepare y k (Eq. 38) directly comparable to the outputs of our numerical geodynamo model.

Numerical experiments
We examine the data assimilation scheme using geomagnetic data from 2004.50 to 2014.25 from the MCM3 model, covering up to 2019.25, as observations, and evaluate the resulting forecasts against data from 2014.25 to 2019.50 using the latest version of the MCM model, MCM6, covering up to 2019.50. The difference between the MCM3 and MCM6 models is as small as ~ 5 nT at the Earth's surface on average. This hindcast can be a rigorous test of our forecast scheme since a geomagnetic jerk occurred in 2014 (Torta et al. 2015). We set the ensemble size as 960 and prepare ensemble members from a single free dynamo run for a long time corresponding to 320,000 years, where the ensemble members are taken randomly to be longer than one turn over time of core convection apart from each other. We found that iterative renewals of the ensemble members after obtaining x 0,m+1 in Eq. (9) is remarkably effective, compared with the original 4DEnVar (Liu et al. 2008). This approach enables us to deal with weak nonlinearity of the physical process in concern by shrinking the matrix ⌣ X 0,m through iterations, which meets the requirement (II) for the linear approximation we discussed after Eq. (6). We refer to Nakano (2020) for the efficiency of the iterative approach in 4DEnVar. In this section, we focus on the results after five iterations ('five' was found enough to minimize the cost function of Eq. (12) from our experiments). We sequentially reduce α S and α UW (Eq. 39) through iterations, and the values in the final iteration are listed in Table 1. By the numerical experiments, we investigate dependences of our data assimilation and forecast performance on several parameters: length of the assimilation window, data weights ( α S , α UW ), and types of prediction approaches for the following 5 years.
In this section, we investigate three types of forecasts: (i) Ensemble-weighted sum, shorten to "Ens. wei. sum" hereafter, is a forecast based on Eq. (13). After the five iterations of our data assimilation procedure, the final MHD ensemble members g k x (n) 0,M for n = 1, . . . , N and g k x 0,M for k = 1, . . . , K are available with the optimized weight ŵ M . We further run MHD dynamo simulations to obtain g k x (n) 0,M and g k x 0,M for k > K and calculate the field in the forecast period by Eq. (13). (ii) The MHD forecast is the forecast by an MHD dynamo simulation starting from a state vector given by data assimilation at time t k ′ (Eq. 16). The MHD forecast differs from the Ens. wei. sum in the sense that it involves only a single MHD dynamo simulation after the data assimilation procedure. Since t k ′ is arbitrary within the assimilation window, we test k ′ = 0 (at the beginning of the window) and k ′ = K (at the end of the window). If the adopted linear approximation (Eqs. 3 or 17) holds precisely, the MHD forecast is equivalent to the Ens. wei. sum by Eq. (13). (iii) The KD forecast is the forecast by a KD simulation starting from a state vector given by assimilation at time t k ′ (Eq. 16), similar to the MHD forecast. The only difference is that the KD forecast fixes the flow field in the outer core through the simulation. As for the KD forecast, we test only the case with k ′ = K.
We test the three types of forecasts to choose the best way to generate an SV candidate model for IGRF-13SV. For comparison of our results with IGRF-12, we assume the release time of our hindcast to be 2014.75, taking into consideration that one can use magnetic data until roughly 0.50 year before the release epoch for generation of IGRF candidate models.
To assess results of our data assimilation and forecasts, we use the square root of the data misfit dP (e.g., Whaler and Beggan 2015) defined by where the subscripts "model" and "data" mean that the Gauss coefficients come from our data assimilation-forecast procedure and from the observation (the MCM6 model in this study), respectively. To obtain g m l model and h m l model in the unit of nT, we use S 0 1 ref in Eq. (21). At the beginning of the data assimilation, S 0 1 ref is set for both the observation and the numerical geodynamo, here named S 0 1 ref_data and S 0 1 ref_model . To multiply S m l from the numerical geodynamo by the ratio of S 0 1 ref_data /S 0 1 ref_model generates dimensional numerical magnetic fields available in Eq. (40). To assess the performance of SV forecast, we further define (40) where �g m l release = g m l model (t release ) − g m l data (t release ) and a similar equation for h m l release , and t release denotes the release time of the forecast. dP wo (t) measures dP without the main field offset at the release time and purely evaluate the accuracy of SV from the release time. Table 1 summarizes all the assimilation and forecast settings with the results of √ dP and √ dP wo 4.5 years after the release time.

Bias and trend terms
We first explain the effect of the bias and trend terms in Eqs. (10)-(13) on our data assimilation. We introduced those trend and bias terms, only in the final iteration and only for the magnetic field variations to reduce the data misfits. We found that the bias and trend terms represent model errors of the dynamo simulation well. Figure 5 shows three examples of the MHD forecasts with the initial condition at k ′ = 0 (Eq. 16): Case O1 for an MHD dynamo forecast without bias and trend terms, Case O2 (41) dP wo (t) =

Table 1 Summary of the numerical experiments
In calculations of √ dP (Eq. 40) and √ dP wo (Eq. 41), the MCM6 model, which covers up to 2019.50, is used as g m l data and h m l data . Release times for all cases are assumed to be 2014.75, except 2015.00 for IGRF-12. In the "Forecast type" column, k ′ = 0 and k ′ = K mean the MHD dynamo or KD simulations running from the optimized state vector at 2004.25 and 2014.25, respectively. See text for details of the three forecast types. In the α S and α UW columns, we specify those for the final 5th step in Eq. 39, where α ′ = 14 and α ′′ = 7 Italic-Case A4 was found to be the best setting in our numerical experiments for that containing the effect of the bias term, −b , and Case A2 for that containing the effect of the combination of the bias and trend terms, −b − ak , when minimizing the cost function V (x 0 ) in the final iteration. One can see that the data misfits became smaller successively by introducing the bias and trend terms. We here briefly discuss the accuracy of the linear approximation Eq. 3 from the results in Fig. 5. We found from all of our experiments that "Ens. wei. sum" (Eq. 13) results in symmetric √ dP within the assimilation window with respect to the central time of the assimilation window (see Case A1 in Fig. 5). This is likely because we adopt the time-independent data covariance R k (Eq. 39) through the window. On the other hand, the MHD forecasts with k ′ = 0 (Eq. 16) often show asymmetric √ dP lines in the assimilation window (see Cases O1 and O2 in Fig. 5). The departures of the MHD forecast lines from the symmetry, especially at the end of the assimilation window, imply break-down of the linear approximation. This break-down probably arises from large deviations, x (n) 0,m −x 0,m (n = 1, . . . , N ) , which breaks the requirement (II) specified after Eq. (6), because the requirement (I) is satisfied by the 10-year assimilation window, much shorter than the error growth e-folding time τ e of our numerical dynamo (~ 130 years). In Fig. 5, however, the MHD forecast line approaches symmetry by introducing the bias and trend terms. The MHD forecast with bias and trend terms (Case A2) fits the Ens. wei. sum (Case A1) from 2004.50 to 2012.00 very well, which indicates that the linear approximation Eq. 3 holds over 75% of the assimilation window. This implies that the introduction of the bias and trend terms is effective not only for reducing the data misfit √ dP but also for shrinking ⌣ X 0,m (Eq. 4) and keeping the linear approximation accurate. From the small departure of Case A2 from Case A1 in Fig. 5, we conclude that both requirements for the linear approximation (I) and (II) are satisfied by the 10-year window and by the introduction of the bias and trend terms in our iterative 4DEnVar scheme, and that the linear approximation Eq. 3 is acceptable in our assimilation scheme. Hereafter, we only discuss results including both bias and trend terms in the final iteration.

Recovery of core surface flow by data assimilation
We next check how our data assimilation recovers the core surface flow. To evaluate this, we use a data misfit for the velocity field based on Holme (2015) by the expression, where the U m l and W m l are the Schmidt quasi-normalized poloidal and toroidal functions for the velocity field defined by Eq. (22). The subscripts "model" and "data" mean functions from the results of data assimilation and the observation, respectively. Note that functions with the subscript "model" are generated from outputs of the numerical geodynamo at the grid r 89 , while those with "data" are prepared by the estimation of core surface flow at r b2 = 3428.5 km (see "Data 2: core surface flow" section). Figure 4b shows the data misfit for the velocity field at r = r b2 from the assimilation result of Case A. The top panel shows the total misfit, √ dP UW , while the each line in the bottom panel indicates the data misfit for each degree, dP UW ,l . The total misfit of the core surface flow, √ dP UW , is relatively large compared to the misfit of the magnetic field, √ dP (Eq. 40). This is caused by the two facts that both the bias and trend terms are not introduced in the final iteration for the velocity data, and that the factor for the velocity covariance matrices, α UW , is larger than that for the magnetic data, α S , in construction of R in Eq. 39 (see Table 1). Because our primary purpose is to predict the magnetic field over 5 years after the assimilation window, the role of the velocity data is complement to the magnetic data in our data assimilation.  Table 1 for details of each calculation Figure 4c shows the global maps of the velocity field for the observation (the MCM3 model) and the optimized one by the data assimilation (Case A1) within the assimilation window. The spatial patterns of the horizontal divergence and the arrows for horizontal velocity are consistent to some degree. However, the intensities of the horizontal velocity at r = r b2 and their horizontal divergence for Case A1 are remarkably smaller than those for the MCM3 model. This may be associated with the large Elsasser number of our numerical geodynamo ( Λ * ≈ 1.778 ). It seems difficult to optimize intensities of the magnetic and velocity field in the outer core simultaneously with our present dynamo parameters far from those of the actual Earth. Alternatively, this may be ascribed to the effect of magnetic diffusion. The horizontal velocity at r = r b2 is derived from the MCM3 model, on the assumption that the contribution of magnetic diffusion to time variations in the magnetic field can be neglected. In the numerical geodynamo model, however, the induction equation containing the magnetic diffusion term as well as the motional induction and advection terms is solved. Therefore, the velocity field at r = r b2 for Case A is probably smaller than that for the MCM3 model by the contribution of the magnetic diffusion. The discrepancy of the intensity of the velocity between the model and data remains the problem that should be solved in the future.

Length of assimilation window
We tried several lengths of the assimilation window with the maximum length of 10 years, and investigated dependence of the MHD forecasts on t ′ k , when the initial condition for the MHD dynamo simulation is given by the data assimilation (Eq. 16). Figure 6a and 6b shows comparison of the assimilation and forecast results between 10-and 5-year assimilation windows. As for the 10-year assimilation window (Fig. 6a), when the initial condition is given at the beginning of the assimilation window, i.e., k ′ = 0 (at 2004.25), the MHD forecast (Case A2; the cyan line) fits the Ens. wei. sum via Eq. 13 (Case A1; the blue dashed line) only up to 2012.0 and slightly deviates from it during the prediction period after the assumed release epoch, 2014.75 (the red downward arrow). When the initial condition is given at the end of the assimilation window, i.e., k ′ = K (at 2014.25), the MHD forecast (Case A3; the green line) fits the Ens. wei. sum well even after the release epoch of 2014.75. As for the 5-year window, both the MHD forecasts starting with initial conditions at k ′ = 0 and k ′ = K generate almost the same results as the Ens. wei. sum even from 2014.50 to 2019.50.
These results imply that the MHD forecast with the initial condition at the end of the assimilation window ( k ′ = K in Eq. 16) fit the corresponding Ens. wei. sum (Eq. 13) well over the forecast period regardless of the length of assimilation window. This comes from the fact that, when k ′ = K , g k (x 0 ) is linearly approximated with respect to the deviation from x K ,m (set k ′ = K in Eq. 17); the accuracy of the linear approximation after the assimilation window is independent of the length of the assimilation window. When k ′ = 0 in the MHD forecast, on the other hand, the shorter assimilation window results in better accuracy of the linear approximation over the forecast period. We can conclude that the adoption of k ′ = K in the MHD forecasts enables us to achieve high accuracy of the linear approximation (here Eq. 17) over the forecast period even when the assimilation window is relatively long.

Dependence on weight, α S
In our numerical experiments, we varied only α S by fixing α UW as listed in Table 1. Figure 6c and 6d shows the dependence of √ dP on α S in Eq. (39). The two panels show the results for a given α S (the blue dashed lines) and for the corresponding 0.1α S (the red dashed lines). We found that longer assimilation window is necessary to effectively reduce the misfit in the forecast period using smaller α S . On the other hand, when the initial conditions were obtained with too small α S , the MHD forecast did not produce reasonable outputs; the forecast was degraded when α S was smaller than shown in Fig. 6 because the problem got unstable.
Note here that the red dashed lines in Fig. 6 (Cases A4 and B4) correspond to results with the almost smallest α S that allows reasonable MHD forecasts for each assimilation window. We decided to adopt 10-year assimilation window for estimation of SV candidate model, because (1) the linear approximation (Eq. 3 and resulting Eq. 6) is valid over most of the assimilation window shown in Figs. 5 and 6a, and (2) smaller √ dP can be achieved by decreasing α S compared with shorter windows.

MHD and KD forecasts
We investigated the three types of forecasts for the forecast period from 2014.50 to 2019.25: the Ens. wei. sum by Eq. (13), the MHD forecast, and the KD forecast. As already shown in Figs. 5 and 6, the Ens. wei. sum is a good candidate SV, which is theoretically identical to the MHD forecast with the optimized initial condition under the linear approximation assumption (Eqs. 3 or 17). Figure 6e and 6f shows comparison of the three types of forecasts for the 10-and 5-year windows, respectively. The MHD and KD forecasts were conducted with the same initial condition given at the end of the assimilation window, 2014.25 ( k ′ = K in Eq. 16). For the 10-year window, the MHD and KD forecasts are almost the same and slightly worse than the Ens. wei. sum. For the 5-year window, the MHD forecast is almost the same as the Ens. wei. sum, while the KD forecast is remarkably worse than the other two. As far as the 10-year and 5-year assimilation windows are concerned, the following three seem to hold: (1) the Ens. wei. sum is always one of the best forecasts among the three, (2) the MHD forecast is the same as or slightly worse than the Ens. wei. sum, and (3) the KD forecast is often similar to MHD forecasts, although it is sometimes remarkably worse than the other two. We speculate that the MHD forecast is always slightly worse than the Ens. wei. sum because of a slight departure from the linear approximation over the forecast period.
We further examined the difference between the MHD and KD forecasts. Figure 7a shows the comparison of secular acceleration (SA) of the radial component of the geomagnetic field at the Earth's surface between the MHD and KD forecasts for both 10-year and 5-year assimilation windows. SA was calculated by the equation, SA(t) = (F (t + 0.25) + F (t − 0.25) − 2 · F (t))/(0.25) 2 , where F (t) represents the geomagnetic main field and t is  Table 1 for the summary). Left and right columns show results of 10-year and 5-year assimilation windows, respectively. All ordinates are √ dP , i.e., the metric for deviation from the MCM6 model. The 10-year window spans from 2004.50 to 2014.25, while the 5-year one from 2009.50 to 2014.25. α ′ (= 14) in the left and α ′′ (= 7) in the right indicate α S (Eq. 39) adopted in the final 5th iteration. a, b The weighted sum of ensemble members by Eq. (13) (the Ens. wei. sum) and MHD dynamo simulation starting from optimized initial conditions at the beginning ( k ′ = 0 ) or the end ( k ′ = 0 ) of the assimilation window (see Eqs. 14-16) (the MHD forecasts). c, d The dependence on the weight coefficient, α S (Eq. 39), where the red dashed lines show the results with 0.1 α ′ and 0.1 α ′′ . e, f The ensemble-weighted sum, MHD dynamo, and KD simulations starting at the end of assimilation window (the MHD and KD forecasts with k ′ = K ) for 0.1 α ′ and 0.1 α ′′ . In all the panels, the black solid lines indicate √ dP calculated between the MCM6 model and IGRF-12. The release time for our forecast, assumed 0.5 year after the end of the assimilation window, is 2014.75, while 2015.00 for IGRF-12 time in year. After the calculation of SA, we applied 0.75year (three-point) running average only to SA from the MCM6 model in order to reduce the noise. The results in Fig. 7a are the same as those in Fig. 6e, f. It is noticeable in Fig. 7a that SAs from the MHD and KD forecasts are totally different from each other; SAs from the MHD forecasts (Cases A5 and B5) generate synthetic jerks, i.e., a flip of spatial patterns of SA in sign from 2015.0 to 2017.0, and the phenomenon is more remarkable for the 5-year window. As for the KD forecasts (Cases A6 and B6), the spatial patterns of SAs do not change significantly from 2015.0 to 2017.0.
We show in Fig. 7b, c the comparisons of SV and SA for two example Gauss coefficients, h 2 3 and h 3 3 , which correspond to Fig. 6e, f, respectively, since we found that the large √ dP over the forecast period by the KD forecast with 5-year window (Case B6 in Fig. 6f ) arises mainly from the misfit in degree l = 3 in comparison to the MHD forecast (Case B5). Top two panels in Fig. 7c, SV of h 2 3 and h 3 3 for the 5-year window, illustrate the reason for  Fig. 6f; the KD forecast (Case B6) shows rapid monotonic changes in the SV components after the initial conditions are imposed at the end of assimilation window. These are attributed to the steady positive or negative SA through the forecast period as shown in the bottom panels in Fig. 7c, which is associated with the fixed flow velocity in the KD forecast. In contrast, since the dynamic processes that allow oscillation of SA in time are present in the MHD forecasts, oscillatory SA variation of the MHD forecast results in moderate variation of SV and small departures from the observation (see Case B5 in Fig. 7c). The rapid monotonic change in SV in the KD forecast does not occur in the 10-year assimilation window case (see Case A6 in Fig. 7b). We speculate that shorter assimilation window is likely to cause large SA and SV intensities in the initial conditions at k ′ = K (Eq. 16) as shown in Fig. 7a, whereas it is easy to reduce √ dP over the assimilation window. When the short assimilation window of 5 years is adopted, the MHD forecast that allows dynamic process seems more stable rather than the KD forecast in the use of our assimilation scheme.

Performance of 5-year forecast in comparison to IGRF-12
At the end of this section, we compare the performance of our data assimilation-forecast scheme with that of IGRF-12 in terms of √ dP (Eq. 40) and √ dP wo (Eq. 41) from the release time. Figure 8 shows √ dP and √ dP wo from our forecasts and IGRF-12 with respect to the time from the release time, where we assume that the release time is 2014.75 for our numerical experiments and 2015.00 for IGRF-12, respectively. Table 1 summarizes the values of √ dP and √ dP wo 4.50 years after the release time. In Fig. 8, we added two new results compared to Fig. 6, "extrapolation of SV of MCM6" (Case C2) and "no SV" (Case C3). In the forecast by "extrapolation of SV of MCM6", we took instantaneous SV at 2014.25 from the MCM6 model and extrapolated it to the forecast period. The "no SV" forecast in turn provides the result literally with "no SV" after the assimilation window. The other four lines (Case C1, A4, A5, and A6) are the forecast by IGRF-12 and our best forecasts with 10-year assimilation window in Fig. 6e.
In both panels of Fig. 8, the orange lines is the worst while the blue line is the best. For the orange line, one can easily understand that any forecast schemes are much better than the result with "no SV". As for the blue line, we should note that the SV at 2014.25 from the MCM6 model is very accurate compared to, for example, IGRF-12SV. This is because the MCM6 model was constructed by using the geomagnetic data up to 2019.50 and applying backward smoothing to determine SV at 2014.25 (Ropp et al. 2020), especially including the information on the following rapid change of SV due to the geomagnetic jerk in 2014. Strictly speaking, our numerical experiments also used the accurate SV at 2014.25 from the MCM3 model in calculation of the core surface flow. However, the effect of the accuracy of SV data at 2014.25, i.e., at the end of the assimilation window, is very small because of the homogeneous data covariance matrix R through the assimilation window (Eq. 39) and of the larger weight to the geomagnetic main field (compare α S with α UW in Table 1).
Our three forecasts in Fig. 8, Cases A4, A5, and A6, result in slightly worse forecast performance in terms of √ dP than that of IGRF-12 and almost the same in terms  Table 1 for the values of √ dP and √ dP wo 4.50 years after the release time of √ dP wo . Since the pure performance of SV forecast for the 5 years after the release time can be measured by √ dP wo , Fig. 8b illustrates that the performance of our data assimilation-forecast scheme is comparable with that of IGRF-12SV. It is also noticeable that our numerical experiments used only data up to 2014.25, almost at the center of the geomagnetic jerk, around 2014.0 (Torta et al. 2015), while the candidate models for IGRF-12SV used data up to around 2014.50, or even up to around 2014.75 (e.g., Lesur et al. 2015, including data slightly after the jerk center. This fact implies that our data-assimilation forecast strategy has the possibility of retrieving part of information about the coming jerks from the data just before the jerks and weakening their effect on the 5-year SV forecast. Figure 8b also reveals that the IGRF-12 is clearly better than our SV forecast during the initial 2 years in the forecast period. We, therefore, speculate that our scheme may not reproduce the future magnetic jerks accurately, but can possibly provide a SV trend including jerks over the 5-year forecast period, compared to conventional SV candidate models for IGRF. Since the Ens. wei. sum by Eq. (13) always provides the best forecast among the three forecasts through our numerical experiments, we decided to adopt it for estimation of our candidate SV model for IGRF-13, which is theoretically, and in most cases practically, equivalent to the forecast by MHD dynamo simulation starting at the end of assimilation window (the MHD forecast with k ′ = K).

Summary of numerical experiments
We found in the application of our data assimilation and to the geomagnetic data from 2004.50 to 2014.25 that: • A longer assimilation window works better than a shorter window, up to the assimilation window of 10 years, although we could not try assimilation windows longer than 10 years in this study. • MHD forecast is not systematically superior to KD forecast. However, KD forecast is sometimes much worse than MHD forecast. Forecast by the Ens. wei. sum by Eq. (13) is always one of the best forecasts among the three. • KD forecasts sometimes give worse results due to the monotonic increase/decrease in the intensity of geomagnetic field arising from the almost fixed pattern of SA since the release time, especially with the short assimilation window in the use of our data assimilation scheme. On the other hand, the MHD forecast seems stable, possibly because it can suppress large SVs by dynamic processes allowed.
• The performance of our data assimilation-forecast scheme turns out to be comparable with that of IGRF-12 for the purpose of 5-year SV forecast, even though the tested data assimilation window ends around the center of the geomagnetic jerk in 2014.0.

The SV candidate model for IGRF-13
We generated a final SV candidate model valid for the period from 2020.0 to 2025.0 based on the numerical experiments described in the previous section. We conducted the data assimilation with the 10-year assimilation window from 2009.50 to 2019.50 using the provided MCM6 model. The ensemble size was set to 960 in the same manner as the numerical experiments. We adopted values of α S and α UW for five iterations as listed in Table 2. In the 5th iteration, the adopted weights for S m l and (U m l , W m l ), i.e., α −1 S and α −1 UW , are approximately 4:1 respecting each observation error.
Before the final iteration, we prepared the final ensemble simulation members by MHD dynamo simulations for both assimilation window and further 5.50 years starting from the end of assimilation window, 2019.50. We then forecast the future magnetic field by the weighted sum of all the ensemble members, i.e., the Ens. wei. sum by Eq. (13), for each epoch with 0.25 year interval from 2019.50 to 2025.00, which yielded time-series of Gauss coefficients, p m l (t) , for each mode from (l, m) = (1, 0) through (8, 8). We generated a SV candidate model not by an averaged SV from 2020.0 to 2025.0 but by fitting of linear lines to our forecasts from 2019.50 to 2025.0. Linefitting was applied to each time-series so as to estimate our candidate secular variation, a m l , during the 5.50 years: where the time origin, t 0 , was set to the end of our data assimilation window, i.e., 2019.50. Figure 9 shows the fitting results for several modes, where we can confirm that the estimated SV (the blue lines) fits to the original SV forecast by the weighted sum of MHD simulation members (the red dashed lines).
(43) p m l (t) = a m l · (t − t 0 ) + p m l (t 0 ), Whether or not the line-fit model in Eq. (43) gives a good secular variation estimate over the 5.50 years depends on each spherical harmonic mode. For some modes, it may be necessary to include higher time derivatives such as secular acceleration in order to describe the time variation in concern properly. It, therefore, is necessary to diagnose the validity of the line-fit model by examining the distribution of the resulting ensemble members for each mode, which will be described in the next section.

Error estimation of the candidate SV model
After carrying out the data assimilation, we estimated a 5-year SV candidate model by fitting Eq. (43) to the prediction using the Ens. wei. sum by Eq. (13) over the period from 2019.5 to 2025.0. We applied the same SV estimation procedure to every ensemble member and obtained 960 SV models from the ensemble. We used the standard deviation of 960 SV values as the error of our SV candidate model. Figure 10 shows two examples of the distribution of SV values including all the ensemble members. In each panel of Fig. 9, the estimated error is displayed by a gray area behind the estimated SV candidate line (the blue line).

Discussion and conclusions
We derived an SV candidate model for IGRF-13 (listed in Table 3) using the 4DEnVar data assimilation technique and MHD dynamo simulations. This candidate model is the first 'model' contribution to the IGRF community from Japanese research groups. In this paper, we illustrated a potential of the 4DEnVar data assimilation in the field of short-term forecast of the geomagnetic field, through the numerical experiments with the existing geomagnetic data and construction of our SV candidate model for IGRF-13 SV.
Our numerical experiments using the iterative 4DEn-Var approach in the geomagnetic field prediction provided us important information. Introduction of the bias and trend terms in Eq. (10) effectively reduces the data misfit in the prediction of the geomagnetic field as shown in Fig. 5. We also confirmed that the bias and trend terms enhance validity of the linear approximation (Eq. 3) over the 10-year assimilation windows. In the numerical experiments, we tested three types of forecasts, the weighted sum of the ensemble members of MHD simulations extended to the forecast period by Eq. (13) (the Ens. wei. sum), the forecast by a single MHD dynamo simulation starting from the optimized initial state vector at t ′ k (Eq. 16) (the MHD forecast), and the forecast by a single KD simulation from the optimized initial state vector at Fig. 10 Histograms of the final distribution of ensemble members for SV ( ġ 0 1 and ġ 0 2 ). The black and red dashed lines denote the range of one sigma (standard deviation) and the mean of the ensemble, respectively. The red dashed means correspond to components of the final SV candidate model t ′ k (the KD forecast). We found that the Ens. wei. sum is always one of the best forecasts among the tested three forecast types, and that longer assimilation window (up to 10 years) is superior to enhance the performance of SV forecast for 5 years after the assimilation window. It is also found that the KD forecasts are sometimes much worse than the MHD forecasts with short assimilation window, probably due to large SA and SV intensities in the initial condition and the lack of dynamic process that suppresses monotonic growth of the SV amplitude. This feature is just a characteristic of our current data assimilation and forecast scheme. If more stable estimation of SA at the end of assimilation windows becomes feasible in the future, the KD forecast may be the best forecast in conjunction with our data assimilation scheme.
We found that the SV forecast performance of our scheme is comparable with that of IGRF-12 in terms of √ dP wo (Eq. 41) from the release time (see Fig. 8). However, there is still plenty of room of further improvement in our iterative 4DEnVar data assimilation. For example, we found synthetic rapid SV changes at the end of the assimilation window in the Ens. wei. sum forecast, e.g., the rapid SV change in h 2 3 around 2014.0 (red lines) in Fig. 7b and c, and the significant SV change in g 0 1 in Fig. 9 before and after the release epoch of 2020.0. These are possibly associated with the fact that the data misfit close to the edge of the assimilation window is likely larger than the center of the assimilation window (see cases of the Ens wei. sum in Fig. 6). Because this arises mainly from the time-independent data covariance matrix R (Eq. 39), one possible solution is to adopt relatively large weight, i.e., small α S and α UW , close to the end of the assimilation window. Suppression of large SV and SA discrepancies from the observations at the end of assimilation window is one of the urgent issues that should be solved in the future.   and η = (µ 0 σ ) −1 is the magnetic diffusivity with the magnetic permeability of vacuum, µ 0 = 4π × 10 −7 H m −1 and the electrical conductivity of core, σ . Then, the typical time scale given as τ ν = (r oc − r ic ) 2 /ν = 13, 049 years leads to ν = 12.467m 2 s −1 , which in turn gives Ω = 4.050 × 10 −8 s −1 and σ = 1.277 × 10 5 S m −1 . Fluid flow can be scaled by (r oc − r ic )/τ ν = 5.504 × 10 −6 m s −1 = 0.1736km year −1 .
The Ekman number used in Matsushima (2015) is defined as E ′ = ν/Ωr 2 oc = E · 2(r oc − r ic ) 2 /r 2 oc = 2.534 × 10 −5 . Thickness of the Ekman boundary layer is, therefore, δ E = r o E ′ 1/2 = 17.544km . It should be noted here that the effect of magnetic field reduces the boundary layer thickness. In our numerical dynamo model, grid points in the radial direction inside the outer core are given as where N r = 99 and k = 0, 1, . . . , N r . As shown in Fig. 4a, r b1 = r oc − ξ 1 and r b2 = r oc − ξ 2 should be defined to be ξ 1 < δ E and ξ 2 ≫ δ E . In this study, we adopted r b1 = 3475.9km and r b2 = 3428.5km corresponding to r * 95 and r * 89 , respectively.