# Broadband a priori temporal cross-covariances for the geomagnetic inverse problem: application to the satellite era

Romain Claveau, Nicolas Gillet, Pierre-Olivier Amblard, Christopher C Finlay

► To cite this version:

Romain Claveau, Nicolas Gillet, Pierre-Olivier Amblard, Christopher C Finlay. Broadband a priori temporal cross-covariances for the geomagnetic inverse problem: application to the satellite era. Earth Planets and Space, 2025, 77, pp.126. 10.1186/s40623-025-02248-z. hal-05389879

HAL Id: hal-05389879

https://hal.science/hal-05389879v1

Submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
FULL PAPER

Open Access

<figure><img src="image_2.png" /></figure>

# Broadband a priori temporal cross-covariances for the geomagnetic inverse problem: application to the satellite era

Romain Claveau¹, Nicolas Gillet¹*, Pierre-Olivier Amblard² and Christopher C. Finlay³

**Abstract**

We present COV-SAT, a magnetic field model covering the era of continuous satellite monitoring from space, constrained by ground-based and virtual observatory data. It incorporates as a priori information temporal cross-covariances associated with auto-regressive processes of order 3. These are derived analytically, and rely on a small number of free parameters (variances and time-scales) deduced from spectral properties extracted from long magnetic records and geodynamo simulations. The new a priori information proposed here encompasses power spectra stemming from paleomagnetic, observatory and satellite records, statistically replicating magnetic variations over a broad range of time-scales. It, furthermore, allows the cutoff in the spectral density expected for periods shorter than ≈ 2 years, the Alfvén time in Earth's core, to be mimicked. Field model coefficients are projected in time based on the a priori cross-correlation functions, avoiding the use of arbitrary temporal basis functions. This formalism is exploited to forecast the field evolution, and to provide uncertainties for the estimated main field and its secular variation over the upcoming 5 years.

**Keywords** Geomagnetic inverse problem, Satellite data, Correlation functions

**Graphical Abstract**

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence:
Nicolas Gillet
nicolas.gillet@univ-grenoble-alpes.fr
Full list of author information is available at the end of the article</footnote>
# 1 Introduction

A series of satellite missions has enabled the continuous monitoring of the geomagnetic field from space since about 2000. Combined with ground-based stations, these provide high quality, global observational constraints on the Earth's magnetic field evolution (e.g. Lesur et al 2022). Meanwhile, two main difficulties limit the modeling of the dynamo field originating in Earth's core (or main field, MF). First, the era covered by space records is short in comparison with the typical time-scale of the core dynamics. Indeed, the eccentric gyre first imaged from satellite data (Pais and Jault 2008), which is responsible for the observed westward drift of the magnetic field (Bullard et al 1950), is present since at least the late 19th century (Aubert 2014; Gillet et al 2019). Knowledge about possible core state evolutions must then take inspiration from magnetic records spanning epochs prior to the modern satellite era (observatory series, historical logs and paleomagnetic samples, e.g. Suttie et al 2025), as well as from the increasingly realistic numerical simulations of the geodynamo (Aubert et al 2022). Second, the core signal is hindered by external sources toward short periods. An accurate description of the rapid core variations thus relies on improved reconstructions of ionospheric sources (in particular toward auroral regions, see Kloss et al 2023) as well as magnetospheric fields patterns (e.g. Ou et al 2024).

Inverting for core field models requires incorporating additional information, to reduce the non-uniqueness inherent to the geomagnetic inverse problem. Under the assumption of Gaussian stochastic processes, this information is stored into a priori cross-covariance matrices. It may take the form of a smoothing (regularization, or damping) if one wishes to infer only an estimated model as free as possible from unresolved small length and time-scale patterns. Conversely, considering a realistic prior, allowing for the physically expected spatio-temporal variations, is necessary when one also targets realistic estimates of model uncertainties from the posterior covariance matrix (for a general formalism in the context of geomagnetism, see Gubbins and Bloxham 1985; Backus 1988).

It is with this spirit in mind that the COV-OBS suite of models has been derived (Gillet et al 2013; Huder et al 2020). These consider a priori information in time based on stochastic processes, namely, second-order auto-regressive or AR-2 processes. The rationale behind this choice is their link, in terms of differentiability properties, with geomagnetic series presenting jerks—defined as abrupt changes in the slope of the secular variation or SV, the rate of change of MF series. Such series could indeed be described as the realization of continuous and once differentiable processes, as are AR-2 processes. This

temporal prior information is incorporated in the COV-OBS inversion scheme via associated temporal cross-correlation functions. An alternative approach is taken in the KALMAG model by Baerenzung et al (2022), which is based on a sequential scheme, where AR-2 stochastic equations are embedded within a Kalman filter. Although AR-2 processes seem coherent with the temporal power spectral density (PSD) of the observed field on periods longer than a few years (Lesur et al 2018), they possibly allow too large variability on shorter time-scales, which is problematic in the presence of noise from imperfectly modelled external and related induced fields (Finlay et al 2017). Analysis of advanced geodynamo simulations shows the temporal PSD of the simulated field severely drops on periods shorter than the Alfvén time, associated with magnetic diffusion (Aubert and Gillet 2021). Scaled to the Earth's core conditions, this cutoff would indeed correspond to ≈ 2 yr (Gillet et al 2010).

In this context, we propose with the present study an alternative family of a priori temporal cross-covariances, based on third-order auto-regressive (or AR-3) processes. The PSDs from AR-3 processes have been considered previously by Sadhasivan and Constable (2022) to represent the behavior of the axial dipole field over a broad range of time-scales. This framework allows us to mimic the rapid decay in the PSD observed on short time-scales in dynamo simulations, mitigating possible spurious variability related to external field leakage when imposing a too loose constraint on rapid fluctuations. This is potentially important in the period band around 1–10 years that presents ambiguities between internal and external sources (e.g. Fig. 5 in Constable and Constable 2023). We use these new a priori cross-covariances to construct a geomagnetic field model over the period 1999.5–2024.5, together with its associated uncertainties. Constrained by satellite and ground-based magnetic records, this model named COV-SAT has been the basis for the calculation of ISTerre's candidates to the 14th edition of the International Geomagnetic Reference Field (for the previous 13th edition, see Alken et al 2021). Satellite data are incorporated by means of geomagnetic virtual observatories (Hammer et al 2021a). In comparison with the COV-OBS models, we not only incorporate a higher order temporal prior, but also follow Hellio and Gillet (2018) and project the model parameters on the auto-covariance functions when interpolating and extrapolating Gauss coefficients. By doing so, our temporal representation is based on geophysically motivated functions, avoiding a projection onto B-spline functions that may alter the high frequency behavior.

We begin by describing in §2 the geomagnetic observations used to build the field model and its associated uncertainties. Next in §3 we describe the model
parametrization, the derivation of the AR-3 a priori information, the inversion procedure and the interpolation/extrapolation method based on AR-3 correlation functions. In §4 we illustrate the fit to magnetic data, and present the time evolution of the field model, before we discuss in §5 some implications of our developments.

## 2 Data

We focus here on the satellite era (1999.5–2024.5), and consider ground observatories (GO) and geomagnetic virtual observatories (GVO). Three-components GO observations are incorporated in the form of annual differences of 4-monthly robust means, to avoid dealing with the crustal biases. The GO 4-monthly SV time series are derived from version 0140 (up to January 2024) of the ground observatory hourly means between January 1997 and April 2024 from the database AUX_OBS prepared by the British Geological Survey (Macmillan and Olsen 2013). Robust 4-monthly means of MF series are computed from hourly means, after subtracting predictions from the CM4 ionospheric and associated induced field (Sabaka et al 2004) and the CHAOS-7 magnetospheric and associated induced field (Finlay et al 2020). This latter uses as input the RC index (Olsen et al 2014). 4-monthly GO SV series are then computed as annual differences of GO MF time series.

Our model also incorporates three-components GVO field estimates from the CHAMP, Ørsted, GRACE, CryoSat-2 and Swarm missions each provided on a grid of 300 locations on an equal area grid at the mean height of the satellite. We refer to Hammer et al (2021a) for details about the GVO data processing algorithm. For all missions but Swarm, we consider 12-monthly means GVO, to avoid possible aliasing effects related to the slow drift in local time of polar orbiting satellites and the existence of remaining ionospheric fields that depend on local

time. Because the Swarm constellation covers local times at a higher cadence, with the Swarm Bravo satellite usually sampling a different local time to the Swarm Alpha and Charlie satellites, Swarm GVO MF data are incorporated as 4-monthly means. They are derived from Swarm L1B Mag-L OPER data versions 0602–0605, as available in June 2024. For each GVO data set we consider both the 4-monthly observed field estimates (built from night-side data with quiet time selection criteria applied but without any further corrections or cleaning, see Hammer et al 2021a) and the core field estimates (for which estimates of crustal, magnetospheric and ionospheric fields have been removed). For core field GVO data, estimates of the CIY4 ionospheric and induced fields have been removed (Sabaka et al 2018), as well as estimates of CHAOS-7 magnetospheric and induced fields. These data are also cleaned from the lithospheric field contributions for harmonic degrees 14–120, as estimated with the LCS-1 model (Olsen et al 2017). Furthermore, a denoising has been applied to remove remaining external and toroidal terms, via a spherical harmonic analysis. To mitigate the impact of field-aligned currents we consider, instead of three-components observations, intensity data $F = \sqrt{B_r^2 + B_\theta^2 + B_\phi^2}$ for GVO located at dipole latitudes $|\lambda_D| > 55^\circ$. We note throughout $(r, \theta, \phi)$ the radial, orthodradial and azimuthal spherical coordinates.

