Geomagnetic core field models and secular variation forecasts for the 13th International Geomagnetic Reference Field (IGRF-13)

Observations of the geomagnetic field taken at Earth’s surface and at satellite altitude are combined to construct continuous models of the geomagnetic field and its secular variation from 1957 to 2020. From these parent models, we derive candidate main field models for the epochs 2015 and 2020 to the 13th generation of the International Geomagnetic Reference Field (IGRF). The secular variation candidate model for the period 2020–2025 is derived from a forecast of the secular variation in 2022.5, which results from a multi-variate singular spectrum analysis of the secular variation from 1957 to 2020.


Introduction
The International Association of Geomagnetism and Aeronomy (IAGA) regularly releases the International Geomagnetic Reference Field (IGRF), which is a mathematical description of Earth's main magnetic field and its secular variation. Previous versions (e.g., Finlay et al. 2010;Thébault et al. 2015) are widely used in many disciplines of Earth sciences and applied for navigational purposes (Jiménez et al. 2012;Canciani and Raquet 2016) and in satellite orientation (Slavinskis et al. 2014).
In this study, we combine geomagnetic field observations taken at Earth's surface and at satellite altitude to construct continuous models of the geomagnetic core field and its secular variation between 1957 and 2020. From these models, candidate models for the IGRF (Alken et al. 2020), i.e., main field models for the epochs 2015 and 2020 and a secular variation model for the period 2020 to 2025 centered at 2022.5 are derived. We apply two modeling techniques to derive these models.
First, a method that descends from the time-dependent modeling technique developed by Bloxham and Jackson (1992); we refer to this as the classical model. The second technique is based on a method for constructing core field models that satisfy the frozen-flux radial magnetic induction equation on the core-mantle boundary (CMB) by imposing the field evolution to be entirely due to advection of the magnetic field at the core surface Wardinski and Lesur 2012), which we refer to as the kinematic field model. The latter method could be understood as a simple data assimilation approach, where the diffusion-less induction equation and assumptions about the dynamical regime of the core flow form the priors, and observations define their likelihood.
Methods of forecasting the future geomagnetic field evolution range from simple linear extrapolation to data assimilation into numerical dynamo simulations Fournier et al. 2015). Here, we devise two strategies to forecast the geomagnetic secular variation. First, a direct forecast based on a multi-variate singular spectrum analysis (MSSA) (Broomhead and King 1986;Plaut and Vautard 1994) of the magnetic field variability of past decades. Second, a kinematic forecast scheme is applied that is also based on the MSSA, but of the core flow variability of past decades. The reconstruction of the past flow variability and its forecast are used to predict the future geomagnetic field by forwardmodeling the diffusion-less radial induction equation on the CMB. This approach is somewhat similar to geomagnetic field forecasts using steady and time invariant flows (Beggan and Whaler 2010;Hamilton et al. 2015;Whaler and Beggan 2015). However, such forecasts are expected to fail at the occurrence of geomagnetic jerks that are sudden changes in the secular variation. Such events occurred in the past decades Brown et al. 2013), most recently in 2014 (Torta et al. 2015). The cause of these events is not fully understood. Their occurrences have been related to different types of rapid wave motion within Earth's liquid core (Bloxham et al. 2002;Aubert and Finlay 2019), temporal changes of the core flow (Wardinski et al. 2008) and Earth's rotation variation de Viron 2005, 2013).
This paper is organized as follows. "Geomagnetic field modeling" section outlines the two techniques to derive the parent geomagnetic field model for the IGRF candidates. In the third section, we develop the methodology to predict future geomagnetic secular variation. "Results and discussion" section provides results of the geomagnetic field modeling, secular variation forecasts and the derivation of the candidate models. The last section discusses the results and concludes the study.

Geomagnetic field modeling
In this section, we summarize the derivation of a parent geomagnetic main field model from which we deduce an IGRF candidate model. The parent model, hereafter C 3 FM3, covers the period from 1957 to 2020. The model derivation follows that of Wardinski and Lesur (2012) and consists of two branches, a classical model without the kinematic constraint applied (see "Classical modeling" section), and a kinematic field model based on the tangential geostrophic flow assumption (see "Kinematic field modeling" section). Like in the previous model, C 3 FM2, we use order 6 B-splines to parameterize field and flow coefficients in time. The spline knot spacing is set to be roughly 1.5 years. Both model branches are constrained to fit a main field model for the epoch 2015. This main field model is based on magnetic measurements taken by the Swarm satellite mission (Lesur, priv. comm.). We choose 2015, as it is the epoch of the last IGRF, with a good data coverage provided by geomagnetic observatories and the Swarm satellite mission. This data coverage decreases towards the model endpoint.

Data
In this work, we use two types of data, measurements taken at a network of ground-based geomagnetic observatories and satellite data taken at satellite virtual geomagnetic observatories (Mandea and Olsen 2006). The idea of combining ground-based and virtual observatory data to perform a geomagnetic field modeling was already carried out by Barrois et al. (2018). However, here we derive secular variation estimates to avoid leakage of sub-annual external field variations into the description of the core field and to obtain a sufficient representation of the short-term secular variation, that may be used to forecast geomagnetic field changes.

