Sequential modelling of the Earth’s core magnetic field

We describe a new, original approach to the modelling of the Earth’s magnetic field. The overall objective of this study is to reliably render fast variations of the core field and its secular variation. This method combines a sequential modelling approach, a Kalman filter, and a correlation-based modelling step. Sources that most significantly contribute to the field measured at the surface of the Earth are modelled. Their separation is based on strong prior information on their spatial and temporal behaviours. We obtain a time series of model distributions which display behaviours similar to those of recent models based on more classic approaches, particularly at large temporal and spatial scales. Interesting new features and periodicities are visible in our models at smaller time and spatial scales. An important aspect of our method is to yield reliable error bars for all model parameters. These errors, however, are only as reliable as the description of the different sources and the prior information used are realistic. Finally, we used a slightly different version of our method to produce candidate models for the thirteenth edition of the International Geomagnetic Reference Field.


Introduction
The magnetic field surrounding the Earth is sustained by-and constantly evolving due to-the motions in the Earth's liquid outer core. It has been observed and studied for centuries, and recent evolution in both technological and mathematical methods has allowed us to understand more and more of its dynamics. Although studies of ancient variations of the geomagnetic field must rely on archeomagnetic and paleomagnetic data which are sparse and thus give only access to long time scales, the study of recent changes relies on a worldwide, dense data distribution, and efficient data acquisition platforms. With the significant increase of the number of magnetic observatories since the 1950s and the launch of several magnetic satellite missions in the last two decades, a continuous set of high quality data has allowed a deeper insight into the magnetic field evolution.
The field measured at the surface of the Earth results from the contributions of numerous sources, that must be separated in order to access their individual variations. We are particularly interested in those of the core field, i.e. the field generated in the liquid outer core of the Earth. One obstacle to progress in this domain is the presence of fields generated by phenomena such as the interactions between the core field and charged particles in the thermosphere, where ionospheric currents flow, and in the magnetosphere. These generate currents that in turn produce signals-the so-called external fields-that contribute to the measured magnetic field. They evolve on time scales ranging from seconds to years and these variations induce currents in the Earth's core, mantle, lithosphere and oceans that also generate magnetic signals. Separating all these contributions from the core field requires an adequate handling of data and has long been an important obstacle to the development of high-resolution core field models. In the past decades, however, the geomagnetic modelling community has been able to build better, Ropp et al. Earth, Planets and Space (2020) 72:153 more accurate models of the core field and its secular variation (SV). The separation of the different sources is now well controlled, with some remaining difficulties, especially at high latitudes. Modern models-e.g. the GRIMM  or, more recently, the Chaos-6 model (Finlay et al. 2016), both of which use splines of order 6 for their time evolution-separate well contributions from the core and external sources. They are typically able to render short core field time scales, although with a precision depending on the spatial scales-e.g. small spherical harmonics degrees have a resolution of the order of 2 years. Besides using splines, sequential approaches have also been used. For example, the POMME model (Chulliat and Maus 2014) uses a 3-year sliding time window, but its time resolution remains of the same order.
In this paper, we describe an original method for modelling the Earth's magnetic field. Our aim is to build a high-resolution times series of field models over the satellite era by separating most external and internal sources at small time scales. In order to achieve this, all modelled sources must be described as reliably as possible. This implies that the models contain a large amount of parameters, leading to very large models and considerable computation times, as soon as the field is modelled for more than a few months. To overcome this problem, we use a Kalman filter process to generate a time series of snapshot models. Each snapshot model covers only a given time period, of about 3 months (or 12 months for the parent model of our IGRF-13 candidate). This reduces the amount of data analysed for a single inversion, so we need to add strong and reliable prior information on the spatial behaviour of the modelled sources, in order to further constrain our problem. The Kalman filter is a 3-step process. In the first step-the analysis, an a priori information on the field is updated through a correlation-based data assimilation method described in Holschneider et al. (2016). In this latter work, the spatial correlations of the different signals are described in the spherical harmonics (SH) domain and then used in the spatial domain to constrain the inversion process. Here, we have a similar approach, but we remain in the SH domain. Then, the information gained by the analysis at a given time step is used to predict an estimate of the field at the next time step, which serves as a new prior. When the whole time period is covered, a backward smoothing is applied. The evolution in time is not directly parameterised-e.g. through splines-but is a consequence of the analysis, prediction and smoothing steps. The Kalman filter has already been used, noticeably by Lesur et al. (2017), but no smoothing was applied there. Gillet et al. (2015) and Beggan and Whaler (2009) also proposed an implementation of the Kalman filter for geomagnetic modelling. The Kalman filter has also been used to compute core field and core surface flow models in geomagnetic modelling. Barrois et al. (2018) describe an application of the Kalman filter to the estimation of the surface core flow. Baerenzung et al. (2018) use a full Ensemble Kalman filter to simulate the evolution of an ensemble of flow models, along with its statistical properties. None of these approaches have attempted to apply the Kalman filter to the modelling of all major sources as presented here.
The remaining of the paper is organised as follows: the first section describes the data sets used and the selection criteria applied to filter them. Then, the composition of the model and the parameterisation of each included source is detailed. Next, each of the three steps of the Kalman filter is described. In the Results section, we present our model time series, spanning 2000.0 to 2019.75. Via different representations, we assess the reliability of our method by comparing the resulting models to the Chaos-6 models, and discuss the features that characterise our models. We show how our approach can bring new insight on the core field secular variation. Before concluding, the application of this methodology in deriving candidates for the thirteenth edition of the International Geomagnetic Reference Field (IGRF) is detailed (see Thébault et al. (2015) for the preceding edition of the IGRF).