The observation operator for the cleaned GO SV and GVO MF (“core”) data is only related to our internal field model coefficients. The observation operator for the difference between observed and core fields GVO data is associated with the external field model coefficients. For each site separately, uncertainty estimates $(\sigma_r, \sigma_\theta, \sigma_\phi)$ for the GO series are obtained from the misfit to the CHAOS-7 field model. Uncertainties associated with the GVO series have been estimated as described in Hammer et al (2021a). Uncertainties on *F* data are obtained by

<table><caption>**Table 1** Top lines: median values of the observation uncertainties $\sigma_{r,\theta,\phi}$ considered in this study for “core” GO SV data (in nT/yr) and for core field GVO MF data (in nT)</caption><thead><tr><th>Data set<br>Units</th><th>GO SV<br>nT/yr</th><th>Ørsted<br>nT</th><th>CHAMP<br>nT</th><th>CryoSat-2<br>nT</th><th>GRACE<br>nT</th><th>Swarm<br>nT</th></tr></thead><tbody><tr><td>σ<sub>r</sub></td><td>2.36</td><td>2.99</td><td>2.30</td><td>2.71</td><td>3.13</td><td>2.39</td></tr><tr><td>σ<sub>θ</sub></td><td>1.94</td><td>2.55</td><td>1.90</td><td>2.09</td><td>2.28</td><td>2.01</td></tr><tr><td>σ<sub>φ</sub></td><td>2.11</td><td>1.58</td><td>1.28</td><td>1.29</td><td>1.95</td><td>1.16</td></tr><tr><td>N<sup>O</sup><sub>r</sub></td><td>9040</td><td>1078</td><td>2368</td><td>2026</td><td>2284</td><td>7333</td></tr><tr><td>N<sup>O</sup><sub>θ</sub></td><td>8982</td><td>1076</td><td>2366</td><td>2024</td><td>2279</td><td>7316</td></tr><tr><td>N<sup>O</sup><sub>φ</sub></td><td>9014</td><td>1077</td><td>2368</td><td>2026</td><td>2283</td><td>7331</td></tr><tr><td>N<sup>O</sup><sub>F</sub></td><td>0</td><td>222</td><td>436</td><td>440</td><td>490</td><td>1501</td></tr></tbody></table>

Bottom lines: number of data $N_{r,\theta,\phi,F}^O$ used to solve the inverse problem. The slight deviation in the number of data for (X, Y, Z) data is the result of our mild rejection criterion at 10σ (see §3.4). There are less *F* data, because these are used only in polar regions, while away from the poles, we only consider (X, Y, Z) data
propagating the three-components errors as $\sigma_F = \frac{|B_r|}{F}\sigma_r + \frac{|B_\theta|}{F}\sigma_\theta + \frac{|B_\phi|}{F}\sigma_\phi$. Since we get rid of the projection onto splines with 2 yr knots spacing, we do not add on top of these observation errors the modeling errors considered by Huder et al (2020) for the construction of COV-OBS.x2. We show in Table 1, for all components, the number of data used and the median uncertainties for core field GVO MF data from all missions, as well as for GO SV data. Uncertainties are generally lower on the east-west ($\phi$) component, which is less affected by external signals. These are as low as $\approx 1.3$ nT for CHAMP or CryoSat-2, and 1.15 nT for Swarm. This is to be compared with an uncertainty level closer to 2 nT in the north-south ($\theta$) component, and larger than 2.3 nT for $B_r$. Uncertainties for GVO based on platform magnetometers (Olsen et al 2020; Olsen 2021; Hammer et al 2021b) are significantly larger for GRACE, while for CryoSat-2, they are closer to error levels from dedicated magnetic missions.

# 3 Methods

## 3.1 Spatial parameterization of the magnetic model

The parameterization and methodology for COV-SAT only partly rely on the procedures described for earlier COV-OBS models (Gillet et al 2013, 2015; Huder et al 2020). We recall here the overall underlying framework and then discuss in more detail the key additions to this new model. We assume the mantle to be insulating. The magnetic field $\mathbf{B} = -\nabla V$ then derives from a potential for $r \ge r_c = 3485$ km, the Earth's core radius. The potential $V = V_i + V_e$ is separated between internal and external sources. The internal (core) field is expanded onto spherical harmonics as

$$V_i(r, \theta, \phi) = r_E \sum_{l=1}^{L} \left(\frac{r_E}{r}\right)^{l+1} \sum_{m=0}^{l} (g_l^m \cos(m\phi) + h_l^m \sin(m\phi)) P_l^m(\cos\theta), \quad (1)$$

with a truncation at spherical harmonic degree $L = 14$. Here $r_E = 6371.2$ km is the Earth's reference radius, $(g_l^m, h_l^m)$ are the Schmidt semi-normalized internal Gauss coefficients, and $P_l^m$ the associated Legendre polynomials. The external field potential

$$V_e(r, \theta, \phi) = r_E \sum_{l=1}^{L_e} \left(\frac{r}{r_E}\right)^l \sum_{m=0}^{l} (q_l^m \cos(m\phi) + s_l^m \sin(m\phi)) P_l^m(\cos\theta) \quad (2)$$

is restricted to a dipole aligned to the internal dipole (so $L_e = 1$), with $(q_l^m, s_l^m)$ the Schmidt semi-normalized

external Gauss coefficients. It is then described by a single coefficient $q_d$, such that

$$\begin{pmatrix} q_1^0 \\ q_1^1 \\ s_1^1 \end{pmatrix} = \frac{-q_d}{\sqrt{(g_1^0)^2 + (g_1^1)^2 + (h_1^1)^2}} \begin{pmatrix} g_1^0 \\ g_1^1 \\ h_1^1 \end{pmatrix}. \quad (3)$$

We consider the induced response of the core to external field variations as performed for COV-OBS-x2 (see Eqs. 5 and 15 in Huder et al 2020), by assuming that on long periods the fluid core behaves as a perfect conductor (see Olsen et al 2005). This assumption is supported by numerical modeling of the hydro-magnetic response of a rapidly rotating fluid core to an external excitation (Thea Lepage, pers. comm.)

One important modification w.r.t. COV-OBS models concerns the temporal projection of the model coefficients. Here we move away from the spline functions commonly used since Bloxham and Jackson (1992) to parameterize the time-dependence. We instead perform an interpolation at any observation epoch by projecting onto physically motivated correlation functions, following Hellio and Gillet (2018). These functions are described in §3.2, and the temporal projection in §3.3. The model inverted here is based on Gauss coefficient values specified every year over $[t_i, t_e] = [1999.5, 2024.5]$, so here a total of $N = 26$ epochs. The spherical harmonic model coefficients for all epochs and all degrees are stored into a vector $\mathbf{m}$, of size $P = N(L(L+2)+1) = 5850$.

## 3.2 Stochastic priors for the Gauss coefficients

### 3.2.1 A priori information for the core field

