LABEL:firstpage–LABEL:lastpage

# Co-estimation of core and lithospheric magnetic fields
by a maximum entropy method

Mikkel Otzen  Christopher C. Finlay and Clemens Kloss
Division of Geomagnetism and Geospace  Department of Space Research and
Technology
Technical University of Denmark  Kgs. Lyngby  Denmark
E-mail<EMAIL_ADDRESS>

(Received XXX; in original form XXX; xxxx)

###### keywords:

Satellite magnetics, Rapid time variations, Magnetic anomalies: modelling and
interpretation, Inverse theory.

††volume: XXX

Satellite observations of the geomagnetic field contain signals generated in
Earth’s interior by electrical currents in the core and by magnetized rocks in
the lithosphere. At short wavelengths the lithospheric signal dominates,
obscuring the signal from the core. Here we present details of a method to co-
estimate separate models for the core and lithospheric fields, which are
allowed to overlap in spherical harmonic degree, that makes use of prior
information to aid the separation. Using a maximum entropy method we estimate
probabilistic models for the time-dependent core field and the static
lithospheric field that satisfy constraints provided by satellite observations
while being consistent with prior knowledge of the spatial covariance and
expected magnitude of each field at its source surface.

For the core field, we find that between spherical harmonic degree 13 and 22
power adds coherently to the established structures, and present a synthetic
test that illustrates the aspects of the small scale core field that can
reliably be retrieved. For the large scale lithospheric field we also find
encouraging results, with the strongest signatures below spherical harmonic
degree 13 occurring at locations of known prominent lithospheric field
anomalies in north-Eastern Europe, Australia and eastern North America.
Although the amplitudes of the small scale core field and large scale
lithospheric field are underestimated we find no evidence that obvious
artefacts are introduced. Compared with conventional maps of the core-mantle
boundary field our results suggest more localized normal flux concentrations
close to the tangent cylinder, and that low latitude flux concentrations occur
in pairs of opposite polarity. Future improvements in the recovery of the
small scale core field and large scale lithospheric field will depend on
whether more detailed prior information can be reliably extracted from core
dynamo and lithospheric magnetisation simulations.

##  1 Introduction

Earth’s magnetic field is a result of sources located both within the Earth
and above its surface in the upper atmosphere (Chapman & Bartels, 1940; Langel
& Hinze, 1998; Olsen & Stolle, 2012). Spherical harmonic analysis indicates
that internal sources are responsible for the majority of the field (Gauss,
1839; Hulot et al., 2015). The spatial power spectrum of the internal field at
Earth’s surface is steep at low degree (up to approximately spherical harmonic
degree 13) and essentially flat at higher degree (Lowes, 1974), indicating two
sources: one deep within the planet and one located near to the surface
(Voorhies et al., 2002; Voorhies, 2004). These deep and shallow sources are
thought to correspond to the core dynamo and lithospheric magnetisation.

Measurements of Earth’s magnetic field from space have provided an
increasingly detailed picture of the magnetic field due to internal sources.
The MAGSAT mission (Langel et al., 1982) delivered the first set of vector
measurements with global coverage allowing the change in the slope of the
spatial power spectra, between the wavelengths where core and respectively
lithospheric sources dominate, to be definitively observed (Langel & Estes,
1982). With more recent satellite missions, in particular the CHAMP (Reigber
et al., 2002) and Swarm (Friis-Christensen et al., 2006; Olsen & Floberghagen,
2018) missions, it is possible to determine the internal field spectrum out to
beyond degree 130 (Maus, 2010; Olsen et al., 2017). Knowledge of the small
scale core field (which we define here as above spherical harmonic degree 13)
has on the other hand advanced rather little since the time of MAGSAT due to
it being obscured by the lithospheric field. The conventional approach is to
estimate a single internal spherical harmonic field model and then to truncate
at degree 13 to study the core field, for example when plotting maps of the
radial field at the core-mantle boundary (CMB) (e.g. Cain et al., 1989; Olsen
et al., 2014; Sabaka et al., 2020).

Truncation at a fixed spherical harmonic degree has however limitations when
seeking to isolate the core field. Abrupt truncation in spectral space may
cause ringing in physical space (Whaler & Gubbins, 1981; Gubbins, 2007).
Furthermore, the lithospheric field does, of course, not stop at degree 14,
and there will be some contribution to the internal field below 13 from
lithospheric sources. Most seriously all information on the small scale core
field above degree 13 is lost.

In the 1980s it was suggested that a better way to estimate the core field
would be to minimize suitable norms of the field complexity at the CMB (Shure
et al., 1985; Gubbins & Bloxham, 1985). This approach, known as spatial
regularization of the field, has been widely adopted for studying the core
field over historical (Bloxham et al., 1989; Jackson et al., 2000) and
paleomagnetic timescales (Korte et al., 2011; Panovska et al., 2018) when data
coverage is sparse; it makes use of prior information from seismology on the
depth of the CMB, along with asking for a field that is simple in a specific
way (as measured by a chosen regularization norm) at the source radius. A
drawback is that traditional regularizations norms, such as the squared value
of the radial field or the horizontal gradient of the radial field integrated
over CMB, or Ohmic heating norms (Gubbins, 1976; Jackson et al., 2000),
strongly penalize small length scales and typically cause the spatial power
spectrum to decay in an unphysical fashion above degree 13 (Backus, 1988;
Buffett & Christensen, 2007).

Geodynamo simulations for which the magnetic Reynolds number is of order 1000,
as expected in Earth’s core (Christensen & Tilgner, 2004; Lhuillier et al.,
2011), involve localized, high amplitude, flux features and spatial spectra at
the CMB that are rather flat, decreasing only very slowly at spherical
harmonic degrees 10 to 30 (see, for example, Schaeffer et al., 2017; Aubert et
al., 2017; Sheyko et al., 2018). Jackson (2003) and Jackson et al. (2007)
showed that regularization norms based on the entropy of the radial field at
the CMB allowed the estimation of core fields with flatter spatial spectra and
localized, high amplitude flux features. The entropy regularization technique
was adapted to time-dependent spherical harmonic field models by Gillet et al.
(2007) and applied to satellite observations from the Ørsted, SAC-C and CHAMP
missions by Finlay et al. (2012). A drawback in these studies was the need to
abritrarily pick a value for the so-called default parameter (the magnitude of
the radial field expected in the absence of data constraints) that controlled
the width of the entropy distribution and hence the sharpness of the field
structures (Maisinger et al., 2004; Jackson et al., 2007). Jackson (2003) and
Jackson et al. (2007) focused on default values around 10 $\mu$T for the core
field, while Gillet et al. (2007) used 30 $\mu$T. Despite many desirable
features, entropy-based field reconstruction techniques have been little
exploited in subsequent years in part due to doubts as to how to pick the
troublesome default parameter.

An important conceptual step forward in co-estimating core and lithospheric
field sources was made by Holschneider et al. (2016). They suggested how
various field sources (including the core and lithospheric fields) could be
co-estimated within a Bayesian framework making use of prior information, for
example on the expected source depth and its correlation (or covariance)
structure. This approach has been used to develop temporal sequences of field
models using a Kalman filter algorithm, being applied to the modelling of
satellite and ground magnetic observations by Ropp et al. (2020), Ropp & Lesur
(2023) and Baerenzung et al. (2020, 2022). Using a simple correlation function
and treating the source depth as a free parameter Baerenzung et al. (2020)
were able to construct stable maps of the field at the CMB up to spherical
harmonic degree 20, although it was found to be difficult to reliably separate
the large scale lithospheric field. More detailed prior information on the
covariance between spherical harmonic coefficients in dynamo simulations has
also been used in combination with observation-based internal field models up
to degree 13 to infer the core field up to degree 30 (Aubert, 2015, 2020). On
the theoretical side Baratchart & Gerhards (2017) have shown that core and
lithospheric fields can be formally separated if the lithospheric field
sources are localized to a sub-region of the spherical surface. Non-Gaussian
field distributions thus seem to aid the separation of fields from different
sources, as is well known in other contexts such as independent component
analysis (e.g. Hyvärinen & Oja, 2000).

Here we build on the above studies and seek to estimate separate models for
the core and lithospheric fields within a Bayesian framework using a maximum
entropy method that accounts for spatial covariances found in first principles
simulations of the core dynamo and the lithospheric magnetisation. Similar
maximum entropy based techniques have previously been applied to signal
separation problems in cosmology (Hobson et al., 2010). Section 2 sets out
details of our Bayesian model estimation scheme and specifies the prior
information used. Section 3 describes the satellite and ground magnetic
observations employed. Section 4 presents our results, with Appendix A
collecting findings from a synthetic test based on a similar data and
modelling setup. We conclude in Section 5 with a discussion of what has been
achieved and suggestions for future improvements of the method.

##  2 Methodology

###  2.1 Geomagnetic field model

We model Earth’s magnetic field $\mathbf{B}$ as a potential field,
representing it by the gradient of a scalar potential $V$ such that

$\mathbf{B}=-\nabla V\qquad\mbox{where}\quad V=V^{int}+V^{ext},$ (1)

with $V^{int}$ the potential due to internal sources and $V^{ext}$ that due to
external sources. Both core and lithospheric sources contribute to $V^{int}$,
we represent each by a separate spherical harmonic expansion

$\displaystyle V^{int}(r,\theta,\phi)$ $\displaystyle=$ $\displaystyle
a\sum^{N_{C}}_{n=1}\bigg{(}\frac{a}{r}\bigg{)}^{n+1}\sum^{n}_{m=0}\Big{[}g^{C}_{n,m}(t)\cos
m\phi+h^{C}_{n,m}(t)\sin m\phi\Big{]}P^{m}_{n}(\cos\theta)$ (2)
$\displaystyle+$ $\displaystyle
a\sum^{N_{L}}_{n=1}\bigg{(}\frac{a}{r}\bigg{)}^{n+1}\sum^{n}_{m=0}\Big{[}g^{L}_{n,m}\cos
m\phi+h^{L}_{n,m}\sin m\phi\Big{]}P^{m}_{n}(\cos\theta)$ (3)

where $(r,\theta,\phi)$ are geocentric spherical polar coordinates, $a$ is the
Earth’s mean spherical reference radius, $n$ is the degree of the spherical
harmonic, $m$ the order of the spherical harmonic and $P^{m}_{n}(\cos\theta)$
are associated Legendre functions. $g^{L}_{n,m}$ and $h^{L}_{n,m}$ are
spherical harmonic coefficients describing the lithospheric field, assumed
here to be static, and here considered up to a maximum degree of $N_{L}=120$.
$g^{C}_{n,m}(t)$ and $h^{C}_{n,m}(t)$ are time-dependent spherical harmonic
coefficients for the core field, considered up to a maximum degree $N_{C}=30$.
These are expanded in time using a B-spline basis, of order 6 and with a 0.5
year knot spacing,