Data
The data set used for our modelling is compiled from ground observatory and satellite data, and is made up exclusively of vector magnetic data. It covers a period of 20 years, from 2000.0 to 2020.0. The observatory data set is built using hourly means from all available ground observatories from 2000.0 to 2020.0, reprocessed according to Macmillan and Olsen (2013), so that it covers this time period continuously. Satellite data are compiled from the Champ and Swarm missions data sets. The Champ data covers a time period ranging from Sept. 2000 to Aug. 2010. A gap of about 3 years in satellite data separates the Champ and Swarm missions. The Swarm data set spans Nov. 2013 to the end of 2020. It includes the latest available versions (0505, 0506 and 0507) of level-1b vectorial data files in December 2019 from the Swarm A (Alpha) and B (Bravo) satellites. Swarm C (Charlie) data were not used, as their information content on the core field is very similar to that of satellite A. We distinguish between "high latitude" (HL) data, with absolute magnetic latitudes above 55 • , and "medium-to-low latitude" (ML) data, of absolute magnetic latitudes below 55 • . HL data are handled in the usual North, East, Centre (NEC) reference frame, whereas ML data are used in a Solar Magnetic (SM) reference frame, reducing this way the correlations between vector data component errors (Lesur et al. 2008). All data are originally taken in the NEC reference frame.

Data selection
We apply an overall light selection on the data set [see Thomson and Lesur (2007)]. The selection criteria for all types of data are detailed in Table 1. Different criteria apply for HL and ML data. ML data are taken only inside the 23:00-5:00 local time window, while HL data are selected for all local times. Different time samplings are set for ML and HL data, to compensate for the higher data density at high latitudes due to the nearly sunsynchronous satellite orbits (see Table 1). Data are also selected for a limited range of values of the D st index, and for positive values of the z component of the interplanetary magnetic field (IMF). Finally, observatory and ML satellite data are selected only if they are located in a nonsunlit area.
We define N v and N d , the number of selected vector data, and the total number of selected data, respectively

Data weights
The variances attributed to each type of data are given in Table 2. X,Y and Z are the coordinates in the NEC reference frame for HL data, and the SM frame for the ML data. In this study, data errors are assumed to be uncorrelated, so the time and spatial covariance is set to zero between different data samples. Similarly, the errors of the different data vector components are assumed to be uncorrelated. The inverse of the variances are used to weight the data.

