A secular variation candidate model for IGRF-13 based on Swarm data and ensemble inverse geodynamo modelling

This paper describes the design of a candidate secular variation model for the 13th generation of the International Geomagnetic Reference Field. This candidate is based upon the integration of an ensemble of 100 numerical models of the geodynamo between epochs 2019.0 and 2025.0. The only difference between each ensemble member lies in the initial condition that is used for the numerical integration, all other control parameters being fixed. An initial condition is defined as follows: an estimate of the magnetic field and its rate-of-change at the core surface for 2019.0 is obtained from a year (2018.5–2019.5) of vector Swarm data. This estimate (common to all ensemble members) is subject to prior constraints: the statistical properties of the numerical dynamo model for the main geomagnetic field and its secular variation, and prescribed covariances for the other sources. One next considers 100 three-dimensional core states (in terms of flow, buoyancy and magnetic fields) extracted at different discrete times from a dynamo simulation that is not constrained by observations, with the time distance between each state exceeding the dynamo decorrelation time. Each state is adjusted (in three dimensions) in order to take the estimate of the geomagnetic field and its rate-of-change for 2019.0 into account. This methodology provides 100 different initial conditions for subsequent numerical integration of the dynamo model up to epoch 2025.0. Focussing on the 2020.0–2025.0 time window, we use the median average rate-of-change of each Gauss coefficient of the ensemble and its statistics to define the geomagnetic secular variation over that time frame and its uncertainties.


Background
The International Geomagnetic Reference Field (IGRF henceforth) encapsulates three standard mathematical models of Earth's magnetic field and is released by the International Association of Geomagnetism and Aeronomy (IAGA) every 5 years. One of these three mathematical models represents a forecast of the evolution of the main geomagnetic field for the 5 years to come, and it is referred to as the secular variation (SV) model. For the current release of the IGRF described in this special issue (Alken et al. 2020a), the 5-year window starts on January 1st, 2020 and ends on December 31st, 2024. The secular variation model is truncated at spherical harmonic degree 8; it rests on the assumption that each Gauss coefficient undergoes an independent linear evolution over the 5 years of interest. The linear interpolation of each Gauss coefficient between its adopted value in 2020.0 and any instant over that time interval will allow any end-user to compute the value of any geomagnetic element he or she is interested in during that period, at any location on Earth.
The latest version of the IGRF, hereafter referred to as the 13th generation of IGRF (IGRF-13), was released in December 2019 (Alken et al. 2020a). IPGP has been involved in a number of candidate models, notably the Open Access *Correspondence: fournier@ipgp.fr Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris, France DGRF candidate by Vigneron et al. (2019), the IGRF candidate by Ropp et al. (2020), the SV candidate described here, in addition to the IGRF candidate by Yang et al. (2020) and the SV candidate proposed by Minami et al. (2020), all of which are presented in this special issue. A total of 14 international teams proposed candidate secular variation models, a substantial increase with regard to the 12th generation, to which eight teams contributed . A task force appointed by IAGA assessed the quality of the candidate models and proposed a composite model based on those candidates (consult Alken et al. 2020b, in this issue to see how the evaluation has been carried out for IGRF-13).
The goal of this paper is to describe the secular variation candidate proposed by IPGP for IGRF-13. To forecast the evolution of the geomagnetic field until 2025, we combine geomagnetic data provided by the Swarm constellation (Olsen et al. 2013) with a physics-based model of the geomagnetic field, in the form of a numerical model of the geodynamo. This methodology, classically referred to as data assimilation, has been investigated for more than a decade in geomagnetism (e.g. Fournier et al. 2010;Wicht and Sanchez 2019). The applications of data assimilation to geomagnetism are fundamental (e.g. Gillet et al. 2010;Aubert 2015) and practical, a recent example of the latter being the forecast of the evolution of the low-altitude proton environment for the next few decades (Bourdarie et al. 2019). Another practical application of geomagnetic data assimilation is the one of interest here: it builds on our previous contribution to IGRF-12  and includes a few novelties, most notably an ensemble framework to approximate the uncertainties impacting the numerical dynamo model.