$g^{C}_{n,m}(t)=\sum_{k}g^{C}_{n,m,k}\mathcal{B}_{k}(t).$ (4)

where $\mathcal{B}_{k}$ is the kth basis function of the order 6 B-splines. We
collect the coefficients describing the core field in a vector
$\mathbf{m}^{C}=\left\\{g^{C}_{n,m,k},h^{C}_{n,m,k}\right\\}$ and the
coefficients describing the lithospheric field in a vector
$\mathbf{m}^{L}=\left\\{g^{L}_{n,m},h^{L}_{n,m}\right\\}$. Note that the core
and lithospheric field representations overlap between spherical harmonic
degrees 1 and 30, additional prior information is therefore needed in order to
separate them.

As in the CHAOS-7 model (Finlay et al., 2020) these internal field
coefficients are supplemented by model coefficients $\mathbf{m}^{ext}$
describing the external field, and coefficients $\mathbf{m}^{q}$ describing
the in-flight alignment of the vector magnetometers on each satellite, to give
the full model vector
$\mathbf{m}=\left[\mathbf{m}^{C},\mathbf{m}^{L},\mathbf{m}^{ext},\mathbf{m}^{q}\right]^{T}$.

###  2.2 Bayesian model estimation

####  2.2.1 Entropic priors for the core and lithospheric fields

We make use of prior information regarding the radial component of the core
and lithospheric fields at their respective source surfaces, the CMB and
Earth’s surface. The radial field is evaluated at each source surface on an
approximately equal area grid, and values are collected into vectors
$\mathbf{b}^{C}$ and $\mathbf{b}^{L}$ for the core and lithospheric fields
respectively. Such knowledge of the radial field at the source surface
completely defines the potential due to an internal source. These are related
to the spherical harmonic model coefficients discussed in the previous section
by

$\mathbf{b}^{C}(t_{p})=\mathbf{G}^{C,t_{p}}\mathbf{m}^{C}\quad\mbox{and}\quad\mathbf{b}^{L}=\mathbf{G}^{L}\mathbf{m}^{L}$
(5)

with $\mathbf{G}^{C,t_{p}}$ and $\mathbf{G}^{L}$ being matrices that
synthesize the radial field from the relevant spherical harmonic model
coefficients, for the core field at some epoch $t_{p}$.

Knowledge regarding the spatial covariance of each field at its source surface
is provided in the form of a-priori model covariance matrices $\mathbf{C}_{C}$
and $\mathbf{C}_{L}$, with lower triangular Cholesky factors $\mathbf{L}_{C}$
and $\mathbf{L}_{L}$ that can be used to transform $\mathbf{b}^{C}$ and
$\mathbf{b}^{L}$ to latent variables $\mathbf{x}^{C}$ and $\mathbf{x}^{L}$
such that

$\mathbf{b}^{C}(t_{p})=\mathbf{L}_{C}\mathbf{x}^{C}(t_{p})\quad\mbox{and}\quad\mathbf{b}^{L}=\mathbf{L}_{L}\mathbf{x}^{L}.$
(6)

The latent variables $\mathbf{x}^{C}$ and $\mathbf{x}^{L}$ therefore describe
the core and lithospheric radial fields at their source surfaces in a space
where their elements are normalized and decorrelated, as is appropriate for
the application of maximum entropy methods (Maisinger et al., 2004).

It is assumed that $\mathbf{x}^{C}$ and $\mathbf{x}^{L}$ are each described by
an entropic probability density function

$P(\mathbf{x}^{C})\propto\exp\left[\lambda_{C}S(\mathbf{x}^{C})\right]\quad\mbox{and}\quad
P(\mathbf{x}^{L})\propto\exp\left[\lambda_{L}S(\mathbf{x}^{L})\right],$ (7)

where $S$ is the information entropy for variables that can take both positive
and negative values (Gull & Skilling, 1990; Hobson & Lasenby, 1998)

$S[\mathbf{x},\omega]=\sum\limits_{i=1}^{M}\left[\psi_{i}-2\omega-
x_{i}\log\left(\frac{\psi_{i}+x_{i}}{2\omega}\right)\right]$ (8)

where $M$ is in our case the number of grid points on the spherical surface,
$\psi_{i}=\sqrt{x_{i}^{2}+4\omega^{2}}$, and $\omega$ is a so-called ’default’
parameter that defines the width of the entropy function. The information
entropy function $S$ is a measure of the amount of uncertainty inherent in the
distribution of values $\mathbf{x}$ (Shannon, 1948; Jaynes, 2003), the form we
use follows from requirements of subset independence, coordinate invariance,
system independence and scaling (Skilling, 1988). In the geomagnetic context
it can be thought of as measuring the number of ways a given distribution of
radial field on the source surface can be arranged from elementary flux
bundles (Jackson, 2003; Jackson et al., 2007); fields with larger entropy are
simpler in the sense that they can be arranged in more ways.

Assignment of an entropic prior is argued to be an appropriate choice in the
absence of precise information as to the form of a prior pdf (Skilling, 1989;
Hobson et al., 1998), and it is compatible with possibly non-Gaussian
distributions of $\mathbf{x}$. Maximizing the entropy essentially broadens the
distribution of $\mathbf{x}$ as much as possible without violating the
available observational constraints. The resulting distribution therefore
agrees with what is known but expresses maximum uncertainty with respect to
all other matters (Jaynes, 1968). Maximizing the entropy does not introduce
additional correlations amongst the latent variables (Gull & Skilling, 1984).

The factors $\lambda_{C}$ and $\lambda_{L}$ appearing in the entropic pdfs are
scaling factors. We are able to set these equal to 1 because $\mathbf{x}^{C}$
and $\mathbf{x}^{L}$ have already been normalized via the transform of
$\mathbf{b}^{C}$ and $\mathbf{b}^{L}$ to the latent space (Hobson et al.,
2010). The transform to latent space using the a-priori covariance functions
also ensures the entropy is computed from uncorrelated variables, an important
condition for the maximum entropy method.

To put into practice the above scheme we require prior information concerning
the covariance structure of the radial fields on the source surfaces and the
expected widths of the distributions of $\mathbf{x}$ for each source. We
obtain these from first principles simulations of the core dynamo and the
lithospheric magnetisation. Full details are provided in Otzen (2022) only a
short summary is given here.

For the core field prior, we use an ensemble of realizations of the core field
produced by versions of the coupled-Earth dynamo of Aubert et al. (2013). This
numerical dynamo is known for producing field structures and patterns of
secular variation similar to those observed over the past centuries. To start
with we used a collection of radial fields realizations, well separated in
time, generated by the original version of the coupled-Earth dynamo (Aubert et
al., 2013) that has been used in previous field modelling and data
assimilation studies (Barrois et al., 2017; Ropp et al., 2020). To this we
added radial field realizations from a long run of an updated version (71% of
path) of the coupled-Earth dynamo (Aubert & Gillet, 2021). Although these two
cases involve different control parameters they lie on a path through control
parameter space along which the field morphology is essentially invariant
(Aubert et al., 2017). We finally augmented our set of realizations by
carrying out rotations of the simulated core fields by an arbitrary amount in
longitude, this was possible because the covariance functions we use do not
depend on longitude and this enabled us to work with a larger ensemble. In all
this resulted in an ensemble of 5688 core field realizations up to spherical
harmonic degree 30. To be more consistent with the observed field we also
adjusted the dipole fields from the dynamo simulations replacing the $n=1$
coefficients with random samples from normal distributions with mean values of
$\bar{g}^{C}_{1,0}=-29000$ nT, $\bar{g}^{C}_{1,1}=0$ nT, and
$\bar{h}^{C}_{1,1}=0$ nT and standard deviations of $5000$nT, $3000$nT, and
$3000$nT respectively, the latter being similar to those seen in the dynamo
realizations. The spread of the power spectra from the resulting ensemble
encompasses the observed internal field (e.g. Finlay et al., 2020) up to
spherical harmonic degree 13 (Otzen, 2022).

(a) (b)
(c) (d)
(e) (f)

Figure 1: Statistical properties of the prior ensemble of core and
lithospheric fields. Top row: Empirical covariance functions (a) for the CMB
radial field, (b) for the lithospheric radial field at Earth’s surface showing
ensemble members (grey), and estimated covariance functions (black). For
computational reasons empirical covariances are shown only out to 18 degrees
for the lithospheric field. Middle row: power spectra for 5000 realizations
generated from the estimated covariance functions (c) at the CMB, and (d) at
Earth’s surface with the CHAOS-7 model (up to degree 13) and from the LCS-1
model (above degree 16) for reference. Bottom row: empirical probability
density functions for radial magnetic field values on approximately equal area
grids at (e) the CMB and (f) Earth’s surface, derived from the prior ensembles
of core and lithospheric magnetic fields respectively, after transformation to
an uncorrelated and normalized latent space.

Our prior for lithospheric field comes from simulations of the lithospheric
magnetization based on the forward modelling scheme developed by Hemant & Maus
(2005), with revised oceanic magnetisation according to Masterton et al.
(2012) and subduction zone magnetisations following Williams & Gubbins (2019).
We produced an ensemble of $6000$ realizations of the lithospheric field by
(i) varying the crustal thickness within a range given by published crustal
thickness models (Nataf & Ricard, 1996; Reguzzoni & Sampietro, 2015), (ii)
varying parameters of the remanent vertically integrated magnetisation model
(Masterton et al., 2012; Williams & Gubbins, 2019), and (iii) using stochastic
perturbations generated using a Gaussian random field approach. The observed
power spectra for the lithospheric field, for example from the LCS-1 model
(Olsen et al., 2017), lies within the spread of the power spectra of this
ensemble (Otzen, 2022). The simulated fields were generated up to spherical
harmonic degree $255$, but here we only used them up to degree $120$.

From each ensemble member, we evaluated the radial field on an approximately
equal area grid on the source surface. We used HEALPix (Górski et al., 2005)
grids with $3072$ points for the core field, and $49152$ points for the
lithospheric field, which are suitable for representing fields up to the
spherical harmonic truncation level chosen for each source. Based on these
gridded values we computed empirical semi-variograms as a function of angular
distance on the spherical surface and fit covariance functions to these. We
used a multi-quadratic covariance function for the lithospheric field
(Gneiting, 2013) and a combination of a multi-quadratic covariance function
and a spline function for large distances for the core field. The resulting
covariance functions, along with the empirical covariances of the ensemble
members, are shown in the top row of Fig. 1 along with their corresponding
power spectra at the source surfaces generated using these covariance
functions are shown in the middle row. These covariance functions provide us
with the a-priori expected covariance structure for our core and lithospheric
fields. Fig. 2 shows example realizations of the core and lithospheric field
generated using these covariance functions.

We also make use of the distribution of the radial field at the sources
surfaces from our core and lithospheric field ensembles, after transformation
to the latent space (see equation (6) and the related discussion). These
distributions and relevant statistics are presented in the bottom row of Fig.
1. In particular we use expected absolute values, $<|x|>=\int|x|\,p(x)\,dx$,
calculated using the mean empirical pdfs shown in Fig. 1, to define the latent
space default parameters i.e. we set $\omega=<|x|>$, separately for the core
and lithospheric fields.