Ground-based geomagnetic observatories
A large portion of the data used in this study comes from ground-based observatories. Like in previous studies (Wardinski and Holme 2006;Wardinski and Lesur 2012;Lesur et al. 2018), we derive estimates of secular variation by annual differences from observatory monthly means, where these monthly means are averages of observatory hourly means. Also, annual means are used for observatories for which hourly mean values are not available from the World Data Centre for Geomagnetism -Edinburgh (2019). These observatory annual means are part of a compilation that is provided by the British Geological Survey -Edinburgh (2020). Over the period 1957-2018 the number of geomagnetic observatories simultaneously in operation that have been providing vectorial hourly means of North, East and downward components ranges between 72 and 155. Data errors were removed when encountered and data gaps were not filled by interpolations. Figure 1 maps locations of the ground-based geomagnetic observatories used in this study.

Satellite virtual geomagnetic observatories
We use vector magnetic field measurements from Swarm Level-1b data product, version 0505 (0506 for some data files). All three Swarm satellites are considered for the period between January 2014 to June 2019. Data is screened for quality flags defined in the Level-1b Product Definition Document (Tøffner-Clausen and Nielsen 2018). We keep only measurements identified as nominal, and also Swarm C vector measurements after 4th November 2014.
We select only data where the Sun is below the horizon (Chambodut et al. 2002). Additionally, we retain only data showing moderate magnetic activity. Sectorial magnetic activity index a σ (Chambodut et al. 2013), provided by the International Service of Geomagnetic indices (2020), were used and we only select data corresponding to a σ < 25 nT.
With the selected data we then construct a global mesh of virtual geomagnetic observatories (VO) following Saturnino et al. (2018), with some small changes. An approximately equal area mesh is obtained with the VO centers separated in latitude by 12.8 • and defining 14 latitudinal bands, with θ vo = ±6.40 • , ±19.20 • , ±32.0 • , . . . , ±83.20 • . In each band, the longitude φ vo of each VO and the number of longitudinal divisions, N φ VO (rounded up to the nearest integer), are chosen so that: The resulting mesh contains 258 VOs. Figure 1 displays the locations of all VOs in the mesh. The data set of each VO consists of selected data acquired inside a cylinder of 3.0 • radius centered around each VO and during a 30-day period, i.e., leading to nearly monthly values. The Equivalent Source Dipole (ESD) technique is then used, following closely Saturnino et al. (2018). For each month (30-day period) the ESD inversion is applied to each VO vector data for the equivalent magnetization of dipoles placed at 2900 km depth inside Earth's interior, by a least squares fit in an iterative, conjugate gradient, inversion scheme (Purucker et al. 1996). Then, the forward calculation is used to estimate a magnetic field value at the VO center location and for a given time period. In this way, time series of magnetic field values at the center of each VO and at a constant altitude of 500 km, are obtained.
The distributions of ground-based and virtual observatories differ as can be seen in Fig. 1. Overall, geomagnetic observatory locations cluster, which may lead to a higher spatial resolution of the model in some parts of the world, whereas in other parts the resolution may be lower than that of the virtual observatories.

Secular variation estimates
We derive secular variation estimates as input for the geomagnetic field modeling. The technique is applied to monthly means of VOs and to ground-based geomagnetic monthly means. The secular variation of the X-component at a given observatory is estimated as where τ denotes a particular month. These are annual differences of observatory monthly means. Likewise, observatory annual means are treated using where t is in calendar years and dt is 1 year. Then, secular variation estimates derived by (2) and (3) are given in nT/year.

Classical modeling
Conventional geomagnetic field modeling approaches rely on the assumption that Earth's magnetic field B(r, θ , φ, t) is a potential field without magnetic sources in the mantle and in the vicinity of satellite virtual observatories. Because of this, the geomagnetic field is determined as a gradient of a scalar potential, i.e., (2) dX/dt| τ = X(τ + 6) − X(τ − 6) /dt, The distribution of ground-based and satellite virtual geomagnetic observatories. The ground-based observatory data (blue circles) cover the period 1957 to 2018, whereas the VO data (red stars) are derived from Swarm and cover the period 2015 to 2019.5 Then in a spherical geometry, the scalar potential of the geomagnetic field can be represented as where a is Earth's radius (6371.2 km), (r, θ, φ) the geocentric spherical radial, co-latitude and longitude coordinates and P m l (cos θ) are the Schmidt quasi-normalized associated Legendre functions, with their normalization defined by with degree l and order m. The maximum spherical harmonic degree in (5) is chosen to be l max = 14 , to minimize the contamination by the crustal field. The Gauss coefficients {g m l , h m l } are expanded in time using order six B-splines M n (t): The objective function �(m) to be minimized in the inversion is: where A is an operator which relates the model vector m containing the Gauss coefficients to the data y . y − Am is the misfit between data and model, subject to the regularization. C e and C m are the error and the prior model covariance matrix, respectively. C m is an expression of the model priors.
Here, we report solutions that are adjusted to minimize the integral of B 2 r over the core surface to obtain a spatial smooth model: The matrix S has the diagonal elements Like in some previous studies ; Wardinski and Lesur 2012), we seek a reliable estimate of the secular acceleration. Therefore, the temporal model 2π φ=0 π θ=0 (P m l (cos θ) cos(mφ)) 2 sin θ dθ dφ = 4π 2l + 1 , constraint is to minimize the integral of the third time derivative of the radial field component over the core surface and in the model period between t S and t E where S(c) is the spherical surface of the core at radius c = 3485 km . The diagonal elements of the matrix T are the same as for S , and the time integral is computed using a Newton-Cotes formula of a closed type, e.g., Bode's rule (Abramowitz and Stegun 1973). Minimization of the third time derivative requires placing further conditions on the second time derivatives of the radial field at the model end-points; best results are obtained when these are set to zero. The model prior covariance matrix C m is then given by: The last term of (8) is the constraint to fit a given satellite field in 2015, which is where 0 g m l and 0 h m l are the Gauss coefficients of the satellite geomagnetic main field model. This constraint is necessary, as our model is based on secular variation data and it needs a main field model at a given epoch in order to provide also description of the main field at all times. Traditionally, we use a main field model that is derived in an independent study from satellite data. The constraint is then written with the damping parameter f .
Solutions are sought iteratively in a very similar manner as for the previous model, C 3 FM2, by deriving an initial model to re-weight the observatory data by their residuals to this initial model. Then, the strength of external field variation is reduced by a noise-removal scheme (Wardinski and Holme 2011). From this data set the final model is derived.