Methods
Here we describe first how we construct an "initial" model for year 2019.0 from Swarm data, and second how this initial model is used to initialize an ensemble of dynamo-based forecasts of the evolution of the main geomagnetic field from 2019.0 to 2025.0. This ensemble allows us to define a candidate secular variation model for IGRF-13 and its uncertainties.

Construction of the initial model for 2019.0 Data selection
We consider only vector data for Swarm satellite A (Alpha) from epoch 2018.5 to epoch 2019.5. We do not utilize magnetic intensity data. We use systematically the latest version of the data and select first the processed data (SW_OPER_MAGA_LR), then, depending on the epoch, either version 0505 (for the most part) or version 0506.
Data selection differs depending on the geomagnetic latitude. We distinguish high-latitude (HL henceforth) data from mid-latitude (ML henceforth) data. HL (resp. ML) data correspond to a geomagnetic latitude of absolute value larger (resp. smaller) than 55 • .
The following selection criteria apply: • ML data are selected for local times between 11:00 pm and 5:00 am, and rejected for sunlit ionosphere. The size of the corresponding data vector d is 493, 620 (= 3 × 164, 540) . The locations of the corresponding measurements are shown in Fig. 1. Note that 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, in order to lower the correlations between vector data component errors (Lesur et al. 2008).

Data weighting
The variances attributed to each type of data were chosen based on experience and they are given in Table 1. The inverse of the variance is used to weigh the data when constructing the initial model.

Modelling method and initial model parameterization
The approach used to build this model is to fit data through a robust re-weighted iterative least-squares process, using Huber weights. The first iteration is a standard L 2 -norm least-squares inversion. The magnetic field potential V is represented using a spherical harmonic expansion of the form (1) V (r, θ, ϕ, t) =V i (r, θ , ϕ, t) + V e (r, θ , ϕ, t) in which (r, θ , ϕ) denote the standard spherical coordinates, t is time, a is the mean radius of the Earth ( a = 6371.2 km), L i is the truncation of the spherical harmonic expansion of the internal sources (g and h), and L e the truncation of the expansion of the external sources (q and s). The P m ℓ are the associated Legendre functions of degree ℓ and order m, whose normalization is subject to the Schmidt convention.
The description of internal and external sources is done according to the guidelines provided by Holschneider et al. (2016). The idea is essentially to over-parameterize the contribution of each source and to constrain the resulting parameters using physical prior information.
Prior information on the model components is provided through the description of a Gaussian model distribution characterized by a mean and a covariance matrix. For all model sources the mean is zero. For the main field and its secular variation, the covariance matrix is built based on the statistics of the numerical dynamo model (more on this below). For all other sources, as in Holschneider et al. (2016), the covariance matrix is defined in terms of a radius R and a scaling factor γ . The radius R is the one at which the source spectrum is flat, and the scaling factor γ is set such that the energy level corresponds to the estimated contribution of the source to the measured signal. Depending on the source, the value of γ is set by optimization, or empirically, or by relying on the energy spectrum of the source if it happens to be known. In spectral space, this methodology amounts to defining covariance matrices C = C mm ′ ℓℓ ′ of the form and where δ denotes the Kronecker delta. The pairs of (γ , R) listed below for the various sources are slightly modified with respect to those appearing in the analysis of Lesur et al. (2018). (2) for external sources,

Internal sources
The truncation applied for internal sources is L i = 30 . The main (dynamo) field is described up to spherical harmonic degree 18, and so is its secular variation. The secular variation is assumed to be constant between 2018.5 and 2019.5, and consequently the Gauss coefficients to vary linearly with time over this time frame. The prior information used to constrain these Gauss coefficients is based on the multivariate statistics of a 70,000-year-long free run integration of the coupled Earth dynamo model by Aubert et al. (2013). This choice is motivated by the will to ensure consistency between the initial model and the forecast that is produced using the same numerical model (more on the forecast below). The crustal field is expanded in spherical harmonic degree from degree 15 to degree L i = 30 (prior covariance properties: R = 6280.0 km , γ = 2.7 10 −2 nT 2 ). A known crustal field (Lesur et al. 2013) is subtracted from the data for degrees 30 to 110.
The mantle field induced by the magnetospheric ring current is described up to spherical harmonic degree 3, and each of its spherical harmonic component is assumed to be proportional to the internal component of the D st index (prior covariance properties: R = 2537 km, γ = 1).