(a) Core field prior realizations

(b) Lithospheric field prior realizations

Figure 2: Example realizations of radial magnetic fields generated using the
estimated covariance functions for the core and lithospheric fields. The core
fields are shown at the core surface (top row) and the lithospheric fields are
shown at Earth’s surface (bottom row). These illustrate the a-prior
correlation structures assumed for the core and lithospheric fields. Note the
covariance models used to generate these prior fields are rotationally
invariant.

The time-dependence of the core field is represented by a 6th order B-spline
representation smoothed by third time derivative regularization, a standard
choice in time-dependent geomagnetic field models when one wishes to study
field accelerations (see e.g. Olsen et al., 2014). This is formally equivalent
within the Bayesian framework to assuming a-priori that the spherical harmonic
coefficients are realizations of a continuous process of the form (Wahba,
1990)

$g^{C}_{n,m}(t)=K^{C,0}_{n,m}+K^{C,1}_{n,m}\,t+K^{C,2}_{n,m}\,t^{2}+f_{n,m}(t),$
(9)

where $f_{n,m}(t)$ is a zero mean Gaussian process (see e.g. Rasmussen &
Williams, 2006) specified by the covariance function

$C_{n}(t_{a},t_{b})=\sigma_{n}^{2}\int{(t_{a}-u)_{+}^{2}\over
2}{(t_{b}-u)_{+}^{2}\over 2}\,du,$ (10)

where $t_{a}$ and $t_{b}$ are arbitrary times, $u$ is a dummy integration
variable, and we have used the notation $(z)_{+}=z$ for $z\geq 0$ and
$(z)_{+}=0$ otherwise. $K^{C,0}_{n,m}$ , $K^{C,1}_{n,m}$ and $K^{C,2}_{n,m}$
are constants associated with constant, linear and quadratic time-dependences
that can be different for each spherical harmonic coefficient.
$\sigma^{2}_{n}$ is the variance of the process that depends on the spherical
harmonic degree $n$ and is related to the choice of regularization parameter
in the standard field modelling framework.

####  2.2.2 Likelihood of geomagnetic observations

Turning to the observations, we assume the geomagnetic measurements are
contaminated by unmodelled signals that follow a long-tailed Huber error
distribution. The appropriate likelihood function is then

$P(\mathbf{d}\,|\,\mathbf{m})\propto\exp\left[-\frac{1}{2}\chi^{2}(\mathbf{m})\right]$
(11)

where $\chi^{2}(\mathbf{m})=\mathbf{e}^{T}\,\mathbf{W}\,\mathbf{e}$ is a
robust (Huber weighted) misfit norm.
$\mathbf{e}=\mathbf{d}-\mathbf{g}(\mathbf{m)}$ are the residuals between the
ground and satellite magnetic observations $\mathbf{d}$ and the associated
model predictions $\mathbf{g}(\mathbf{m)}$.
$\mathbf{W}=\mathbf{C_{e}}^{-1/2}\mathbf{W_{h}}\mathbf{C_{e}}^{-1/2}$ where
$\mathbf{C_{e}}$ is the a-priori data error covariance matrix and
$\mathbf{W_{h}}$ a diagonal matrix that implements robust (Huber) weights and
has elements
$W^{i,i}_{h}=\mbox{min}[1,\left(c\,\sigma^{d}_{i}/|e_{i}|\right)]$ where
$\sigma^{d}_{i}$ is the a-priori expected error for the $i$th datum and
$c=1.5$ is a constant (Constable, 1988; Olsen, 2002; Sabaka et al., 2004).

####  2.2.3 Estimation of the model posterior probability density function

Applying Bayes theorem the posterior probability function is

$P(\mathbf{m}\,|\,\mathbf{d})\propto\exp\left[-\frac{1}{2}\chi^{2}(\mathbf{m})+S_{tav}(\mathbf{x}^{C})+S(\mathbf{x}^{L})-\frac{1}{2}\mathbf{m}^{T}\mathbf{C}_{T}^{-1}\mathbf{m}\right]$
(12)

This can be maximized by minimizing the loss function

$\Phi(\mathbf{m})=\frac{1}{2}\chi^{2}(\mathbf{m})-S_{tav}(\mathbf{x}^{C})-S(\mathbf{x}^{L})+\frac{1}{2}\mathbf{m}^{T}\mathbf{C}_{T}^{-1}\mathbf{m}.$
(13)

Here $S_{tav}(\mathbf{x}^{C})={1\over
N_{P}}\sum\limits_{p=1}^{N_{P}}S[\mathbf{x}^{C}(t_{p})]$ approximates the
information entropy of the (spatially decorrelated) CMB radial field averaged
over time (Gillet et al., 2007) by summing the entropy at each discrete epoch
$t_{p}$ over $N_{P}=1000$ epochs.
$\mathbf{x}^{C}(t_{p})=\mathbf{L}_{C}^{-1}\mathbf{G}^{C,t_{p}}\mathbf{m}^{C}$
changes with $t_{p}$,
$\mathbf{x}^{L}=\mathbf{L}_{L}^{-1}\mathbf{G}^{L}\mathbf{m}^{L}$ is assumed to
be static. $\mathbf{C}_{T}$ is the a-prior model (temporal) covariance matrix,
which includes the temporal prior information from Eqn. 10 for the core part
of the model.

The minimization is achieved iteratively using a Newton-type descent
algorithm. The $(k+1)$th estimate of the posterior mean model is obtained
based on the model at the previous $k$th iteration by

$\mathbf{m}_{k+1}=\mathbf{m}_{k}+\left[\mathbf{A}_{k}^{T}\mathbf{W}_{k}\mathbf{A}_{k}+\bm{\alpha}^{C}_{k}+\bm{\alpha}^{L}_{k}+\mathbf{C}_{T}^{-1}\right]^{-1}\cdot\left[\mathbf{A}_{k}^{T}\mathbf{W}_{k}\mathbf{e}_{k}-\mathbf{C}_{T}^{-1}\mathbf{m}_{k}-\bm{\beta}^{C}_{k}-\bm{\beta}^{L}_{k}\right]$
(14)

where $\mathbf{A}_{k}=(A_{ij})_{k}={\partial
g_{i}(\mathbf{m}_{k})\over\partial m_{j}}$ is the Jacobian matrix of partial
derivatives of the forward model for each datum with respect to the model
parameters, evaluated using the model parameters $\mathbf{m}_{k}$. We have
used a notation similar to that of Stockmann et al. (2009) to define matrices,
related to the Hessian matrix of the second-order partial derivatives of the
entropy function, of the form

$\bm{\alpha}=\left(\mathbf{L}^{-1}\mathbf{G}\right)^{T}\mathbf{S}_{H}\left(\mathbf{L}^{-1}\mathbf{G}\right)\quad\mbox{where}\quad\mathbf{S}_{H}=\,\mbox{diag}\left[{4\omega\over\psi_{1}},{4\omega\over\psi_{2}},...,{4\omega\over\psi_{M}}\right]^{T},$
(15)

and vectors, related to the first order partial derivatives of the entropy
function, of the form

$\bm{\beta}=\left(\mathbf{L}^{-1}\mathbf{G}\right)^{T}\mathbf{S}_{G}\quad\mbox{where}\quad\mathbf{S}_{G}=4\omega\left[\log\left({\psi_{1}+x_{1}\over
2\omega}\right),\log\left({\psi_{2}+x_{2}\over
2\omega}\right),...,\log\left({\psi_{M}+x_{M}\over
2\omega}\right)\right]^{T}.$ (16)

In $\bm{\alpha}^{C}_{k}$, $\bm{\beta}^{C}_{k}$ and $\bm{\alpha}^{L}_{k}$ ,
$\bm{\beta}^{L}_{k}$ in Eq. (14) the superscripts $C$ and $L$ refer to the
core and lithospheric fields respectively, while subscript $k$ denotes that
the computation is carried out using model parameters from the previous $k$ th
iteration. In $\bm{\alpha}^{C}$ and $\bm{\beta}^{C}$ the expressions for
$\bm{\alpha}$ and $\bm{\beta}$ at each epoch $t_{p}$ must be averaged over
time in the same way as $S_{tav}$ is defined above.

Minimizing measures of the data misfit and the temporal complexity while
maximizing the entropy of the (decorrelated) core and lithospheric radial
fields at their source surfaces results in internal fields that fit the
observations and are compatible with the temporal prior but allow high dynamic
ranges of $\mathbf{x}$ at the source surfaces while satisfying the spatial
covariance properties of the core and lithospheric priors.

After convergence of the above scheme we describe the dispersion of the
posterior distribution using an approximate model covariance matrix defined
about the maximum of the posterior pdf, computed from the Hessian of the loss
function $\Phi(\mathbf{m})$ by (e.g. Tarantola, 2005; Hobson et al., 1998)

$\mathbf{C_{m}}\approx(\mathbf{A}^{T}\mathbf{W}\mathbf{A}+\mathbf{\bm{\alpha}^{C}}+\mathbf{\bm{\alpha}^{L}}+\mathbf{C}_{T}^{-1})^{-1}$
(17)

where $\mathbf{A}$, $\mathbf{W}$, $\mathbf{\bm{\alpha}^{C}}$ and
$\mathbf{\bm{\alpha}^{L}}$ are the values of $\mathbf{A}_{k}$
$\mathbf{W}_{k}$, $\mathbf{\bm{\alpha}^{C}_{k}}$ and
$\mathbf{\bm{\alpha}^{L}_{k}}$ from the final iteration, when the scheme is
considered to have converged.

##  3 Observations

The models reported here are built from a dataset of geomagnetic observations
similar to that used to construct the CHAOS-7 geomagnetic field model (Finlay
et al., 2020), but restricted to the period between 2005 and 2020.

Observations from the CHAMP, Cryosat-2 and Swarm A and B satellites are used,
three-component vector measurement for quasi-dipole latitudes up to 55 degrees
and scalar intensity data at higher latitudes. Level 3 CHAMP magnetic field
data, Cryosat-2 L1b magnetic field data (FGM 1, the dataset used in the
CHAOS-7 model, here pre-calibrated using CHAOS-7) and Swarm L1b magnetic field
data (version 0601) are used with 1 minute sampling for CHAMP and Cryosat-2
and 2 minute sampling for each Swarm satellite. We also used along-track
gradients from CHAMP, Swarm A and Swarm B; gradients are particularly useful
for constraining the high degree lithospheric field. Geomagnetic quiet-time
selection criteria were employed such that Kp $\leq 2^{0}$, d$|$RC$|$/dt
$\leq$ 2 nT/hr (Olsen et al., 2014), averaging over the previous 2 hours the
merging electric field at the magnetopause $E_{m}\leq 0.8$ mV/m, the
interplanetary magnetic field (IMF) $B_{Z}>0$ , and IMF $B_{Y}$ is less than 3
nT in the northern hemisphere and greater than -3 nT in the southern
hemisphere. Only data from dark conditions (sun at least 10 degrees below the
horizon) were used. A more detailed description of these data selection
criteria is found in Finlay et al. (2020). In addition to satellite
observations, as in CHAOS-7 we used annual differences of revised monthly
means from ground observatories are used, as in the CHAOS-7 model, based on
hourly mean values from the BGS AUX$\\_$OBS database, version 0129. A stacked
histogram of the number of vector field observations used versus time is
presented in Fig. 3.