Kinematic field modeling
In this section, we describe our method to invert geomagnetic observations for field and flow at the core surface, which extends our previous study (Wardinski and (11) Lesur 2012), where we imposed the core flow to be purely toroidal. The inversion of secular variation data for field and flow at the core surface formulates to a Bayesian inference. Assuming a Gaussian distribution, then this leads to an objective function similar to (8) that is minimal for the preferred solution m.
In the following, we provide details of the prior information (constraints) used to derive a preferred solution of the non-unique and non-linear inverse problem. The constraints are applied on the portion of the model vector that represents the core surface flow.
Here, we use the data set with the reduced external field noise obtained in "Classical modeling" section to jointly invert for the field and flow at the core surface. This eases the joint-inversion process, as it avoids the iterative solving scheme to re-weight the data for the non-linear problem. Different assumptions of the flow dynamic could be applied (Holme 2007). Among them: purely toroidal (Whaler 1980), tangential geostrophic (LeMouël 1984) and quasi-geostrophic flow (Pais and Jault 2008). However, we focus on the tangential geostrophic flow assumption, as it is more comprehensive than a purely toroidal flow, but less restrictive than a quasi-geostrophic flow. (Note that the term flow refers to its horizontal part only, as the radial part vanishes at the core surface).
The flow is decomposed into toroidal and poloidal components: T and P are scalars which are expanded in Schmidt-normalized real spherical harmonics in space and B-splines of order 6 in time, represented by t m l (n), s m l (n). Following Lesur et al. (2010), the objective function of the joint inversion for the field and the flow at the core surface reads where the model vector m now contains the sets of Gauss and flow coefficients. The functional 1 is related to the kinematic constraint and defined bẏ �(t) is a design matrix based on the radial induction equation in the kinematic assumption, i.e., According to (18), the secular variation of the radial field component at the core surface, ∂ t B r , can be expressed in terms of a core field, B r , advected by the core fluid flow, u , where ∇ h is the horizontal divergence. Then where g is a vector that contains the Gauss coefficients and u contains the flow coefficients. The elements of the diagonal weight matrix C g are defined as w g l = 4π(l+1) 2 (2l+1) . Minimizing the mean square difference, integrated over the core surface S(c) and time, between the observed secular variation and the secular variation generated by the flow, then the functional 1 is equivalent to the integral and we can write this (similarly to (9) and (11)) as where m is now the model vector that holds the flow and Gauss coefficients. The diagonal elements of the field part of K are given by C g . The parameter 1 controls the conformity of the model to the kinematic constraint.
Because A g (t) in the functional 1 depends on the Gauss coefficients g m l (t), h m l (t) and is multiplied by u , this optimization problem (inversion) is clearly nonlinear and has to be solved iteratively. However, the iterative process is unlikely to converge unless some constraints are applied on the flow model, as finding a flow model is an ill-posed inverse problem (Holme 2007). In order to obtain the optimal field model and simultaneously reduce the null space for the flow, two types of constraints are considered. The flow model is forced to have a convergent spectrum, i.e., to be large scale, and to minimize Bloxham's "strong norm" (Bloxham 1988;Jackson 1997), N has the diagonal elements (l 3 (l + 1) 3 )/(2l + 1) . The damping parameter S controls to what extent the flow follows this constraint. Minimizing (21) constrains efficiently the secular variation. Secondly the flow model is chosen such that it varies smoothly in time, with T as the associated control parameter of the flow temporal evolution that efficiently regularizes the inverse problem. The constraints (21) and (22) are similar to the temporal constraint of the classical model, i.e., (11), as they involve a temporal integration to be minimized. Finally, it is required that the flow acceleration at starting and ending epochs is minimized by t 1 and t 2 are the epochs 1957 and 2020.0, respectively. This becomes necessary, as, if the flow acceleration is unconstrained at the endpoint, it may exceed realistic values. The factor E controls this constraint. These are the basic settings for the joint inversion for the magnetic core field and the core flow. We impose a further constraint, in order to derive models that are based on different dynamical assumptions of the flow. One possible assumption commonly used is a tangential geostrophic flow (Hills 1979;LeMouël 1984) which is established by minimizing where θ is the co-latitude. Elements of G are given by Pais et al. (2004). The constraint is controlled by setting TG .
We apply two measures to discriminate models: the rms secular variation misfit is measured by differences between model and the observed secular variation, i.e., where N obs is the number of observations. In addition, for the kinematic field models we derive the ratio that specifies to what extent the frozen-flux radial induction equation is satisfied. Values similar to 1 and larger mean that this constraint is not fulfilled by the model, and conversely values significantly smaller than 1 indicate that the frozen-flux radial induction equation is approximately satisfied.