External sources
The external sources considered are the following: • A static external field in the Geocentric Solar Magnetospheric (GSM) coordinate system for the remote magnetosphere (prior covariance properties: R = 16000 km, γ = 5400 nT 2 ). • A static external field in the Solar magnetic (SM) coordinate system for the near magnetosphere (prior covariance properties: R = 6900 km, γ = 3.56 nT 2 ). • A time-varying external field dependent upon the external D st index (prior covariance properties: R = 16000 km, γ = 5.4). Each of these external sources is described up to spherical harmonic degree L e = 3.

The initial model
The prior information on Gauss coefficients and their sources described above is encapsulated in a covariance matrix C y . If y i denotes the vector of unknown Gauss coefficients after iteration i, iterations are governed by where A and W i are the design and weight matrices, respectively. While the former is constant, the latter, W i , is updated at each iteration according to the re-weighting (5)   (1998). Since the initial vector of Gauss coefficients is filled with zeros, y 0 = 0 , and the previous equation simplifies to The converged mean initial model, centred on 2019.0, is obtained after 3 iterations. The corresponding posterior covariance matrix C obs reads where we understand that W contains the final weights. We can sketch its block-structure according to Matrices on the diagonal express the confidence one should place in the sets of Gauss coefficients associated with each (internal or external) source. Of particular interest for what follows are C obs mf and C obs sv , which quantify the uncertainties of the main field (MF) and the secular variation (SV) of the mean model at epoch 2019.0. We shall refer to this mean model for 2019.0 as our "initial" model henceforth. Off-diagonal matrices describe the covariances between different sets of coefficients. They are not explicitly shown in our sketch, save for matrices C obs mf↔sv and C obs sv↔mf which express the covariances between the main field and secular variation coefficients. We will not consider these in the following, given the sequential nature of the inverse modelling framework to be described below.
Note that the fit to the data (weighted by the Huber weight to data ratio) for the initial model is 1.31 nT. The axisymmetric field lines of the initial field model are shown in Fig. 2.