Figure 3: Stacked histogram showing number of vector field observations
employed per month between 2005 and 2020. Colours indicate the data source.

##  4 Results

###  4.1 Implementation and model diagnostics

We now briefly document here some details regarding the practical
implementation before moving on to the results. In order to make our modelling
setup as close as possible to the CHAOS field modelling scheme, in
$\mathbf{C}_{T}^{-1}$ we consider sub-matrices associated with calculating
quadratic norms of (i) the 3rd time derivative of the internal radial field
integrated over the core surface and and throughout the model time span
(formally equivalent to the the prior defined in Eqn. 9) (ii) the acceleration
of the core field at the model endpoints and (iii) temporal first differences
of the estimated offsets of the external dipole in solar-magnetic coordinates
(related to imperfections in the RC index). The related hyperparameters
denoted by $\lambda_{i3}$, $\lambda_{i2e}$, and $\lambda_{sm}$ are fixed
throughout this study and implicitly included within $\mathbf{C}_{T}^{-1}$.
They were chosen so as to produce time-variations similar to the CHAOS-7 field
model although for simplicity we did not use degree-dependent tapering or
treat zonal terms differently. The relevant hyperparameters related to the
entropy default parameters for the dimensionless latent variables,
$\omega^{C}$ and $\omega^{L}$, were set to values 0.412 and 0.422 based on the
expected absolute values of $\mathbf{x}^{C}$ and $\mathbf{x}^{L}$ from the
distributions of the dimensionless latent variables found in the prior
ensembles (see Fig. 1, bottom). The adopted hyperparameters are collected in
table 1.

Table 1: Modelling hyperparameters Parameter |  Value
---|---
$\omega^{C}$ |  $0.412$
$\omega^{L}$ |  $0.422$
$\lambda_{i3}$ |  $1$ $(\text{nT yr}^{-3})^{-2}$
$\lambda_{i2e}$ |  $100$ $(\text{nT yr}^{-2})^{-2}$
$\lambda_{sm}$ |  $1200$ $(\text{nT yr}^{-1})^{-2}$

The full estimated model, including the time-dependent core field to degree
30, static lithospheric field to degree 120, external field parameters and
alignment parameters for each satellite, consists of $49495$ parameters in
all. We started the iterative model estimation scheme with model parameters
for the core taken from CHAOS-6.9 up to degree $n=12$, and for the lithosphere
from the LCS-1 model at degree $n=16$ and above. The small-scale core field
and large-scale lithospheric field were otherwise initialized with zeros.
After $24$ iterations the largest change in a model parameter relative to its
amplitude was $0.0037\%$ and no further change was seen in CMB maps of the
posterior mean core field. Below we refer to the resulting model, including
co-estimated core and lithospheric parts, as model CL.

For comparison, we also built a more traditional CHAOS-type field model, with
a single time-dependent internal field up to spherical harmonic degree 20 and
a static internal field for degrees 21 to 120, using the same external field
parameterization and hyperparameters, covering the same period, and based on
the same dataset as used to build the model CL. For this model we considered
the estimated core field to be the time-dependent internal field up to degree
13, as has been the standard practise when interpreting the CHAOS model (e.g.
Olsen et al., 2014).

Table 2 collects Huber-weighted means and RMS residuals between the vector
field data and the model predictions (in nT), comparing model CL and our
CHAOS-type reference model. Model CL fits the satellite and ground data
overall to a similar level as the CHAOS-type model, while simultaneously
minimizing the information entropy of the (spatially-decorrelated) time-
dependent core and lithospheric fields. The two models are found to have very
similar temporal regularization norms, which is not surprising at they were
built using the same temporal hyperparameters. The non-dimensional information
entropy norms for the decorrelated core and lithospheric fields, $S^{C}_{tav}$
and $S^{L}$, after 24 iterations were respectively $0.94$ and $5.31$ for the
CL model.

Table 2: Misfit statistics for vector field data in the non-polar region and scalar data in the polar region. Note that Huber weights were employed when calculating the reported statistics. (Misfits for gradient data not shown). |  |  | CL | CHAOS-type
---|---|---|---|---
Source |  QD lat |  N |  Mean (nT) |  RMS (nT) |  Mean (nT) |  RMS (nT)
CHAMP |  non-polar |  $2,080,146$ |  $-0.003$ |  $2.397$ |  $-0.003$ |  $2.397$
|  polar |  $98,751$ |  $-0.001$ |  $3.498$ |  $-0.004$ |  $3.496$
CryoSat-2 |  non-polar |  $1,017,960$ |  $-0.005$ |  $6.160$ |  $-0.012$ |  $6.160$
|  polar |  $59,509$ |  $1.371$ |  $6.556$ |  $1.538$ |  $6.604$
Swarm-A |  non-polar |  $1,085,118$ |  $0.007$ |  $2.182$ |  $0.007$ |  $2.182$
|  polar |  $47,936$ |  $0.039$ |  $3.117$ |  $0.045$ |  $3.117$
Swarm-B |  non-polar |  $1,079,508$ |  $-0.016$ |  $2.176$ |  $-0.016$ |  $2.176$
|  polar |  $48,433$ |  $0.194$ |  $2.905$ |  $0.203$ |  $2.906$
|  |  |  (nT/yr) |  (nT/yr) |  (nT/yr) |  (nT/yr)
Ground observatory |  |  $67,632$ |  $-0.027$ |  $3.304$ |  $-0.026$ |  $3.306$

###  4.2 The core-mantle boundary field

The spatial power spectra (Lowes, 1966; Mauersberger, 1956) of model CL at the
CMB in 2020 is presented in Fig. 4. The obtained posterior mean model closely
follows the internal field from CHAOS-7 up to degree 11 but contains slightly
less power at degree 12 and 13. Here and below comparisons to CHAOS-7 used
version CHAOS-7.9. Internal field models (such as the CHAOS model) are usually
truncated at degree 13 when carrying out interpretations at the CMB since
above this degree their CMB spectra diverge as they also include signals from
the lithospheric field. The posterior realizations from the model CL core
field have power spectra that are approximately flat out to degree 30 and do
not diverge. The posterior mean model shows a gradual drop in power for
degrees 15 to 19 and a slight increase again for degrees 20 to 22. Above
degree 22 the power in the posterior mean model drops to much lower levels
indicating that field realizations essentially average to zero at higher
degrees where almost all the observed signal comes from the lithospheric
field.

The power in the lithospheric field realizations and mean models mapped down
to the CMB is also shown in Fig. 4, these do diverge. The core and
lithospheric field spectra cross between degree 14 and 16 for realizations of
model CL, and at degree 15 for the mean models. Note that the estimated
lithospheric field is presented only for degree 2 and above, at degree 1 it is
not well separated from the core field, which we believe is a consequence of
the entropy function of the decorrelated latent variables not being greatly
affected by changes in the dipole field.

Figure 4: Lowes-Mauersberger spectra showing power as a function of spherical
harmonic degree at the core-mantle boundary in 2020.0 for model CL. Solid
black lines with circles show the estimated posterior mean core and
lithospheric field models which cross at degree 15. Grey lines show
realizations based on the adopted prior covariance model, cyan lines show
posterior realizations of the core field and purple lines posterior
realizations of the lithospheric field. Degree 1 is not shown for the
lithospheric field as it is not well separated from the core field. CHAOS-7
(up to degree 13) and LCS-1 (at degree 16 and above) are shown for reference.

Maps of the radial magnetic field at the CMB in 2020.0 from model CL, for the
posterior mean core field model and four example posterior realizations all
truncated at degree 22, are presented in Fig. 5. A similar map from the CHAOS-
type reference model, truncated in the conventional fashion at degree 13, is
shown for reference. The posterior realizations contain more power at small
length scales, but all realizations agree on the larger-scale structure as
represented by the mean model.

Figure 5: Radial magnetic field maps at the core-mantle boundary in 2020 for
model CL. Estimated posterior mean model (top left), a CHAOS-type model
constructed from the same dataset (top right), and four example posterior
realizations. All models are truncated at degree 22 except for the CHAOS-type
model which is truncated at degree 13. The common features across all
posterior models, as captured in the mean model, provide information on the
core field to beyond degree 13.

Figure 6: Time sequence of maps of radial field at the core-mantle boundary in
2006, 2010, 2014 and 2019 from the model CL posterior mean truncated at degree
22. Changes are due both to the true evolution of the core field and changes
in observational constraints across the epochs.

More details of the CMB field structures are evident in the model CL posterior
mean compared to the traditional CHAOS-type model truncated at degree 13. Some
low latitude flux concentrations are split, see for example the two strong
positive radial field features near the equator between Africa and South
America which are usually interpreted as single feature (as in the CHAOS-type
model). A strong high latitude flux feature under Siberia, located under the
Taymyr peninsula in central northern Siberia in 2020 is found to be more
localized and stronger than in models truncated at degree 13. This feature has
moved north-westwards between 2005 and 2020, as seen in Fig. 6 which shows the
posterior mean map up to spherical harmonic degree 22 at a sequence of times.

Concerning reversed flux patches in the South Atlantic, we find evidence for
two reversed patches under South Africa adjacent to strong norm flux patches
under central Africa (see also their time evolution in Fig. 6). Regarding the
reverse flux region under the Southern Atlantic ocean, there are several
distinct reversed flux concentrations visible within this region, which are
observed to evolve separately.

It is also evident in Fig 6 that there is more power in the mean model at
degrees 16 to 22 for the first five years and last six years of the model,
when CHAMP and Swarm observations respectively were available. Fig. 7 further
illustrates the time-dependence of the core field, focusing on coherent east-
west motions of flux features in time-longitude plots of the CMB radial field
at the equator and at 55 degrees south. Intense equatorial features are
observed to drift coherently westwards under the mid-Atlantic over the 15
years studied. There is also evidence for eastward drift of a reversed flux
feature under the Southern mid-Atlantic ocean at 55 degrees south heading from
South America toward Africa. Such coherently drifting flux features involving
power above spherical harmonic degree 13 cannot be ascribed to lithospheric
sources.

(a) (b)

Figure 7: Time-longitude plots of the posterior mean from model CL truncated
at degree 22, showing the evolution of the radial magnetic field at core-
mantle boundary (a) on the equator and (b) at latitude 55 degrees south
(right).