Forecasting schemes
In this study, we aim to obtain a description of (a) the temporal dynamics of the secular variation and (b) the temporal variability of the advective motion within the liquid outer core. These should not differ largely, but remaining external signals in the secular variation estimates may pose problems to identify clearly signals from the core. To forecast states of the physical system that lead to the secular variation it is necessary to analyze the observed time series. Our strategy relies on the derivation of multi-variate time-series models, where individual secular variation and flow coefficients are treated as time series. (23) Our development of time-series models and forecasts is based on the multi-variate singular spectrum analysis (MSSA) (Broomhead and King 1986;Plaut and Vautard 1994, and see "Appendix A" section) of field and core flow variability. Other methods could be chosen to find a sufficient description of the field and flow variability, such as ARIMA-models, or vector auto-regressions models; however, these model types are based on assumptions like stationarity of the time series and normality of the residuals (Box and Jenkins 1976;Brockwell and Davis 2002, and references therein), which are unlikely for the temporal variability of the geomagnetic field . Furthermore, the singular spectrum analysis does not rely on assumptions about the linearity or nonlinearity of the process that is generating the time series, whereas auto-regression models implicitly assume a linear behavior of the observed data. If 1, . . . , T serves as the time range for a training set, then the model parameters are estimated from these observations. The (outof-sample) prediction is generated over the time range T + 1, . . . , T + h according to the generation mechanism of the model. In this study we consider 10 years as the time range of the prediction. This sets h = 120 months. Prediction models are derived using the first 22 eigenmodes obtained by the MSSA of the temporal flow variability (see "Appendix A" section).

Direct secular variation forecast
The time series of the Gauss coefficients are taken from the classical model branch of C 3 FM3. The analysis considers discrete series of secular variation coefficients consisting in N observations regularly sampled in time. The sampling time τ s is chosen to be one month.
Time-series models are derived using a multi-variate singular spectrum analysis, where the temporal variability of the series is decomposed in eigenmodes (Vautard et al. 1992;Golyandina et al. 2001;Ghil et al. 2002). The time-series models are then reconstructions based on superposition of these eigenmodes, and similarly forecasts. For the technical details of the model selection, we refer the reader to "Appendix A" section.

Kinematic secular variation forecast
The kinematic forecasting scheme to predict future secular variation is motivated by results that the observed secular variation can be mostly explained by the advection of magnetic field (Roberts and Scott 1965;Wardinski and Lesur 2012). Additionally, results of a previous kinematic field model show that the spatial complexity imposed by a satellite magnetic field model of 2004 is maintained backward in time over decades (Wardinski and Lesur 2012). Here, we (27) y t :=ġ m l (n), n = 1, . . . , N want to explore to what extent the approach can be used to forecast geomagnetic field changes. The scheme consists of two steps. First, we derive timeseries models that give a sufficient accurate description of the temporal variability of the fluid motion in the liquid outer core. Therefore, we consider toroidal and poloidal flow terms, which are treated as discrete multi-variate time series from which we derive time-series models for a given 'learning' phase. The learning phase is set to cover the interval 1957 to 2020.
The second step, the kinematic forecasting step, employs forward computation of the secular variation by using the diffusion-less induction equation and is initialized by where ũ(t 0 ) represents the flow ( t m l (n), p m l (n) ), and B r (t 0 ) the radial magnetic field at an initial epoch t 0 . Here, we use the flow coefficients of the kinematic branch of C 3 FM3. The forecasts of ũ(t i ) for i = 0 are obtained by deriving forecasts of the flow coefficients from their timeseries models, similar to the operation in "Direct secular variation forecastDirect secular variation forecast" section. Generally, the forecasts of the secular variation are computed by where δt represents the sampling interval of 1 month, and ũ(t n ) is the flow forecast vector that contains forecasts of the individual flow coefficients for the first month after 2020.0. Equation (28) is computed on a spherical grid, and then transformed to spherical harmonic secular variation coefficients ( ġ m l ,ḣ m l ) similar to the scheme of Lloyd and Gubbins (1990). The forecast period is 120 months, and ends in 2030.0.

Assessing the accuracy of forecasts
Choosing an appropriate error measure in forecasting is problematic, because no single measure gives an unambiguous indication of forecasting performance, while the use of multiple measures makes comparisons between forecasting methods difficult and unwieldy (Mathews and Diamantopoulos 1994). We apply a derivative of the widely used mean absolute percentage error, i.e., the symmetric absolute percentage error (sAPE) to assess the difference between the observation, i.e., the actual field and flow coefficient y(t) and its forecast ŷ(t) as: This follows the definition by Chen and Yang (2004), which differs from those of Armstrong (1985), Flores (1986), as it does not become singular and has a maximum value of two when either y(t) or ŷ(t) is zero, but is undefined when both are zero. The range of (30) is (0, 2), i.e., the maximum value corresponds to 200% percentage error. We define a prediction length, for which the sAPE(t) becomes larger than 10%. This is a rather cautious definition and other limits could also be considered. Another measure that is widely used (e.g., Aubert 2015; Whaler and Beggan 2015) is the rms-difference between models where A g m l , A h m l and B g m l , B h m l are Gauss coefficients of the compared field models A and B.
√ dP represents the total (difference) field integrated over the surface of the Earth, and similarly is the secular variation. We compute it in order to allow comparisons to other studies which use this measures as primary diagnostic. However, we must note, that √ dP is biased towards large-scale contributions, when computed at Earth's surface where the spatial energy spectrum of the field is decreasing, and is biased towards small scales at the core surface, where the spectrum is increasing. This is not the case for the sAPE(t) , as it is relative to the amplitude of the actual coefficient.