Model parametrisation
To parameterise the model in time, we introduced a grid over the time period 2000.0-2020.0, composed from N t + 1 knots denoted t k , such that t k = t 0 + k t , for k = 0 . . . , N t , with �t = 365.25/4 days (i.e. roughly 3 months) and t 0 corresponding to decimal year 2000.0. It results that N t = 80 . A discrete time series of N t snapshot models is computed over the whole time period.
One snapshot model is constructed for time t k using data from the time interval [t k ; t k+1 ] . In this time interval, the model includes multiple internal and external sources, nearly all of which are parameterised through spherical harmonics (SH). Internal sources include the static core field and its secular variation (up to SH degree 18) and the lithospheric field (for SH degrees 15 to 30). They also include internal fields induced by magnetospheric currents that evolve on time scales up to a few months, along with their variations (up to SH degree 6), and the internal part of the D st indexed field, which is denoted I st (up to SH degree 3). Including this dependance to I st allows to track as fast as hourly variations of the induced fields. A known lithospheric field model (Lesur et al. 2013), computed from SH coefficients of degrees 30 to 120, is subtracted from the data. External sources include the outer magnetospheric field in Geocentric Solar Magnetic (GSM) coordinates, the inner magnetospheric field in Solar Magnetic (SM) coordinates, a time varying field indexed on the E stthe external part of the D st (in SM coordinates)-and another one indexed on hourly mean values of the Y component of the interplanetary magnetic field (IMF), in SM coordinates. All these sources are modelled for where the various g m ℓ and q m ℓ are the model parameters for internal and external sources, respectively. The sources and coefficients are listed in Table 3.
are the coordinates in the SM (resp. GSM) system of coordinates. The reference radius for the SH development, denoted a, is set to the usual Earth's surface radius, a = 6371.2 km . The symbol ℓ,m stands for the double sum L ℓ=1 ℓ m=−ℓ , with L the maximal SH degree considered. The Ŷ m ℓ,ℓ+1 and Ŷ m ℓ,ℓ−1 are the vector spherical harmonics defined by:

Table 3 Summary of all modelled sources with relevant parameters and prior information values
From left to right: name of the source, minimum ( L i ) and maximum ( L a ) harmonic degrees (when relevant), time scale, radius and scale for the HS prior (see "Analysis step" section), total number of parameters, and coefficient notation used in Eq. (1). t is set to 3 months in this study. No R and S values are given for the induced field and its variation. See "Analysis step" section for the construction of the statistics for these latter contributions

Source
L i , L a Time scale Prior Nb. of par. Coef.