###  4.3 The lithospheric field at Earth’s surface

Although not our main focus here, we present for reference details of our co-
estimated lithospheric field, which includes an estimate of the large scale of
the lithospheric field which is usually neglected.

Fig.8 shows the spatial power spectra of the co-estimated core and
lithospheric fields from model CL at the Earth’s mean spherical reference
surface. The spread in the posterior realizations is very small (almost
invisible in the plot) at degrees 17 up to 70, indicating that the
lithospheric field is very well constrained by the observations for these
degrees. At lower degrees the posterior spread increases, with the mean model
containing lower power than most of the individual realizations. The posterior
spread also increases above degree 70, becoming as large as the spread in the
prior above degree 110 by which point the mean model contains less power than
any of the posterior realizations. The cross-over between the mean core and
lithospheric field models also occurs at degree 15 at Earth’s surface.

Figure 8: Lowes-Mauersberger spectra showing power as a function of spherical
harmonic degree at Earth’s mean spherical reference radius in 2020.0 for model
CL. Solid black lines with circles show the estimated posterior mean core and
lithospheric field models which cross at degree 15. Grey lines show prior
realizations based on the adopted spatial covariance model, cyan lines show
posterior realizations of the core field and purple lines posterior
realizations of the lithospheric field. CHAOS-7 (up to degree 13) and LCS-1
(at degree 16 and above) are shown for reference.

In Fig. 9 we present a map of the posterior mean lithospheric field from model
CL at Earth’s surface for degree 2 to 120.

Figure 9: Map of estimated posterior mean lithospheric field at Earth’s
spherical reference radius, for degrees 2 to 120.

###  4.4 Modelled secular variation and comparisons with ground observatories

To document the time dependence in model CL and show that this is also
reasonable, Fig. 10 presents example time-series of the first time derivative
(or secular variation, SV) of spherical harmonic coefficients from model CL’s
core field model. It shows the posterior mean, 5000 posterior realizations,
and the CHAOS-7 model for reference. The time-dependent SV of $g^{0}_{1}$ is
slightly smoother in model CL than in CHAOS-7, agrees well at intermediate
degrees, following the same trends and showing some (expected) differences
close to the endpoints. Model CL generally shows lower amplitude changes in SV
at high degree. This behaviour is expected because the CHAOS-7 model tapered
its regularization to lower values at high degrees. The dispersion of the
posterior realizations also increases with spherical harmonic degree. We see
no evidence for unrealistic features in the SV coefficients of model CL.

Figure 10: First time derivative (secular variation) of selected spherical
harmonic coefficients of the core field over the model timespan from 2005 to
2020. Top left: $dg^{0}_{1}/dt$, top right: $dh^{2}_{3}/dt$, bottom left:
$dg^{1}_{7}/dt$ and bottom right: $dh^{14}_{14}/dt$. The estimated posterior
mean model is shown as the black solid line, posterior realizations in grey
and for reference CHAOS-7 in red.

This conclusion is supported by comparisons of SV from model CL with observed
SV data at ground observatories, for example as shown in Fig. 11. Model CL
predictions agree well with the annual differences of monthly means for
stations from Africa (M’Bour), Europe (Niemegk) and in the Pacific (Honolulu),
with the fit being similar to that of CHAOS-7 although with somewhat smoother
time variations. The posterior realizations give an indication of the formal
uncertainties in the SV predictions of model CL.

Figure 11: Comparison of secular variation predicted by model CL for the
posterior mean model (black solid line) and posterior realizations (blue
lines) with annual differences of revised monthly means at selected ground
observatories (thin black line with dots). CHAOS-7 is shown for reference in
red. Left column shows the time derivative of the radial field, middle column
the time derivative of the southward field component and right column the time
derivative of the eastward field component. Top row: M’bour observatory (MBO)
from low latitudes in west Africa, middle row: Niemegk observatory (NGK) from
mid-latitudes in Europe and bottom row: Honolulu observatory (HON) from mid-
latitudes in the Pacific.

##  5 Discussion and Conclusions

###  5.1 Insights from a synthetic test

A key question is the extent to which our method is able to retrieve the core
field above degree 13. To investigate this issue we carried out a synthetic
test where the true core and lithospheric fields were known. We took as input
a time-dependent core field up to spherical harmonic degree 30 from a dynamo
simulation, along with a synthetic lithospheric field based on a simulation of
the induced and remanent lithospheric magnetisation up to degree 120. Magnetic
data were synthesized at the same locations, and for the same field
components, as in the observed dataset and this was inverted using the maximum
entropy co-estimation scheme described in 2. Full details and the results from
this synthetic test are collected in Appendix A. The synthetic data used in
the test is consistent with our prior, however the utilized prior information
is rather weak (involving only the source radii and isotropic spatial
covariance functions); the main purpose of the test is to investigate what
level of separation of the core and lithospheric fields is possible using our
approach.

Fig. 14 presents the resulting spatial spectra at the CMB and at Earth’s
surface while Fig. 15 compares side-by-side the radial field at the CMB from
the estimated posterior mean core field model and the input dynamo field for
increasing truncation degrees of 13, 16, 19 and 22. The posterior mean model
obviously has less power than the dynamo synthetic truth from degree 16 to 22
with the missing power at small scales particularly obvious at low latitudes
where the dynamo solution is most complex.

Despite under-estimation of the power at degrees 16 to 22, Fig. 15 shows the
estimated posterior mean model does contain useful information on CMB field
structures above spherical harmonic degree 13. Field structures remain
coherent as power is added from degree 13 to 22, with some important details
recovered. For example, in the southern polar region, under the Australian-
Antarctic basin (near latitude 60∘S, longitude 135∘E) there is a localized
intense flux feature present in both the dynamo model and in the posterior
mean; this is weaker and smeared when the CMB field is truncated at degree 13.
Similarly in the northern polar region under Siberia and Alaska, both intense
normal polarity features and reversed flux features are better retrieved in
the posterior mean model to degree 22 compared to a more conventional model
truncated at degree 13.

At low latitudes a number of intense features are better retrieved in the
estimated posterior mean model to degree 22 than in the model truncated at
degree 13, for example the flux concentration under central Africa and the
positive-negative pair of flux patches arranged north-south across the equator
under India and under the equator south of Mexico. However some smaller scale
features, for example a positive flux concentration in the dynamo model under
the equator at 15 degrees west, are still poorly retrieved. Isolated small
scale flux features with no imprint on larger scales, and which change rapidly
in time, are not well recovered.

Flux features recovered in the posterior mean core field model up to degree 22
are however always related to features in the true dynamo field. We do not
find any evidence for artefacts due to leakage of the lithospheric field. The
reconstructed fields are by construction as simple as possible (in terms of
maximizing the information entropy) while satisfying the observational
constraints. The adopted temporal prior (spline smoothing) also involves
strong time-averaging over small length scales that is absent in the dynamo,
this will contribute to the loss of detail on small scales. We note that
interpretations in terms of the model resolution matrix (e.g. Bloxham et al.,
1989) are not straightforward here. With the available satellite data coverage
we are able to well resolve the internal field up to high spherical harmonic
degree. It is however constraints from prior information that allow us to
partially separate the core and lithospheric fields.

Variations are found in the synthetic test results regarding the amount of
small scale power in the posterior mean model as the data constraints change.
For example there is less power in the small scale field in the gap between
the CHAMP and Swarm missions when only high altitude and lower quality
satellite data from Cryosat-2 was available. This is an expected consequence
of the maximum entropy method - at times when the data constraints are weaker
the posterior mean model becomes simpler and the model uncertainty at high
degree larger.

###  5.2 Interpretations of the inferred CMB field

Returning to the model CL derived from the real observations, in Fig. 12 we
present maps of the radial magnetic field at the CMB from the posterior mean
model, for increasing truncation degrees of 13, 16, 19 and 22. As in the
synthetic test, power adds coherently with increasing degree with some
structures becoming more localized and intense. This is also the case for the
SV at the CMB, which we have examined in manner similar to Holme et al. (2011)
(not shown); SV features seen at lower degree remain coherent although more
small scale SV features begin to appear by degree 22 indicating this is at the
edge of what we are able to reliably study.

Figure 12: Maps of radial magnetic field at the core-mantle boundary in 2020.0
from the model CL posterior mean, as a function of the truncation degree, up
to degree 13 (top left), to degree 16 (top right), to degree 19 (bottom left)
and to degree 22 (bottom right). Note how the power adds coherently to
existing features already evident at lower truncation degrees.

At high latitudes strong normal flux features are found close to the where the
inner core tangent cylinder intersects the CMB (at latitudes 69.6 degrees
north and south). There is a particularly intense flux concentration under the
Taymyr Peninsula in Siberia in 2020, and a number of high amplitude features
arranged near latitude 60 degrees north under Greenland, Canada and eastern
Siberia. In the Southern polar region there are strong normal flux features
under Wilkes land in eastern Antarctica and under western Antartica near the
Antarctic peninsula. Normal flux patches localized close to the tangent
cylinder are consistent with the poloidal dipole field being produced by an
alpha-effect in energetic eddies originating in vigorous convection close to
the inner core boundary. Another possibility is that such eddies are spawned
by powerful azimuthal flows inside the tangent cylinder occasionally ejected
across the tangent cylinder (e.g. Schaeffer et al., 2017). The different
locations of the flux concentrations in the northern and southern hemisphere
seem difficult to explain in terms of purely columnar flows.

Reversed flux features are found in the northern polar region under the
Canadian Arctic, centred under the Queen Elizabeth Islands, and also under the
Nansen Basin (between Spitzbergen and the new Siberian Islands). In the
Southern polar region there is an extended but weak reverse flux feature under
Eastern Antarctica southwards from Africa, this feature, also present in
models truncated at degree 13, persists to higher degree and is interesting as
it crosses the tangent cylinder. Reverse flux features inside the tangent
cylinder could be related to a strong omega effect driven by strong azimuthal
flows inside the tangent cylinder and associated flux expulsion, similar
features have been seen in turbulent dynamos (Schaeffer et al., 2017; Sheyko
et al., 2018).

At low latitudes the strong flux feature found under the equator between
Africa and south America in models truncated at degree 13 is split into two
features, as is a normal feature under central Africa. Such splitting of low
latitude flux features has been suggested in previous attempts to retrieve the
core field above degree 13 (Baerenzung et al., 2020; Otzen et al., 2022), and
it is consistent with the patterns of core surface SV retrieved above degree
13 (Finlay et al., 2020). Many of the field concentrations at low latitudes
seem to occur in oppositely signed pairs and they are observed to move
westwards. A possible explanation could be that they are the signature of
toroidal flux being expelled from the core at low latitudes and subsequently
propagating as a wave (e.g. Aubert et al., 2013, 2022). Such features are
ubiquitous at low latitudes in strongly forced geodynamo simulations when the
viscosity is sufficiently low (Sheyko, 2014; Schaeffer et al., 2017; Aubert &
Finlay, 2019).