Results and discussion
In this section, we present results of the core field modeling and the secular variation forecasts. At the end of this section, we also present and discuss our candidates for the definitive geomagnetic reference field model (DGRF) in 2015, and IGRF candidate model for 2020.
We focus on two models: Model 1 was derived with the classical modeling approach and Model 2 is a kinematic field model based on the tangential geostrophic flow assumption. Table 1 compiles the set of model parameters and related model characteristics. Both models provide almost equal fits to the data in terms of their residual standard deviation.
We note that, by strongly imposing the tangential geostrophic constraint, the fit to the radial induction equation deteriorates, when compared to solutions with a weaker tangential geostrophic constraint, see Appendix Table 3. The reason for this is not clear, but, perhaps, could be explained by the penalization of the poloidal p 0 n terms when the tangential geostrophic constraint is imposed.

Comparison of modeled and observed secular variation
In Fig. 2, we show the rms differences (31) between the two models during the model period. The rms differences of the main field (Fig. 2, left) and secular variation coefficients (Fig. 2, right) are continuously decreasing during the period, except close to the model endpoint in 2020, where they increase. These increases can be explained by the constraint on the flow acceleration in Model 2. We argue that the decrease in differences during the model period is not directly related to the growth of the available ground-based geomagnetic observatory data, because this affects both models. Rather, the decrease in differences could be explained by the different modeling techniques. As it was found in a previous study (Wardinski and Lesur 2012), the kinematic field modeling seems to project the field morphology of a spatially high-resolved satellite field model imposed in 2010 backward in time. Therefore, the kinematic model is less subject to a varying data distribution, whereas the classical field model is. Figure 3 compares the observed and modeled secular variation at different observatory sites. The model curves are derived from our preferred classical and kinematic field models. The observed secular variation series from ground-based geomagnetic observatories stop in 2017.5, as more recent data were not available at the time of this study.
The overall impression is that the classical and kinematic field models fit the observed secular variation equally well. An exception is, at the model endpoints where the two curves deviate. Most likely, this is due to the endpoint constraint of the kinematic field models. The known jerk occurrences (A-J) are all clearly visible in the Y-component of Chambon-la-Forêt, whereas in other observatory data only some jerks are detectable, e.g., the series of Sitka shows only events in 1969 and 2011. Data of the Sitka observatory show large changes of the secular variation of Y and Z components in the past decade, that may be related to the rapid movement of the magnetic North Pole.
The same large similarity between the models can be observed in Fig. 4, which shows six secular variation coefficients with the largest amplitude of the two different models. These coefficients largely determine the morphology of the secular variation at Earth's surface. Particularly, the coefficient h 1 2 is closely related to the prominent patch of reverse magnetic flux in the southern hemisphere (Terra-Nova et al. 2017). Overall, differences between the model coefficients rarely exceed 2 nT/year, except close to model endpoints, where models disperse largely. Figure 4 also shows markers of known geomagnetic jerks, to allow their identification in the evolution of the secular variation coefficients. The coefficients ġ 0 1 and ḣ1 2 represent equatorial anti-symmetric contributions of the secular variation. These coefficients carry most of the known geomagnetic jerks, apart from the event in 1969 which is either not clearly visible in ġ 0 1 , or appears later in 1970. The identification of geomagnetic jerks in the equatorial symmetric parts of the secular variation described by ġ 1 1 ,ḣ 1 1 ,ġ 0 2 is less clear. The temporally averaged secular variation spectra of the two models are shown in Fig. 5. On large and small scales, differences between spectra remain small at Earth's surface. The spectra also indicate a very high temporal variability of the kinematic field secular variation at degrees 9 and 10, suggesting that at some epochs the  . 2 Temporal evolution of the rms difference ( √ dP ) of the main field (left) and secular variation (right) between models 1 and 2 power of the secular variation of these models is weaker by a few orders of magnitude than the secular variation of the classical field model. At the core surface (Fig. 5, right) the spectral power of the secular variation grows with spherical harmonic degree. The spectra of the two models start to deviate from degree 11, where the spectrum of Model 2 flattens, whereas the spectrum of Model 1 continues to increase. We note the same high temporal variability at spherical harmonic degree 9 of the kinematic secular variation model. Maps of the radial component of the magnetic field and its secular variation derived from Model 1 at the top of the core are shown in Fig. 6. The radial magnetic field component shows a dipole-dominated morphology, but with considerably small-scale features, such as reversed flux patches in both hemispheres. The secular variation at the core surface tends to be dominated by small scales. Particularly notable features are low latitude pairs of opposite polarity secular variation which are typical to advection (Amit 2014). Also note larger patches of the radial secular variation in the vicinity of the magnetic North Pole, e.g., under Eastern Siberia. In the context of equatorial symmetry, they argued that such a zonal jet at high latitudes of the southern hemisphere would not produce detectable secular variation, because the field is oriented in the east-west direction there. Indeed in Fig. 6 in 2020 the secular variation is very weak around the South Pole (and quite strong under the North Pole). All this makes the anomalously strong secular variation around the South Pole in 2020 as seen in Fig. 7 quite suspicious. However, secular variation maps of other epochs (1969 and 2010) are very similar to maps of Fig. 6, with dominant small-scale features at mid-and low latitudes.
We conclude that classical and kinematic field models largely agree, both in spatial morphology as well as in the temporal evolution, but significantly deviate in their secular variation towards the endpoints. The cause for this deviation is the constraint of the kinematic field model onto the flow acceleration, while the classical model gets along without such constraint. is the steep decrease of ġ 0 1 by almost 15nT/year during a very short interval of about 3 years. While this decrease is determined by the statistical properties of the timeseries model, it may be difficult to accept this when considering the past secular variation. In fact, this decrease would be caused by a large secular acceleration. Advective sources and sinks of ġ0 1 exhibit a lot of cancellations and it is the rather small residual that gives the historical dipole decay (Olson and Amit 2006;Finlay et al. 2016), i.e., subtle change in the core flow pattern and its interaction with the field may yield considerable changes in g 0 1 . However, apparently such a change did not happen since 1840. The other terms vary within ranges of previous oscillations (see Fig. 4).