360
Cġm ℓ Lithospheric field 15, 30 τ lt = 10 6 year R = 6280 km S = 2.7 · 10 −2 nT 2 SH degrees 1 to 3. For each ground observatory, local static contributions are modelled by 3 constant values, one for each vector direction. All modelled sources, along with their characteristics, are listed in Table 3. Note that the separation of the D st in I st and E st , its internal and external components, respectively, is part of the data preprocessing [see Maus and Weidelt (2004)].
The equation linking the model parameters to the value of the field B at a point (r, θ, ϕ, t) in spherical coordinates (i.e. radius, colatitude and longitude), and time t = t k + δt , with δt < �t , is given in Eq. (1), The error e is a zero mean multivariate random variable with a distribution described by a pdf N (0, ) . is a diagonal matrix, its elements being the data variances given in Table 2. Equation (4) is solved for the mean model m k and covariance C k via a re-weighted least-square (RLS) process using Huber weights (Huber 1981). The solution is where the weight matrix is W j = − 1 2 U j − 1 2 , and j denotes the index of the iterations. W j is updated at each of the 3 iterative steps of the RLS process. U 0 = I d , and U j for j > 0 is the diagonal matrix for Huber weights. The super-script t denotes the transpose.
In Eqs. (5) and (6), the prior mean model m k and its prior covariance C k are updated according to the information extracted from the data, to give the posterior mean model m k and covariance C k . This prior model distribution N (m k ,C k ) describes what we know of the model before assimilation of the data. It is, for all time steps but the first, derived from the previous time step posterior model distribution through the prediction process that is described below. However, the prior distribution N (m 0 ,C 0 ) for time interval [t 0 : t 1 ] needs to be defined.
The initial mean prior model is null: m 0 = 0 . This means we allow for each parameter to vary around zero, in a range characterised by the variances and covariances given in the matrix C 0 .
Regarding its structure, C 0 is a block diagonal matrix, with one block for the core field and SV, and one block for each of the other sources.
Regarding the construction of the different matrix blocks (designated as covariance blocks in the following), two possibilities are investigated. The first possibility is to use the Holschneider et al. (2016) type of prior that consists in information on the energy spectrum of each modelled source. For each contribution, this spectrum is defined by a scaling S, and a radius R. The value of S is set either empirically, or by optimisation, or alternatively by using the actual spectrum of the source if it is known. R is the radius at which the spectrum of the source is flat. The values of R and S are given for all modelled sources in Table 3. The resulting covariance block is diagonal, containing the respective a priori variances of each model parameter. The variances v n s ℓ , where n s refers to the source type (i.e. core, where ∇ is the gradient operator. The Y m ℓ (θ, ϕ) are the Schmidt semi-normalised real spherical harmonics usually employed in geomagnetism. Positive orders ( m ≥ 0 ) are associated with cos(mθ) terms, and negative orders ( m < 0 ) are associated with sin(|m|θ) terms. The vector harmonics Ŷ m ℓ,ℓ−1 (θ SM , ϕ SM ) (resp. Ŷ m ℓ,ℓ−1 (θ GSM , ϕ GSM ) ) are vectors in the SM (resp. GSM) system of coordinates.

Modelling method
The models we compute are model distributions that we assume to be Gaussian, a property which is required for the Kalman filter framework. Each distribution is described by the normal probability density function , m k is the mean model, and C k is the covariance of the model distribution. The mean model m k is a N m -sized vector containing the mean values of the model parameters ( N m = 2227 , see Table 3) and C k is a N m × N m matrix. The Gaussian nature of the distribution is preserved through the whole process, since all operations applied are linear.
Our modelling approach relies on a Kalman filter (Kalman 1960), a 3-step process where data assimilation relies on a correlation-based technique (Holschneider et al. 2016). The first step of the process is to model the magnetic field from a subset of data spanning the time interval [t k ; t k+1 ] , through a re-weighted least square process (hereinafter the analysis step). The model thus obtained is used to predict the model for the next time interval, through an extrapolation of its mean and covariance (hereinafter the prediction step). When the full time series is built, a backward smoothing is applied, in order to constrain each model but the last with information from posterior time intervals (hereinafter the smoothing step). Each one of these three steps is described below.

Analysis step
The data available in the time interval [t k , t k+1 ] are gathered in a vector d k . The model parameters and the data are linked by a linear operator A k , and the uncertainty on the data is accounted for by an error vector e k . The elements of A k are directly derived from Eq. (1).
To lighten notations, the index k is dropped for d , A and e as it does not impede the understanding of the equations. Therefore, the relation between the data and mean model is SV, etc.) and ℓ is the spherical harmonics degree, are defined as for internal and external sources, respectively. Note that for a given SH degree ℓ , variances do not depend on the SH order m. The initial variance for observatory offsets is set to 1000 nT 2 for all observatories.
The second possibility is to derive covariance information from a range of parameter samples. We used this approach for the core and induced fields. For the core field and SV, Gauss coefficient samples were obtained from numerical dynamo runs, of the Coupled Earth model described by Aubert (2013). Nine thousand samples were taken for each Gauss coefficient, and variance and covariances derived from them. This yields a block containing variances and covariances for both the core field and SV, as well as the cross-covariances between the two sources. For the induced field, we built the prior through an empirical approach. We first built a model where all I g m ℓ and Iġ m ℓ were imposed to take zero values. This resulted in very noisy core and SV time series of Gauss coefficients. These time series were smoothed using a 2-years averaging window, sliding with time. The residuals of this smoothing process were used to derive variances for the I g m ℓ and Iġ m ℓ Gauss coefficients. Given the simplicity of this ad hoc process, no covariance terms were calculated. The corresponding covariance block is therefore diagonal. It should be noted that these matrices concern only the spatial correlations of the core and induced fields, and do not affect their time dependance, which is handled separately (via the prediction and smoothing).
In this work, two series of model distributions were derived. Both series are parameterised in the same way and are using the same data and the same prior, except for the core field and SV. The HS series uses the Holschneider et al. (2016) type of prior information for the core field and its secular variation, whereas the CE series use the statistics derived from the Coupled Earth model for these two contributions.