###  5.3 Features of the inferred large-scale lithospheric field

In the synthetic test reported in Appendix A we were also able to test the
retrieval of the lithospheric field. The lithospheric field at the Earth’s
surface generally compares well with the input synthetic truth lithospheric
field, albeit with less power below degree 15 and above degree 90 (Fig.14 and
Fig. 16). Of particular interest is whether or not any details of the
synthetic truth large scale lithospheric field were retrieved. Fig. 17 shows
that, somewhat surprisingly, some details of the large scale lithospheric
field can be recovered, albeit with reduced amplitude. The largest anomalies
in the recovered large scale lithospheric field, for example between the North
Pole and the Bering strait, near Australia and in north-Eastern Europe are
also present in the synthetic truth model. On the other hand some prominent
structures are missing or incomplete, for example in the Atlantic-Indian-
Antarctic basin or under North America, and the recovered amplitudes are too
low. Fig. 13 shows a similar map of the large scale lithospheric field from
the model CL posterior mean derived from the real data. It shows strong
anomalies in the northern part of Eastern Europe, Australia, around eastern
Antarctica and under eastern North America, which are known locations of
strong continental magnetic anomalies.

Figure 13: Estimated posterior mean large-scale lithospheric field (degrees 2
to 14) plotted on a spherical surface at Earth’s mean reference radius.

###  5.4 Limitations and future prospects

The aim of this study was to better separate core and lithospheric magnetic
fields; it has only been partially successful. The recovered small scale core
field and large scale lithospheric field lack power. Use of more informative
priors, if these can be justified, would certainly bring improvements. For
example, we used a single spatial covariance function for each source, which
assumes field structures that are statistically the same for all locations on
the sphere. The ensembles of prior fields from the dynamo and magnetisation
models could instead be used to build full covariance matrices characterising
the statistical covariances between all locations on the sphere for each
source. Such dense spatial covariance matrices have already been used by other
authors (Gillet et al., 2019; Ropp et al., 2020; Istas et al., 2023), in the
context of lower resolution core field and flow modelling. Once sufficiently
large prior ensembles are available this will be a relatively simple extension
of the method presented above. A concern with this approach is that imperfect
aspects of the dynamo and magnetisation simulations might be mapped into the
estimated field models, it was for this reason we started here by using rather
simple information from the prior ensembles.

Regarding the small scale core field, the spline-based temporal prior employed
here prevents the recovery of rapid changes on small length scales. This
limitation could be particularly serious at low latitudes. This situation
could be remedied by adopting temporal priors that better reflect the expected
physics, for example based on AR2 or AR3 processes (Gillet et al., 2013;
Sadhasivan & Constable, 2022), or perhaps using temporal statistics from high
resolution simulations (Aubert, 2023).

A variant of the approach presented here is to take the spherical harmonic
(Gauss) coefficients of the internal potential as input data for the
separation into core and lithospheric fields rather than satellite data. This
has some computational advantages and may prove useful in future applications,
further details and an example are presented in Appendix B. Our scheme could
also be applied to dedicated studies of the lithospheric field. This would
require improving the lithospheric prior to allow for more power at small
length scales, use of higher data sampling rates, and use of Swarm east-west
gradients that were not included in this study.

Much is still to be learnt regarding the small scale core field. The maps of
the posterior mean core field presented here show how information on the core
dynamo is lost when CMB fields are truncated at degree 13. On the other hand
our ability to retrieve the small scale core field, and avoid lithospheric
field contamination, depends on correctly formulating and utilizing prior
information regarding the sources. Further effort is needed on how best to
extract reliable prior information from a variety of simulations of the core
dynamo and the lithospheric magnetisation. Improved Bayesian field modelling
requires prior ensembles that are both informative and broadly representative
of the diversity of possible core and lithospheric fields.

## Availability of datasets and material

The data used in this article and the derived field models are archived at
https://doi.org/10.11583/DTU.24968763.v1. Magnetic field data from the Swarm
mission are freely available from https://earth.esa.int/web/guest/swarm/data-
access; CHAMP data are available from https://isdc.gfz-potsdam.de/champ-isdc;
Ground observatory data are available from ftp://ftp.nerc-
murchison.ac.uk/geomag/Swarm/AUX_OBS/hour/; The RC-index is available from
http://www.spacecenter.dk/files/magnetic-models/RC/; The CHAOS-7 model and its
updates are available at http://www.spacecenter.dk/files/magnetic-
models/CHAOS-7/; solar wind speed, interplanetary magnetic field, and Kp-index
are available from https://omniweb.gsfc.nasa.gov/ow.html.

###### Acknowledgements.

The authors thank Julien Aubert for providing data from his dynamo simulations
and Simon Williams for help with the lithospheric magnetisation simulations.
Nils Olsen and Lars Tøffner Clausen are thanked for numerous helpful
discussions. Andy Jackson is thanked for discussions regarding maximum entropy
and correlation functions. Nicolas Gillet and Phil Livermore provided
constructive comments on an earlier version of this work that appeared in MO’s
PhD thesis. The editor Richard Holme and two anonymous reviewers are thanked
for comments that helped us improve the manuscript. We thank the GFZ German
Research Centre for Geoscience for providing access to the CHAMP MAG-L3 data
and the European Space Agency (ESA) for providing prompt access to the Swarm
L1b data. The high resolution 1-min OMNI data were provided by the Space
Physics Data Facility (SPDF), NASA Goddard Space Flight Centre. We thank the
staff of the geomagnetic observatories and the INTERMAGNET for providing high-
quality observatory data. This work was funded by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 772561).