We represent internal field coefficients as a perturbation w.r.t. a background, or $g_l^m(t) = \bar{g}_l^m + \tilde{g}_l^m(t)$. We consider a background $\bar{g}_l^m = 0$ for all coefficients except the axial dipole. For this latter we choose from paleomagnetic reconstructions $\bar{g}_l^0 = -24000$ nT, the same value as used for building COV-OBS.x2 (see Huder et al 2020). In the construction of COV-OBS models, AR-2 processes have been considered to build the a priori temporal auto-covariance functions of Gauss coefficients, $C_l^m(\tau) = E[\tilde{g}_l^m(t)\tilde{g}_l^m(t+\tau)]$. This was justified based on the spectral index $p = 4$ found for the PSD $S(f) \propto f^{-p}$ of the magnetic field in the range $2 \lesssim 1/f \lesssim 70$ yr, with $f$ the frequency (Lesur et al 2018). However, advanced numerical simulations of the geodynamo (Aubert and Gillet 2021) indicate that $S(f)$ should steepen for periods shorter than a cutoff time-scale which, one scaled to Earth-like conditions, is found to be commensurate with the Alfvén time $\tau_a (\approx 2$ yr in Earth's core, see Gillet et al 2010). The cutoff for the 71p dynamo of Aubert and Gillet
<figure><img src="image_5.png" /><figcaption>**Fig. 1** Temporal PSDs for the Gauss coefficients from the 71p dynamo of Aubert and Gillet (2021). For the various degrees *l* shown, the PSDs have been averaged over all orders *m* ∈ [0, *l*], and normalized to the value at the period *T* = *τ*<sub>A</sub>. Periods on the *x*-axis are in unit of Alfvén time, which correspond to 5.8 yr in the 71p dynamo. Each series is 10 kyr long. The PSDs have been obtained with a multitaper method, applying a Hanning window on each of the 10 considered tapers</figcaption></figure>

(2021) is illustrated in Fig. 1, where we show the normal-
ized PSD for the magnetic Gauss coefficients series. It
clearly shows a more severe decay of the power toward
short periods, evolving approximately as S(f) ∝ f⁻⁶ and
independent of the spherical harmonic degree *l*. One may
wonder if this could be the consequence of the enhanced
dissipation of the velocity and codensity fields on short
length-scales in the 71p dynamo. Although this may
affect the subgrid induction, we recall that the magnetic
field itself is not subject to hyperdiffusion, motivated by
a magnetic field diffusing faster than momentum in liq-
uid metals (i.e. low magnetic Prandtl number). The cutoff
witnessed for periods less than τ<sub>a</sub> is, therefore, unlikely
the consequence of the enhanced dissipations.

The core field evolution on time-scales shorter than τ_a is difficult to constrain, due to ambiguities between internal and external sources (e.g. Lesur et al 2022), and the projection onto spline functions with 2 yr knot spacing in practice constrained rapid changes in previous COV-OBS models. Shortening the knot spacing with an AR-2 prior led to unsatisfying leakage from unmodelled external contributions on short periods, in part because the

AR-2 PSD allows considerable freedom on short time-
scales (Finlay et al 2017). Here we build a correlation
function that accounts for the steep decay of the mag-
netic field PSD at high frequencies by relying on higher
order (AR-3) stochastic processes (recalling that the PSD
for AR-n processes behave as S(f) ∝ f⁻²n toward short
time-scales). A linear AR-3 process X(t) is governed by a
stochastic equation of the form

$$
X''' + a_2 X'' + a_1 X' + a_0 X = a_z Z, \quad (4)
$$

where Z represents here a centered white noise of unit
variance. The determination of the 4 parameters a₀,₁,₂,₃
comes down to finding one variance and 3 time-scales, as
detailed below. Our approach is based on the shape of the
magnetic field PSD.

To mimic the spectrum of the core magnetic field, we
wish to build a correlation function, as simple as pos-
sible, such that S(f) describes 4 spectral ranges, with
indices *p* = 0, 2, 4 and 6 indicating low to high frequen-
cies. The cutoff between *p* = 0 and 2 is seen only for the
axial dipole from paleomagnetic records, for a period
τ<sub>d</sub> = 2π/ω<sub>d</sub> ≈ 50 – 100 kyr (Constable and Johnson
2005). Parameter estimations using stochastic models
indicate ω<sub>d</sub><sup>-1</sup> ≈ 10 kyr, linked to magnetic diffusion (Buf-
fett and Puranam 2017; Morzfeld and Buffett 2019). The
cutoff between *p* = 4 and 6 (in practice a gradual transi-
tion toward even steeper slopes) is deduced from its scal-
ing from advanced geodynamo simulations, at a period
about the Alfvén time τ<sub>a</sub> = 2π/ω<sub>a</sub> ≈ 2 yr and independ-
ent of the spherical harmonic degree *l* and order *m* (see
Fig. 1). The cutoff between *p* = 0 and 4 is seen for all
Gauss coefficients except for the axial dipole, at a period
commensurate with the SV time-scale

$$
\tau_{SV}(l) = \frac{\sum_{m=0}^{l} \sigma_{g_l^m}^2 + \sigma_{h_l^m}^2}{\sum_{m=0}^{l} \sigma_{\tilde{g}_l^m}^2 + \sigma_{\tilde{h}_l^m}^2}, \qquad (5)
$$

which can be related to the vortex turn-over time in the
core (Gillet et al 2013; Bouligand et al 2016). We use here
the notation
$$
\sigma_{g_l^m}^2 = \mathbb{E}[(g_l^m - \mathbb{E}[g_l^m])^2]
$$
and
$$
\sigma_{\tilde{g}_l^m}^2 = \mathbb{E}[(\partial_t g_l^m - \mathbb{E}[\partial_t g_l^m])^2],
$$
where $\mathbb{E}[\dots]$ denotes the
mathematical expectation, with a similar notation for the
$h_l^m$. Lhuillier et al (2011) shows that $\tau_{SV}(l)$ (in yr) evolves
approximately as $\approx 400/l$ for $l \ge 2$.

The above situation can be represented by a PSD of the
form
<figure><img src="image_6.png" /><figcaption>**Fig. 2** PSD associated with the a priori cross-correlation function considered for the axial dipole for the construction of COV-SAT, compared to previous editions</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3** PSD associated with the COV-SAT model a priori cross-correlation function considered for all coefficients but the axial dipole. The color indicates the coefficient degree, where dark blue is associated with the lowest degrees (starting at $l = 1$) and yellow the highest ones (up to $l = 14$)</figcaption></figure>

$$S(\omega) = \frac{K}{(\omega_d^2 + \omega^2)(\omega_u^2 + \omega^2)(\omega_a^2 + \omega^2)}, \quad (6)$$

with K a constant and $\omega = 2\pi f$ the angular frequency, and where the parameters $\omega_{d,u,a}$ have to be adjusted for all Gauss coefficients. We show in Appendix A how to derive the correlation functions associated with such a PSD. They take the simple analytical form (A5). We search for functions depending on as few parameters as possible. For all non-dipole coefficients, the parameters entering (6) are assumed to be independent of the order $m$, and vary only with the degree $l$. The axial dipole is treated separately because of its particular behavior presenting an intermediate spectral range with index $p = 2$.

<table><caption>**Table 2** Parameters for the correlation function (7) used as a prior on $q_d$ in the construction of COV-SAT</caption><thead><tr><th>σ<sub>2</sub><sup>2</sup> (nT<sup>2</sup>)</th><th>1/α<sub>2</sub> (years)</th><th>2π/ω<sub>2</sub> (years)</th><th>σ<sub>1</sub><sup>2</sup> (nT<sup>2</sup>)</th><th>1/α<sub>1</sub> (years)</th></tr></thead><tbody><tr><td>25</td><td>25</td><td>11.5</td><td>65</td><td>2.5</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 4** Auto-covariance function for the external dipole coefficient $\tilde{q}_d(t)$, as inferred from the RC index lowpass filtered for periods longer than 2 yr (see text for details)</figcaption></figure>

We present in Appendix B details of how we determine in practice the frequencies $\omega_{d,u,a}$ and the variances entering our correlation functions. There we also compare the parameters considered in the present study with those proposed by Sadhasivan and Constable (2022), that have been determined by adjusting the AR-3 PSD to a composite spectrum from models of the main field and of the dipole moment.

The associated PSD are illustrated in Fig. 2 for the axial dipole, and in Fig. 3 for the other coefficients. The a priori PSD for the axial dipole is close to that considered for COV-OBS-x2 on periods longer than $\approx \tau_a$. In comparison with Huder et al (2020), the update of the cutoff frequencies that define the transitions between period ranges with a spectral index $p = 0$ and $p = 2$, and then $p = 4$, slightly reduces the intermediate range of frequencies, where $S(f) \propto f^{-2}$, with a slight decrease in the power. We clearly see in Fig. 3 the decrease with the degree $l$ of the variance, as well as for the cutoff period between frequency ranges presenting a 0 and 4 spectral index. At all degrees, the cutoff on short periods toward a spectral index $p = 6$ is characteristic of AR-3 processes.

### 3.2.2 A priori correlation function for the external dipole field

We represent the external field coefficient as a perturbation w.r.t. a background, or $q_d(t) = \bar{q}_d + \tilde{q}_d(t)$. We consider $\bar{q}_d = 20$ nT, which approximates the sum of the Geocentric Solar Magnetospheric (GSM) and Solar Magnetic (SM) parts at zero frequency. The former is
almost constant on the time-scales considered here (see Fig. 1 of Lühr and Maus 2010), at a level of ≈ 8 nT. The latter can be represented via the RC index for the magnetospheric activity (Olsen et al 2014), and averages out to about 11 nT over the period 1997–2024. To estimate the auto-covariance function for $\tilde{q}_d(t)$ we consider, as suggested by Huder et al (2020), the combination of an AR-1 (i.e. Laplacian) and AR-2 (or damped oscillator) correlation functions:

$$C_{qd}(\tau) = \sigma_1^2 \exp(-\alpha_1 |\tau|) + \sigma_2^2 \exp(-\alpha_2 |\tau|) \left[ \cos(\omega_2 \tau) + \frac{\alpha_2}{\omega_2} \sin(\omega_2 \tau) \right]. \quad (7)$$

To determine the 5 parameters (2 variances and 3 times) entering Eq. (7), we rely on the RC index. Because we use satellite and observatory data built from annual or 4-monthly means, our estimate of $q_d$ cannot be considered as instantaneous, but instead as a weighted time-average. We build our prior based on this weighted average of the series $RC(t)$ over the period 1997–2024, after applying a low-pass filter with cutoff period 2 yr (the Nyquist period for the 1 yr sampling of our inverted model coefficients). In practice we use a Butterworth filter of order 2. The fit of Eq. (7) to the correlation for the filtered RC (using the *curve_fit* python function based on Vugrin et al 2007) gives $\sigma_1^2 = 64.9$ nT$^2$, $\sigma_2^2 = 23.1$ nT$^2$, $1/\alpha_1 = 2.6$ yr, $1/\alpha_2 = 23.6$ yr, and $2\pi/\omega_2 = 11.7$ yr, that we approximate to the values given in Table 2. The period $2\pi/\omega_2$ of the damped oscillator is reminiscent of the solar cycle. The amplitude of the AR-1 process is larger by a factor ≈ 1.6 than that of the AR-2 process. The fit to the RC correlation function is shown in Fig. 4.

### 3.3 Temporal representation of the Gauss coefficients

The model vector **m** contains Gauss coefficients ($^k g_l^m$, $^k h_l^m$, $^k q_d$) at the discrete set of N epochs $t_k$ considered to build the field model. The $t_k$ are stored in a vector **t**. We need a continuous representation for the model coefficients, e.g. $g_l^m(t)$ with t any epoch of interest, to estimate the model prediction at an observation epoch, or for interpolating and forecasting the field model. To this purpose we perform a Gaussian process regression (Rasmussen and Williams 2006), based on the auto-covariance functions defined in §3.2. We follow the approach of Hellio and Gillet (2018), considering each Gauss coefficient individually. We store in a vector **g**$_x$ its values to be interpolated and/or extrapolated at a set of epochs **t**$_x$, and in a vector **g**$_m$ its inverted values at the set of epochs **t**$_m$. Based on the auto-covariance functions defined above, we build for each coefficient the covariance matrices

$$\begin{gathered} P_{mm} = \mathbb{E}(\tilde{g}_m \tilde{g}_m^T), \quad P_{xx} = \mathbb{E}(\tilde{g}_x \tilde{g}_x^T) \\ \text{and } P_{xm} = \mathbb{E}(\tilde{g}_x \tilde{g}_m^T), \end{gathered} \quad (8)$$

which contain the a priori cross-covariances for all pairs of epochs within the vectors **t**$_x$ and **t**$_m$. We use the notation $\tilde{g}$ to represent the anomaly w.r.t. the background model $\bar{g}$. Within our setup, this latter is non-zero only for the internal axial dipole and the external dipole. Under Gaussian assumptions for the model distribution and uncertainties, the most likely estimate at epochs **t**$_x$ is then

$$\mathbf{g}_x = \bar{\mathbf{g}} + P_{xm}[P_{mm} + \mathbf{R}]^{-1}(\mathbf{g}_m - \bar{\mathbf{g}}), \quad (9)$$

where $\mathbf{R}$ is the error covariance matrix for the model parameters $\mathbf{g}_m$. Equation (9) is used to calculate the model predictions $H_i(\mathbf{m})$ to any datum $y_i$ at epoch $t_i$ (see §3.4). In this case we consider within the iterative inversion procedure $\mathbf{R} = 0$. We use this approach as an alternative to the spline expansion introduced by Bloxham and Jackson (1992).

Once the final model parameters have been derived, the same Eq. (9) is used to provide our model estimates and their associated uncertainties at any epochs, both for interpolating and for forecasting. In this latter case $\mathbf{R}$ is obtained from diagonal elements of the posterior model covariance matrix (the inverse of the Hessian, as performed by Gillet et al 2013). The above regression is applied separately for each coefficient, so that we avoid manipulating too large covariance matrices, thus reducing the risk of numerical instabilities associated with too wide ranges of eigenvalues.

On top of this estimate for the model coefficient series $\mathbf{g}_x$ at epochs $\mathbf{t}_x$, we also provide the associated uncertainties based on the posterior model error covariance matrix

$$P_{xx}^{\text{post}} = P_{xx} - P_{xm}[P_{mm} + \mathbf{R}]^{-1}P_{xm}^T. \quad (10)$$

To derive numerical values of the uncertainties on $\mathbf{g}_x$, we first perform the Cholesky factorization $P_{xx}^{\text{post}} = \mathbf{LL}^T$. From the the lower triangle matrix $\mathbf{L}$ we then generate an ensemble of $N^e = 500$ realizations

$$\mathbf{g}_{x,i} = \mathbf{L}\mathbf{w}_i \text{ for } i \in [1, N^e], \quad (11)$$

where $\mathbf{w}_i$ is a centered and unit variance normal random vector. The dispersion within this ensemble of series provides a useful measure of the model uncertainties through time.
### 3.4 Inversion procedure

As for previous editions of the COV-OBS models, COV-SAT Gauss coefficients have been obtained iteratively via a Newton-type method minimizing a cost function $J$ that consists of a sum of a misfit term $M$ and a term $N$ measuring the model complexity:

$$J(\mathbf{m}) = M(\mathbf{m}) + N(\mathbf{m}). \qquad (12)$$

One evolution w.r.t. earlier releases is the consideration of Huber weights for the measure of the residuals (data misfit) instead of a L2 norm previously. This avoids the use of a strict rejection criterion (previously $3\sigma$, now at $10\sigma$), and we now rely on an iterative weighted least-square method (Farquharson and Oldenburg 1998). The use of the Huber norm is a common way for dealing with outliers. The misfit to the data is then given by

$$M(\mathbf{m}) = \sum_{i=1}^{N^o} \mathcal{L}_{\delta}(\epsilon_i), \quad \text{with } \epsilon_i = |y_i - H_i(\mathbf{m})|/\sigma_i, \qquad (13)$$

where $H_i$ is the forward operator associated with the ith datum $y_i$, with estimated uncertainties $\sigma_i$. $N^o$ is the total number of observations. $\mathcal{L}_\delta$ is the Huber loss function defined as (e.g. Olsen 2002)

$$\mathcal{L}_{\delta}(\epsilon) = \begin{cases} \epsilon^2/2 & \text{for } |\epsilon| \le \delta \\ \delta|\epsilon| - \delta^2/2 & \text{otherwise,} \end{cases} \qquad (14)$$

with $\delta = 1.5$ the parameter controlling the normalized residual value above which $\mathcal{L}_\delta$ morphs from a L2 measure

(for low $|\epsilon|$) to a L1 measure. Observation errors are supposed independent.

The model complexity is measured as

$$\mathcal{N}(\mathbf{m}) = (\mathbf{m} - \bar{\mathbf{m}})^T \mathbf{P}^{-1} (\mathbf{m} - \bar{\mathbf{m}}). \qquad (15)$$

The matrix $\mathbf{P}$ contains, for all coefficients, the temporal cross-covariance between all epochs, as described in §3.2.1 and 3.2.2. We assume a priori that the coefficients of different orders or degrees are independent. $\bar{\mathbf{m}}$ is the background model (see above). The iterative process is stopped after 10 iterations, after which the relative change in the radial SV norm, averaged over the CMB of surface $\Sigma$ and over the full era $[t_i, t_e]$,

$$S(\mathbf{m}) = \frac{1}{\Sigma(t_e - t_i)} \int_{t_i}^{t_s} \oint_{\Sigma} |\partial_t B_r(c, \theta, \phi, t)|^2 d\Sigma dt, \qquad (16)$$

is less than $10^{-3}$ from one iteration to the next.

Note that the projection onto wide correlation functions (instead local functions such as of B-splines) comes with a numerical cost, as the filling of the matrices entering the iterative algorithm does not rely anymore on a banded storage. With global functions, the update of the Frechet derivatives for each datum requires one to fill all entries of the Hessian; this is the most time consuming part of the algorithm. Conversely, when using local functions a banded storage can be employed. Then Frechet derivatives have to be calculated only w.r.t. model parameters at knots entering a limited bandwidth $N^{bw}$ around the observation epoch. In practice one has $N^{bw} = j + 1$ for

<table><caption>**Table 3** Dimensional and dimensionless misfit and bias statistics for GO SV data and GVO ("core") MF data used when fitting the COV-SAT model. Stars indicate dimensionless values (see equ. 17)</caption><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Kind</th><th colspan="4">σ*</th><th colspan="4">μ*</th></tr><tr><th>X</th><th>Y</th><th>Z</th><th>F</th><th>X</th><th>Y</th><th>Z</th><th>F</th></tr></thead><tbody><tr><td>GO</td><td>SV</td><td>0.89</td><td>0.59</td><td>0.49</td><td>-</td><td>0.02</td><td>-0.02</td><td>-0.04</td><td>-</td></tr><tr><td>Ørsted</td><td>MF</td><td>0.98</td><td>0.69</td><td>0.69</td><td>0.48</td><td>0.18</td><td>0.02</td><td>0.07</td><td>0.11</td></tr><tr><td>CHAMP</td><td>MF</td><td>0.35</td><td>0.41</td><td>0.48</td><td>0.28</td><td>-0.01</td><td>0.01</td><td>0.01</td><td>-0.05</td></tr><tr><td>Cryosat-2</td><td>MF</td><td>0.54</td><td>0.52</td><td>0.77</td><td>0.51</td><td>0.09</td><td>0.00</td><td>0.12</td><td>0.29</td></tr><tr><td>GRACE</td><td>MF</td><td>0.61</td><td>0.72</td><td>0.71</td><td>0.32</td><td>0.19</td><td>0.02</td><td>0.03</td><td>0.04</td></tr><tr><td>Swarm</td><td>MF</td><td>0.43</td><td>0.36</td><td>0.53</td><td>0.30</td><td>-0.05</td><td>-0.02</td><td>-0.01</td><td>-0.09</td></tr></tbody></table><table><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Kind</th><th colspan="4">σ</th><th colspan="4">μ</th></tr><tr><th>X</th><th>Y</th><th>Z</th><th>F</th><th>X</th><th>Y</th><th>Z</th><th>F</th></tr></thead><tbody><tr><td>GO</td><td>SV</td><td>1.37</td><td>1.27</td><td>1.58</td><td>-</td><td>0.04</td><td>-0.04</td><td>-0.12</td><td>-</td></tr><tr><td>Ørsted</td><td>MF</td><td>1.92</td><td>1.63</td><td>1.68</td><td>1.47</td><td>0.37</td><td>0.08</td><td>0.14</td><td>0.33</td></tr><tr><td>CHAMP</td><td>MF</td><td>0.69</td><td>0.83</td><td>1.01</td><td>0.81</td><td>-0.16</td><td>0.02</td><td>0.02</td><td>-0.13</td></tr><tr><td>Cryosat-2</td><td>MF</td><td>1.18</td><td>1.12</td><td>1.58</td><td>1.41</td><td>0.24</td><td>0.03</td><td>0.33</td><td>0.92</td></tr><tr><td>GRACE</td><td>MF</td><td>1.51</td><td>1.81</td><td>1.88</td><td>1.20</td><td>0.50</td><td>0.07</td><td>0.03</td><td>0.19</td></tr><tr><td>Swarm</td><td>MF</td><td>0.81</td><td>0.77</td><td>1.00</td><td>0.80</td><td>-0.11</td><td>-0.04</td><td>-0.01</td><td>-0.24</td></tr></tbody></table>

See text for details. All dimensional values for $\sigma$ and $\mu$ are given in nT (resp. nT/yr) for the MF (resp. SV) data
B-splines of order j, reducing the number of operations by a factor ~ (N<sup>bw</sup>/N)<sup>2</sup> when filling in the matrices (neglecting the extra knots at endpoints, so in practice a bit less). In our configuration with N = 26, projecting onto order 6 B-splines (for the SA to be continuous, as for AR-3 processes) would reduce the computational cost by a factor about 10. Note that this reduction is valid when using the same knot spacing for both local and global functions, and that the projection of field changes with time-scales close to the knot spacing is sensitive to the choice of basis functions. This issue could be mitigated by reducing the knot spacing when using B-splines (6 months is used for the CHAOS-7 model by Finlay et al 2020), reducing the relative advantage of local basis functions in terms of computer cost, as one then needs to invert a Hessian of larger dimension (another potentially costly part of the algorithm). With our global correlation functions, the computational cost may become significant if using a much higher sampling rate of the model, or when covering a much longer time-span. On the other hand, when considering periods with lower observational coverage, the choice of geophysically motivated global functions may be beneficial, since model parameters at sparsely sampled eras are then influenced by data at other epochs (Hellio and Gillet 2018).

# 4 Results

## 4.1 Misfits and predictions to GO and GVO series

We provide in Table 3 some statistics regarding the COV-SAT predictions for the GO SV and GVO MF data. These are produced separately for all components and for all satellite missions. The diagnostics are provided for the X (northward), Y (eastward) and Z (downward) components

of the field, as well as for its intensity F. We introduce for each subset of data both dimensional and normalized measures of the misfits and biases. For this purpose we fit a Huber distribution (see Eq. 14)

$$p_{\delta}(\epsilon) = \frac{1}{k_{\delta}\sigma^{*}} \exp\left(-\mathcal{L}_{\delta}\left(\frac{\epsilon - \mu^{*}}{\sigma^{*}}\right)\right) \quad (17)$$

to the normalized histograms of the dimensionless residuals $\epsilon$, with the dimensionless bias $\mu^*$ and misfit $\sigma^*$ being the two adjustable parameters. The normalisation constant $k_{\delta} = \sqrt{2\pi}\operatorname{erf}(\delta/\sqrt{2}) + \frac{2}{\delta}\exp(-\delta^2/2) \approx 0.3839$

ensures the condition $\int_{-\infty}^{+\infty} p_{\delta}(\epsilon)d\epsilon = 1$ for a probability density function (PDF). The dimensional diagnostics are obtained similarly from histograms of the dimensional residuals, by adjusting a dimensional bias $\mu$ and misfit $\sigma$. The adjusted misfit values are given in Table 3 for “core” field data and in Table 4 for “observed” minus “core” field data.

For “core” data, apart from Ørsted (X and F), GRACE (X) and Cryosat-2 (Z and F) data, dimensionless biases are less than 0.1 in absolute value. For all components, GO as well as CHAMP and Swarm data, that bring the most important observational constraint, show almost centered distributions for normalized residuals. Dimensionless residuals are generally less than unity, denoting a possible over-estimation of the a priori data uncertainties, in particular for CHAMP and Swarm (recall that these uncertainties do not refer directly to the satellite measurements, but to a local potential field model fit to these for a given time window). The lower dimensional

<table><caption>**Table 4** Same as Table 3 for GVO “observed” minus “core” MF data, associated with our external field model and its induced counterpart</caption><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Kind</th><th colspan="4">σ*</th><th colspan="4">μ*</th></tr><tr><th>X</th><th>Y</th><th>Z</th><th>F</th><th>X</th><th>Y</th><th>Z</th><th>F</th></tr></thead><tbody><tr><td>Ørsted</td><td>MF</td><td>0.67</td><td>0.38</td><td>0.42</td><td>0.52</td><td>0.18</td><td>-0.06</td><td>0.04</td><td>0.15</td></tr><tr><td>CHAMP</td><td>MF</td><td>0.90</td><td>0.51</td><td>0.75</td><td>0.60</td><td>0.17</td><td>0.00</td><td>0.04</td><td>0.30</td></tr><tr><td>Cryosat-2</td><td>MF</td><td>0.47</td><td>0.40</td><td>0.35</td><td>0.37</td><td>-0.05</td><td>-0.10</td><td>-0.02</td><td>0.29</td></tr><tr><td>GRACE</td><td>MF</td><td>0.30</td><td>0.41</td><td>0.30</td><td>0.19</td><td>-0.06</td><td>0.04</td><td>0.01</td><td>0.13</td></tr><tr><td>Swarm</td><td>MF</td><td>1.28</td><td>0.59</td><td>0.77</td><td>0.84</td><td>0.23</td><td>0.03</td><td>0.00</td><td>0.23</td></tr></tbody><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Kind</th><th colspan="4">σ</th><th colspan="4">μ</th></tr><tr><th>X</th><th>Y</th><th>Z</th><th>F</th><th>X</th><th>Y</th><th>Z</th><th>F</th></tr></thead><tbody><tr><td>Ørsted</td><td>MF</td><td>3.12</td><td>2.23</td><td>2.75</td><td>3.13</td><td>0.84</td><td>-0.44</td><td>0.00</td><td>1.01</td></tr><tr><td>CHAMP</td><td>MF</td><td>3.20</td><td>2.12</td><td>3.07</td><td>2.70</td><td>0.62</td><td>0.02</td><td>-0.12</td><td>1.60</td></tr><tr><td>Cryosat-2</td><td>MF</td><td>2.65</td><td>2.51</td><td>2.23</td><td>2.42</td><td>-0.32</td><td>-0.66</td><td>-0.20</td><td>2.32</td></tr><tr><td>GRACE</td><td>MF</td><td>2.54</td><td>4.25</td><td>3.04</td><td>2.10</td><td>-0.57</td><td>0.46</td><td>0.06</td><td>1.54</td></tr><tr><td>Swarm</td><td>MF</td><td>3.82</td><td>2.03</td><td>2.48</td><td>3.10</td><td>0.71</td><td>0.12</td><td>-0.11</td><td>0.94</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Fig. 5** Histogram of the normalized misfits $\epsilon$ for the various GO and GVO data sets, split into spatial components. The cyan curve represents the Huber distribution fitted to each histogram, with parameters $\mu^*$ and $\sigma^*$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Comparison between the SV observed at Honolulu (HON, Eastern Pacific) and Hermanus (HER, South Africa), and the predictions from COV-SAT and CHAOS-7, for the three components. Shaded areas represents the spread within ±1σ</figcaption></figure>

<figure><img src="image_11.png" /></figure>

misfits are found for CHAMP and Swarm, where it is less than 1 nT for X, Y and F, and only slightly above 1 nT for Z data. All components from Cryosat-2 data are fitted within 1.6 nT. The larger misfits values are found for Ørsted and GRACE GVO vector data. Ground SV data misfits are less than 1.6 nT/yr, the weakest for Y and the largest for Z data.

Dimensional misfits for “observed” minus “core” field data, used to estimate our external field model, show relatively larger values, ranging from 2 to 4.4 nT depending on the mission and component. The distinction between the missions is less obvious than for “core”

field data, due to the crude external model considered in this study. Normalized misfits are nevertheless less than unity (except, for Swarm X data), in particular for the less accurate Ørsted, CryoSat-2 and GRACE data, suggesting that a priori error uncertainties considered in our inverse problem are possibly over-estimated. We note also for all satellite missions significant biases for F data used in the polar region (and to a lesser extent on the north-south X component used at nonpolar latitudes).

In addition, we represent in Fig. 5 for all missions and components the distribution of normalized residuals ε, superimposed with the best fitting Huber distribution.
<figure><img src="image_12.png" /><figcaption>**Fig. 7** Comparisons between the MF observed at two Swarm GVO and the predictions from COV-SAT and CHAOS-7, for the three components. Shaded areas represent the spread within ±1σ</figcaption></figure>

Most histograms behave according to our assumptions on the residuals distribution. In some cases the PDF behaves closer to a Laplacian (e.g. for the *X* component for CHAMP and Swarm), with a sharper peak around 0. Likewise, some asymmetries are seen for some PDFs for *F* data.

We further illustrate the ability of our model to fit the observations in Fig. 6, where we present SV series at the Honolulu (158°W, 21°N, Eastern Pacific) and Hermanus (19°E, 34°S, South Africa) ground stations, and in Fig. 7, where we present MF series at nearby GVO positions. Our model closely fits the observed series, with predictions comparable to those from the CHAOS-7 field model.

## 4.2 Time evolution of internal field model coefficients

We next illustrate in Fig. 8 the evolution of some SV coefficients obtained from COV-SAT, in comparison

with alternative models CHAOS-7, KALMAG and COV-OBS.x2. The estimation is made using $t_x$ sampling 2000–2030 every 0.2 yr, using the method described in §3.3. We note some significant changes in the rapid evolution of the axial dipole between COV-SAT and COV-OBS.x2. Our new model is closer to CHAOS-7, and generally to KALMAG (with the exception of the period from 2010 to 2014, where CHAOS-7 and KALMAG differ the most, a consequence of CHAOS-7 using Cryosat-2 data and KALMAG lacking satellite data). We believe this change is due to our revised setup now being based on “core” MF GVO and SV GO data that have been cleaned as much as possible from external sources (contrary to COV-OBS.x2 for which magnetospheric corrections had not been applied). For other coefficients at low to moderate degrees *l*, the series of the several considered models largely overlap. Larger discrepancies are seen toward
<figure><img src="image_13.png" /><figcaption>**Fig. 8** SV time series for several Gauss coefficients showing COV-SAT, compared with COV-OBS.x2, CHAOS-7 and KALMAG. For models provided with uncertainties, shaded areas represent the spread within ±1σ</figcaption></figure>

small length-scales (see for instance $\dot{g}_{12}^5$, although the trend for CHAOS-7, KALMAG and COV-SAT are overall similar, with differences most often within the estimated errorbars.

Models based on stochastic priors can be used for forecasting possible evolutions after the date of the

latest observation. We observe for COV-SAT that the spread of the SV forecast behaves smoothly once the observational constraint is relaxed after $t_e = 2024.5$. This is characteristic of AR-3 processes that are the basis of the temporal variations considered in this study, for which the SV is continuous and differentiable,
<figure><img src="image_14.png" /><figcaption>**Fig. 9** Top: MF and SV spectra for models COV-SAT, COV-OBS.x2, CHAOS-7 and KALMAG in 2018. Dashed lines correspond to the spectra for the model uncertainties (spread within ±1σ). Bottom: MF and SV spectra for the difference between field models, with spectra for the model uncertainties again shown as dashed lines</figcaption></figure>

with a spread evolving ∝ (t - t_e)³/². It contrasts with AR-2 processes for which the SV is continuous but not differentiable, which translates into a spread that increases sharply (∝√(t - t_e)), as observed for instance for the KALMAG model in Fig. 8.

We further compare the proposed uncertainties and the model differences in Fig. 9 (top), with MF and SV spectra at the Earth's surface for the several models in 2018 (qualitatively similar behaviors would be witnessed at other epochs). Spectra for the mean models almost superimpose. Our estimate of the MF errors is lower than those provided for KALMAG, in particular toward high degrees. SV uncertainties are on the contrary very similar between the two models (only slightly less for COV-SAT at low degrees). Toward low (resp. high) degrees

our estimates of the SV uncertainties are revised upward (resp. downward) compared with COV-OBS.x2. This is partly the consequence of removing the projection onto splines with 2 yr knot spacing. Indeed, high degree coefficients are hardly accessible on short time-scales, while there is potentially some observational constraint on low degrees (e.g. Ropp et al 2020). Because the projection onto splines in COV-OBS models filters period shorter than 3–4 years, the estimated uncertainty level was characteristic of a low-pass version of the model coefficients, which is by essence less than for instantaneous Gauss coefficients.

Spectra for the proposed SV model uncertainties and for the differences within field models are comparable (see Fig. 9, bottom), which suggests our SV error estimate
correlation functions characteristic of order 2 autoregressive processes and projected onto B-spline functions (Gillet et al 2013; Huder et al 2020).

The stochastic framework considered here allows us to propose errorbars on the model parameters. The magnitude of the proposed SV errors is comparable to that of the difference between various models proposed by the community, which is encouraging for the purpose of geomagnetic data assimilation. Our measure of the SV model uncertainties is comparable to that provided with the KALMAG model (a bit less at the largest lengthscales) despite the use of different input data and modelling strategies. The KALMAG model is developed by time-stepping AR-2 stochastic equations in a sequential algorithm, whereas in our inversion, the prior information enters via the correlation function. There is no barrier to the incorporation of the AR-3 prior that we propose into a sequential scheme such as employed by Ropp et al (2020) or Baerenzung et al (2022). We indeed propose in Appendix A the corresponding stochastic equation. It will be worthwhile in the future to consider the AR-3 prior that we propose here for building a field model from observatory and satellite observations sampled at a higher cadence, to see if the transition from −4 to −6 spectral at a period about 2 yr is also found in the posterior PSD for the Gauss coefficients.

In the meantime the mean model we derived, estimated from ground-based and virtual observatory data cleaned as much as possible for external sources, is close to alternative models such as the well-accepted CHAOS-7 (Finlay et al 2020) or KALMAG (Baerenzung et al 2022) models, where data along satellite tracks are considered, and that involve the co-estimation of a complex external field model. This indicates that GVO data produced following Hammer et al (2021a) provide useful observational constraints that allows one to infer the evolution of the core magnetic field from a reduced number of higher level post-processed magnetic data. The sampling at which we can consider GVO data is limited by the drift rate of the satellite orbits in local time. The Swarm constellation configuration with two well-separated polar tracks reduces it to ≈ 4 months, which is less than that needed with a single satellite such as CHAMP. The inclined orbit mission MSS-1 reduces it further down to 2 months, covering latitudes less than 41° (Jiang et al 2024). We can thus hope to get insight on more rapid core field changes with the current configuration (MSS-1 plus Swarm), and the up-coming NanoMag-Sat mission from ESA with 3 spacecrafts (one polar orbit and two inclined at 60°, see Deconinck et al 2025). The

production of GVO data from these new combinations of satellites will then be helpful to isolate core processes on periods shorter than a few years (Brown et al 2023).

This concerns first the presence of Magneto-Coriolis waves with short periods, as well as long period inertial waves. Such transient dynamics is expected from studies of eigenvalue and initial value problems in a rapidly rotating sphere in the presence of a non-axisymmetric imposed magnetic field (Gerick and Livermore 2024; Barrois and Aubert 2024). It also concerns identifying the actual high frequency cutoff for the magnetic signal, currently hindered by the dominant external sources on short time-scales, and only witnessed in geodynamo simulations. Extending to millenial timescales dynamo simulation series at Earth-like parameters, such as that by Aubert (2023), will be helpful to more accurately constrain the shape of the temporal spectrum at periods close to the Alfvén time. Detecting a decay in the PSD of the core signal is challenging but may help constrain physical properties of the Earth’s deep interior. We think here of bounds on the magnetic diffusivity η of the core material, to which the cutoff is related (Aubert and Gillet 2021). We recall that η (and the associated thermal conductivity) is a subject of debate within the community of mineral physics (e.g. Davies et al 2015; Zhou et al 2022). We also have in mind bounds on the electrical conductivity σm of the deep mantle, a quantity inaccessible using electro-magnetic sounding with external field variations as the source (Grayver 2024). Indeed, while σm ≠ 0 over a thick layer will act as a non-causal low-pass filter on the magnetic signal originating from the core (Jault 2015), σm ≠ 0 at the base of the mantle will alter the transient dynamics within the core through the electromagnetic boundary condition (Schaeffer and Jault 2016; Firsov et al 2023).

Appendices

From the PSD to the correlation function

We derive here an AR-3 correlation function that satisfies the observational and numerical constraints from §3.2. Below we rely on the properties of AR-1 processes with variance σ02 and decay frequency ω0, which are characterized by a auto-covariance function

$$
C(\tau) = \sigma_0^2 \exp(-\omega_0 |\tau|) \quad (A1)
$$
and a PSD (e.g. Yaglom 2004)

$$S(\omega) = \frac{\sigma_0^2 \omega_0 / \pi}{\omega_0^2 + \omega^2}. \qquad (A2)$$

First we perform a decomposition in simple fractions of a PSD of the form (6). Consider three constants (a, b, c), and three frequencies (α, β, γ). Lets choose a PSD such that

$$S(\omega) = \frac{1}{\pi} \left( \frac{a\alpha}{\alpha^2 + \omega^2} + \frac{b\beta}{\beta^2 + \omega^2} + \frac{c\gamma}{\gamma^2 + \omega^2} \right) \qquad (A3) \\ = \frac{k_0 + k_2\omega^2 + k_4\omega^4}{(\alpha^2 + \omega^2)(\beta^2 + \omega^2)(\gamma^2 + \omega^2)},$$

where

$$\left\{ \begin{array}{l} \pi k_4 = a\alpha + b\beta + c\gamma \\ \pi k_2 = a\alpha(\beta^2 + \gamma^2) + b\beta(\alpha^2 + \gamma^2) + c\gamma(\alpha^2 + \beta^2) \\ \pi k_0 = a\alpha\beta^2\gamma^2 + b\beta\alpha^2\gamma^2 + c\gamma\alpha^2\beta^2 \end{array} \right. . \qquad (A4)$$

Using Eqs. (A1 and A2), the correlation function associated with the PSD (A3) takes the simple form

$$C(\tau) = a \exp(-\alpha|\tau|) + b \exp(-\beta|\tau|) + c \exp(-\gamma|\tau|). \qquad (A5)$$

For $S(\omega)$ to evolve as $\omega^{-6}$ toward high frequencies, we deduce the two conditions $k_2 = k_4 = 0$. Imposing the variance $\sigma^2 = C(0)$ of the signal leads to a third condition $\sigma^2 = a + b + c$. Knowing the values for the three frequencies (α, β, γ) and $\sigma^2$, we then have to solve the linear problem $Mx = y$ where

$$M = \begin{bmatrix} 1 & 1 & 1 \\ \alpha & \beta & \gamma \\ \alpha(\beta^2 + \gamma^2) & \beta(\alpha^2 + \gamma^2) & \gamma(\alpha^2 + \beta^2) \end{bmatrix}, \qquad (A6) \\ x = \begin{bmatrix} a \\ b \\ c \end{bmatrix} \quad \text{and} \quad y = \begin{bmatrix} \sigma^2 \\ 0 \\ 0 \end{bmatrix}.$$

We have

$$\det(M) = \gamma\beta(\beta^2 - \gamma^2) - \alpha\gamma(\alpha^2 - \gamma^2) + \alpha\beta(\alpha^2 - \beta^2), \qquad (A7)$$

and the solution is

$$\left\{ \begin{array}{l} a/\sigma^2 = \frac{\gamma\beta(\beta^2 - \gamma^2)}{\det(M)} \\ b/\sigma^2 = \frac{-\alpha\gamma(\alpha^2 - \gamma^2)}{\det(M)} \\ c/\sigma^2 = \frac{\alpha\beta(\alpha^2 - \beta^2)}{\det(M)} \end{array} \right. . \qquad (A8)$$

From these values we deduce the auto-covariance functions (A5), that we use to build the prior covariance matrix when building our field model.

For all coefficients but the axial dipole, we could resort to only 3 spectral ranges of indices $p = 0, -4$ and $-6$. However, we cannot use $\beta = \alpha$, since in this case $c = 0$ and then the PSD will behave as $\omega^{-4}$ on short time-scales. However, this pathological case can be avoided by having $\alpha$ close to $\beta$, a bit larger.

The PSD and auto-covariance functions obtained this way are characteristic of the AR-3 process represented as

$$X''' + (\alpha + \beta + \gamma)X'' + (\alpha\beta + \beta\gamma + \gamma\alpha)X' + \alpha\beta\gamma X = Z, \qquad (A9)$$

where Z is a centered white noise of variance

$$k_0 = \frac{\alpha\beta\gamma\sigma^2}{\pi \det(M)} (\beta^2\gamma^2(\beta^2 - \gamma^2) - \alpha^2\gamma^2(\alpha^2 - \gamma^2) + \alpha^2\beta^2(\alpha^2 - \beta^2)). \qquad (A10)$$

For the parameters chosen here, the roots of the denominator of (A3) being complex conjugates, a process such as (A9) is stationary. Note that we may simplify Eq. (A8) in the limit $\gamma \gg \beta \gg \alpha$ (as for the axial dipole), in which $\det(M) \simeq -\beta\gamma^3$, leading to

$$\left\{ \begin{array}{l} a/\sigma^2 \simeq 1 \\ b/\sigma^2 \simeq -\frac{\alpha}{\beta} \\ c/\sigma^2 \simeq \frac{\alpha\beta^2}{\gamma^3} \end{array} \right. . \qquad (A11)$$

Then one has $k_0 \simeq \sigma^2\alpha\beta^2\gamma^2/\pi$.

### Choice of the input parameters

The AR-3 correlation functions designed above rely on the choice of some parameters. We consider here $\alpha^{AD} = \omega_d = 10^{-4} \text{ yr}^{-1}$ for the axial dipole. We also fix $\gamma = \omega_a = 2\pi/\tau_a = \pi \text{ yr}^{-1}$ for all coefficients. We now wish to choose the frequency $\beta(n) = \omega_u$. We define $\omega_u(n)$ as performed for previous COV-OBS models, i.e. based on the SV time-scale (5). The reason is that $\tau_{SV}(n)$ is a quantity accessible from geophysical observations via magnetic field models. To this purpose we use the relationship (e.g. Hulot and Le Mouël 1994)

$$\mathbb{E}(\dot{X}^2) = - \frac{d^2 C(\tau)}{d\tau^2} \Big|_{\tau=0}, \qquad (B1)$$

for a process $X(t)$ with auto-covariance function $C(\tau)$. For a function such as that defined by Eq. (A5), we then have
$$ \sigma_{g_n}^2 = -(a\alpha^2 + b\beta^2 + c\gamma^2). \qquad (B2) $$

Using (A8) in the limit where $γ \gg α, β$ we get

$$ \alpha\beta(n) = 1/\tau_{SV}^2(n). \qquad (B3) $$

For all coefficients but the axial dipole, if we assume as above the condition $c\alpha = \beta$ with c a bit larger than one (to reduce the frequency range presenting a -2 spectral index), this gives

$$ \beta(n) = c\alpha(n) = \omega_u(n) = \frac{\sqrt{c}}{\tau_{SV}(n)}. \qquad (B4) $$

This condition is close to that introduced previously by Gillet et al (2013) for the cutoff frequency of Matern AR-2 processes. For the axial dipole (noted with a subscript 'AD') we obtain instead

$$ \omega_d \beta^{\text{AD}} = 1/\tau_{\text{SV}}^{\text{AD}} \qquad (B5) $$

where $τ_{SV}^{AD} = \sqrt{\sigma_{g_1}^2/\sigma_{g_1^0}^2}$. A numerical application with $1/\omega_d = 10$ kyr, $\sigma_{g_1^0} = 10$ nT/yr and $\sigma_{g_1} = 7700$ nT (or $τ_{SV}^{AD} = 7700$ yr, see Huder et al 2020) leads to $1/\beta^{AD} \approx 60$ yr $\gg \tau_a$. We check a posteriori that $\omega_d \ll \beta^{AD} \ll \omega_a$.

A similar AR-3 formalism was considered by Sadhasivan and Constable (2022) for the dipole moment. By fitting the analytical AR-3 spectrum to a composite spectrum of historical and paleomagnetic models derived from magnetic observations, they obtained $1/\omega_d \simeq 15$ kyr, $1/\beta^{AD} \simeq 30$ yr and $1/\omega_a \simeq 2$ yr. Their estimate of $\omega_d$ (resp. $\beta^{AD}$) is comparable to ours, only 1.5 (resp. 2) times smaller. The main difference concerns their value for $\omega_a$ that is smaller by a factor about 6. This is because the high frequency cutoff is constrained in their approach by the PSD from a historical field model that does not capture variations on periods as short as the Alfvén time $\tau_a$.

### Acknowledgements

We thank the two anonymous referees for their useful comments that helped improve the quality of our manuscript. We thank ESA for the prompt availability of *Swarm* L1b data and for providing access to the CryoSat-2 platform magnetometer data and related engineering information. The support of the CHAMP mission by the German Aerospace Center (DLR) and the Federal Ministry of Education and Research for the CHAMP mission is gratefully acknowledged. The Ørsted Project was made possible by extensive support from the Danish Government, NASA, ESA, CNES, DARA and the Thomas B. Thriges Foundation. We thank JPL and GFZ for providing GRACE L1b data. The staff of the geomagnetic observatories and INTERMAGNET are thanked for supplying high-quality observatory data.

### Author Contributions

The numerical developments have been implemented and run by RC. The theoretical part about the AR-3 prior has been performed by NG, with the help of POA. VO data have been operated and provided by CF. The manuscript has been primarily written by NG, with inputs from all co-authors.

### Funding

This work has been funded by ESA in the framework of EO Science for Society, through contract 4000127193/19/NL/IA (Swarm + 4D Deep Earth: Core). This project was financially supported by CNES as an application of the Core_Swarm mission. RC and NG are part of Labex OSUG@2020 (ANR10 LABX56).

### Data availability

● GO 4-monthly SV time series are available at http://www.spacecenter.dk/files/magnetic-models/GO/● The database AUX OBS is accessible from https://auxobs-api-uat.bgs.ac.uk/● MF GVO data for all satellite missions are available at http://www.spacecenter.dk/files/magnetic-models/GVO/● The field model generated within the current study are accessible from https://geodyn.univ-grenoble-alpes.fr/

### Declarations

### Competing interests

The authors declare having no Conflict of interest.

### Author details

¹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTER, 38000 Grenoble, France. ²CNRS, Univ. Grenoble Alpes, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. ³DTU Space, Technical University of Denmark, Centrifugevej 356, Kongens Lyngby, Denmark.

Received: 11 February 2025 Accepted: 26 June 2025

Published online: 30 July 2025

### References

Alken P, Thébault E, Beggan CD, Amit H, Aubert J, Baerenzung J, Bondar T, Brown W, Calif S, Chambodut A et al (2021) International geomagnetic reference field: the thirteenth generation. Earth, Planets and Space 73:1–25

Aubert J (2014) Earth's core internal dynamics 1840–2010 imaged by inverse geodynamo modelling. Geophys J Int 197(3):1321–1334

Aubert J (2023) State and evolution of the geodynamo from numerical models reaching the physical conditions of earth's core. Geophys J Int 235(1):468–487

Aubert J, Gillet N (2021) The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions. Geophys J Int 225(3):1854–1873

Aubert J, Livermore PW, Finlay CC, Fournier A, Gillet N (2022) A taxonomy of simulated geomagnetic jerks. Geophys J Int 231(1):650–672

Backus GE (1988) Bayesian inference in geomagnetism. Geophys J Int 92(1):125–142

Baerenzung J, Holschneider M, Wicht J, Sanchez S, Lesur V (2018) Modeling and predicting the short-term evolution of the geomagnetic field. J Geophys Res Solid Earth 123(6):4539–4560

Baerenzung J, Holschneider M, Saynisch-Wagner J, Maik T (2022) Kalmag: a high spatio-temporal model of the geomagnetic field. Earth, Planets and Space 74(1):139

Barrois O, Aubert J (2024) Characterization of hydromagnetic waves propagating over a steady, non-axisymmetric background magnetic field. In: Proc. R. Soc. A. The Royal Society, vol 480, p 20240067

Bloхam J, Jackson A (1992) Time-dependent mapping of the magnetic field at the core-mantle boundary. J Geophys Res Solid Earth 97(B13):19,537–19,563

Bouligand C, Gillet N, Jault D, Schaeffer N, Fournier A, Aubert J (2016) Frequency spectrum of the geomagnetic field harmonic coefficients from dynamo simulations. Geophys J Int 207(2):1142–1157

Brown W, Gomez Perez N, Beggan C, Hammer M, Finlay C, Hulot G (2023) The potential of geomagnetic virtual observatories with nanomagsat. In: XXVIII General Assembly of the International Union of Geodesy and Geophysics (IUGG), GFZ German Research Centre for Geosciences

Buffett B, Puranam A (2017) Constructing stochastic models for dipole fluctuations from paleomagnetic observations. PhysEarth Planet Int 272:68–77
Bullard EC, Freedman C, Gellman H, Nixon J (1950) The westward drift of the Earth's magnetic field. Phil Trans R Soc A 243:67–92

Constable C, Constable S (2023) A grand spectrum of the geomagnetic field. Phys Earth Planet Int 344(107):090

Constable C, Johnson C (2005) A paleomagnetic power spectrum. Phys-Earth Planet Int 153(1–3):61–73

Davies C, Pozzo M, Gubbins D, Alfe D (2015) Constraints from material properties on the dynamics and evolution of earth's core. Nat Geosci 8(9):678–685

Deconinck F, Hulot G, Léger JM, Clausen LB, Pastena M, Lejault JP, Wilkinson J, Caballero DM, Tostado MV, Coísson P, et al (2025) Nanomagsat status: a 3x16u low-earth orbit constellation to monitor the earth magnetic field and the ionospheric environment. In: Small Satellites Systems and Services Symposium (4S 2024), SPIE, vol 13546, pp 1086–1100

Farquharson CG, Oldenburg DW (1998) Non-linear inversion using general measures of data misfit and model structure. Geophys J Int 134(1):213–227. https://doi.org/10.1046/j.1365-246x.1998.00555.x

Finlay C, Lesur V, Thébault E, Vervelidou F, Morschhauser A, Shore R (2017) Challenges handling magnetospheric and ionospheric signals in internal geomagnetic field modelling. Space Sci Rev 206(1):157–189

Finlay C, Kloss C, Olsen N, Hammer M, Tøffner-Clausen L, Grayver A, Kuvs-hinov A (2020) The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic Anomaly. Earth Planets Space 72(1):1–31

Firsov I, Jault D, Gillet N, Aubert J, Mandea M (2023) Radial shear in the flow at the Earth's core surface. Geophys J Int 235(3):2524–2539

Gerick F, Livermore P (2024) Interannual Magneto-Coriolis modes and their sensitivity on the magnetic field within the Earth's core. Proc R Soc A. https://doi.org/10.1098/rspa.2024.0184

Gillet N, Jault D, Canet E, Fournier A (2010) Fast torsional waves and strong magnetic field within the Earth's core. Nature 465(7294):74–77

Gillet N, Jault D, Finlay C, Olsen N (2013) Stochastic modeling of the earth's magnetic field: inversion for covariances over the observatory era. Geochem Geophys Geosyst 14(4):766–786

Gillet N, Barrois O, Finlay CC (2015) Stochastic forecasting of the geomagnetic field from the COV-OBS.x1 geomagnetic field model, and candidate models for IGRF-12. Earth, Planets and Space 67:1–14

Gillet N, Huder L, Aubert J (2019) A reduced stochastic model of core surface dynamics based on geodynamo simulations. Geophys J Int 219(1):522–539

Gillet N, Dall'Asta F, Amblard PO, Claveau R, Aubert J (2024) Waves in earth's core and geomagnetic field forecast. Phys Earth Planet Int 357(107):284

Grayver A (2024) Unravelling the electrical conductivity of earth and planets. Surv Geophys 45(1):187–238

Gubbins D, Bloxham J (1985) Geomagnetic field analysis-III. Magnetic fields on the core-mantle boundary. Geophys J R Astr Soc 80(2):695–713

Hammer MD, Cox GA, Brown WJ, Beggan CD, Finlay CC (2021) Geomagnetic virtual observatories: monitoring geomagnetic secular variation with the swarm satellites. Earth Planets Space 73(1):1–22

Hammer MD, Finlay CC, Olsen N (2021) Applications for cryosat-2 satellite magnetic data in studies of earth's core field variations. Earth, Planets and Space 73:1–22

Hellio G, Gillet N (2018) Time-correlation-based regression of the geomagnetic field from archeological and sediment records. Geophys J Int 214(3):1585–1607

Huder L, Gillet N, Finlay CC, Hammer MD, Tchoungui H (2020) COV-OBS.x2: 180 years of geomagnetic field evolution from ground-based and satellite observations. Earth, Planets and Space 72(1):160

Hulot G, Le Mouël J (1994) A statistical approach to the earth's main magnetic field. Phys Earth Planet Int 82(3–4):167–183

Jault D (2015) Illuminating the electrical conductivity of the lowermost mantle from below. Geophys J Int 202(1):482–496

Jiang Y, Finlay CC, Olsen N, Tøffner-Clausen L, Yan Q, Zhang K (2024) Macau scientific satellite-1 initial magnetic field model. Geophys Res Lett 51(22):e2024GL112,305

Kloss C, Finlay CC, Laundal KM, Olsen N (2023) Polar ionospheric currents and high temporal resolution geomagnetic field models. Geophys J Int 235(2):1736–1760

Lesur V, Wardinski I, Baerenzung J, Holschneider M (2018) On the frequency spectra of the core magnetic field gauss coefficients. Phys Earth Planet Int 276:145–158

Lesur V, Gillet N, Hammer M, Mandea M (2022) Rapid variations of Earth's core magnetic field. Surv Geophys 43:1–29

Lhuillier F, Fournier A, Hulot G, Aubert J (2011) The geomagnetic secular-variation timescale in observations and numerical dynamo models. Geophys Res Lett 38(9):1

Lühr H, Maus S (2010) Solar cycle dependence of quiet-time magnetospheric currents and a model of their near-earth magnetic fields. Earth, Planets and Space 62:843–848

Macmillan S, Olsen N (2013) Observatory data and the swarm mission. Earth, Planets and Space 65:1355–1362

Morzfeld M, Buffett BA (2019) A comprehensive model for the kyr and myr timescales of earth's axial magnetic dipole field. Nonlinear Process Geophys 26(3):123–142

Olsen N (2002) A model of the geomagnetic field and its secular variation for epoch 2000 estimated from orsted data. Geophys J Int 149(2):454–462

Olsen N (2021) Magnetometer data from the grace satellite duo. Earth, Planets and Space 73:1–20

Olsen N, Sabaka TJ, Lowes F (2005) New parameterization of external and induced fields in geomagnetic field modeling, and a candidate model for igrf 2005. Earth, Planets and Space 57(12):1141–1149

Olsen N, Lühr H, Finlay CC, Sabaka TJ, Michaelis I, Rauberg J, Tøffner-Clausen L (2014) The CHAOS-4 geomagnetic field model. Geophys J Int 197(2):815–827

Olsen N, Ravat D, Finlay CC, Kother LK (2017) LCS-1: a high-resolution global model of the lithospheric magnetic field derived from CHAMP and Swarm satellite observations. Geophys J Int 211(3):1461–1477

Olsen N, Albini G, Bouffard J, Parrinello T, Tøffner-Clausen L (2020) Magnetic observations from cryosat-2: calibration and processing of satellite platform magnetometer data. Earth, Planets and Space 72:1–18

Ou J, Jiang Y, Yan Q (2024) Study on the local time asymmetry of the earth's ring currents based on multiple low-earth-orbit satellite magnetic observations. Sci China Technol Sci 67(11):3354–3364

Pais M, Jault D (2008) Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field. Geophys J Int 173(2):421–443

Rasmussen CE, Williams C (2006) Gaussian processes for machine learning the mit press. Cambridge vol. 32. pp 68

Ropp G, Lesur V, Baerenzung J, Holschneider M (2020) Sequential modelling of the earth's core magnetic field. Earth, Planets and Space 72(1):153

Sabaka TJ, Olsen N, Purucker ME (2004) Extending comprehensive models of the earth's magnetic field with ørsted and champ data. Geophys J Int 159(2):521–547

Sabaka TJ, Tøffner-Clausen L, Olsen N, Finlay CC (2018) A comprehensive model of earth's magnetic field determined from 4 years of swarm satellite observations. Earth, Planets and Space 70:1–26

Sadhasivan M, Constable C (2022) A new power spectrum and stochastic representation for the geomagnetic axial dipole. Geophys J Int 231(1):15–26

Sanchez S, Wicht J, Bärenzung J (2020) Predictions of the geomagnetic secular variation based on the ensemble sequential assimilation of geomagnetic field models by dynamo simulations. Earth, Planets and Space 72:1–20

Schaeffer N, Jault D (2016) Electrical conductivity of the lowermost mantle explains absorption of core torsional waves at the equator. Geophys Res Lett 43(10):4922–4928

Suttie N, Nilsson A, Gillet N, Dumberry M (2025) Large-scale palaeoflow at the top of earth's core. Earth Planet Sci Lett 652(119):185

Vugrin KW, Swiler LP, Roberts RM, Stucky-Mack NJ, Sullivan SP (2007) Confidence region estimation techniques for nonlinear regression in groundwater flow: three case studies. Water Resour Res 43(3):1

Yaglom AM (2004) An introduction to the theory of stationary random functions. Courier Corporation

Zhou Y, Dong ZY, Hsieh WP, Goncharov AF, Chen XJ (2022) Thermal conductivity of materials under pressure. Nature Rev Phys 4(5):319–335

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.