Prediction step
The prediction step defines the prior model distribution N (m k+1 ,C k+1 ) at time t k+1 from the posterior model The prediction is based on the assumption that the model is a multivariate random variable with stationary second order statistics. It evolves according to where P is the prediction operator. Since this prediction step is not exact, an error term is introduced. This term is also a multivariate random variable normally distributed: N (w, C w ) . The mean of the error is null: w = 0 , but its covariance matrix C w depends, as the operator P , on the physics governing the evolution of the different sources.
In particular, for all the sources, but the core field and its secular variation, it is assumed that the parameters evolve in time as auto-regressive processes of order 1 (AR1). For such processes a single parameter g with a timescale τ evolves following: where α is defined by: The timescale τ for the different sources is specified in Table 3. The prediction error ω is a random variable with zero mean. The requirement that the statistics of the parameter g are stationary over time defines the variance v ω of ω . Assuming the parameter g has a variance v g constant in time, i.e. v g is the corresponding diagonal element on C 0 , then the variance of ω is v ω = v g (1 − α 2 ) . For external sources, the timescales are smaller than the time step t . Therefore, α = 0 and the predicted parameter has, as ω , a zero mean and a variance v ω = v g . On the contrary, for the lithospheric field, the prediction is strongly dependent on the precedent state, α ≃ 1 and the variance of ω is v ω ≃ 0. For the core field and secular variation, the evolution of a given coefficient is given through the coupled set of equations: where ω c and ω c are the errors. These errors have, to the first order, amplitudes proportional to the coefficient accelerations g m ℓ . As in the case of AR1 processes, we assume that the errors have zero means, and we set ω c (resp. ω c ) variances to vω c = vġm ℓ ( �t τ sa ) 2 (resp. (9) m k+1 = Pm k + w, where τ sa is the timescale for acceleration. The value of τ sa is known to be around 11 years up to spherical harmonics degree 13 (Christensen et al. 2012), and we used this same value up to degree 18. vġm ℓ is the known variance of the secular variation Gauss coefficient ġ m ℓ -i.e. a diagonal element of the C 0 matrix. From there, the stationary hypothesis defines the α ℓ and α ℓ : where v g m ℓ is the variance of the core field Gauss coefficient g m ℓ as estimated in C 0 . τ sv is the timescale for the secular variation. Note that we used in this work the definition of α ℓ that does not involve τ sv . As an indication, τ sv ≃ 415/ℓ is in the range of acceptable values for Chaos [see Christensen et al. (2012)]. In Additional file 1: Appendix S1, more details are given on the way the α , α ℓ and α ℓ values are derived. The exact construction of the matrices P and C w is also detailed.

Smoothing
The result of the process presented above is a time series of model distributions N (m k , C k ) for k = 0, . . . , N t . The final model distributions are smoothed versions of this time series. The smoothing consists in re-computing the model at a time step k by using the information provided by the data analysed at every time step k ′ > k . The smoothed model distributions are identified with a upper script s : N (m s k , C s k ) . This smoothing is achieved through the following equations ( Anderson and Moore (1979), Rauch et al. (1965)) for the mean and covariance of the smoothed model at time t k , respectively: where the matrix G k is defined by and P and C w are, respectively, the prediction operator and covariance matrix introduced in the previous section. These equations are similar to Eqs. (5) and (6), as they give the solution m s k of the inverse problem set by Eq. (9), where m s k+1 defines the data. (15)

Results
We present the results of the process described above applied to our data set. We recall that the full output is described by a series of mean models and covariances, each pair defining a normal distribution of models. Strictly speaking, a single model would have to be randomly drawn from this ensemble. Here, we systematically use the series of mean models for our representations. When specified, the mean models are presented with a ±2σ wide error range. We recall also that two series of model distributions have been derived, the HS and CE series, that differ only by the prior information used for the core field and secular variation statistics. Through this section, our results are compared with the Chaos-6 model (Finlay et al. 2016) in its version Chaos-6-×9.
For each time interval [t k , t k+1 ] the normalised misfit to the data is defined by: where W i is defined in "Analysis step" section. The values of R k are in the range [0.98, 1.92] nT for both CE and HS models, depending essentially on the type of data available during this time interval. Over the whole time period, the mean of the misfits for both models is 1.44 nT . It shows that our models achieve a fairly good fit to the data. These values are estimated after the analysis process, but before the smoothing steps. The misfit after smoothing is not reported here as its estimation requires adjusting other contributions to the geomagnetic field such as external and induced fields.
In Fig. 1, the CE and HS model core field radial components are plotted at the core-mantle boundary (CMB) for year 2019.5. The field models are truncated at degree 14. Figure 2 displays the radial components of the SV of the same two models at the CMB. The Chaos-6 SV model is also displayed. All three are truncated at degree 14. Both CE and HS models are quite consistent with Chaos-6, although some differences appear in small-scale features, over the Pacific or Indian oceans, for example. Some noticeable and interesting differences can be highlighted at other time periods. Concerning the core field, the largest differences occur at the poles, where they rarely exceed ±10 nT at the Earth's surface, when satellite data are available. Secular variation differences occur mostly during years 2010 to 2014. They are characterised by localised maxima in the Indian Ocean, as shown in Fig. 3. These maxima change signs between 2011.75 and 2013.75 suggesting a spike of acceleration in the HS/CE models during this period of time. Otherwise, SV radial component differences are mostly contained inside a ±5nT/year interval. The maps shown in Fig. 3 do not include the HS model, for which the results described above are also valid. All models are truncated at degree 12. Figure 4 displays the spherical harmonics power spectra of the calculated core field and SV (HS on the left, CE on the right) for year 2019.5. The crosses picture the Chaos-6 power spectra, and the spectra of the prior are displayed in dotted lines. These spectra show the The green and red crosses represent the power spectra of the core field and SV, respectively, for Chaos-6. The spectra of the lithospheric field, at SH degree greater than 14, is not displayed general agreement of both HS and CE model core field with those obtained by other modellers, but they have too much energy at small scales compared to the prior (from degree 14 for the SV, degree 15 for the core field). As a simple, first analysis, this suggests either that our priors are poor-e.g. too small for some Gauss coefficients of the lithospheric field, or that there are some un-modelled signals in the data at these wavelengths. This excess of energy probably results from a combination of the two hypotheses. To discriminate between them would require a more thorough analysis of the results, focusing on the output distributions of the core field Gauss coefficients, their variances and possible covariances with other contributions. The conclusions of this analysis are not yet completely clear. However, it suggests that sources (e.g. small-scale induced signals and high latitude ionospheric contributions) are missing in our description of the observed signal. In favour of the first hypothesis, the CE distribution of models is closer to the prior at earlier epochs (not shown) and keeps closer to the prior at later epochs. This highlights the importance of cross-correlations, which are absent from the HS prior information.
In Fig. 5, we compare the time series for several SH coefficients of the secular variation models before and after the Kalman smoothing step (see the process description in "Modelling method" section). Results are shown for both CE and HS secular variation models. The non-smoothed time series present a time lag, as compared to the smoothed and Chaos-6 time series, which agree on the general trend of the variation for most large-scale coefficients. This lag results from the difficulty to resolve the secular variation within a single time step, using mainly information from the past-thanks to the prediction step, but very little information from the future, as only 3 months worth of data are used. It leads to an underestimation of the secular variation when the acceleration is positive, and an overestimation when the acceleration is negative. The smoothing step clearly alleviates this time lag. Figure 6 displays the time series of CE and HS secular variations compared with Chaos-6 for some low SH degrees. The variations are overall compatible with Chaos-6. Major differences can appear during the time period 2010-2014, because of the lack of satellite data combined with radically different ways of handling ground observatory data. The bottom part of the figure displays the Fourier Transform (FT) spectra of the time series to which a linear trend has been first removed to avoid spectral leakage. Coefficients time series (and their Fourier spectra) are displayed for higher SH degrees (i.e. lower spatial scales), in Figs. 7 and 8. At these spatial scales, the SV models sometimes present very different behaviours. The Fourier spectra of several Gauss coefficients show significant periodicities of 4 to 6 years (e.g. for ġ 1 6 and ġ 6 6 at 5 years). Higher degree Gauss coefficients show various peaks in their Fourier spectra, such as ġ 0 7 , which presents very strong variations between 4 and 10 years, or ġ 1 8 with peaks at 3 and 4 years. Chaos-6 is systematically and strongly smoothed, with little energy for periods under 5 years. This most certainly results from the temporal smoothing used in the process leading to this latter model. It contrasts with the approach used here where the time scales controlling the temporal behaviour of the core field are set to realistic values. However, although our HS and CE models present more energy at periods under 5 years, they are not free of anomalous features. The high variability of the Fourier spectra at small time scales is probably a signature of noise (or unidentified contributions) affecting the models. For example, the periodicity observed in some coefficients (e.g. ġ 0 7 in Fig. 8) might be explained by an incorrect separation of some ionospheric fields with the core field. This is a significant risk, considering that satellites see ionospheric fields as internal sources. These ionospheric fields could be co-estimated in the future, using a suitable parametrisation and proper prior information. Furthermore the ad hoc way in which we handled the induced fields has smoothed out short time variations of the core field Gauss coefficients, but did not allow a clear separation of core and induced field signals for periodicity longer than 2 years. To progress in this matter, it should be assessed whether the amplitude of the estimated induced fields is physically plausible, considering the mantle conductivity and external fields amplitudes. In particular, harmonics of the 11 years solar cycle (e.g. 5 years periodicity) are often more pronounced in our models than, e.g. in CHAOS-6. This is visible, for example, in the Cġ 0 1 time series of the HS model, but not in the CE model or the Chaos-6 model.
To test the coherency of our models, we compared the modelled SV time series (HS and CE SV) to the variations of the modelled core field time series. For this purpose, we computed a SV series by finite differences of the mean core field model series (referred to as the FD SV). The coherency between both the CE and HS SV series and the FD SV series was verified by evaluating the power In green: the k −2 slope, where k is the frequency, i.e. one over the period spectra of their respective differences. The residuals of this comparison have a total energy of less than 0.03% of the total energy of the SV spectrum, for both CE and HS models. The difference between the FD SV and the CE or HS SV is most important where the secular variation presents steep slopes.
Finally, we would like to point out that, because there is no smooth parametrisation in time-such as splines would provide-in our model, the core field secular acceleration cannot be derived in a robust way directly from our mean SV Gauss coefficient series. If, for further analysis, a core field acceleration model is required, it is recommended to pick a SV Gauss coefficient series in the model distribution such that the derived acceleration varies smoothly in time, rather than using the mean SV.

Discussion
The results presented in the previous section show that the core field and secular variation general features are well recovered through the process presented. Some difficulties have been identified, e.g. extracting information on the secular acceleration requires further processing of the derived series of snapshots. There are nonetheless technical advantages in using a sequential modelling approach over more classical methods, such as, for example, the relatively light computing power requirements. However, the main motivation for applying this method is first to extract from the data better models through the information provided by the prior, and second to have, as an output, reliable estimates of the model uncertainties. We discuss in the following these two points.
The spatial prior information for a time interval [t k , t k+1 ] is provided to the model through the covariance matrix C k in Eq. (5). This matrix cannot be singular because it would imply that some parts of the model are perfectly known, an hypothesis that we reject by imposing a minimum variance to all eigenvectors of the matrix. However, in the same equation the matrix A t W j A can be singular when the data set cannot resolve the model and the spatial prior information has an important role in this case. The prior information on the temporal behaviour of the model components intervene in the prediction part of the Kalman filter through the definition of the α ℓ and α ℓ in Eq. (9). They influence the output of the model as soon as they differ from unity. It follows that in our core field model, the prior information plays an important We point out that the HS and CE models differ only by the spatial prior information. It is therefore the spatial prior information that leads to different temporal behaviour for the HS and CE Gauss coefficients of the secular variation, for intermediate-to-high SH degrees as shown in Figs. 7 and 8. In contrast, both the spatial and temporal prior information influence the deviations of the Chaos-6 SV model from our model series. At these spatial scales, and for a given SH degree, the HS and CE prior variances are very similar. The observed differences are therefore mainly due to the covariance between coefficients of different degree and order. These covariances are present only in the initial prior of the CE model. In particular, covariances between well resolved large scales and poorly defined small scales can be responsible for observed differences. As long as the information carried by the Coupled Earth dynamo model outputs is correct, there is no reason to believe less in the variations of the SV Gauss coefficient distributions of our CE models than in those of our HS or Chaos-6 models. Of course, the spread of the model distributions, characterised by the variance values, plays an important role, and it is worth studying what influences these values. The example of the induced fields is here particularly instructive.
The first point to notice is that if the contribution of the induced field is neglected, the output variances of the core field distributions are particularly small. In that case, of course, the induced field signals are partly described by the core field model, the noise, and possibly spread over other components of the model. The prior information we can provide to separate them from the core field is extremely limited, particularly regarding their low frequency components. We used here the fact that they are potential fields of internal origins, with small amplitudes. This information is not sufficient to separate them from the core contributions, and when co-estimated, the output variances of the latter increase considerably. However, we also assumed that the induced fields over a 3-months period are uncorrelated to the induced field of the next time period (whereas the core field components are strongly correlated in time). It is this characteristic, although not entirely valid, that allows an acceptable separation of the two contributions 75. The blue dotted line represents the Chaos-6-×9 model. The purple and cyan blue lines, respectively, represent the CE and HS model. Both CE and HS series are presented within a 2σ interval. Bottom: time spectra of the corresponding SV Gauss coefficients with the same colour code. In green: the k −2 slope, where k is the frequency, i.e. one over the period in our model. In principle, it should be possible to use our knowledge of the deep mantle conductivity and the external fields behaviours to improve the separation.
Through this example, it is clear that a contribution which is not described in the model "leaks" inside other modelled contributions. These model components may then present spurious behaviours and small variances. As another example, the lithospheric field contribution at SH degrees lower than 15 cannot be distinguished from the static core field, and trying to model it would again increase the variance of the core field model. In contrast, the tidal signals separate well from the other components because of their well defined periodicities in time. When un-modelled they are likely to remain in the noise, whereas trying to co-estimate them should not increase the variances of the core field model. Overall, the core field model we obtained is probably not completely free of induced field contributions, it also certainly describes the small SH degree lithospheric field, some ionospheric fields, and possibly unexpected other minor, more or less static contributions. Nonetheless, while improvements are possible, our field model provides a good description of the main field, together with an estimation of its variance, that is probably slightly under-estimated.

Candidate models to the IGRF
The models presented here have contributed to set two independent contributions to the IGRF-13: the IPGP candidate to the IGRF main field 2020, and the Japanese team candidate to the IGRF predictive secular variation (Minami et al. 2020). In both cases, the data selection criteria were those described in the first section, but Swarm-B satellite data were not used.
Regarding the contribution of the Japanese candidate, the model is nearly the same as the HS model presented here with a data set including data only up to 2019.5. The scaling and reference radius used for the different covariance matrices were those of Table 3. The obtained core field and secular variation components of the model were then used as input data to the En4dVar assimilation process that led ultimately to a prediction of the mean SV over 2020-2025-see Minami et al. (2020) for details on this assimilation process.
For the IPGP main field model, the general scheme for deriving the candidate is the same as the CE model presented here, but with a data set reduced to 2013.5-2019.5 and a Kalman time step set to a year. In total only 6 time steps had to be done. This length of the Kalman step insures robust estimations of the secular variation and reduces significantly the contribution of the induced field components that required to be modelled. The scaling and radius parameter were the same as in Table 3 except for the I g m ℓ and Iġ m ℓ . Their radii R were set to 2200 km and their S values to 1 · 10 −3 . From the main field model and its SV derived this way for 2019.0, the main field was linearly extrapolated to 2020.0 giving a main field candidate not much different from other candidates, although with generally slightly lesser energy.

Conclusion
In this paper, we have presented a sequential approach to core field modelling based on a combination of a Kalman filter and a correlation-based modelling method. We aim at modelling separately most major sources contributing to the observed Geomagnetic field. The separation of these different contributions relies on a strong prior information on their spatial and temporal behaviours. We built a sequence of snapshot models constituting a time series that spans 2000.0 to 2019.75, using data from 2000.0 to 2020.0. Our model time series present mean values that are generally in agreement with recent, reliable models such as the Chaos-6 model (Finlay et al. 2016).
Nonetheless, the results suggest that more temporal variability exists in the small spatial scales of the core field compared to what is shown by classic modelling techniques. In particular, several Gauss coefficient time series present significant periodicities at time scales ranging from 3 to 10 years that are absent in the Chaos-6 time series. This technique offers the further advantage to give reasonably good estimates of the Gauss coefficients variances, provided that the separation of the different sources is appropriately handled. Our main field and SV models, as it is, probably yield slightly underestimated variances for their Gauss coefficients.
The idea behind the method presented in this paper is that, starting from a prior information on the behaviour of the different sources contributing to the magnetic field, we seek to improve this information through the analysis of data. This knowledge comes in the posterior model distribution, via the Gauss coefficients mean values and variances, that are reduced through the process. This can come, however, at the cost of increasing covariances between model coefficients, due to the incomplete separation of sources. The parametrisation of the modelled sources, and the tuning of their associated prior information is therefore the backbone of this technique. As further work is achieved for this purpose, the produced models should provide more precise and reliable information on the dynamics of the core field. The setup described in this paper, used with a different time step, has allowed the production of a candidate model for the IGRF-13.