## References

  * Aubert (2015) Aubert, J., 2015.  Geomagnetic forecasts driven by thermal wind dynamics in the Earth’s core, Geophys. J. Int., 203(3), 1738–1751.
  * Aubert (2020) Aubert, J., 2020.  Recent geomagnetic variations and the force balance in Earth’s core, Geophys. J. Int., 221, 378–393.
  * Aubert (2023) Aubert, J., 2023.  State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth’s core, Geophysical Journal International, 235(1), 468–487.
  * Aubert & Finlay (2019) Aubert, J. & Finlay, C. C., 2019.  Geomagnetic jerks and rapid hydromagnetic waves focusing at Earth’s core surface, Nature Geosci., 12(5), 393–398.
  * Aubert & Gillet (2021) 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, 1854–1873.
  * Aubert et al. (2013) Aubert, J., Finlay, C. C., & Fournier, A., 2013.  Bottom-up control of geomagnetic secular variation by the Earth’s inner core, Nature, 502, 219–223.
  * Aubert et al. (2017) Aubert, J., Gastine, T., & Fournier, A., 2017.  Spherical convective dynamos in the rapidly rotating asymptotic regime, Journal of Fluid Mechanics, 813, 558–593.
  * Aubert et al. (2022) Aubert, J., Livermore, P., Finlay, C., Fournier, A., & Gillet, N., 2022.  A taxonomy of simulated geomagnetic jerks, Geophysical Journal International, p. (in press).
  * Backus (1988) Backus, G. E., 1988.  Bayesian inference in geomagnetism, Geophysical Journal International, 92(1), 125–142.
  * Baerenzung et al. (2020) Baerenzung, J., Holschneider, M., Wicht, J., Lesur, V., & Sanchez, S., 2020.  The kalmag model as a candidate for IGRF-13, Earth, Planets and Space, 72(1).
  * Baerenzung et al. (2022) Baerenzung, J., Holschneider, M., Saynisch-Wagner, J., & Thomas, M., 2022.  Kalmag: a high spatio-temporal model of the geomagnetic field, Earth, Planets and Space, 74(1).
  * Baratchart & Gerhards (2017) Baratchart, L. & Gerhards, C., 2017.  On the recovery of core and crustal components of geomagnetic potential fields, SIAM Journal on Applied Mathematics, 77(5), 1756–1780.
  * Barrois et al. (2017) Barrois, O., Gillet, N., & Aubert, J., 2017.  Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics, Geophysical Journal International, 211(1), 50–68.
  * Bloxham et al. (1989) Bloxham, J., Gubbins, D., & Jackson, A., 1989.  Geomagnetic secular variation, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 329(1606), 415–502.
  * Buffett & Christensen (2007) Buffett, B. A. & Christensen, U. R., 2007.  Magnetic and viscous coupling at the core—mantle boundary: inferences from observations of the Earth’s nutations, Geophysical Journal International, 171(1), 145–152.
  * Cain et al. (1989) Cain, J. C., Wang, Z., Schmitz, D. R., & Meyer, J., 1989.  The geomagnetic spectrum for 1980 and core-crustal separation, Geophysical Journal International, 97(3), 443–447.
  * Chapman & Bartels (1940) Chapman, S. & Bartels, J., 1940.  Geomagnetism, vol. I+II, Clarendon Press, Oxford, London.
  * Christensen & Tilgner (2004) Christensen, U. & Tilgner, A., 2004.  Power requirement of the geodynamo from ohmic losses in numerical and laboratory dynamos, Nature, 429, 169–171, doi: 10.1038/nature02508.
  * Constable (1988) Constable, C. G., 1988.  Parameter estimation in non-gaussian noise, Geophysical Journal International, 94(1), 131–142.
  * Finlay et al. (2012) Finlay, C. C., Jackson, A., Gillet, N., & Olsen, N., 2012.  Core surface magnetic field evolution 2000-2010, Geophysical Journal International, 189(2), 761–781.
  * Finlay et al. (2020) Finlay, C. C., Kloss, C., Olsen, N., Hammer, M. D., Tøffner-Clausen, L., Grayver, A., & Kuvshinov, A., 2020.  The CHAOS-7 geomagnetic field model and observed changes in the south atlantic anomaly, Earth, Planets and Space, 72, 156.
  * Friis-Christensen et al. (2006) Friis-Christensen, E., Lühr, H., & Hulot, G., 2006.  Swarm: A constellation to study the Earth’s magnetic field, Earth, planets and space, 58(4), 351–358.
  * Gauss (1839) Gauss, C. F., 1839.  Allgemeine theorie des Erdmagnetismus, in: Resultate aus den Beobachtungen des Magnetischen Vereins im Jahre 1838, edited by: Gauss, C. F. and Weber, W., pp. 1–57.
  * Gillet et al. (2007) Gillet, N., Jackson, A., & Finlay, C. C., 2007.  Maximum entropy regularization of time-dependent geomagnetic field models, Geophys. J. Int., 171, 1005–1016.
  * Gillet et al. (2013) Gillet, N., Jault, D., Finlay, C. C., & Olsen, N., 2013.  Stochastic modeling of the earth’s magnetic field: Inversion for covariances over the observatory era, Geochemistry, Geophysics, Geosystems, 14(4), 766–786.
  * Gillet et al. (2019) Gillet, N., Huder, L., & Aubert, J., 2019.  A reduced stochastic model of core surface dynamics based on geodynamo simulations, Geophysical Journal International, 219(1), 522–539.
  * Gneiting (2013) Gneiting, T., 2013.  Strictly and non-strictly positive definite functions on spheres, Bernoulli, 19(4), 1327–1349.
  * Gubbins (1976) Gubbins, D., 1976.  Observational constraints on the generation process of the Earth’s magnetic field, Geophysical Journal International, 47(1), 19–39.
  * Gubbins (2007) Gubbins, D., 2007.  Core‐based inversions for the main geomagnetic field, in Encyclopedia of Geomagnetism and Paleomagnetism, pp. 122–124, eds Gubbins, D. & Herrero-Bervera, E., Springer Netherlands, Dordrecht.
  * Gubbins & Bloxham (1985) Gubbins, D. & Bloxham, J., 1985.  Geomagnetic field analysis- III. Magnetic fields on the core-mantle boundary, Geophysical Journal of the Royal Astronomical Society, 80, 695–713.
  * Gull & Skilling (1984) Gull, S. & Skilling, J., 1984.  Maximum entropy method in image processing, IEE Proc. - F Radar Signal Proc., 131, 646–659.
  * Gull & Skilling (1990) Gull, S. & Skilling, J., 1990.  The MEMSYS5 User’s Manual, Maximum Entropy Data Consultants Ltd, Royston.
  * Górski et al. (2005) Górski, M., Hivon, E., Banday, A. J., Wandelt, B. D., Hansen, F. K., Reinecke, M., & Bartelmann, M., 2005.  Healpix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere, Astrophysical Journal, 622(2 I), 759–771.
  * Hemant & Maus (2005) Hemant, K. & Maus, S., 2005.  Geological modeling of the new CHAMP magnetic anomaly maps using a geographical information system technique, Journal of Geophysical Research: Solid Earth, 110(B12).
  * Hobson et al. (2010) Hobson, M., Ashdown, M., & Stolyarov, V., 2010.  Signal separation in Cosmology, in Bayesian Methods in Cosmology, pp. 45–52, eds Hobson, M. P., Jaffe, A., Liddle, A., Mukherjee, P., & Parkinson, D., Cambridge University Press.
  * Hobson & Lasenby (1998) Hobson, M. P. & Lasenby, A. N., 1998.  The entropic prior for distributions with positive and negative values, Monthly Notices of the Royal Astronomical Society, 298(3), 905–908.
  * Hobson et al. (1998) Hobson, M. P., Jones, A. W., Lasenby, A. N., & Bouchet, F., 1998.  Foreground separation methods for satellite observations of the cosmic microwave background, Monthly Notices of the Royal Astronomical Society, 300, 1–29.
  * Holme et al. (2011) Holme, R., Olsen, N., & Bairstow, F. L., 2011.  Mapping geomagnetic secular variation at the core–mantle boundary, Geophysical Journal International, 186(2), 521–528.
  * Holschneider et al. (2016) Holschneider, M., Lesur, V., Mauerberger, S., & Baerenzung, J., 2016.  Correlation-based modeling and separation of geomagnetic field components, Journal of Geophysical Research: Solid Earth, 121(5), 3142–3160.
  * Hulot et al. (2015) Hulot, G., Sabaka, T. J., Olsen, N., & Fournier, A., 2015.  The present and future geomagnetic field, in Treatise on Geophysics (Second Edition), vol. 5 - Geomagnetism, pp. 33–78, Elsevier.
  * Hyvärinen & Oja (2000) Hyvärinen, A. & Oja, E., 2000.  Independent component analysis: algorithms and applications, Neural Networks, 13(4), 411–430.
  * Istas et al. (2023) Istas, M., Gillet, N., Finlay, C., Hammer, M., & Huder, L., 2023.  Transient core surface dynamics from ground and satellite geomagnetic data, Geophysical Journal International, 233, 1890–1915.
  * Jackson (2003) Jackson, A., 2003.  Intense equatorial flux spots on the surface of Earth’s core, Nature, 464, 760–763.
  * Jackson et al. (2000) Jackson, A., Jonkers, A. R. T., & Walker, M. R., 2000.  Four centuries of geomagnetic secular variation from historical records, Philosophical Transactions of the Royal Society of London A., 358, 957–990.
  * Jackson et al. (2007) Jackson, A., Constable, C., & Gillet, N., 2007.  Maximum entropy regularization of the geomagnetic core field inverse problem, Geophys. J. Int., 171, 995–1004.
  * Jaynes (1968) Jaynes, E., 1968.  Prior probabilities, IEEE Transactions On Systems Science and Cybernetics, 4, 227–241.
  * Jaynes (2003) Jaynes, E., 2003.  Probability Theory: The Logic of Science, Cambridge University Press.
  * Korte et al. (2011) Korte, M., Constable, C., Donadini, F., & Holme, R., 2011.  Reconstructing the holocene geomagnetic field, Earth and Planetary Science Letters, 312(3), 497–505.
  * Langel & Estes (1982) Langel, R. A. & Estes, R. H., 1982.  A geomagnetic field spectrum, Geophysical Research Letters, 9(4), 250–253.
  * Langel & Hinze (1998) Langel, R. A. & Hinze, W. J., 1998.  The magnetic field of the Earth’s lithosphere: The satellite perspective, Cambridge University Press, Cambridge.
  * Langel et al. (1982) Langel, R. A., Ousley, G., Berbert, J., Murphy, J., & Settle, M., 1982.  The MAGSAT mission, Geophys. Res. Letters, 9, 243–245.
  * Lhuillier et al. (2011) Lhuillier, F., Fournier, A., Hulot, G., & Aubert, J., 2011.  The geomagnetic secular-variation timescale in observations and numerical dynamo models, Geophysical Research Letters, 38(9).
  * Lowes (1966) Lowes, F. J., 1966.  Mean-square values on sphere of spherical harmonic vector fields, J. Geophys. Res., 71, 2179.
  * Lowes (1974) Lowes, F. J., 1974.  Spatial power spectrum of the main geomagnetic field, and extrapolation to the core, Geophys. J. Roy. Astron. Soc., 36, 717–730.
  * Maisinger et al. (2004) Maisinger, K., Hobson, M., & Lasenby, A., 2004.  Maximum-entropy image reconstruction using wavelets, Mon. Not. R. Astron. Soc., 347, 339–354.
  * Masterton et al. (2012) Masterton, S. M., Gubbins, D., Müller, R. D., & Singh, K. H., 2012.  Forward modelling of oceanic lithospheric magnetization, Geophysical Journal International, 192(3), 951–962.
  * Mauersberger (1956) Mauersberger, P., 1956.  Das Mittel der Energiedichte des geomagnetischen Hauptfeldes an der Erdoberfläche und seine säkulare Änderung, Gerl. Beitr. Geophys., 65, 207–215.
  * Maus (2010) Maus, S., 2010.  Magnetic field model MF7, www.geomag.us/models/MF7.html.
  * Nataf & Ricard (1996) Nataf, H.-C. & Ricard, Y., 1996.  3smac: an a priori tomographic model of the upper mantle based on geophysical modeling, Physics of the Earth and Planetary Interiors, 95(1-2), 101–122.
  * Olsen (2002) Olsen, N., 2002.  A model of the geomagnetic field and its secular variation for epoch 2000 estimated from Ørsted data, Geophysical Journal International, pp. 454–462.
  * Olsen & Floberghagen (2018) Olsen, N. & Floberghagen, R., 2018.  Exploring Geospace from Space: the Swarm Satellite Constellation Mission, Space Research Today, 203, 61–71.
  * Olsen & Stolle (2012) Olsen, N. & Stolle, C., 2012.  Satellite geomagnetism, Annual Review of Earth and Planetary Sciences, 40, 441–465.
  * Olsen et al. (2014) Olsen, N., Lühr, H., Finlay, C. C., Sabaka, T. J., Michaelis, I., Rauberg, J., & Tøffner-Clausen, L., 2014.  The CHAOS-4 geomagnetic field model, Geophysical Journal International, pp. 817–819.
  * Olsen et al. (2017) Olsen, N., Ravat, D., Finlay, C. C., & Kother, L. K., 2017.  LCS-1: a high-resolution global model of the lithospheric magnetic field derived from CHAMP and Swarm satellite observations, Geophysical Journal International, 211(3), 1461–1477.
  * Otzen (2022) Otzen, M., 2022.  Geostatistical simulation and deep learning in geomagnetism, Ph.D. thesis, Technical University of Denmark.
  * Otzen et al. (2022) Otzen, M., Finlay, C. C., & Hansen, T. M., 2022.  Direct sequential simulation for spherical linear inverse problems, Computers and Geosciences, 160, 105026.
  * Panovska et al. (2018) Panovska, S., Constable, C. G., & Korte, M., 2018.  Extending Global Continuous Geomagnetic Field Reconstructions on Timescales Beyond Human Civilization, Geochemistry, Geophysics, Geosystems, 19(12), 4757–4772.
  * Rasmussen & Williams (2006) Rasmussen, C. E. & Williams, C. K. I., 2006.  Gaussian processes for machine learning, MIT press.
  * Reguzzoni & Sampietro (2015) Reguzzoni, M. & Sampietro, D., 2015.  GEMMA: An Earth crustal model based on GOCE satellite data, International Journal of Applied Earth Observation and Geoinformation, 35, 31–43, GOCE earth science applications and models(Based on the ESA GOCE solid earth workshop, 16-17 October 2012).
  * Reigber et al. (2002) Reigber, C., Lühr, H., & Schwintzer, P., 2002.  CHAMP mission status, Advances in Space Research, 30(2), 129–134.
  * Ropp & Lesur (2023) Ropp, G. & Lesur, V., 2023.  Mid-latitude and equatorial core surface flow variations derived from observatory and satellite magnetic data, Geophysical Journal International, 234(2), 1191–1204.
  * Ropp et al. (2020) Ropp, G., Lesur, V., Baerenzung, J., & Holschneider, M., 2020.  Sequential modelling of the Earth’s core magnetic field, Earth, Planets and Space, 72(1).
  * Sabaka et al. (2004) Sabaka, T. J., Olsen, N., & Purucker, M. E., 2004.  Extending comprehensive models of the earth’s magnetic field with ørsted and champ data, Geophysical Journal International, 159(2), 521–547.
  * Sabaka et al. (2020) Sabaka, T. J., Tøffner-Clausen, L., Olsen, N., & Finlay, C. C., 2020.  CM6: A comprehensive geomagnetic field model derived from both champ and swarm satellite observations, Earth, Planets, Space, 72, 80.
  * Sadhasivan & Constable (2022) Sadhasivan, M. & Constable, C., 2022.  A new power spectrum and stochastic representation for the geomagnetic axial dipole, Geophysical Journal International, 231(1), 15–26.
  * Schaeffer et al. (2017) Schaeffer, N., Jault, D., Nataf, H.-C., & Fournier, A., 2017.  Turbulent geodynamo simulations: a leap towards Earth’s core, Geophysical Journal International, 211(1), 1–29.
  * Shannon (1948) Shannon, C. E., 1948.  A mathematical theory of communication, The Bell System Technical Journal, 27(3), 379–423.
  * Sheyko (2014) Sheyko, A., 2014.  Numerical investigations of rotating MHD in a spherical shell, Ph.D. thesis, ETH Zurich.
  * Sheyko et al. (2018) Sheyko, A., Finlay, C., Favre, J., & Jackson, A., 2018.  Scale separated low viscosity dynamos and dissipation within the Earth’s core, Sci. Rep., 8(1), 12566.
  * Shure et al. (1985) Shure, L., Parker, R. L., & Langel, R. A., 1985.  A preliminary harmonic spline model from Magsat data, J. Geophys. Res., 90(B13), 11505–11512.
  * Skilling (1988) Skilling, J., 1988.  The axioms of maximum entropy, in Maximum-Entropy and Bayesian Methods in Science and Engineering, vol. 1, pp. 173–187, eds Erickson, G. & Smith, C., Dordrecht: Kluwer.
  * Skilling (1989) Skilling, J., 1989.  Classical Maxent, in Maximum Entropy and Bayesian Methods, pp. 45–52, ed. Skilling, J., Dordrecht: Kluwer.
  * Stockmann et al. (2009) Stockmann, R., Finlay, C., & Jackson, A., 2009.  Imaging Earth’s crustal magnetic field with satellite data: a regularized spherical triangle tessellation approach, Geophys. J. Int., 179, 929–944.
  * Tarantola (2005) Tarantola, A., 2005.  Inverse problem theory and methods for model parameter estimation, SIAM.
  * Voorhies (2004) Voorhies, C. V., 2004.  Narrow-scale flow and a weak field by the top of Earth’s core: Evidence from Ørsted, Magsat, and secular variation, J. Geophys. Res., 109, B03106.
  * Voorhies et al. (2002) Voorhies, C. V., Sabaka, T. J., & Purucker, M., 2002.  On magnetic spectra of Earth and Mars, J. Geophys. Res., 107(E6), 5034,.
  * Wahba (1990) Wahba, G., 1990.  Spline models for observational data, SIAM.
  * Whaler & Gubbins (1981) Whaler, K. A. & Gubbins, D., 1981.  Spherical harmonic analysis of the geomagnetic field: an example of a linear inverse problem, Geophysical Journal International, 65(3), 645–693.
  * Williams & Gubbins (2019) Williams, S. E. & Gubbins, D., 2019.  Origin of long-wavelength magnetic anomalies at subduction zones, Journal of Geophysical Research: Solid Earth, 124(9), 9457–9473.