Direct secular variation forecast
Most of the secular variation terms of Model 1 indicate possible occurrences of two future geomagnetic jerks in late 2020-early 2021, and in early 2024. In Fig. 8, at these dates the predicted secular variation shows notable changes in slope.

Kinematic secular variation forecast
The kinematic secular variation forecast is based on forecasts of individual flow coefficients derived from a multivariate singular spectrum analysis (MSSA). Figure 9 shows the past and future temporal evolution of two toroidal flow coefficients, t 0 1 , t 0 3 . These coefficients are of a particular interest, as they relate to changes in Earth's rotation (Jault et al. 1988;Jackson et al. 1993), which The MSSA of the temporal flow variability is truncated at the degree where past model and forecast of the flow coefficients becomes almost continuous, which is achieved for a truncation level k = 22 , i.e., the first 22 eigenmodes. Both flow coefficients show spurious accelerations from 2017 on, which is caused the endpoint constraint of the flow acceleration, cf. (23). Therefore, we consider the kinematic field and flow coefficients to be flawed and start the kinematic forecast from 2017. Figure 8 also provides a comparison between the direct and kinematic forecasts of secular variation coefficients. The kinematic forecast starts in 2017 to reduce the influence of the model's faulty temporal behavior close to its endpoint. These forecasts differ by their temporal evolution to those of the direct secular variation forecast. Most apparent are these differences in the forecast of ġ0 1 , where both forecasts have opposite signs and differ by about 20 nT/year. The kinematic forecast predicts a gently varying ġ0 1 . This forecast discrepancy is also clearly seen for ġ0 1 ,ġ 1 1 and ḣ1 2 . However, these forecasts show also common features related to future geomagnetic jerks. The cause for the forecast discrepancy is not understood, and it is not in relation to the anomalous flow acceleration at the kinematic model endpoint, as the kinematic forecast starts a few years prior to the model's endpoint.

Predictability
Previous studies Whaler and Beggan 2015) used rms-based measures (31) to quantify differences between a reference model and forecasts. The reference model is a model for the epoch 2017 of the classical and kinematic model branches and labeled as M1 2017 and M2 2017, respectively. Table 2 lists the rms-differences between different epoch models of Model 1 and Model 2. Generally, the rms-differences of the direct secular variation forecasts (Model 1) is smaller than those of the kinematic forecasts. However, it is not clear whether it allows to conclude upon the performance of the forecasts. Globally, it suggests that the main field and secular variation of Model 1 tend to be more similar over longer time intervals than for Model 2. This may agree with results based on numerical dynamos ) and data assimilation using numerical dynamos . Their studies used rms-based error estimates, and their results suggested a possible predictability of several decades to a century.
We perform several prediction experiments to measure the prediction length of the first 26 secular variation coefficients for different prediction setups. In these setups, the starting of the forecast is set to 1995, 2005 and 2015, respectively. The forecast period is in all cases 10 years, for which we compute the sAPE(t) , cf. (30).
First, it is found for the direct secular variation forecasts of the classical field model that the prediction length is largely independent of the number of eigenvalues considered in the MSSA forecast. Merely, this number has to be larger than the number of significant eigenvalues.
Results of the forecast experiments are shown in Figs. 10 and 11. Generally, the prediction error (sAPE(t)) of the zonal secular variation coefficients  increases with time, except for the prediction of ġ 0 4 of the 1995-experiment which improves with time. It is also noted that prediction error of the 2015-experiment is mostly smaller than for the other experiments (Fig. 10). We focus on the prediction of the large-scale secular variation represented by the first coefficients. For practical purposes, we arbitrarily examine the first 26 coefficients. The mean prediction length for which the prediction error sAPE(t) < 10% of the first 26 secular variation coefficients is around 3 years (Fig. 11). This is substantially shorter than results by , who found a possible predictability of a century.
Similar experiments are performed for the flow forecasts of the kinematic field model (Model 2), with the resulting prediction lengths shown in Fig. 12. The mean prediction length for all experiments is always shorter than 2.5 years, i.e., shorter than the mean prediction length of direct secular variation forecasts (Fig. 11). However, we note that the prediction length of t 0 1 for experiments started in 1995 and 2005 is longer or equal to 5 years, and that only the forecast deteriorates only when the experiment is started closer to the model endpoint.
We avoided performing such experiments for the kinematic secular variation forecast, because of the faulty temporal behavior of Model 2 close to its endpoints. This behavior is caused by the constraint that controls the flow acceleration at the model endpoints by E . Constraining the flow acceleration too strongly leads to a very small flow acceleration, which is as unrealistic as the counter case with large flow acceleration. To this end, conclusions about the flow acceleration at the endpoints are loose, and in the future we may constrain secular variation of kinematic field models to be similar to the secular variation of the classical model.