Ensemble inverse geodynamo modelling
Our approach is similar to the one followed for our contribution to IGRF-12 . Our workflow is tentatively summarized in Fig. 2. With the initial model (and its error covariance matrices C obs mf and C obs sv ) at hand, we apply the inverse geodynamo modelling framework of Aubert (2015), based on the coupled Earth dynamo (Aubert et al. 2013), with a few novel features, namely an ensemble approach to deal with uncertainties and the possibility to enforce a QG-MAC force balance at the core surface for the initial condition to be prescribed (Aubert 2020). The acronym QG-MAC stands for Quasi-Geostrophic Magneto-Archimedean-Coriolis and defines the hierarchical force balance estimated to prevail in Earth's core (Aubert et al. 2017): a dominant geostrophic balance between the pressure gradient and the Coriolis force, followed at the next order of importance by a triple balance between the ageostrophic Coriolis force (the remainder of the Coriolis force not balanced by the pressure gradient), the buoyancy force and the magnetic force; consult Schwaiger et al. (2019) for a systematic study of the force balance in numerical dynamo simulations.
The multistep approach can be summarized as follows: (1) N e = 100 decorrelated independent dynamo states are extracted from a 70,000-year-long free run integration of the coupled Earth Dynamo.
(2) For each ensemble member e, a Kalman filter estimates the three-dimensional structure of the magnetic field B dyn e in the core interior from the field provided by the initial model B obs and the prior information based on the statistics (mean and covariance) of the 70,000-year-long free run integration. This step is referred to as the "Initialization" stage in Fig. 2. For more details regarding this step, and in particular the correlations between the magnetic field at the core surface and in its interior provided by the prior statistics, the interested reader is referred to the Appendix of Fournier et al. (2015). Figure 3 shows how the magnetic variance of the ensemble is impacted by this Kalman step. The variance is defined according to in which · and �·� denote spherical and ensemble averaging, respectively. Inspection of Fig. 3 reveals several important properties: before and after Kalman filtering, the variance is maximum on the fluid side of the interface between solid and liquid core regions, an indication of the active dynamics in this area, where plumes are initiated as this boundary layer is destabilized, and shear occurs between the solid inner core and the fluid outer core, leading to strong, time-dependent azimuthal fields. Also, the variance decreases almost linearly over the top The thick black lines correspond to the Earth's surface, the core-mantle boundary and the inner-core boundary. Middle panel: the information on B obs is used together with the statistics of the coupled Earth dynamo model to estimate an ensemble of dynamo states for 2019.0. Shown is the radial magnetic field at the core-mantle boundary for a few members of the ensemble, with a zoom on ensemble member 43. For this ensemble member, an additional meridional slice comprises the axisymmetric meridional field lines on the right, and the axisymmetric toroidal field inside the core on the left, with a scale ranging from −2 mT (red) to +2 mT (blue). Right panel: each dynamo state is advanced in time until 2025.0. The 100 trajectories so obtained are used to estimate the average secular variation between 2020.0 and 2025.0 300 km of the core, an indication of the decrease of the strength of the toroidal field. Kalman filtering preserves the native trend of the ensemble, while leading to a decrease of the variance. The variance reduction is largest in the topmost 50 km of the core.
(3) For each ensemble member e, the knowledge of the three-dimensional magnetic structure inside the core makes it possible to compute the three-dimensional magnetic diffusion vector inside the core, D mag e . (4) We next solve a diffusion-free core-flow problem at the core surface, seeking the core surface flow u s,e which satisfies where ∇ h · is the horizontal divergence and r is the unit vector in the spherical radial direction. This problem is solved under the weak constraint that the QG-MAC balance is satisfied at the core surface (Aubert 2020), in order for the resulting acceleration to remain in the range spontaneously exhibited (10) by the coupled Earth model while preserving an adequate fit to the data-this is the mid-constraint option described by Aubert (2020). (5) Finally, another Kalman filter converts this estimate of u s,e into a three-dimensional estimate of the flow and buoyancy fields, for this member of the ensemble.
Note that this sequence of operations is performed using a spherical harmonic truncation of degree and order 30, which suffices to take the unmodelled secular variation into account, as B obs and Ḃ obs are both truncated at spherical harmonic degree 13 (Aubert 2015(Aubert , 2020. The three-dimensional estimates of the buoyancy, magnetic field and flow define a unique initial condition for the integration of the numerical dynamo model, details of which can be found in the study of Aubert et al. (2013) (and references therein). This integration is performed for each ensemble member e at the native resolution of the CE dynamo (maximum spherical harmonic degree and order 133). This implies that kinetic and magnetic energy beyond spherical harmonic degree 30 get progressively enriched by virtue of the nonlinearities of the dynamo system. This forecast stage is referred to as the

Results and discussion
Before we introduce our candidate SV model, let us provide the reader with additional information with regard to the error budget and the steady flow assumption we decided to make.

Error budget and steady flow assumption
As noted in our 2015 study , the methodology presented here "leaves in principle no room for tuning parameters, save for the scheme through which non-dimensional numerical dynamo quantities are cast into the dimensional world (for this last point, we use robust physical laws presumed to hold both in the numerical models and the Earth's core (e.g. Aubert et al. 2013;Fournier et al. 2011)). " Yet at that time we had to address two issues: the first was connected with the error budget in the core-flow problem, and our interpretation was that it originated from the small size of the elements of the covariance matrix of the secular variation model, C obs sv . The second was connected with the statistical (as opposed to dynamical) nature of the initial dynamo state, which resulted in instabilities detrimental to the quality of the forecast.

Error budget
Values of the standard deviations for the initial secular variation model were on the order of 0.02 − 0.05 nT/ year for the large-scale components in our IGRF-12 study, a low value which may not reflect the fact that unmodelled external field variations can induce internal errors, and also that induced fields may not be adequately taken into account. In other words, an incomplete representation of non-orthogonal sources can lead to overestimate the accuracy of the modelled secular variation of the main field. In Fournier et al. (2015), we found it preferable to uniformly inflate C obs sv in order to enhance the compatibility of the initial field model with the statistics of the coupled Earth dynamo. That amounted to a tenfold increase of the standard deviations of the secular variation coefficients. If we define the normalized misfit J sv of the core-flow problem as where y sv and y o sv are the vectors containing the N y = 195 Gauss coefficients of the estimated secular variation and that of the initial model, respectively, the inflation allowed the normalized misfit to decrease from an initial value of 3.6 to 1.0.
For our IGRF-13 candidate, the large-scale diagonal coefficients of C obs sv are such that they imply standard deviations on the order of 0.02-0.04 nT/year. These values are similar to the ones discussed above for our IGRF-12 candidate, in spite of the distinct methodology we followed to design the initial model. Again, these low values underestimate internal errors due to unmodelled external sources on the one hand and an inadequate representation of induced fields on the other hand. In fact, these large-scale low values are entirely data-driven, as the large-scale standard deviations of the coupled Earth SV Gauss coefficients (which are used to build the initial field model in order to ensure consistency) are O(1) nT/ year. Yet, we did not inflate C obs sv this time. The normalized misfit we get without inflation is equal to 2.2 (taking the ensemble mean for y sv ), a value we consider acceptable. We propose two explanations to this lower value of J sv : • In 2015, our initial model centred in 2014.3 was not regularized; this led to a very energetic SV at intermediate scales, even more so at the core-mantle boundary, where the core-flow problem is solved. By virtue of its design, our initial model centred in 2019.0 is less energetic at those scales, and more compatible with the Coupled Earth dynamo prior. • Using an ensemble in step 4 above increases the compatibility between the dynamo and the initial SV model. The core-flow problem is solved for each ensemble member e as opposed to the ensemble mean. Each member corresponds to a fully saturated dynamo state, with a meaningful representation of intermediate to subgrid-scales. The ensemble mean, on which our 2015 approach was based, is underenergetic at these scales, which makes it harder to fit the SV model. Figure 4 illustrates the situation from an energetic standpoint, before the core-flow problem is solved. The black envelope is the one Fournier et al. Earth, Planets and Space (2021) 73:43 that would correspond to the approach we followed for IGRF-12, taking the mean and the statistics from the Coupled Earth dynamo. This tight envelope hardly intersects the blue envelope (that of the initial model), especially at large scales. Using an ensemble (thin red spectra) provides a broader range of possibilities which is beneficial to the overall scheme, once the mean (or the median) of the ensemble is considered, as will be demonstrated quantitatively in the next paragraph where the results of forecast-in-thepast experiments are presented.

The steady flow assumption
The statistical (as opposed to dynamical) nature of the initial dynamo state for each ensemble member e leads to instabilities (i.e. unwanted flow accelerations) that are detrimental to the quality of the forecast, even when the QG-MAC balance is mildly enforced, according to the methodology described by Aubert (2020). Considering the evolution equation for the secular acceleration in which σ is the electrical conductivity and µ 0 the magnetic permeability of vacuum, Christensen et al. (2012) demonstrated that in numerical dynamo simulations, the ∇ × (u × B) term is responsible for the low-degree ( ℓ ≤ 10 ) secular acceleration, thereby confirming the analysis of Lesur et al. (2010) based on co-estimated core field and flow models. Errors in this term can therefore have a large impact at the surface of the Earth. Assuming that the flow is steady over the forecast window is a conservative option. We assessed the validity of this assumption by performing forecast-in-the-past experiments, over the following time windows We begin by making some general remarks about the temporal and spectral behaviour of the error in those experiments we conducted using our framework. The error appears to grow quadratically with time over the time window of interest. In other words, it is during the last year of the forecast that most of the error is created (typically a third of the error at the end of the time window is generated during the last year of forecast). Considering now how the error is distributed as a function of spherical harmonic degree ℓ , we find that spherical harmonic degrees ℓ = 2, 3, and 4 are the main contributors to the error budget.
The results of these numerical experiments are summarized in Table 2. Considering lines 3 and 4 in the table, we observe that the steady flow strategy is systematically superior to the fully dynamical strategy, and in two instances out of three, to the linear extrapolation. Superiority means being more accurate by a handful of nT at most. Inspection of the last 4 lines reveals first that the ensemble approach always yields better forecasts, the gain ranging from 1 to 5 nT (all other parameters being equal). Second, enforcing a mild QG-MAC constraint for the force balance at the top of the core has a positive influence in two experiments out of three, with a gain of approximately 5 nT for the most recent time frame (2015.0-2019.3). These series of observations prompted us to define our candidate using a forecasting strategy Fig. 4 Secular variation spectra at the Earth's surface. Each thin red line corresponds to an ensemble member, prior to the ingestion of the initial model for 2019.0. The blue envelope is that defined by the initial model plus or minus one standard deviation. In black, the envelope defined by the mean of the ensemble plus or minus one standard deviation based on the steady flow assumption, an ensemble framework, and the imposition of the QG-MAC constraint at the top of the core.

Secular variation candidate model for IGRF-13
For each ensemble member e, under the assumption of a steady flow, we time-step the induction equation to compute the evolution of B dyn e between 2019.0 and 2025.0. We next simply define the average secular variation between 2020.0 and 2025.0 for this ensemble member as We further describe the restriction of this quantity at the core surface in terms of the time rate of change of the usual Gauss coefficients, ġ m ℓ,e ,ḣ m ℓ,e . For each spherical harmonic degree and order, our candidate Gauss coefficient is the median of either ġ m ℓ,e or ḣ m ℓ,e , e = 1, . . . , N e . The uncertainties affecting each coefficient are determined by means of the 90% credible interval. We decided to resort to the median and the 90% credible interval to define candidate coefficient values and their uncertainties to account for the fact that the distribution of some coefficients deviated substantially from a normal distribution. Coefficients and their uncertainties are listed in Table 3.
The corresponding radial secular variations at Earth's surface and at the core surface are shown in

Comparison with the IGRF-13 SV model and the historical trend
Upon inspection of the 14 SV candidate models it received, the IGRF-13 task force decided to create a composite model, using a strategy based on Huber weighting in space (see Alken et al. 2020b, for details). Information on the weights attributed to our candidate can be found in Appendix 1. The spectra of our candidate, the IGRF-13 SV model (IGRF-13-SV henceforth), and the distance between the two are shown in Fig. 6. In addition, we also show in Fig. 6 the spectra of the difference of each candidate model with IGRF-13-SV in dashed grey lines. Besides the well-known fact that the recent secular variation at Earth's surface peaks at degree 2, we notice that the spectra of IGRF13-SV and our candidate are close to one another, with a root-mean-squared distance between the two models equal to 6.83 nT/year. The difference between the two shows a decreasing trend with harmonic degree, but we note that the difference at degree 8 is on the same order than that of degree 5, a probable cause being that one year of Swarm data are not sufficient to stabilize the secular variation of the initial model for degree 8 (and above). The dashed lines in Fig. 6 indicate that all candidates are within acceptable distance of IGRF-13 SV, in spite of the variety of modelling approaches. We now examine in Fig. 7 the proposed rates of change of some Gauss coefficients of degree 1 and 2 Table 2 The forecast error is expressed as the root-mean-squared difference between the true geomagnetic field, defined by an update of the CHAOS-6 field model (Finlay et al. 2016

), and the forecast at the terminal epoch
Italics values indicate the best strategy for each period Nocast assumes that the field does not change. Linear extrapolation assumes a linear variation whose slope is specified by the secular variation up to degree 13 at the start of the forecast period, as specified by the initial field model we construct; "Coupled Earth, ensemble, QG-MAC" indicates that the full coupled Earth dynamo model is integrated for the 100 ensemble members, with the median defining the forecast, and an initialization where the QG-MAC constraint is mildly enforced. In the last four rows, "'Steady flow" implies that the sole three-dimensional induction equation (with magnetic diffusion) is integrated, with or without an ensemble approach, and with or without the mild imposition of the QG-MAC constraint. If an ensemble approach is adopted, it is the median that defines the forecast  for IGRF-13-SV and our candidate against the statistics (in terms of probability density function) provided on the one hand by the COV-OBS.x1 geomagnetic model by Gillet et al. (2015) for the 1840-2015 time span, and by the Coupled Earth dynamo model on the other hand. The values of the 13 time derivatives of the other SV candidates are shown for completeness. The candidates (and mechanically IGRF-13-SV) are clustered around specific values, with a few noticeable exceptions (one candidate has a negative value for ġ 0 1 for example). Interestingly, these consensual estimates (especially those of degree 2 and order 1) appear to be distinct from the recent historical trend. The statistics of the CE model, which cover a simulated time-span of 70,000 years, provide broader (and centred) distributions, but the recent geomagnetic secular variation seems to challenge these statistics, especially concerning ḣ1 2 . A line of research for future work appears to be the design of a dynamo model capable of accounting for these recent degree-2 features of the geomagnetic secular variation.

Conclusions
In this paper, we have presented the derivation of IPGP's secular variation candidate model for IGRF-13. In comparison to the methodology we resorted to for our previous contribution to IGRF-12, we have essentially added an ensemble formulation to the inverse geodynamo modelling framework. This addition enables two interesting novelties: first, the assessment of the uncertainties impacting our candidate is straightforward. These uncertainties are currently ignored by the IGRF task force, but the situation may evolve in the future. Second, by increasing the range of possible core states, the ensemble made it possible not to inflate the covariance of the secular variation model for 2019.0 (derived from one year of Swarm data) in order to produce decent forecasts, keeping in mind that our initial model was less energetic at small scales than the one we resorted to for IGRF-12, and therefore more compatible with the average behaviour of the coupled Earth dynamo. This is conceptually satisfying. Yet, there remains room for improvement. Limiting unwanted acceleration through the enforcement of the QG-MAC balance within the core was not efficient enough, according to our forecast-in-the-past experiments, to make the full integration of the dynamo model superior to a mere three-dimensional integration of the Coefficients (ġ m ℓ ,ḣ m ℓ ) and their uncertainty range (�ġ m ℓ , �ḣ m ℓ ) are in given nT/ year induction equation, under the assumption of a steady flow. An explanation is that the QG-MAC balance is enforced solely at the core surface, not in its bulk. Also, since the estimate of the magnetic field in the core is a statistical one, the Lorentz force based on this estimate remains of statistical nature, and therefore prone to triggering unwanted flow acceleration.
Dynamo models can now account for interannual features of the geomagnetic field (Aubert and Finlay 2019), which tells us that the acceleration issue can be at least partly resolved from the forward modelling standpoint, if the statistical nature of the estimated core is indeed partly circumvented by the combination of the enforcement of the QG-MAC balance and the ensemble formulation. From the inverse modelling standpoint, the framework can also be improved, a possibility being to merge the two steps corresponding to the estimates of i) the three-dimensional magnetic structure inside the core and ii) the core-surface flow into a single step, that can incorporate the covariances between B obs and Ḃ obs which are present in the initial field model and currently ignored. Regarding the design of the "initial" field model fed to the inverse geodynamo modelling framework, progress is also mandatory, in particular with respect to the proper identification and separation of sources. Key to this last aspect is the continuous and sustained monitoring of Earth's magnetic field with dedicated, low-Earth orbiting satellites. Mauersberger-Lowes spectra of the secular variation of the IGRF-13 SV model, the IPGP candidate and their difference at Earth's surface (thick solid lines). Also shown are the spectra of the differences between the IGRF-13 SV model and the 13 other candidates that were submitted in response to the IGRF call (grey dashed lines)

Fig. 7
Distribution of proposed rates of change for some Gauss coefficients of degrees 1 and 2 against the statistics of the COV-OBS.x1 geomagnetic field model  and the Coupled Earth dynamo model (Aubert et al. 2013). The grey dashed lines correspond to the rates of change proposed by the 13 other SV candidates to IGRF-13 candidate model In this appendix we provide the reader with information on the Huber weights that were attributed to our SV candidate model by the IGRF task force. Those weights range between 0 and 1 and indicate the degree of confidence in the candidate. They are defined on a spatial grid. At every grid point, three weights are specified, one for each component of the geomagnetic secular variation. The spatial grid comprises 10,000 points, which implies that the total number of weights is 30,000. The distribution of weights for our candidate shown in Fig. 8 indicates that in total, the fraction of weights smaller than 0.9 is smaller than 10 %. Inspection of their geographical distribution per component reveals that most of those weights markedly different from unity are given to the vertical component of the secular variation, essentially underneath Siberia. In practice, from a core-mantle boundary perspective, this points to the fact that our candidate predicts an evolution of the Siberian normal flux patch over the next 5 years sightly at odds with that predicted by most candidates.