##  Appendix A An application to synthetic data

Here we report results from a synthetic test designed to test the extent to
which we can retrieve the small scale core field above spherical harmonic
degree 13 and the large scale lithospheric field below degree 13 using the
scheme described in Section 2.

We use as input a time-dependent core field taken from a preliminary version
of the dynamo model assimilation runs described by (Aubert, 2023). This
simulation contained realistic core field structures and time dependence up to
degree 30 and was not contained in the prior ensemble. For the synthetic
lithospheric field we used one realizatation from our set of prior
magnetisation models, based on the forward models by Hemant & Maus (2005),
Masterton et al. (2012) and Williams & Gubbins (2019) with the perturbations
described by Otzen (2022), but generated separately and not included in the
ensemble used to construct the prior statistics. From these models we
synthesized data at the same times and positions, and with the same measured
components, as the real satellite and ground data described in Section 3.
Gaussian noise was added using a standard deviation of $3.0$nT for CHAMP,
$5.0$nT for CryoSat-2, $2.0$nT for Swarm and $3.0$nT/yr for ground observatory
data.

Inversions were then carried out in exactly the same way as for the real data.
The same priors, and same fixed regularization settings the same starting
models and the same number of iterations (24) were performed.

In Fig. 14 we present the Lowes-Mauersberger spectrum at the core surface
(top) and at the Earth’s surface (bottom), with the synthetic truth marked in
the dashed line. In Fig. 15 we present how the estimated posterior mean and
the synthetic truth for the core field changes with the truncation degree.
Although some small scale feature are lost and the amplitude is reduced, the
posterior mean model above degree 22 retrieves more details of real features.

Figure 14: Lowes-Mauersberger spectra showing power as a function of spherical
harmonic degree of magnetic fields at the CMB (top) and at Earth’s surface
spherical reference radius, (bottom) in 2020.0. The input synthetic truth
models for the core and lithospheric fields are shown in the black dashed
lines. Solid black lines with circles show the estimated posterior mean core
and lithospheric field models which cross at degree 15. Grey lines show prior
realizations based on the adopted spatial covariance model, cyan lines show
posterior realizations of the core field and purple lines posterior
realizations of the lithospheric field. CHAOS-7 (up to degree 13) and LCS-1
(at degree 16 and above) are shown for reference.

Figure 15: Maps of radial magnetic field at the core-mantle boundary from the
estimated posterior mean model derived from the synthetic test dataset (left
column) and the input truth dynamo model (right column), as a function of the
truncation degree, up to degree 13 (top row), to degree 16 (second row), to
degree 19 (third row) and to degree 22 (bottom row).

Fig. 16 presents a comparison of the estimated posterior mean lithospheric
field and Fig. 17 compares the recovered (posterior mean) large scale
lithospheric field, and the input synthetic truth for the large scale
lithospheric field at Earth’s spherical reference radius, for degrees 2-14.

Figure 16: Map of the radial component of the estimated posterior mean
lithospheric field on a spherical surface at Earth’s reference radius, for
degrees 2 to 120 (top) and the same quantity from the synthetic truth input
model based on lithospheric magnetisation models (bottom).

Figure 17: Estimated posterior mean large-scale lithospheric field (degrees 2
to 14) derived from the synthetic dataset plotted at Earth’s spherical
reference surface (top) and the same quantity from the synthetic truth model
based on a prior realization of lithospheric magnetisation (bottom). Note the
change in the scale used for the two panels.

##  Appendix B Application to Gauss coefficients of the internal field

Instead of co-estimating core and lithospheric fields models directly from
satellite observations it is possible to start with Gauss coefficients for an
internal potential field, $g^{int}_{nm}(t_{i})$ and $,h^{int}_{nm}(t_{i})$,
given at a series of reference times $t_{i}$. In this case the relevant input
data consisting of the Gauss coefficients at all times can be collected in a
vector

$\mathbf{g}=\left[g^{int}_{10}(t_{0}),g^{int}_{11}(t_{0}),h^{int}_{11}(t_{0})...,g^{int}_{10}(t_{1}),...\right]^{T}.$
(18)

If the input field model is already smooth in time, no additional temporal
prior is needed and the loss function to be minimized takes the form

$\Phi(\mathbf{m})=\frac{1}{2}\chi_{g}^{2}(\mathbf{m})-S_{tav}(\mathbf{x}^{C})-S(\mathbf{x}^{L}).$
(19)

where $\chi_{g}^{2}=\mathbf{e_{g}}^{T}\mathbf{C_{g}^{-1}}\mathbf{e_{g}}$ with
$\mathbf{e_{g}}=\mathbf{g}-\widehat{\mathbf{g}}$ and $\widehat{\mathbf{g}}$
are the predicted internal field Gauss coefficients (core plus lithosphere) at
the relevant times and the diagonal values of $\mathbf{C_{g}^{-1}}$ are
$1/\sigma_{g}^{2}$ which define how closely the input Gauss coefficients
should be matched. In the experiments reported here we set
$\sigma_{g}=10^{-5}$nT since we wanted the estimated core and lithospheric
coefficients to closely match those of the input field model. As before
$S_{tav}(\mathbf{x}^{C})$ describes the information entropy of the spatially
decorrelated CMB radial field (averaged over time), and $S(\mathbf{x}^{L})$
describes the information entropy of the spatially decorrelated lithospheric
radial field at Earth’s surface.

Results of applying this procedure to the total internal field (core plus
lithosphere) from the CL model from the main text, and to the CHAOS-7.16
internal field model are shown in Fig. 18. In these tests we started from an
initial lithospheric field set to zero below degree 16 and to the LCS-1 model
at higher degree and an initial core field set to the time-average of the
time-dependent internal field from the CHAOS-6.9 model up to degree 12 and
zero at higher degree, and we found convergence after 5 to 6 iterations.

(a) (b)

(c) (d)

Figure 18:  (a) Radial magnetic field at the core-mantle boundary in 2019.5,
from a model estimated taking as the input combined internal field spherical
harmonic coefficients (core and lithospheric fields) from the CL model
presented in the main text (top left, denoted MaxEnt-Coeffs), (b) from a model
estimated in a similar way but taking spherical harmonic coefficients from
CHAOS-7.16 as input data (top right, denoted MaxEnt-CHAOS), (c) for reference,
the result from the CL model presented in the manuscript (bottom left, denoted
MaxEnt-Sat$\\_$data), (d) Spherical harmonic power spectra for the three
models at the CMB and at Earth’s surface (bottom right).

We find the CMB field retrieved in this way is rather similar starting with
CHAOS-7.16 and the full CL model, with the minor differences in 2019.5 seen in
Fig. 18 being due to (i) the limitation of the CHAOS-7.16 time-dependent field
to degree 20 and possible aliasing related to this, and (ii) because
CHAOS-7.16 is derived using data up to 2023 while the CL model used data only
up to 2020. On the other hand we find that the large-scale lithospheric fields
estimated here from input Gauss coefficients contain slightly lower power than
that estimated from satellite observations, perhaps because there is less
freedom here when fitting the input data. We find that time-dependence of the
input internal field is crucial to the separation.

The power spectra of the combined internal field from these estimated models
is presented and compared to that of CHAOS-7.16 in Fig. 19. The spectra of the
combined internal field model derived from the CHAOS model coefficients
(labelled MaxEnt-CHAOS) matches the CHAOS spectra to machine precision - this
was enforced by construction. The CL model from the main text (labelled
MaxEnt-Sat$\\_$data) is, as expected, not identical to CHAOS, but it
nevertheless agrees well at all degrees, the difference being less than 0.1
nT2 in the power (i.e. around 0.3 nT) which is less than the likely errors in
the model coefficients. The combined CL model is thus compatible with the
observational constraints, in that it fits the satellite data from which it
was constructed, and its combined internal field agrees well with established
field models such as the CHAOS model.

Figure 19: Comparison of spatial power spectra of the combined (core plus
lithosphere) internal field for the models CHAOS-7.16 (green), the CL model
from the main text labelled MaxEnt-Sat$\\_$data (blue) and the model from this
Appendix, labelled MaxEnt-CHAOS, derived from the spherical harmonic
coefficients of CHAOS-7.16. Difference to CHAOS-7.16 are marked as dashed
lines.

The scheme outlined in this Appendix provides a means of separating a given
internal field model into core and lithospheric parts as a post-processing
procedure; this may prove useful for future applications.