Derivation of candidate models
Our candidate models for the DGRF 2015, IGRF 2020, and the secular variation for the period 2020 to 2025 are derived from Model 1 and its respective forecast. This is justified by the large resemblance of the main field description of two models. The candidates for the DGRF 2015 and IGRF 2020 are given by the main field model in 2015.0 and 2020.0, respectively, truncated at spherical harmonic degree 13, in nT with two decimal places. The candidate of the secular variation model is given by the forecast for the epoch 2022.5 truncated to spherical harmonic degree 8.

Conclusion
We derive field models for the period 1957 to 2020 from ground-based geomagnetic observatory data as well as satellite-based virtual observatory data. These models are constructed using two different techniques, a classical modeling (Model 1) that provides an optimal fit to observations, and a data assimilation to a dynamical assumption of Earth's core flow (Model 2). These strategies provide similar results for the core field and its temporal variations over the past decades.
In this study we set up two forecasting schemes that rely on analyses of multi-variate time series of secular variation coefficients. We derived time-series models of the field variability, from which forecasts of individual secular variation coefficients are obtained. These serve as direct secular variation predictions. Prediction experiments indicate a robust forecast of the secular variation of about three years. This might be a lower boundary and is determined by our cautious definition of the prediction length; the time interval when sAPE(t) < 10%.
Forecasts of the flow coefficients are derived in the same manner and used in forward calculations of the secular variation due to advection in the core, where contributions from magnetic diffusion are neglected and a tangential geostrophy assumption couples the toroidal and poloidal flows (kinematic secular variation forecast). This approach extends beyond the approach of Beggan and Whaler (2010); Whaler and Beggan (2015) which used steady and constantly accelerated flow to predict future secular variation. However, our kinematic secular variation forecast suffers from faulty estimates of the core flow at the model endpoints due to a possibly inaccurate restriction of the flow acceleration. Therefore, the latter approach may not provide robust forecasts of the secular variation, unless the endpoint constraint can be dropped. Interestingly, both strategies consistently indicate the occurrence of future geomagnetic jerks in 2021 and 2024. The uncertainty in the exact timing of these events is related to the original temporal resolution of our data, which is optimistically smaller than ±1 year and realistically larger than ±6 months, as well as to the strength of the temporal constraint. Of course, our strategies of forecasting are limited by the lack of observations from within the (geodynamo) system (which are not available), but not by inferences made upon the geodynamo. Similar limitations may occur to approaches using a full description of the geodynamo in a data assimilation framework to forecast geomagnetic secular variation, like , Fournier et al. (2015). We hypothesize that when longer time series are considered, a longer behavior of the field can be modeled; however, this will not improve the predictability of the short-term (decadal) field variations, as they may be chaotic.
To this end, results of the classical modeling (Model 1) are used to provide candidates for the definitive geomagnetic reference field model (DGRF) in 2015, and IGRF candidate model for 2020. The candidate of the secular variation model is given by the forecast for the epoch 2022.5.
We briefly recall some aspects of the singular spectrum analysis (SSA). The SSA is a non-parametric spectral estimation method that incorporates aspects of classical timeseries analysis, signal processing, multi-variate statistics and geometry (Vautard et al. 1992;Golyandina et al. 2001). We first formulate the theoretical concepts for uni-variate time series, and then extend them to multi-variate time series. For a deeper review and also its application to a variety of phenomena in meteorology, oceanography and climatology we refer to Ghil et al. (2002) and references therein.

Uni-variate analysis
The analysis is based on the embedding of a time series y t in a vector space of dimension M that determines the longest periodicity captured by the analysis. The spectral information on a time series are obtained by diagonalizing the M × M covariance matrix C y of y t . The covariance matrix C y can be estimated directly from the data, i.e., its entries c ij depend only on the lag δ = |i − j|: Usually, the decomposition of the covariance matrix C y in q orthogonal eigenvectors E q that are called temporal empirical orthogonal functions (EOF) is done by singular value decomposition. The eigenvalues ǫ q of C y account for the partial variance in the direction E q and the sum of the eigenvalues, i.e., the trace of C y , gives the total variance of the original time series y t . By a projection of the time series onto the EOF the temporal principal component A k t can be obtained: In fact, this method decomposes the time series in parts that correspond to trends, oscillatory modes or noise. Therefore, it allows the time series to be reconstructed and to be forecasted by using linear combinations of the temporal principal components and EOFs. The reconstructed components R k t are given by where K is the set of k EOFs and temporal principal components on which the reconstruction is based. We refer to k as the truncation level of the reconstruction and forecast. This truncation level is k ≤ q . The values of the normalization factor N t , as well as of the lower and upper bound of summation L t and U t , differ between the central part of the time series and in the vicinity of its endpoints (Ghil et al. 2002). An oscillatory mode is characterized by a pair of nearly equal eigenvalues and by associated principal components that are in approximate phase quadrature (Ghil et al. 2002). Such a pair can represent a non-linear, non-harmonic oscillation, because a single pair of eigenmodes are more sensitive to the basic periodicity of an oscillatory mode than methods with fixed basis functions, such as the sines and cosines used in the Fourier transform.

Multi-variate analysis
The multi-variate (or multi-channel) singular spectrum analysis (MSSA) is a generalization of the SSA from univariate to multi-variate time series, such as time series of individual Gauss coefficients. Its use was proposed theoretically in the context of non-linear dynamics (Broomhead and King 1986) and numerous examples of successful application of this methods can be found, e.g., Plaut and Vautard 1994.
The MSSA allows the identification of dynamically relevant temporal patterns that are coherent in series that form a multi-variate time series. These individual series are often called channels. By analogy to the SSA, each of L-channel data vectors y l,t : l = 1, . . . , L, t = 1, . . . , N is expanded to a state vector for each channel l = 1, . . . , L , and the window length M. Following the approach of Broomhead and King (1986); Allen and Robertson (1996), then the MSSA relies on the construction a grand covariance matrix C x like where each block C ll ′ is a covariance matrix between channels l and l ′ . The blocks C ll ′ have the entries The LM × LM matrix C x is symmetric and by diagonalizing, LM eigenvectors E k : k = 1, . . . , LM can be obtained. The eigenvectors are composed of L consecutive segments with length M. Similarly to (A.2) the associated principal components A k can be computed by projecting Y in the directions of the eigenvectors (i.e. onto the EOFs): where L T and U t are the lower and upper bound of summation, respectively. Figure 13 shows the first 12 eigenvectors of the toroidal flow decomposition. Clearly, the first two eigenvectors capture the slow variation of flow coefficients, whereas other eigenvectors represent the shorter temporal variations of the flow. The first few eigenvectors explain nearly the entire signal variance, which is given on top of each single plot. However, higher indexed eigenvectors also show non-zero partial variance, suggesting that these eigenvectors may carry relevant temporal information.
An important aspect of the analysis is how well components of the time series are separated from each other. The components generally group in two disjunctive parts, one corresponds to the signal, the other to the noise. The signal could be composed of slowly varying, periodic and/or quasi periodic components. A way to measure the separation between components, is to calculate the weighted correlation or w-correlations, as given by Golyandina et al. (2001).
In Fig. 14, the so-called, w-correlation matrix is displayed. It shows the weighted correlations for principal components, A k , of the temporal flow variability. Accordingly, large values of w-correlation exist for diagonal elements of the matrix, where individual modes correlate with themselves (dark red color). Whereas, small values of the w-correlation (lighter colors) indicate orthogonality between components, as the scalar product between perpendicular vectors vanishes. Therefore, small correlation values indicate a good separation of components or noise. In particular, noise components do not correlate with each other and other components, and therefore, show zero correlation.
The w-correlation matrix (Fig. 14) is found to be diagonal symmetric, and we mainly identify two different regions. The first region relates to the first 11 eigenvalues and their modes, where correlations quickly diminish away from the central diagonal. In the second region, from the 11th eigenvector onward, off-diagonal correlations emerge. This suggests that principal components do not fully separate. However, these correlations are minor and those principal components may still carry important information. This is almost identical for toroidal and poloidal flow coefficients.
The forecast of the flow variability requires to identify significant and non-noisy components. There are several possibilities to quantify the statistical significance of the temporal principal components. In the next section we develop a rule that provides a criterion for the selection of a sufficient correct time-series model of the flow variability.

Model selection criterion
Generally, to what extend the variability of the time series can be explained strongly depends on k, the 1 (96.39%) 2 (2.84%) 3 (0.64%) 4 (0.13%) 5 (0.01%) 6 (0%) 7 (0%) 8 (0%) 9 (0%) 10 (0%) 11 (0%) 12 (0%) Fig. 13 The first 12 eigenvectors of the toroidal flow temporal variability. The x-axes of each subplot represent the model interval . The percentages of the partial variance of each eigenvector are given in the top bar of the individual plots number of model parameter. Therefore, a selection criterion is applied which ascertains the statistical significance of the eigenvalues and their related EOF detected by the MSSA. The recovery of the temporal variability of the time series can be considered to be most complete, when the maximum number k of EOF's (truncation level) are used for the reconstruction of the time series. However, this is for several reasons not desirable, for instance some modes may represent noise that is part of the time series (and we want a noise-free representation of the flow temporal variability). Figure 15 shows the eigenvalue spectra of different analyses of the secular variation and the toroidal and poloidal flow variability. Broadly, the first 7 eigenvalues are ≥ 1 . A first break in the spectral slopes occurs at k ∼ 10 , which is related to features of the w-correlation matrices shown in Fig. 14. A second break in slope occurs at k ∼ 15 . Thereafter, spectra of the kinematic secular variation and the flow variability continuously decrease. Based on these results a MSSA-model truncated at a level k > 15 should provide a sufficient reconstruction of the temporal flow variability. However, in the kinematic forecast, we use k = 22 , at this value past and predicted flow variability become continuous. This may indicate significant portion of the signal to be carried by eigenmodes for 15 < k < 23.
The eigenvalue spectrum of the secular variation of Model 1 breaks at eigenvalue 11. Thereafter, it is flat, but eigenvalues remain an order of magnitude larger, than those of the kinematic secular variation and the flow coefficients.