Polar ionospheric currents and high temporal resolution geomagnetic field
models–Polar ionospheric currents and high temporal resolution geomagnetic
field models

# Polar ionospheric currents and high temporal resolution geomagnetic field
models

Clemens Kloss1 Christopher C. Finlay1 Karl M. Laundal2 and Nils Olsen1
1 DTU Space  Technical University of Denmark  Centrifugevej 356  2800 Kongens
Lyngby  Denmark
2 Birkeland Centre for Space Science  Department of Physics and Technology
University of Bergen  5007 Bergen  Norway

(Accepted 2023 August 13. Received 2023 August 10; in original form 2023 March
8; 2023)

###### keywords:

Core, Magnetic field variations through time, Satellite magnetics, Inverse
theory, Polar ionospheric currents

††volume: 00

Estimating high resolution models of the Earth’s core magnetic field and its
time variation in the polar regions requires that one can adequately account
for magnetic signals produced by polar ionospheric currents, which vary on a
wide range of time and length scales. Limitations of existing ionospheric
field models in the challenging polar regions can adversely affect core field
models, which in turn has important implications for studies of the core flow
dynamics in those regions. Here we implement a new approach to co-estimate a
climatological model of the ionospheric field together with a model of the
internal and magnetospheric fields within the CHAOS geomagnetic field
modelling framework. The parametrization of the ionospheric field exploits
non-orthogonal magnetic coordinates to efficiently account for the geometry of
the Earth’s magnetic field and scales linearly with external driving
parameters related to the solar wind and the interplanetary magnetic field.
Using this approach we derive a new geomagnetic field model from measurements
of the magnetic field collected by low Earth orbit satellites, which in
addition to the internal field provides estimates of the typical current
system in the polar ionosphere and successfully accounts for previously
unmodelled ionospheric signals in field model residuals. To resolve the
ambiguity between the internal and ionospheric fields when using satellite
data alone, we impose regularisation. We find that the time derivative of the
estimated internal field is less contaminated by the polar currents, which is
mostly visible in the zonal and near-zonal terms at high spherical harmonic
degrees. Distinctive patches of strong secular variation at the core-mantle
boundary, which have important implications for core dynamics, persist.
Relaxing the temporal regularisation reveals annual oscillations, which could
indicate remaining ionospheric field or related induced signals in the
internal field model. Using principal component analysis we find that the
annual oscillations mostly affect the zonal low-degree spherical harmonics of
the internal field.

##  1 Introduction

The ionospheric magnetic field is generated by electrical currents that
circulate in the Earth’s ionosphere, the electrically conducting layer of the
atmosphere from about $90\text{\,}\mathrm{km}$ to $1000\text{\,}\mathrm{km}$
altitude. The ionospheric field undergoes daily, seasonal and solar cycle
variations, which depend on solar activity and illumination [66]. In
particular in the polar regions, where the ionospheric field is highly dynamic
and very sensitive to changes in the solar wind and the Interplanetary
Magnetic Field (IMF) thanks to field-aligned currents that facilitate the
coupling to the magnetosphere, it is a difficult task to estimate accurate
ionospheric field models. Imperfect modelling and the fact that the involved
time and length scales overlap with those of the large-scale time-varying
internal field, which originates in the Earth’s core, makes the separation
between the two fields a major challenge in geomagnetic field modelling [12].
New strategies to deal with the ionospheric field are therefore crucial for
studies of the core field and its time variation. Such core field models are
used to infer flow patterns in Earth’s core and to study the geodynamo process
and its related waves and oscillations. The high latitude regions of the outer
core shell, within what is known as the inner core tangent cylinder (a
cylinder aligned with the rotation axis just touching the inner core in the
equatorial plane), play a special role in core processes because they are
dynamically separated from the remainder of the shell and so can be an
important source of equatorial symmetry breaking. Recent studies indicate
strong jet like flows near to the tangent cylinder [32], while highly time-
dependent turbulent polar vortices are expected within the tangent cylinder
[2, 58, 60]. Detailed study of such processes requires that ionospheric
signals in the polar regions are adequately separated.

Earlier studies have developed different techniques to model the ionospheric
field or have tried to reduce its effect on the recovered internal field
during geomagnetic field modelling. One common technique to minimise the
ionospheric disturbance field in geomagnetic field modelling is data
selection. By focusing on data under geomagnetic quiet conditions and by
choosing suitable magnetic components, one seeks to reduce ionospheric signals
that are not well parametrized in the model. For example, the CHAOS model [46,
44, 45, 47, 13], which is a model of the recent geomagnetic field and provides
estimates of the time-dependent and static internal fields and the quiet-time
magnetospheric field, is derived from magnetic vector and total intensity
observations using, among other criteria, a dark-time selection criterion
based on the sun elevation angle to remove the strong ionospheric disturbances
that are present under sunlit conditions. At polar latitudes only the scalar
magnitude of the field is used in an effort to minimise the effect of field-
aligned currents, which mainly disturb the field direction but not its scalar
magnitude. Similarly, in the sequential approach of [49] for modelling the
internal, quiet-time magnetospheric and associated internally-induced fields,
the magnetic vector data at mid and low latitudes are selected according to
local night-time and dark conditions, but in the polar regions vector data
from all local times are used. Although data selection is effective,
significant ionospheric signals often remain in the data, especially at polar
latitudes, where strong horizontal currents in the E-layer of the ionosphere
continue to disturb the scalar magnitude of the field at all local times
including during dark conditions [[, e.g.]]FriisChristensen2017.

In the comprehensive modelling approach [56, 55, 53, 52, 51], all major
sources of the geomagnetic field are parametrized and co-estimated in a single
step, including the magnetic field produced by ionospheric currents. In the
CM6 model [56], the latest in the series of Comprehensive Models (CM), the
magnetic field due to the currents in the E-layer of the ionosphere are
parametrized in space using special basis functions that involve projecting
spherical harmonics in the quasi-dipole coordinate system [48], which is
believed appropriate for describing these currents, whose geometry is
organised by the main magnetic field. Temporal variations are expressed in
terms of specific daily and sub-daily harmonics with periods of
$24\text{\,}\mathrm{h}$, $12\text{\,}\mathrm{h}$, $8\text{\,}\mathrm{h}$, and
$6\text{\,}\mathrm{h}$, which are further modulated with annual and semiannual
harmonics and scaled by a three-monthly average of $F_{10.7}$ solar radiation
index, which tracks long term variations in solar activity. In addition, the
model takes into account the Earth-induced field using a model of the
electrical conductivity of the Earth’s surface. Thanks to the sophisticated
parametrization, the CM models are very successful at describing the slowly
varying averaged ionospheric magnetic field at mid and low latitudes. This is
why the same parametrization has been adopted in the dedicated ionospheric
field inversion chain [6] to produce spherical harmonic models of the
ionospheric magnetic field at low-to-mid latitudes using the magnetic data
collected by the satellites of the European Space Agency’s (ESA) Swarm mission
[15]. However, the approach of using a finite set of specific harmonics may
not be as suitable for describing the polar ionospheric field, which varies on
a much wider range of frequencies in response to changes in the solar wind
speed and the IMF. In addition, it is not clear whether the basis functions
for parametrizing the ionospheric magnetic field are also appropriate at high
latitudes.

In the Kalmag geomagnetic field models [4, 3] the magnetic field associated
with ionospheric currents including field-aligned currents are also co-
estimated. These models are sequentially derived using a Kalman filter
approach after applying data selection to reduce the dayside ionospheric field
signal in the input magnetic data. For the parametrization of the ionospheric
sources they use poloidal and toroidal potentials in magnetic coordinate
systems and represent the evolution in time through random processes based on
a-priori spatio-temporal statistics. [31] estimate ionospheric currents in the
polar ionosphere as part of the first generation of the GFZ Reference Internal
Magnetic Models (GRIMM). However, they did not co-estimate the ionospheric
field but had to use a two-step procedure whereby they first derived the model
part corresponding to the internal, large-scale external and associated
induced fields from the data and then used its residuals to build the
ionospheric part of the model.

Apart from geomagnetic field models, there are dedicated ionospheric field
models such as the Average Magnetic field and Polar current System (AMPS)
model [28] that seek to better model the fields and currents in the polar
regions. Instead of using specific periodicities to model the variability of
the ionospheric disturbance field explicitly in time as in the CM models, the
AMPS model focuses on its climatological aspects, i.e. it seeks to model the
long term average of the field as a function of external driving parameters.
The AMPS model expresses the ionospheric field in terms of poloidal and
toroidal potentials, which are expanded into a global basis of spherical
harmonics. It exploits magnetic apex coordinates [48] to efficiently take into
account the geometry of the main magnetic field, which organises the large-
scale spatial structure of the ionospheric field. To express the variability
of the average ionospheric field in time, the model uses a combination of
external driving parameters related to the solar wind speed and IMF
components. It is, however, derived using vector residuals, i.e. observations
of the magnetic vector taken by the CHAllenging Minisatellite Payload (CHAMP)
and Swarm satellites after the removal of estimates of the internal and
magnetospheric fields given by the CHAOS model.

In this study we combine the climatological approach of the AMPS model for
modelling the ionospheric field with the CHAOS framework for modelling the
internal and magnetospheric fields. More specifically, we implement a co-
estimation approach, where an AMPS-type ionospheric field model is derived at
the same time as a geomagnetic field model similar to the CHAOS model. Making
use of satellite magnetic observations made by the CHAMP, CryoSat-2 and Swarm
satellites during geomagnetic quiet conditions, we derive a new model of the
geomagnetic field. Using this model, we study the quiet-time ionospheric field
and the associated electrical currents in the polar regions and go on to
investigate the effect on the time-variation of the internal field at polar
latitudes when ionospheric currents are co-estimated. In addition, we explore
cases when the temporal smoothness imposed on the internal field model is
considerably relaxed. Note that the goal is not to derive an all-purpose model
of the ionospheric field but rather to improve the time-dependent internal
field model in the CHAOS modelling framework for geomagnetic quiet conditions
in the challenging polar regions.

The paper is organised as follows. In Sect. 2 we describe the satellite
magnetic data used and the applied data selection. In Sect. 3 we provide
details about the model parametrization, giving special attention to the
ionospheric part taken from the AMPS model. There, we also give the equations
for the model estimation and the applied regularisation, and list the chosen
regularisation parameters. In Sect. 4 we evaluate the performance of the
estimated model in terms of the fit to the magnetic data, study the polar
ionospheric currents during geomagnetic quiet conditions and investigate the
recovered core field and its time variations at polar latitudes. In the last
part of that section we study the variations in time of the internal field as
given by a test model where we apply weaker temporal smoothing. We finish with
a discussion of the obtained results in Sect. 5 and the conclusions in Sect.
6.

##  2 Magnetic observations and data selection

We used vector observations of the magnetic field made by the CHAMP and
CryoSat2 satellites, and the three satellites of the Swarm constellation,
Swarm-A, Swarm-B and Swarm-C, from 2001 to the end of 2021.

From the CHAMP mission, we used the Level 3 $1\text{\,}\mathrm{Hz}$ magnetic
data, product CH-ME-3-MAG [50], between January 2001 and August 2010, which we
downsampled to $1\text{\,}\mathrm{min}$ values. We selected data according to
the recommended quality flags that are provided in the distributed CHAMP data
product files [16]. However, we did not require that both star camera heads on
the boom close to the vector magnetometer were active and provided attitude
information at the time of measurement since this created gaps in the global
distribution of magnetic data at low and mid latitudes during dusk and dawn.
More specifically, we allowed data if at least one of the two star camera
heads was available. To account for the corresponding increase in the
uncertainty of the attitude information, we chose larger a-priori attitude
errors compared to when both star camera heads were active (see Sect. 3.4.1
for details).

Concerning CryoSat2, we used fully calibrated $4\text{\,}\mathrm{s}$ magnetic
vector data from the onboard fluxgate magnetometer FGM1 (CryoSat2-1), version
0103, from August 2010 to the end of 2013. These data have been calibrated as
described in [42]. We reduced the dataset to $1\text{\,}\mathrm{min}$ values
through the following steps. First, we used estimates of the time-dependent
internal field and the CryoSat2-1 Euler angles from CHAOS-7.9 to compute
residuals in the calibrated magnetometer frame. Then, we performed a Huber-
weighted linear regression of the residuals within $20\text{\,}\mathrm{s}$
intervals and kept one fit value from each interval. Finally, we added back
the previously subtracted model estimates but retained only every third value
to obtain a reduced time series of approximately $1\text{\,}\mathrm{min}$
resolution. By using $20\text{\,}\mathrm{s}$ intervals for the linear fit, we
followed [42], who recommends averaging over five successive values to reduce
the intrinsic noise. In addition, we removed data if the attitude uncertainty
$q_{\mathrm{error}}$, which is provided in the CryoSat2 data product files,
was larger than
$40\text{\,}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{s}$.

From the Swarm mission, we made use of the Level 1b $1\text{\,}\mathrm{Hz}$
magnetic vector data from all three satellites (Swarm-A, Swarm-B and Swarm-C),
versions 0505-0508 as available, from November 2013 to the end of 2021. We
downsampled the magnetic data from each satellite to $3\text{\,}\mathrm{min}$
values to have a similar amount of data per time interval as for CHAMP and
CryoSat2.

On the entire dataset of magnetic observations, we applied several selection
criteria to focus on geomagnetic quiet-time conditions. First, we removed
gross outliers for which vector residuals with respect to the CHAOS-7.9 field
model were greater than $1000\text{\,}\mathrm{nT}$. We note that this approach
also removed magnetic signals in the data associated with field-aligned
currents, which can reach several thousands of $\mathrm{nT}$ in the polar
regions also during geomagnetic quiet-time conditions. Similarly, the
averaging of the CryoSat-2 data, as described above, removed high-frequency
ionospheric magnetic signals, i.e., signals that varied along the satellite
orbit on timescales much shorter than the 20-second interval used for
averaging. Nevertheless, since we do not expect that our approach of modelling
the average ionospheric field is able to capture intermittent high-amplitude
events, we preferred to remove these data and to process the CryoSat-2 data in
this way to improve the overall quality of the model. Next, to focus on
geomagnetically quiet conditions, we selected data for which the $\mathit{Kp}$
index [34, 35] was below 2o and the absolute rate of change of the
$\mathit{RC}$ index [47], a quantitative measure of the magnetic disturbance
at equatorial and mid-latitudes similar to the $\mathit{Dst}$ index [62], was
below $2\text{\,}\mathrm{nT}\text{\,}{\mathrm{h}}^{-1}$. Furthermore, we
selected data if, on average over $2\text{\,}\mathrm{h}$ prior to the time of
measurement, the Newell coupling function [[[]for the exact definition used in
this study, see Eq. (17)]Newell2007], measuring the rate of magnetic flux
opened at the magnetopause, was below $2.4$ and the IMF at the magnetopause
was pointing northward, i.e. having a positive z-component in the Geocentric
Solar Magnetic (GSM) frame.

The data processing and selection resulted in
$N_{d}=$2\text{,}472\text{,}746$$ vector observations, which we used for
estimating models of the geomagnetic field. To illustrate the data
distribution in time, we show in Fig. 1 a stacked histogram of the amount of
data in 3-month intervals for each satellite dataset.

Figure 1: Number of selected vector data every 3 months for each satellite
shown as stacked histogram.

We did not treat data differently depending on dark and sunlit conditions
during the model estimation to avoid seasonal variations in the data
distribution. Otherwise using, for example, only dark data for the estimation
of the internal field would adversely affect the time-dependence of the
associated model parameters, unless sufficiently smoothed in time through
regularisation. This is due to an annual variation in the data distribution,
which is created by the periodic exclusion of the data in the polar region on
the summer hemisphere. Similarly, we did not select data based on magnetic
local time for the estimation of the internal field to uniformly sample the
polar electrojets. Note that we did not use ground-based magnetic field
observations as input data for the modelling both to allow comparisons of the
model predictions with independent data and to be sure not to bias the
geographical distribution of the input data.

##  3 Model parametrization and estimation

In this paper, we largely follow the modelling approach of the CHAOS
geomagnetic field model series [46, 44, 45, 47], version CHAOS-7.9 [13].
However, a significant difference to the CHAOS models is that we also co-
estimate a model of the ionospheric currents based on the AMPS model [28]. The
following summarises the parametrization of the geomagnetic field sources that
are represented in our model and gives the equations used for the model
parameter estimation.

###  3.1 Internal magnetic field

Satellites in low Earth orbit take magnetic measurements in a region that is
free of electrical currents associated with the internal sources. In the
quasi-static approximation, the internal magnetic field can therefore be
represented by an internal scalar potential, $V^{\mathrm{int}}$, such that
$\mathbf{B}^{\mathrm{int}}=-\nabla V^{\mathrm{int}}$. In spherical coordinates
$V^{\mathrm{int}}$ is given by

$V^{\mathrm{int}}(t,r,\theta,\phi)=a\sum_{n=1}^{N^{\mathrm{int}}}\sum_{m=-n}^{n}\bigg{(}\frac{a}{r}\bigg{)}^{n+1}g_{n}^{m}(t)Y_{n}^{m}(\theta,\phi),$
(1)

where $a=$6371.2\text{\,}\mathrm{km}$$ is the mean surface radius of the
Earth, $g_{n}^{m}$ are the spherical harmonic coefficients of degree $n$ and
order $m$, $Y_{n}^{m}$ are the spherical harmonic functions, and
$N^{\mathrm{int}}=55$ is the chosen truncation degree to limit the spatial
resolution of the model. The spherical harmonic functions are defined as

$Y_{n}^{m}(\theta,\phi)\equiv\left\\{\begin{aligned}
&\cos(m\phi)P_{n}^{m}(\cos\theta),\,m\geq 0\\\ &\sin(\lvert
m\rvert\phi)P_{n}^{\lvert m\rvert}(\cos\theta),\,m<0,\end{aligned}\right.$ (2)

where $\theta$ and $\phi$ are respectively the geocentric colatitude and
longitude, and $P_{n}^{m}$ are the associated Legendre functions using the
Schmidt semi-normalization. We allow the spherical harmonic coefficients for
$n\leq 20$ to be time-dependent using a basis of 6th-order B-splines to
account for the slow time changes of the internal field

$g_{n}^{m}(t)=\sum_{k=1}^{K}g_{n,k}^{m}\mathcal{B}_{6,k}(t),$ (3)

where $\mathcal{B}_{6,k}$ ($k=1,\dots,K$) are the B-spline basis functions
defined on the model interval using a sequence of knots with a
$0.5\text{\,}\mathrm{yr}$ knot spacing and a 6-fold knot multiplicity at the
model endpoints. The coefficients for $21\leq n\leq N^{\mathrm{int}}$ are kept
constant to represent the high-degree part of the assumed static lithospheric
field.

###  3.2 External magnetic field

The sources of the external field are located in the space above the Earth’s
surface. In our model we distinguish between the magnetospheric field and the
ionospheric field. The parametrization of the magnetospheric field is
identical to the CHAOS model, whereas the one for the ionospheric field is
basically taken from the AMPS model. In this section, we will also introduce
magnetic apex coordinate systems, which are important for an efficient
parametrization of the ionospheric magnetic field.

####  3.2.1 Ionospheric field

Following the approach of [29] we write the ionospheric magnetic field as

$\mathbf{B}^{\mathrm{ion}}=\mathbf{B}^{\mathrm{pol}}+\mathbf{B}^{\mathrm{tor}}=-\nabla
V^{\mathrm{ion}}+\hat{\mathbf{r}}\times\nabla T^{\mathrm{ion}},$ (4)

where the poloidal magnetic field $\mathbf{B}^{\mathrm{pol}}$, written in
terms of the scalar potential $V^{\mathrm{ion}}$, is associated with the
currents in the ionospheric E-layer, which flow entirely below the measurement
shell of the satellites, whereas the toroidal magnetic field
$\mathbf{B}^{\mathrm{tor}}$, written in terms of the potential
$T^{\mathrm{ion}}$, is associated with the field-aligned currents that couple
the polar ionosphere to the magnetosphere (Birkeland currents). To take
advantage of the fact that the currents in the ionosphere are highly organised
with respect to the geomagnetic field, we specify the potentials in magnetic
apex coordinate systems defined by [48]. There are two of these systems:
Quasi-Dipole (QD) and Modified-Apex (MA). In QD coordinates the latitude is
defined as

$\lambda_{\mathrm{QD}}=\pm\arccos\sqrt{\frac{a+h}{a+h_{\mathrm{A}}}},$ (5)

where positive (negative) values refer to the northern (southern) magnetic
hemisphere, $h$ is the geodetic height of the point of interest, and
$h_{\mathrm{A}}$ is the geodetic height of the apex, which is the highest
point above the Earth’s ellipsoidal surface along the magnetic field line, as
given by a geomagnetic field model, that passes through the point of interest.
The longitude of the QD coordinate system is defined as the longitude of the
apex in centred dipole coordinates, a coordinate system where the z-axis
points along Earth’s magnetic dipole axis towards the northern hemisphere, the
y-axis is perpendicular to both the dipole axis and the rotation axis, and the
x-axis completes the right-handed system [30]. In Modified-Apex (MA)
coordinates the latitude is defined as

$\lambda_{\mathrm{MA}}=\pm\arccos\sqrt{\frac{a+h_{\mathrm{R}}}{a+h_{\mathrm{A}}}},$
(6)

where $h_{\mathrm{R}}$ is a chosen reference height for the mapping, which we
set to $h_{\mathrm{R}}=$110\text{\,}\mathrm{km}$$. The MA latitude is positive
for points that map to the northern magnetic hemisphere and negative
otherwise. The longitude of the MA coordinate system is identical to the QD
longitude. Since both are equal, they can be used interchangeably. Coordinates
and base vectors of the two magnetic apex coordinate systems can be
conveniently computed with the Python software package Apexpy [38], which is a
wrapper of the Fortran library by [10]. As the reference model for the field
line tracing, we used the 13th generation of the International Geomagnetic
Reference Field [[, IGRF;]]Alken2021 at epoch 2015.0 throughout the entire
model time interval.

Using a combination of the apex coordinate systems, again following [29], we
express the ionospheric potentials in terms of spherical harmonic functions

$\displaystyle V^{\mathrm{ion}}(h,\theta_{\mathrm{QD}},\phi_{\mathrm{MLT}})$
$\displaystyle=a\sum_{n=1}^{N^{\mathrm{ion}}}\sum_{\begin{subarray}{c}m=-n\\\
\lvert m\rvert\leq
M\end{subarray}}^{n}\left(\frac{a}{a+h}\right)^{n+1}g_{n}^{m,\mathrm{ion}}Y_{n}^{m}(\theta_{\mathrm{QD}},\phi_{\mathrm{MLT}})$
(7a) $\displaystyle
T^{\mathrm{ion}}(\theta_{\mathrm{MA}},\phi_{\mathrm{MLT}})$
$\displaystyle=(a+h_{\mathrm{R}})\sum_{n=1}^{N^{\mathrm{tor}}}\sum_{\begin{subarray}{c}m=-n\\\
\lvert m\rvert\leq
M\end{subarray}}^{n}T_{n}^{m,\mathrm{ion}}Y_{n}^{m}(\theta_{\mathrm{MA}},\phi_{\mathrm{MLT}}),$
(7b)

where $\theta_{\mathrm{QD}}=\frac{\uppi}{2}-\lambda_{\mathrm{QD}}$ and
$\theta_{\mathrm{MA}}=\frac{\uppi}{2}-\lambda_{\mathrm{MA}}$ are the QD and MA
colatitudes, respectively. We chose to truncate the spherical harmonic
representations at $N^{\mathrm{ion}}=45$ and $N^{\mathrm{tor}}=65$. In
addition, we used a maximum spherical harmonic order of $M=3$ for both
potentials in agreement with [28]. Instead of the QD and MA longitudes, we
used the Magnetic Local Time (MLT)

$\phi_{\mathrm{MLT}}=\phi_{\mathrm{QD}}-\phi_{\mathrm{noon}}+\uppi,$ (8)

where $\phi_{\mathrm{noon}}$ is the QD longitude of the subsolar point,
computed on a sphere with radius $r\gg a$, in practice $r=50a$. Using MLT
takes account of the fact that the ionospheric field stays fixed with respect
to the sun. By writing $T^{\mathrm{ion}}$ only in dependence of the MA
latitude and MLT, we assume the potential to be constant along the IGRF
magnetic field lines; we do not however impose north-south symmetry.

By inserting QD colatitude and MLT into the spherical harmonic functions in
Eq. 7a, we assume that $V^{ion}$ defines a harmonic potential in the source-
free region. To test this assumption, we performed numerical computations and
found that $-\nabla V^{\mathrm{ion}}$ is approximately but not strictly
divergence-free. The deviations from zero, which are largest in the auroral
regions and along the magnetic dip equator, are usually smaller in absolute
value than $1\text{\,}\mathrm{nT}$. This is smaller than typical errors due to
other unmodelled sources, which remain larger in the polar regions despite co-
estimating a climatological model of the ionospheric magnetic field. Our
conclusions from these tests is therefore that $V^{\mathrm{ion}}$, organised
in magnetic apex coordinates and magnetic local time, satisfactorily
approximates a potential field and is useful for parametrizing the geometry of
the ionospheric magnetic field and current densities.

Inserting the expressions for the potentials into Eq. (4) and evaluating the
gradients yields

$\displaystyle\mathbf{B}^{\mathrm{pol}}$
$\displaystyle=-\frac{1}{(a+h)\sin\theta_{\mathrm{QD}}}\frac{\partial
V^{\mathrm{ion}}}{\partial\phi_{\mathrm{MLT}}}\mathbf{f}_{2}\times\hat{\mathbf{k}}-\frac{1}{a+h}\frac{\partial
V^{\mathrm{ion}}}{\partial\theta_{\mathrm{QD}}}\mathbf{f}_{1}\times\hat{\mathbf{k}}-\sqrt{\lvert\mathbf{f}_{1}\times\mathbf{f}_{2}\rvert}\frac{\partial
V^{\mathrm{ion}}}{\partial h}\hat{\mathbf{k}}$ (9a)
$\displaystyle\mathbf{B}^{\mathrm{tor}}$
$\displaystyle=\frac{1}{(a+h_{\mathrm{R}})\sin\theta_{\mathrm{MA}}}\frac{\partial
T^{\mathrm{ion}}}{\partial\phi_{\mathrm{MLT}}}\hat{\mathbf{k}}\times\mathbf{d}_{1}+\frac{\sqrt{4-3\sin^{2}\theta_{\mathrm{MA}}}}{2(a+h_{\mathrm{R}})\cos\theta_{\mathrm{MA}}}\frac{\partial
T^{\mathrm{ion}}}{\partial\theta_{\mathrm{MA}}}\hat{\mathbf{k}}\times\mathbf{d}_{2}$
(9b)

where $\\{\mathbf{d}_{1},\mathbf{d}_{2},\mathbf{f}_{1},\mathbf{f}_{2}\\}$ are
the non-orthogonal base vectors for the magnetic apex coordinate systems [30],
and $\hat{\mathbf{k}}$ is a unit vector in the geodetic upward direction.
Here, we used that $\hat{\mathbf{r}}\approx\hat{\mathbf{k}}$, which means that
the expressions are best suited for describing the ionospheric field in the
polar regions. Nevertheless we assume that they also approximate the
ionospheric field at low latitudes well. Note that the last term in Eq. (9a)
is multiplied with $\sqrt{\lvert\mathbf{f}_{1}\times\mathbf{f}_{2}\rvert}$
based on the assumption that the vertical component of the ionospheric field
scales with the linear dimension of the horizontal current system [48].

The ionospheric magnetic field can be related to an electric sheet current
density (in units of $\mathrm{A}\text{\,}{\mathrm{m}}^{-1}$) that flows at a
fixed height, chosen to be $h_{\mathrm{R}}$, written in the form of

$\mathbf{J}^{\mathrm{sh}}=\mathbf{J}^{\mathrm{df}}+\mathbf{J}^{\mathrm{cf}}=\hat{\mathbf{k}}\times\nabla\psi^{\mathrm{df}}+\nabla\psi^{\mathrm{cf}},$
(10)

where $\mathbf{J}^{\mathrm{df}}$ is the divergence-free part of the sheet
current density associated with $\mathbf{B}^{\mathrm{pol}}$ and
$\mathbf{J}^{\mathrm{cf}}$ is the curl-free part associated with
$\mathbf{B}^{\mathrm{tor}}$. The potentials of the sheet current density parts
are

$\displaystyle\psi^{\mathrm{df}}(t,\theta_{\mathrm{QD}},\phi_{\mathrm{MLT}})$
$\displaystyle=-\frac{a}{\mu_{0}}\sum_{n=1}^{N^{\mathrm{ion}}}\sum_{\begin{subarray}{c}m=-n\\\
\lvert m\rvert\leq
M\end{subarray}}^{n}\frac{2n+1}{n}\left(\frac{a}{a+h_{\mathrm{R}}}\right)^{n+1}g_{n}^{m,\mathrm{ion}}(t)Y_{n}^{m}(\theta_{\mathrm{QD}},\phi_{\mathrm{MLT}})$
(11a)
$\displaystyle\psi^{\mathrm{cf}}(t,\theta_{\mathrm{MA}},\phi_{\mathrm{MLT}})$
$\displaystyle=-\frac{a+h_{\mathrm{R}}}{\mu_{0}}\sum_{n=1}^{N^{\mathrm{tor}}}\sum_{\begin{subarray}{c}m=-n\\\
\lvert m\rvert\leq
M\end{subarray}}^{n}T_{n}^{m,\mathrm{ion}}(t)Y_{n}^{m}(\theta_{\mathrm{MA}},\phi_{\mathrm{MLT}}),$
(11b)

which were derived by treating the apex coordinates as if they were
orthogonal, following the approach of [28]. The curl-free part of the sheet
current density can be furthermore related to an upward current density
$J_{u}$ (in units of $\mathrm{A}\text{\,}{\mathrm{m}}^{-2}$) through
$J_{u}=-\nabla\cdot\mathbf{J}^{\mathrm{cf}}$, a statement of current
continuity, which yields at the reference height

$J_{u}(t,\theta_{\mathrm{MA}},\phi_{\mathrm{MLT}})=-\frac{1}{\mu_{0}(a+h_{\mathrm{R}})}\sum_{n=1}^{N^{\mathrm{tor}}}\sum_{\begin{subarray}{c}m=-n\\\
\lvert m\rvert\leq
M\end{subarray}}^{n}n(n+1)T_{n}^{m,\mathrm{ion}}(t)Y_{n}^{m}(\theta_{\mathrm{MA}},\phi_{\mathrm{MLT}}).$
(12)

At polar latitudes, where the magnetic field lines are close to vertical,
$J_{u}$ can be interpreted as field-aligned currents and
$\mathbf{J}^{\mathrm{cf}}$ as the horizontal closure of these currents in the
form of a sheet current.

Instead of parametrizing the expansion coefficients $g_{n}^{m,\mathrm{ion}}$
and $T_{n}^{m,\mathrm{ion}}$ in time explicitly, we followed the
climatological approach of the AMPS model [28] and wrote these coefficients as
linear combinations of external driving parameters $X_{i}$ ($i=1,\dots,19$) so
that

$g_{n}^{m,\mathrm{ion}}(t)=g_{n,0}^{m,\mathrm{ion}}+\sum_{i=1}^{19}g_{n,i}^{m,\mathrm{ion}}X_{i}(t)$
(13)

and similarly for $T_{n}^{m,\mathrm{ion}}$. The $X_{i}$ are combinations of
solar wind parameters and IMF components that have been found suitable for
characterising the external driving of the ionospheric current system [28, 65]

$\displaystyle X_{1}=\sin\theta_{\mathrm{c}}$ $\displaystyle
X_{2}=\cos\theta_{\mathrm{c}}$ $\displaystyle X_{3}=\epsilon$ $\displaystyle
X_{4}=\epsilon\sin\theta_{\mathrm{c}}$ (14) $\displaystyle
X_{5}=\epsilon\cos\theta_{\mathrm{c}}$ $\displaystyle
X_{6}=\beta_{\mathrm{tilt}}$ $\displaystyle
X_{7}=\beta_{\mathrm{tilt}}\sin\theta_{\mathrm{c}}$ $\displaystyle
X_{8}=\beta_{\mathrm{tilt}}\cos\theta_{\mathrm{c}}$ $\displaystyle
X_{9}=\epsilon\beta_{\mathrm{tilt}}$ $\displaystyle
X_{10}=\epsilon\beta_{\mathrm{tilt}}\sin\theta_{\mathrm{c}}$ $\displaystyle
X_{11}=\epsilon\beta_{\mathrm{tilt}}\cos\theta_{\mathrm{c}}$ $\displaystyle
X_{12}=\tau$ $\displaystyle X_{13}=\tau\sin\theta_{\mathrm{c}}$ $\displaystyle
X_{14}=\tau\cos\theta_{\mathrm{c}}$ $\displaystyle
X_{15}=\tau\beta_{\mathrm{tilt}}$ $\displaystyle
X_{16}=\tau\beta_{\mathrm{tilt}}\sin\theta_{\mathrm{c}}$ $\displaystyle
X_{17}=\tau\beta_{\mathrm{tilt}}\cos\theta_{\mathrm{c}}$ $\displaystyle
X_{18}=F_{10.7}$ $\displaystyle X_{19}=\mathit{SML},$

which are all functions of time. The terms in Eq. (14) involve the clock angle

$\theta_{\mathrm{c}}=\operatorname{arctan2}(B_{\mathrm{IMF},y},B_{\mathrm{IMF},z}),$
(15)

where the components of the IMF, $B_{\mathrm{IMF},y}$ and
$B_{\mathrm{IMF},z}$, are with respect to the GSM frame; the dipole tilt angle

$\beta_{\mathrm{tilt}}=\arcsin(\hat{\mathbf{s}}\cdot\hat{\mathbf{m}}_{\mathrm{dip}}),$
(16)

where $\hat{\mathbf{s}}$ is a unit vector in the direction of the sun and
$\hat{\mathbf{m}}_{\mathrm{dip}}$ is the dipole moment of the IGRF magnetic
field, parametrizes seasonal effects; the solar wind-magnetospheric coupling
function [40]

$\epsilon=${10}^{-3}$\lvert
v_{\mathrm{sw}}\rvert^{4/3}{B_{t}}^{2/3}\sin^{8/3}\frac{\lvert\theta_{\mathrm{c}}\rvert}{2},$
(17)

where $B_{t}=\sqrt{B^{2}_{\mathrm{IMF},y}+B^{2}_{\mathrm{IMF},z}}$ (given in
$\mathrm{n}\mathrm{T}$) and $v_{\mathrm{sw}}$ (given in
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$) is the solar wind velocity component
antiparallel to the x-axis of the GSM frame, maximises for southward IMF and
measures the rate of reconnection on the dayside magnetopause; and the
coupling function

$\tau=${10}^{-3}$\lvert
v_{\mathrm{sw}}\rvert^{4/3}{B_{t}}^{2/3}\cos^{8/3}\frac{\theta_{\mathrm{c}}}{2}$
(18)

maximises for northward IMF and measures the rate of lobe reconnection in the
magnetotail. To approximate the delay of the near-Earth space environment to
adjust to changes in the external driving, we used $20\text{\,}\mathrm{min}$
moving averages of $\epsilon$, $\tau$ and the clock angle, based on
$1\text{\,}\mathrm{min}$ values propagated to the magnetopause as provided by
the OMNI database [27]. The solar radiation index $F_{10.7}$ in units of solar
flux,
$$\mathrm{sfu}$\equiv${10}^{-22}\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{Hz}}^{-1}$$,
parametrizes solar cycle variations. Expanding on the original parametrization
of AMPS, we included as $X_{19}$ the $\mathit{SML}$ index [39], developed by
the SuperMAG initiative [17, 18], to parametrize indirectly driven currents in
the polar ionosphere. With typically more than 100 contributing ground-based
magnetometer stations, the $\mathit{SML}$ index can be considered an extension
of the traditionally used $\mathit{AL}$ index, which is based on only 12
stations, to monitor nightside auroral activity. Fig. 2 shows stacked
histograms of the number of selected magnetic vector data in dependence of the
external driving parameters at the time of measurement.

Figure 2: Stacked histograms showing the number of selected magnetic vector
data in dependence of the driving parameters for the ionospheric field at the
time of measurement. The colours indicate the different satellite data sets.

####  3.2.2 Magnetospheric field

The magnetic field $\mathbf{B}^{\mathrm{mag}}$ produced by electric currents
in the magnetosphere can be separated into contributions due to the ring
current in the near-Earth magnetosphere, $\mathbf{B}^{\mathrm{near}}$, and the
currents in the remote magnetosphere, $\mathbf{B}^{\mathrm{far}}$, so that
$\mathbf{B}^{\mathrm{mag}}=\mathbf{B}^{\mathrm{near}}+\mathbf{B}^{\mathrm{far}}$.
We present the parametrization of each contribution in detail below.

The magnetic field produced by the ring current in the near-Earth
magnetosphere is written as $\mathbf{B}^{\mathrm{near}}=-\nabla
V^{\mathrm{near}}$ using a scalar potential in Solar Magnetic (SM)
coordinates, where the z-axis is anti-parallel to the dipole axis of the
Earth’s magnetic field, the x-axis is in the plane spanned by the dipole axis
and the Earth-Sun line, and the y-axis completes the right-handed system [30].
The scalar potential is given by

$\displaystyle V^{\mathrm{near}}(t,r,\theta_{\mathrm{SM}},\phi_{\mathrm{SM}})$
$\displaystyle=a\sum_{n=1}^{N^{\mathrm{near}}}\sum_{m=-n}^{n}\left(\frac{r}{a}\right)^{n}q_{n}^{m,\mathrm{SM}}(t)Y_{n}^{m}(\theta_{\mathrm{SM}},\phi_{\mathrm{SM}})$
(19)
$\displaystyle+a\sum_{m=-1}^{1}\hat{q}_{1}^{m,\mathrm{SM}}\bigg{[}\mathit{RC}_{i}(t)\left(\frac{a}{r}\right)^{2}+\mathit{RC}_{e}(t)\left(\frac{r}{a}\right)\bigg{]}Y_{1}^{m}(\theta_{\mathrm{SM}},\phi_{\mathrm{SM}})$
$\displaystyle+\text{Earth-induced counterpart}$

where $N^{\mathrm{near}}=2$ is the chosen truncation degree,
$\hat{q}_{1}^{m,\mathrm{SM}}$ are constant regression parameters multiplying
the $\mathit{RC}$ index, which consists of an internal part,
$\mathit{RC}_{i}$, and an external part, $\mathit{RC}_{e}$, so that
$\mathit{RC}=\mathit{RC}_{i}+\mathit{RC}_{e}$. We estimated the spherical
harmonic coefficients $q_{n}^{m,\mathrm{SM}}$ with $n=1$, called $RC$-baseline
corrections, in bins of 30 days, except for a single bin covering the period
from August 2010 to January 2014, when only platform magnetometer data of
CryoSat2 was available. The coefficients for $n=2$ were treated as constants
over the entire model time interval. The potential of the internally induced
field was not estimated separately but coupled to the external potential by
means of Q-responses, which are based on models of Earth’s electrical
conductivity. These Q-responses have also been used for the decomposition of
the $\mathit{RC}$ index into internal and external parts. The reader is
referred to [13] for details concerning the treatment of induced fields in
CHAOS-7, which is also the approach used here.

The magnetic field produced by the remote sources in the magnetosphere,
assumed to primarily be the magnetopause and magnetotail currents, is written
as $\mathbf{B}^{\mathrm{far}}=-\nabla V^{\mathrm{far}}$ using an axisymmetric,
static scalar potential in GSM coordinates, where the x-axis points sunward
along the Earth-Sun line, the z-axis is contained within the dipole axis and
the Earth-Sun line, while the y-axis completes the right-handed system [30].
The potential is given by

$\displaystyle V^{\mathrm{far}}(r,\theta_{\mathrm{GSM}},\phi_{\mathrm{GSM}})$
$\displaystyle=a\sum_{n=1}^{N^{\mathrm{far}}}\left(\frac{r}{a}\right)^{n}q_{n}^{0,\mathrm{GSM}}Y_{n}^{0}(\theta_{\mathrm{GSM}},\phi_{\mathrm{GSM}})$
(20) $\displaystyle+\text{Earth-induced counterpart}$

where $N^{\mathrm{mag}}=2$ is the chosen truncation degree. The treatment of
the internally induced part is similar to the approach used for
$\mathbf{B}^{\mathrm{near}}$.

###  3.3 Alignment parameters

We estimate alignment parameters in the form of three Euler angles $\alpha$,
$\beta$ and $\gamma$ for each satellite to rotate the magnetic vector
components from the frame of the Vector Field Magnetometer (VFM) to the Common
Reference Frame (CRF), which is defined by the orientation of the onboard star
cameras. The alignment can be written in matrix notation as

$\mathbf{B}_{\mathrm{CRF}}=\mathbf{R}_{3}(\gamma)\mathbf{R}_{2}(\beta)\mathbf{R}_{1}(\alpha)\mathbf{B}_{\mathrm{VFM}},$
(21)

where $\mathbf{B}_{\mathrm{CRF}}$ and $\mathbf{B}_{\mathrm{VFM}}$ are column
vectors that contain the magnetic field components with respect to the CRF and
the VFM frame, respectively, and $\mathbf{R}_{1}$, $\mathbf{R}_{2}$ and
$\mathbf{R}_{3}$ are rotation matrices given by

$\mathbf{R}_{1}(\alpha)=\begin{pmatrix}1&0&0\\\ 0&\cos\alpha&-\sin\alpha\\\
0&\sin\alpha&\cos\alpha\\\
\end{pmatrix},\quad\mathbf{R}_{2}(\beta)=\begin{pmatrix}\cos\beta&0&\sin\beta\\\
0&1&0\\\ -\sin\beta&0&\cos\beta\\\
\end{pmatrix},\quad\mathbf{R}_{3}(\gamma)=\begin{pmatrix}\cos\gamma&-\sin\gamma&0\\\
\sin\gamma&\cos\gamma&0\\\ 0&0&1\\\ \end{pmatrix}.$ (22)

Another rotation based on the quaternions provided in the data product files,
which describe the rotation from CRF to an Earth-fixed frame in dependence on
satellite position and orientation, was then performed to obtain the vector
magnetic field in terms of geocentric spherical components. We estimated the
above Euler angles in bins of 30 days to allow for time variations.

###  3.4 Model estimation

The model parameters are arranged into a column vector
$\mathbf{m}=[\mathbf{p}^{\mathrm{T}},\mathbf{q}^{\mathrm{T}}]^{\mathrm{T}}$,
where the column vector $\mathbf{p}$ contains the parameters of the
geomagnetic field model and the column vector $\mathbf{q}$ contains the Euler
angles for the alignment of the magnetic vector observations. We solved for
the model parameter vector by iteratively minimising the following cost
function using a quasi-Newton scheme

$\Phi(\mathbf{m})=[\mathbf{g}(\mathbf{p})-\mathbf{d}(\mathbf{q})]^{\mathrm{T}}\mathbf{C}_{d}^{-1}[\mathbf{g}(\mathbf{p})-\mathbf{d}(\mathbf{q})]+\mathbf{m}^{\mathrm{T}}\boldsymbol{\Lambda}\mathbf{m}$
(23)

where $\mathbf{g}(\mathbf{p})$ is a column vector containing the model
estimates of the magnetic vector components, $\mathbf{d}(\mathbf{q})$ is a
column vector containing the aligned magnetic vector observations expressed in
terms of spherical geocentric components, $\mathbf{C}_{d}^{-1}$ is the inverse
of the data error covariance matrix, and $\boldsymbol{\Lambda}$ is the model
regularisation matrix. At each iteration $k$ the model parameter vector is
updated through

$\mathbf{m}_{k+1}=\mathbf{m}_{k}+(\mathbf{G}_{k}^{\mathrm{T}}\mathbf{C}_{\mathrm{d}}^{-1}\mathbf{G}_{k}+\mathbf{\Lambda})^{-1}[\mathbf{G}_{k}^{\mathrm{T}}\mathbf{C}_{\mathrm{d}}^{-1}(\mathbf{d}_{k}-\mathbf{g}_{k})-\mathbf{\Lambda}\mathbf{m}_{k}],$
(24)

where $\mathbf{g}_{k}=\mathbf{g}(\mathbf{p}_{k})$,
$\mathbf{d}_{k}=\mathbf{d}(\mathbf{q}_{k})$ and $\mathbf{G}_{k}$ is the matrix
of partial derivatives of the residuals with respect to the model parameter
vector

$\big{(}\mathbf{G}_{k}\big{)}_{ij}=\frac{\partial[\mathbf{g}(\mathbf{p})-\mathbf{d}(\mathbf{q})]_{i}}{\partial(\mathbf{m})_{j}}\bigg{|}_{\mathbf{m}=\mathbf{m}_{k}}.$
(25)

####  3.4.1 Data error covariances

In the data error covariance matrix we account for the instrument error and
the uncertainty in the attitude information provided by the star trackers. The
error contributions are most conveniently described in the B23 frame, which is
defined by unit base vectors in the direction of $\mathbf{B}$,
$\hat{\mathbf{n}}\times\mathbf{B}$ and
$\mathbf{B}\times(\hat{\mathbf{n}}\times\mathbf{B})$, where $\hat{\mathbf{n}}$
is the star camera bore sight assumed not parallel to $\mathbf{B}$. In this
reference frame the data error covariance matrix is diagonal and given by [21]

$\mathbf{C}_{\mathrm{B23}}=\mathrm{diag}[\sigma^{2},\sigma^{2}+B^{2}\chi^{2}-(\chi^{2}-\psi^{2})(\hat{\mathbf{n}}\cdot\mathbf{B})^{2},\sigma^{2}+B^{2}\psi^{2}]$
(26)

where $\sigma$ (in $\mathrm{nT}$) is an isotropic instrument error in the
vector component magnitudes, $\chi$ (in radians) is an error in the attitude
about $\hat{\mathbf{n}}$, $\psi$ (in radians) is an error in the attitude
about the two axes perpendicular to $\hat{\mathbf{n}}$. In the B23 frame, we
multiplied the inverse of the data error covariance matrix, which is diagonal,
by the Huber weights [22, 7], which we recomputed from the residuals at each
iteration. The use of Huber weights allows the robust estimation of model
parameters in the presence of long-tailed error distributions due to sources
of error besides the instrument and attitude errors. We also applied a
$\sin\theta$ weighting to compensate for the larger amount of data near the
poles due to the high-inclination orbits of the satellites. Tab. 1 gives an
overview of the used a-priori instrument and attitude errors, which are based
on results of previous modelling efforts, most notably the CHAOS model series.

Table 1: Adopted instrument and attitude errors for the satellite datasets. |  $\sigma$ ($\mathrm{nT}$) |  $\psi$ (arcsec) |  $\chi$ (arcsec)
---|---|---|---
Dataset |  |  |
CHAMP | 2.5 | 10 | 10†
CryoSat2-1 | 6.0 | 30 | 30
Swarm-A | 2.2 | 5 | 5
Swarm-B | 2.2 | 5 | 5
Swarm-C | 2.2 | 5 | 5
† When both head units of the star camera are active, otherwise
$60\text{\,}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{s}$.

The star camera bore sight $\hat{\mathbf{n}}$, which is aligned with the
z-axis of the CRF frame, was taken from the data product files. Note that star
cameras often consist of several head units, in which case the bore sight
direction is a weighted average of the directions of the individual head units
that were active at the time of measurement. Also note that in our case, the
value of $\hat{\mathbf{n}}$ is in fact arbitrary since we assume $\psi$ and
$\chi$ to be equal. However, $\hat{\mathbf{n}}$ is important in the case of
CHAMP when only one of the two head units was active at a time.

####  3.4.2 Model regularisation

The model regularisation matrix $\boldsymbol{\Lambda}$ aids the convergence of
the model estimation by applying smoothing penalties on the model parameters.
It is a block-diagonal matrix, where each block corresponds to a penalty
measure scaled with an adjustable parameter, the regularisation parameter. To
reduce the temporal variation of the internal field, we used a regularisation
term based on the squared value of the third time-derivative of the radial
internal field at the Core-Mantle Boundary (CMB)
($r=$3485\text{\,}\mathrm{km}$$), averaged over both the entire model time
interval and the CMB, and another regularisation term based on the squared
value of the second time-derivative of the radial internal field at the CMB,
evaluated at the model endpoints, $t_{\mathrm{s}}=2001.0$ and
$t_{\mathrm{e}}=2022.0$ in units of decimal years, and averaged over the CMB.
The corresponding regularisation parameters are $\lambda_{t}$,
$\lambda_{t_{s}}$ and $\lambda_{t_{e}}$, respectively. The time variation of
each $\mathit{RC}$-baseline correction,
$\\{q_{1}^{-1,\mathrm{SM}},q_{1}^{0,\mathrm{SM}},q_{1}^{1,\mathrm{SM}}\\}$,
was minimised using a quadratic norm of the bin-to-bin differences, which is
scaled by the regularisation parameter $\lambda_{\mathrm{mag}}$. For the
ionospheric field, we implemented two regularisation terms. For the first
term, instead of directly applying a regularisation on the poloidal
ionospheric magnetic field $\mathbf{B}^{\mathrm{pol}}$, we designed a
quadratic norm based on the associated divergence-free sheet currents in the
ionospheric E-layer. More specifically, this regularisation term is based on
the squared magnitude of the average divergence-free sheet currents as seen by
an Earth-fixed observer, integrated over the spherical surface, which can be
written as a quadratic form

$\mathbf{m}^{\mathrm{T}}\boldsymbol{\Lambda}^{\mathrm{pol}}\mathbf{m}=\sum_{s\in\\{r,\theta,\phi\\}}\frac{1}{4\uppi}\int_{S(r_{0})}\left[\frac{1}{N_{d}}\sum_{i=1}^{N_{d}}J_{s}^{\mathrm{df}}(t_{i},r_{0},\theta,\phi)\right]^{2}\sin\theta\mathrm{d}\theta\mathrm{d}\phi,$
(27)

where $S(r_{0})$ is the spherical surface of radius $r=r_{0}\equiv
a+h_{\mathrm{R}}$, $J_{s}^{\mathrm{df}}$ with $s\in\\{r,\theta,\phi\\}$ are
the geocentric spherical components of the divergence-free sheet current
density [see Eqs. (10) and (11a)] as given by the model. The surface integral
was implemented by, first, computing the components of the divergence-free
sheet current density on a Gauss-Legendre grid in spherical geocentric
coordinates given the external driving parameter values at the times $t_{i}$
in the input dataset, then, forming the arithmetic mean of each component,
and, finally, integrating the sum of the squared component means over the
sphere using the integration weights. When strongly enforced by choosing a
large value of the associated regularisation parameter
$\lambda_{\mathrm{pol}}$, the regularisation pushes to zero the component
means of the divergence-free sheet current density with respect to an Earth-
fixed frame. Note that the currents can still change in time as required by
the magnetic data but only to the extent that the time average remains small.
We found that this form of regularisation helps to resolve the ambiguity
between the internal field and the poloidal ionospheric field, which is caused
by the fact that both fields have sources that are internal with respect to
the satellites. Without the regularisation, the internal field showed
artefacts in the form of near-zonal patterns that were almost time-invariant
and parallel to lines of constant QD latitude and the divergence-free sheet
currents were organised into a single cell of current encircling the magnetic
poles, very different from the expected configuration consisting of two cells
of current separated by the noon-midnight meridian. Regarding the second
regularisation term, for the toroidal ionospheric field, we followed the AMPS
model by using a regularisation term based on the spatial power spectrum of
the toroidal field to prevent large amplitudes close to the magnetic dip
equator, where the mapping of points at satellite altitude to MA coordinates
leaves a gap in $\theta_{\mathrm{MA}}$. The associated regularisation matrix
$\boldsymbol{\Lambda}^{\mathrm{tor}}$ is diagonal with entries
$\frac{n(n+1)}{2n+1}$, which depend on the degree of the expansion
coefficients $T_{n}^{m,\mathrm{ion}}$. The associated regularisation parameter
is $\lambda_{\mathrm{tor}}$.

We derived three geomagnetic field models: the first model, referred to as
Model-A, accounts for the ionospheric field and is our preferred model,
whereas the second model, denoted Reference, is identical except for omitting
the ionospheric part. The third model, Model-B is identical to Model-A, but we
reduced the temporal regularisation of the internal field part of the model.
Tab. 2 summarises the parametrization of the three models and gives the
numerical values of the regularisation parameters used in this study.

Table 2: Summary of the model parametrization and the chosen numerical values of the regularisation parameters for the three estimated models. The description of the parametrization is further divided into spatial (S) and temporal (T), when applicable. Internal field |  |  |
---|---|---|---
Time-dependent field |  S: Spherical harmonics in geographic coordinates ($n\leq 20$)
|  T: 6th-order B-splines, $0.5\text{\,}\mathrm{yr}$ knot spacing, 6-fold
endpoint multiplicity
Static field |  S: Spherical harmonics in geographic coordinates ($21\leq n\leq 55$)
|  T: Static in geographic coordinates
Ionospheric field |  |  |
|  Reference | Model-A/Model-B
Poloidal field |  n/a |  S: Spherical harmonics in QD/MLT,
|  |  S: $n\leq 45$, $m\leq 3$
|  |  T: 19 external driving parameters + constant
Toroidal field |  n/a |  S: Spherical harmonics in MA/MLT,
|  |  S: $n\leq 65$, $m\leq 3$
|  |  T: 19 external driving parameters + constant
Magnetospheric field |  |  |
Near-magnetospheric field |  S: Spherical harmonics in SM, $n\leq 2$
|  T: Degree-1 spherical harmonic coefficients scaled by hourly $\mathit{RC}$
index,
|  T: degree-2 coefficients static in SM, $\mathit{RC}$-baseline corrections
estimated
|  T: in bins of 30 days, except for a single bin between 2010-08/2014-01
Far-magnetospheric field |  S: Spherical harmonics in GSM, $n\leq 2$, $m=0$
|  T: Static in GSM
Alignment |  |  |
CHAMP | 3 Euler angles estimated in 118 bins of 30 days length
Swarm-A | 3 Euler angles estimated in 99 bins of 30 days length
Swarm-B | 3 Euler angles estimated in 99 bins of 30 days length
Swarm-C | 3 Euler angles estimated in 99 bins of 30 days length
CryoSat2-1 | 3 Euler angles estimated in 42 bins of 30 days length
Regularisation |  |  |
|  Reference |  Model-A |  Model-B
$\lambda_{t}$ ($\mathrm{[}\mathrm{n}\mathrm{T}\mathrm{/}\mathrm{yr}\mathrm{{}^{3}}\mathrm{]}^{-2}$) |  $1.0$ |  $1.0$ |  $0.0125$
$\lambda_{t_{s}}$ ($\mathrm{[}\mathrm{n}\mathrm{T}\mathrm{/}\mathrm{yr}\mathrm{{}^{2}}\mathrm{]}^{-2}$) |  ${10}^{-2}$ |  ${10}^{-2}$ |  $1.25\text{\times}{10}^{-4}$
$\lambda_{t_{e}}$ ($\mathrm{[}\mathrm{n}\mathrm{T}\mathrm{/}\mathrm{yr}\mathrm{{}^{2}}\mathrm{]}^{-2}$) |  ${10}^{-2}$ |  ${10}^{-2}$ |  $1.25\text{\times}{10}^{-4}$
$\lambda_{\mathrm{mag}}$ ($\mathrm{[}\mathrm{n}\mathrm{T}\mathrm{/}\mathrm{yr}\mathrm{]}^{-2}$) |  $5\text{\times}{10}^{3}$ |  $5\text{\times}{10}^{3}$ |  $5\text{\times}{10}^{3}$
$\lambda_{\mathrm{pol}}$ ($\mathrm{[}\mathrm{m}\mathrm{A}\mathrm{/}\mathrm{m}\mathrm{]}^{-2}$) |  n/a |  ${10}^{6}$ |  ${10}^{6}$
$\lambda_{\mathrm{tor}}$ ($\mathrm{n}\mathrm{T}^{-2}$) |  n/a |  ${10}^{5}$ |  ${10}^{5}$

The iterative minimisation for the model estimation was initialised with a
starting model $\mathbf{m}_{0}$, which we chose to be CHAOS-7.9 for the
internal and magnetospheric model parts and, if included, zero-valued
parameters for the ionospheric model part. Concerning the Euler angles, we
used the initial values that have been determined during the pre-flight
calibration on ground for CHAMP [59] and Swarm [63], and during in-flight
calibration for CryoSat2-1 [42]. Convergence was typically reached after 15
iterations.

##  4 Results

In the following, we report on the achieved misfit of Reference and Model-A,
present the estimated ionospheric field of Model-A, and compare the estimated
internal fields of Reference and Model-A. Model-B is presented in the second
half of this section, where we study the effect of relaxing the temporal
regularisation of the internal field model when an ionospheric field is co-
estimated.

###  4.1 Fit to the magnetic data

To illustrate how well the magnetic field estimates of Model-A fit the
magnetic data, we computed vector and scalar residuals, i.e. the component-
wise differences $\Delta B_{r}$, $\Delta B_{\theta}$ and $\Delta B_{\phi}$,
and the difference in scalar magnitude, $\Delta F$, between the vector
estimates of the magnetic field from Model-A and the magnetic observations in
the dataset used for the model estimation. Note that although the scalar
component was not used in constructing the models, it is a useful diagnostic
and included here. Fig. 3 presents histograms of the vector and scalar
residuals for each satellite in bins of $0.5\text{\,}\mathrm{nT}$ width.

Figure 3: Histograms of vector and scalar residuals for each satellite dataset
with respect to Model-A. The bin width is $0.5\text{\,}\mathrm{nT}$ and the
histograms are normalised to integrate to unit area. The residuals outside the
range of $\pm 30\text{\,}\mathrm{nT}$ are not shown but taken into account for
the normalisation.

Irrespective of the residual component and the satellite, the histograms are
fairly symmetric and have a single maximum close to zero. The peaks of the
histograms for the three Swarm satellites are narrow and very similar in
appearance to the extent that they practically overlap. The peaks for CHAMP
are slightly broader with correspondingly lower maxima, especially regarding
the radial and southward components. The histograms for CryoSat2-1 are even
broader with maxima that are at approximately half the value of the Swarm and
CHAMP satellites, reflecting the generally higher noise level of these data.
Model-A clearly fits the radial and scalar components better than the $\theta$
and $\phi$ components, which are more influenced by rapidly varying field-
aligned currents, which our parametrization does not capture.

Tab. 3 presents Huber-weighted mean and Root-Mean-Square (RMS) values for each
satellite and field component, distinguishing between polar
($\lvert\lambda_{\mathrm{QD}}\rvert>$55\text{\,}\mathrm{\SIUnitSymbolDegree}$$)
and non-polar
($\lvert\lambda_{\mathrm{QD}}\rvert\leq$55\text{\,}\mathrm{\SIUnitSymbolDegree}$$)
latitudes.

Table 3: Statistics of vector and scalar residuals with respect to Model-A and Reference for each satellite dataset. In the table $N$ is the number of data, and Mean and RMS refer to the Huber-weighted mean and the Huber-weighted root-mean-square of the deviation from the mean, respectively. Polar QD latitude refers to $\lvert\lambda_{\mathrm{QD}}\rvert>$55\text{\,}\mathrm{\SIUnitSymbolDegree}$$ and non-polar QD latitude to the opposite. |  |  | $N$ | Model-A | Reference
---|---|---|---|---|---
|  |  |  | Mean ($\mathrm{nT}$) | RMS ($\mathrm{nT}$) | Mean ($\mathrm{nT}$) | RMS ($\mathrm{nT}$)
Dataset | QD latitude | Residual |  |  |  |  |
CHAMP | Non-polar | $\Delta B_{r}$ | 661357 | 0.43 | 4.18 | 0.38 | 5.27
$\Delta B_{\theta}$ | 661357 | -0.13 | 4.62 | -0.05 | 4.94
$\Delta B_{\phi}$ | 661357 | 0.05 | 5.92 | 0.07 | 6.48
$\Delta F$ | 661357 | -0.03 | 2.97 | -0.23 | 3.96
Polar | $\Delta B_{r}$ | 410810 | 0.14 | 6.60 | 0.25 | 8.70
$\Delta B_{\theta}$ | 410810 | -0.01 | 17.35 | 0.15 | 21.07
$\Delta B_{\phi}$ | 410810 | 0.03 | 18.89 | 0.31 | 23.14
$\Delta F$ | 410810 | -0.11 | 6.20 | -0.75 | 8.52
CryoSat2-1 | Non-polar | $\Delta B_{r}$ | 285278 | 0.02 | 4.97 | 0.03 | 5.46
$\Delta B_{\theta}$ | 285278 | -0.40 | 6.16 | -0.43 | 6.23
$\Delta B_{\phi}$ | 285278 | -0.02 | 6.32 | 0.11 | 6.66
$\Delta F$ | 285278 | 0.39 | 4.82 | 0.35 | 4.80
Polar | $\Delta B_{r}$ | 178099 | -0.61 | 7.02 | -0.57 | 7.67
$\Delta B_{\theta}$ | 178099 | 1.01 | 18.18 | 0.94 | 20.75
$\Delta B_{\phi}$ | 178099 | -0.15 | 19.45 | 0.08 | 23.52
$\Delta F$ | 178099 | 0.37 | 6.67 | 0.14 | 7.36
Swarm-A | Non-polar | $\Delta B_{r}$ | 190878 | 0.09 | 2.78 | 0.04 | 3.96
$\Delta B_{\theta}$ | 190878 | -0.05 | 3.53 | -0.06 | 3.89
$\Delta B_{\phi}$ | 190878 | 0.03 | 5.32 | 0.06 | 5.84
$\Delta F$ | 190878 | -0.09 | 2.62 | -0.26 | 3.56
Polar | $\Delta B_{r}$ | 118640 | -0.04 | 5.58 | -0.08 | 7.42
$\Delta B_{\theta}$ | 118640 | 0.11 | 15.89 | 0.11 | 19.65
$\Delta B_{\phi}$ | 118640 | 0.05 | 18.83 | 0.39 | 23.28
$\Delta F$ | 118640 | -0.22 | 5.18 | -0.75 | 7.27
Swarm-B | Non-polar | $\Delta B_{r}$ | 192185 | -0.02 | 2.71 | -0.08 | 3.89
$\Delta B_{\theta}$ | 192185 | -0.12 | 3.46 | -0.10 | 3.79
$\Delta B_{\phi}$ | 192185 | -0.01 | 5.27 | 0.04 | 5.78
$\Delta F$ | 192185 | -0.00 | 2.52 | -0.21 | 3.44
Polar | $\Delta B_{r}$ | 119224 | -0.12 | 5.18 | -0.13 | 6.83
$\Delta B_{\theta}$ | 119224 | 0.29 | 15.85 | 0.23 | 19.54
$\Delta B_{\phi}$ | 119224 | 0.11 | 18.61 | 0.44 | 22.99
$\Delta F$ | 119224 | 0.05 | 4.70 | -0.48 | 6.62
Swarm-C | Non-polar | $\Delta B_{r}$ | 194920 | 0.06 | 2.79 | -0.00 | 3.93
$\Delta B_{\theta}$ | 194920 | -0.11 | 3.54 | -0.13 | 3.90
$\Delta B_{\phi}$ | 194920 | 0.03 | 5.29 | 0.07 | 5.82
$\Delta F$ | 194920 | -0.01 | 2.61 | -0.17 | 3.54
Polar | $\Delta B_{r}$ | 121355 | 0.01 | 5.54 | -0.02 | 7.38
$\Delta B_{\theta}$ | 121355 | 0.10 | 15.82 | 0.09 | 19.62
$\Delta B_{\phi}$ | 121355 | 0.05 | 18.77 | 0.41 | 23.19
$\Delta F$ | 121355 | -0.05 | 5.15 | -0.58 | 7.24

The RMS values of Model-A are generally lower than those of Reference,
especially in the case of horizontal and scalar residuals at polar latitudes.
For example, the RMS values of $\Delta B_{\theta}$, $\Delta B_{\phi}$, and
$\Delta F$ for CHAMP in the polar regions are $17.35\text{\,}\mathrm{nT}$,
$18.89\text{\,}\mathrm{nT}$, and $6.20\text{\,}\mathrm{nT}$ for Model-A and
$21.07\text{\,}\mathrm{nT}$, $23.14\text{\,}\mathrm{nT}$, and
$8.52\text{\,}\mathrm{nT}$ for the reference model, respectively, which
corresponds to a reduction of approximately $20\text{\,}\mathrm{\char
37\relax}$. This suggests that co-estimating an ionospheric field improves the
fit to the data, in particular, at high latitudes. To further characterise the
improvement in the data fit, we investigated median scalar residuals as a
function of QD latitude and MLT, which we computed in cells of approximately
equal area using a HEALPix [19] grid in QD/MLT coordinates over the entire
globe. In Fig. 4 we compare these maps of median scalar residuals for Model-A
and the reference model.

Figure 4: Median scalar residuals with respect to the reference model (top)
and Model-A (bottom) in QD/MLT coordinates using an equal-area pixelation. The
global maps on the left are Mollweide projections, where the central vertical
line corresponds to midnight
($\phi_{\mathrm{MLT}}=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$) and the
central horizontal line to the magnetic dip equator
($\lambda_{\mathrm{QD}}=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$). The maps
on the right are orthographic projections of the northern magnetic hemisphere
(North) and, as if looking down on the Earth from the north pole, the southern
magnetic hemisphere (South). The labels indicate noon (12), midnight (00),
dawn (06) and dusk (18). The dashed lines show parallels and meridians at
$30\text{\,}\mathrm{\SIUnitSymbolDegree}$ intervals.

Looking at the polar views for the reference model, median scalar residuals
are organised into a pattern of two crescent-shaped cells of positive and
negative values around each magnetic pole, which reflects the well-known two-
cell current system in the polar ionosphere [[, e.g.]]Dungey1961. Likewise, in
the global view for the reference model, the strongly negative values of the
median scalar residuals on the dayside at the dip-equator and slightly less
negative values at mid-latitudes are associated with the solar quiet current
system and the equatorial electrojet [[, e.g.]]Yamazaki2016. In Model-A, the
median value of the scalar residuals in QD/MLT coordinates are dramatically
reduced not only at polar latitudes, where the AMPS approach is expected to
work best, but also at low and mid-latitudes, in particular, close to the dip
equator on the dayside. The remaining patterns are relatively weak and could
possibly be captured by using a higher truncation degree of the ionospheric
field model. The fact that the patterns found for the reference model are
largely absent for Model-A, in the bottom panel of Fig. 4, shows that our
approach accounts for previously unmodelled signals associated with
ionospheric current systems in the residuals.

Since the models in this study are derived only from satellite magnetic data,
we can test how well time variations are modelled in comparison to independent
observations made by ground-based magnetic observatories of the International
Real‐time Magnetic Observatory Network (INTERMAGNET). These data are available
as hourly mean values at the World Data Centres for Geomagnetism in Edinburgh
and have been quality checked as explained in [33]. In Fig. 5, for example, we
compare monthly means of the ionospheric field given by Model-A compared with
those recorded at the ground-based magnetic observatory in Hornsund (HRN),
located at $77.0\text{\,}\mathrm{\SIUnitSymbolDegree}$N,
$15.5\text{\,}\mathrm{\SIUnitSymbolDegree}$E on Svalbard (Norway) near the
northern edge of the auroral oval.

Figure 5: Monthly mean values of the quiet-time ionospheric field extracted
from the records of the magnetic observatory in Hornsund, Svalbard in Norway
(red) and those given by Model-A at the same location (black). Note that we
ignored the effect of induction during the downward continuation of our
ionospheric field model.

For this comparison, we extracted monthly means of the ionospheric field from
the timeseries of hourly means at HRN by applying the quiet-time selection in
Sect. 2, removing the internal and magnetospheric field estimates given by
Model-A, and centring the corrected timeseries with the component-wise average
to remove the remaining crustal field. Concerning the modelled ionospheric
field, we produced monthly averages using the hourly estimates of the poloidal
ionospheric field of Model-A at the times of the quiet-time HRN timeseries
after downward continuing this part of the model below the reference height to
the Earth’s surface. For this, we assumed that the poloidal ionospheric field
below the reference height can be represented through an external potential
whose radial field component is continuous across the reference height. Note
that the toroidal ionospheric part of the model does not contribute to the
ionospheric field on ground since it does not exist inside the non-conducting
part of the atmosphere. Fig. 5 shows that Model-A closely follows yearly and
slower variations of the HRN timeseries, especially the fit to the southward
component is encouraging. However, Model-A certainly underestimates the peak
values seen in the observatory timeseries, for example, for those in 2003, and
cannot reproduce the more dynamic time periods between 2001 and 2005, and
between 2012 and 2016, which are associated with solar maximum conditions.
This shows that our approach sensibly models the typical variations of the
ionospheric field but is not able to reproduce the dynamic ionospheric field
produced by local currents above the observatory.

To document our estimated ionospheric field at mid and low latitudes, we show
in Fig. 6 the radial component of the ionospheric magnetic field from Model-A
and from CM6 [56], the sum of primary and secondary parts, at satellite
altitude at $450\text{\,}\mathrm{km}$ during noon in Greenwich on March 21,
2018.

Figure 6: Radial component of the ionospheric magnetic field from Model-A
(left) and CM6 (right) at $450\text{\,}\mathrm{k}\mathrm{m}$ altitude during
noon in Greenwich on March 21, 2018. Note that we show the sum of inducing and
induced parts of the ionospheric field from CM6.

The overall pattern of the ionospheric field from Model-A at low and mid
latitudes is similar to the CM6 model. The radial field is stronger around
local noon north and south of the magnetic dip equator and the geometry
broadly follows that of the internal field. However, there are important
differences between Model-A and CM6. At mid and low latitudes the amplitude of
the radial field from Model-A is much weaker than for CM6 on the dayside,
while it is comparable on the nightside but with different signs in the
northern and southern hemispheres. At high latitudes the two cell pattern
around the magnetic poles is clearly visible for Model-A, while the pattern is
weaker and smeared out for CM6, in particular in the southern hemisphere. This
shows that our ionospheric field model captures the basic Solar-quiet pattern
but has limitations on the dayside at non-polar latitudes. Note that the
ionospheric field from Model-A becomes more similar to CM6 at low and mid
latitudes by relaxing the regularisation imposed on the divergence-free sheet
currents associated with this part of the model.

###  4.2 Ionospheric currents during geomagnetic quiet-time conditions

In this section we report the spatial structure of the estimated currents from
Model-A and their response to changes in the external driving.

In Fig. 7 we show polar views of the divergence-free and field-aligned current
densities in QD/MLT coordinates as a function of clock angle for the quiet
solar wind conditions represented in the dataset during winter in the northern
hemisphere.

Figure 7: Divergence-free and field-aligned current densities in the north
polar region for different clock angles. Each panel shows the current
densities in QD/MLT coordinates above
$60\text{\,}\mathrm{\SIUnitSymbolDegree}$ QD latitude with noon at the top,
dawn on the right, midnight at the bottom and dusk on the left. The contours
show the potential of the divergence-free sheet current density in steps of
$5\text{\,}\mathrm{kA}$ (solid for positive and dashed for negative values)
and the colours indicate the field-aligned current density. The location of
the largest ($\bigtriangleup$) and smallest ($\bigtriangledown$) field-aligned
current densities are indicated with the coloured triangles, and their
strength is given in the lower right corner of each panel. The total field-
aligned current ($\parallel$) and the divergence-free part of the sheet
current ($\perp$) flowing between the maximum and minimum of
$\psi^{\mathrm{df}}$, poleward of $60\text{\,}\mathrm{\SIUnitSymbolDegree}$ QD
latitude, are given in the lower left corner. The dotted lines indicate QD
parallels at $10\text{\,}\mathrm{\SIUnitSymbolDegree}$ and MLT meridians at
$2\text{\,}\mathrm{h}$ intervals. This figure is similar in form to those
originally presented by [28] for the AMPS model.

For all clock angles, the divergence-free sheet currents circulate in two
cells, roughly separated by the noon-midnight meridian, whereas the field-
aligned currents form concentric patterns with the centre slightly offset from
the magnetic pole towards midnight, known as R1 and R2 currents [25, 24]. For
northward IMF ($\theta_{\mathrm{c}}=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$)
the cell of the divergence-free sheet currents in the dawn sector is slightly
more pronounced than the one in the dusk sector, and the maxima of the field-
aligned currents are located close to noon. When the IMF rotates southward,
the currents generally gain in strength. However, a clockwise rotation to
$\theta_{\mathrm{c}}=$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$ leads to
stronger currents than a counter-clockwise rotation to
$\theta_{\mathrm{c}}=$-90\text{\,}\mathrm{\SIUnitSymbolDegree}$$, which
corresponds to an asymmetry in the currents that depends on the y-component of
the IMF. The currents are strongest when the IMF is southward
($\theta_{\mathrm{c}}=$180\text{\,}\mathrm{\SIUnitSymbolDegree}$$) with a
maximum of downward field-aligned currents in the noon-dawn sector and a
maximum of upward field-aligned currents in the midnight-dusk sector. The case
of southward oriented IMF, however, should be interpreted with caution since
it is poorly represented in the data due to the chosen data selection (see
panel for the clock angle in the top centre of Fig. 2).

Fig. 8 shows the divergence-free and field-aligned current densities for the
same conditions as in Fig. 7 but in dependence of the $\mathit{SML}$ index and
only for purely northward IMF.

Figure 8: Divergence-free and field-aligned current densities in the north
polar region for different values of the $\mathit{SML}$ index and purely
northward IMF. The figure is otherwise identical to Fig. 7.

With decreasing $\mathit{SML}$ index, corresponding to an increase in auroral
activity, the divergence-free sheet currents and the field-aligned currents
grow in strength. The locations where the field-aligned currents are strongest
move from near noon ($\mathit{SML}=$-40\text{\,}\mathrm{nT}$$) to the midnight
sector ($\mathit{SML}=$-160\text{\,}\mathrm{nT}$$), reaching a configuration
in which an upward directed current dominates the pre-midnight and a downward
current the post-midnight sectors. This pair of strong downward and upward
field-aligned currents centred on midnight agrees well with the current wedge
that is thought to exist during substorms [37, 26]. Again, however, the case
for $\mathit{SML}=$-160\text{\,}\mathrm{nT}$$ should be considered an
extrapolation given that there was a relatively small amount of data with
large negative values of the $\mathit{SML}$ index available during the
modelling thanks to the quiet-time data selection (for the distribution of the
$\mathit{SML}$ see lower right panel in Fig. 2).

For interpreting the variations of the patterns in Figs. 7 and 8 in terms of
physical processes, it is worth mentioning that the terms in the ionospheric
magnetic field parametrization related to the $\mathit{SML}$ index and the
$\epsilon$ coupling function probably compete to some extent for the same
signal. So a part of the structure which may go into the $\epsilon$ terms, if
the $\mathit{SML}$ index was not included in the parametrization, may now be
contained in the $\mathit{SML}$ terms. For example in Fig. 7, the changes due
to variations in the clock angle may be less than in the AMPS model, which did
not include the $\mathit{SML}$ index, because some of the signal for southward
IMF conditions is now contained in the $\mathit{SML}$ terms.

Finally, to illustrate the ability of the climatological approach to allow
solar cycle changes in the ionospheric currents, we produced a sequence of
plots, similar to Figs. 7 and 8, that show the current densities in the
northern polar region averaged in time over successive 2.5 year intervals
between 2008.0 and 2020.5. Specifically, we averaged the currents over the
intervals 2008.0–2010.5 (solar minimum), 2010.5–2013.0 (ascending solar
cycle), 2013.0–2015.5 (solar maximum), 2015.5–2018.0 (descending solar cycle),
and 2018.0–2020.5 (again solar minimum) to illustrate different phases of the
solar cycle. To further emphasise the changes in the currents, we removed the
average current density for the entire solar cycle shown here. Fig. 9 shows
this sequence of plots and the removed average current density along with the
$F_{10.7}$ and $\mathit{SML}$ indices, averaged using the same 2.5 year
intervals.

Figure 9: Average $F_{10.7}$ and $\mathit{SML}$ indices (top row), and
divergence-free and field-aligned current densities (middle row) in the north
polar region successively averaged over 2.5 year periods after removal of the
average current density (bottom row) over the solar cycle from 2008 to 2020.5.
The markers on the curves of the indices are placed at the midpoints of the
2.5 year intervals used for the averaging, and the contours of the potential
of the divergence-free sheet current density are shown in steps of
$5\text{\,}\mathrm{kA}$.

The $F_{10.7}$ index reaches a maximum around 2014, indicating solar maximum,
whereas the $\mathit{SML}$ index has a minimum later, close to 2017,
coinciding with the descending phase of the solar cycle. Turning to the field-
aligned currents, differences with respect to the solar cycle average mostly
occur in the noon sector with the exception of the descending phase, when the
differences are most prominent in the midnight sector. The differences in the
divergence-free sheet currents with respect to the solar cycle average are
generally more complex. However, noteworthy is the two-cell pattern that is
visible during the descending phase of the solar cycle, which is consistent
with the enhancement of this pattern for decreasing values of the
$\mathit{SML}$ index as shown in Fig. 8. Overall, we conclude that our
approach is able to capture at least part of the changes that occur in the
strength and appearance of polar ionospheric currents during the solar cycle.

###  4.3 Core field secular variation at polar latitudes

We demonstrated above that the co-estimation of the polar ionospheric currents
helps with accounting for previously unmodelled signals in the residuals. We
now turn to the impact of co-estimating the ionospheric field on the estimated
core field, by considering differences between the internal field estimates of
Model-A and the reference model.

In Fig. 10 we show the spatial power spectra of the Secular Variation (SV) and
Secular Acceleration (SA) at the CMB in 2019.0 from Model-A, Reference,
CHAOS-7.9 , and the difference between Model-A and Reference.

Figure 10: Spatial power spectrum of the SV (left) and SA (right) at the CMB
in 2019.0 from Model-A (blue), Reference (orange), CHAOS-7.9 (green), and the
difference between Model-A and Reference (black dashed).

The SV spectra are very similar at low spherical harmonic degree, causing the
curves for the models to overlap in the plot. Above degree 13 the spectra
deviate more clearly but continue to stay closely together as the degree
increases. The SA spectra increase with spherical harmonic degree and reach a
maximum at degree 10. Above degree 10 the spectra decrease but much steeper
for Model-A and Reference than for CHAOS-7.9. This difference in the high-
degree SA is the result of the temporal regularisation used in CHAOS-7.9,
which is tapered to allow for more power at high degrees. However, the taper
was not applied in the models of this study. Since the spectra from Model-A
and Reference are very similar, we show in Fig. 11 the sensitivity matrix of
the SV and the SA in 2019.0 from Model-A with respect to Reference.

Figure 11: Sensitivity matrix of the SV (left) and the SA (right) in 2019.0
from Model-A with respect to Reference.

The sensitivity matrix is defined as the coefficient-wise difference between
recovered spherical harmonic coefficients and chosen target coefficients, here
the coefficients from Reference, normalised with the mean amplitude of the
target coefficients at degree $n$ [[, e.g.]]Sabaka2013. The SV coefficients of
Model-A and Reference are very similar at low spherical harmonic degree.
However, the sensitivity increases above degree 14 in particular for near-
zonal ($m\approx 0$) and near-sectorial ($m\approx n$) coefficients. A similar
pattern can be observed in the sensitivity matrix of the SA, but the numerical
values are overall larger. Noteworthy are relatively strong sensitivity values
for the zonal SA coefficients at degree 2 and 3.

The spatial power spectra and the sensitivity matrix for the SV show that
appreciable differences between Model-A and the reference model can be
expected at high spherical harmonic degrees. To examine this we plot in Fig.
12 snapshots of the radial SV for Model-A and Reference up to degree 19 at the
CMB in the north polar region in 2007.0, 2012.0 and 2018.0, after removing the
snapshot average to emphasise changes in time.

Figure 12: Snapshots of the radial SV for $n\leq 19$ in the north polar region
at the CMB as given by Model-A (left column), the reference model (middle
column) and their difference (right column) in 2007.0 (second row), 2012.0
(third row), 2018.0 (fourth row), after removing the snapshot average (first
row). The projection is orthographic. The dashed lines show geographic
parallels and meridians at $30\text{\,}\mathrm{\SIUnitSymbolDegree}$
intervals. Note the change in colour scale for the difference plots.

The average radial SV for Model-A shows small-scale flux patches of positive
and negative polarity, with a distinct pair of patches located over Siberia.
This pair is also visible in the average map for the reference model but is
less prominent due to relatively strong patches west of Greenland, which are
slightly elongated in longitude. For Model-A, the flux patches around that
area are less extended in longitude, showing that the co-estimation of the
ionospheric field leads to better focused patches of SV at the CMB. This can
also be seen in the difference between the average radial SV from Model-A and
Reference, which exhibits near-zonal features of radial SV around the north
pole. Despite this improvement in Model-A, the presence of a weak pattern of
stripes in the SV along geographic parallels over North America may indicate
that the SV at high latitudes is still contaminated by the ionospheric field
at high degree.

The snapshots showing the deviation of the radial SV from the average reveal
that the three non-axisymmetric SV flux patches over Siberia and Alaska
intensify for Model-A over the model time interval, and similarly for
Reference, by an almost linear trend. The lack of structure in the difference
between Model-A and Reference for each of these snapshots shows that the
difference between the radial SV from Model-A and from Reference is mostly
static.

Overall, it seems that the co-estimation of the ionospheric magnetic field
leads to less contaminated models of the SV in the north polar region that are
derived from magnetic vector data at all latitudes and local times. It is
clear that the distinctive non-axisymmetric radial SV flux patches of
alternating sign over Siberia and Alaska, found in earlier geomagnetic field
models and used for inferring accelerating jets of core flow [32], persist
even when polar ionospheric currents are estimated. A clear limitation in
examining the high degree CMB SV in these models is that they are strongly
smoothed in time. The effect of reducing the strong temporal regularisation
will be investigated in the next section.

###  4.4 Relaxing the temporal regularisation of the internal field model

Our approach for modelling the geomagnetic field relies on regularisation to
smooth the spatio-temporal complexity of the model and, thus, ensures the
convergence of the estimation procedure. The regularisation also allows
control over magnetic signals in the data that are not accurately
parametrized. However, in the case of the time-dependent internal field model,
the temporal regularisation severely degrades the resolution in time, in
particular, of the high-degree spherical harmonics, which limits studies of
the dynamics of the Earth’s interior. In this section, we therefore explore a
model in which, building on Model-A, we in addition considerably relaxed the
temporal regularisation of the internal field model. More specifically, we
reduced the temporal regularisation parameters $\lambda_{t}$,
$\lambda_{t_{s}}$ and $\lambda_{t_{e}}$ by a factor of $80$ to $0.0125$,
$1.25\text{\times}{10}^{-4}$ and $1.25\text{\times}{10}^{-4}$, respectively,
to produce a weakly-regularised version of Model-A, called Model-B.

On the left of Fig. 13, we show the spatial power spectrum of the SA at the
CMB in 2019.0 for Model-B, Model-A and CHAOS-7.9.

Figure 13: Spatial power spectrum of the SA at the CMB in 2019.0 (left) and
timeseries of the SV coefficients $\dot{g}_{1}^{0}$ (top right) and
$\dot{h}_{12}^{1}$ (bottom right) for Model-B (blue), Model-A (orange) and
CHAOS-7.9 (green).

The spectra show that there is significantly more power in the SA of Model-B
at all spherical harmonic degrees. Most striking is the increase in SA power
for Model-B above degree 10, while the power for Model-A sharply decreases.
But also noteworthy is the enhanced power below degree 4 for Model-B. For
comparison we also show the spectrum for CHAOS-7.9, which overlaps below
degree 10 with that for Model-A but decreases more slowly as the degree
further increases. On the right of Fig. 13, we show example timeseries of SV
coefficients for the three models. The timeseries of $\dot{g}_{1}^{0}$ reveals
a distinct annual oscillation for Model-B. Similarly, we found annual
oscillations in the timeseries of zonal and low-order coefficients for other
low-degree spherical harmonics. They are investigated further below. We
emphasise that these annual oscillations are not an effect of the applied data
selection, which does not vary with season. The oscillations are not as
apparent in the coefficient timeseries of Model-A and CHAOS-7.9, due to the
relatively strong temporal regularisation. Compared to $\dot{g}_{1}^{0}$, the
$\dot{h}_{12}^{1}$ timeseries is much smoother for all three models and there
is no oscillatory behaviour with a period close to one year visible. So
Model-B is still quite strongly smoothed in time at high degree. This is by
construction due to the chosen regularisation norm. Note that we chose to show
the SA spectrum in 2019.0 in Fig. 13 because the rate of change of the SV
maximises around that time for Model-B. In fact, we find that the spectrum for
Model-B varies significantly with time below degree 8, reflecting the annual
variations found in the low-order and low-degree SV coefficients. These
temporal variations are well known and have been handled by other modelling
efforts in different ways, by using regularisation [e.g., CHAOS models [13],
GRIMM [31], or CM [56]], low resolution basis functions [e.g., POMME [36],
CovObs [23]], or by including other internal sources [e.g., Kalmag models [3],
or sequential models of [49]].

To better characterise the annual oscillations found in the low spherical
harmonic coefficients of the time-dependent internal field model of Model-B,
we performed a principal component analysis (PCA) of the difference between
Model-B and Model-A. The PCA is a data-based tool for extracting spatio-
temporal patterns and has previously been applied to the magnetic data from
ground-based observatories [[, e.g.]]Shore2016 and magnetic observations made
by satellites [[, e.g.]]Domingos2019,Saturnino2021. For the analysis, we
generated timeseries of vector components of the internal time-dependent field
from each model at the centre points of equal-area pixels [19], covering the
entire Earth’s surface at approximately
$2\text{\,}\mathrm{\SIUnitSymbolDegree}$ resolution (10800 pixels), in
spherical geocentric coordinates using a sampling rate of one sample per month
(201 samples in time). Here, we omitted the period during the CryoSat-2 data,
from the middle of 2010 to the end of 2013, when Model-B varies more strongly
in time than during CHAMP and Swarm, and we omitted the first and last 6
months of the model time interval, when Model-B shows sharp time variations
due to the weaker data constraint at the model endpoints. At each pixel, by
subtracting the vector timeseries of Model-A from the one of Model-B, we
obtained timeseries of component-wise differences, which we arranged as
columns in a matrix $\mathbf{X}$ of size $201\times 32400$. Finally, we
centred $\mathbf{X}$ by removing the column-wise mean value,
$\bar{\mathbf{X}}$, such that

$\tilde{\mathbf{X}}=\mathbf{X}-\bar{\mathbf{X}}.$ (28)

For a data matrix such as $\tilde{\mathbf{X}}$, PCA can be used to find a
finite set of modes that maximise the variance of the data in time and are
mutually orthogonal. The $i$th mode obtained through the PCA consists of the
Empirical Orthogonal Function (EOF), $\mathbf{v}_{i}$, which represents the
spatial pattern of the time variation, and the principal component (PC)
$\mathbf{y}_{i}=\tilde{\mathbf{X}}\mathbf{v}_{i}$, which is the timeseries of
variance $\sigma^{2}_{i}$. These modes are typically sorted in decreasing
order of variance and can be used to reconstruct the data matrix through

$\tilde{\mathbf{X}}=\sum_{i}\mathbf{y}_{i}\mathbf{v}^{\mathrm{T}}_{i},$ (29)

or, if a subset of modes is chosen, to perform a partial reconstruction.

We applied the PCA on $\tilde{\mathbf{X}}$ and obtained 201 modes. But only a
small number of modes are needed to explain most of the variance in the data.
In Fig. 14, we show the PCs and the radial part of the EOFs for the first six
modes, which account for $70\text{\,}\mathrm{\char 37\relax}$ of the variance
(13 modes account for $90\text{\,}\mathrm{\char 37\relax}$).

Figure 14: PCs (left) and radial part of the corresponding EOFs (right) for
the first six modes found by the PCA. The EOFs are scaled with the square-root
of the mode variance to indicate the relative importance.

The first PC, PC-1, is a modulated annual oscillation for which the amplitude
maximises around 2002 and 2014, and minimises around 2008 and 2019. The
corresponding EOF, EOF-1, consists of a large-scale pattern that varies mostly
in latitude, similar in appearance to the spherical harmonic function
$Y_{2}^{0}$ except that the zero lines in latitude are shifted slightly
northward. Given the spatio-temporal behaviour, we can assume that the first
mode is responsible for most of the annual oscillations in the coefficient
timeseries, as shown for example by $\dot{g}_{1}^{0}$ in Fig. 13 (phase shift
between $\dot{g}_{1}^{0}$ and PC-1 reflects the difference in time derivative
between the two). The amplitude modulation of PC-1 suggests a dependency on
the solar cycle since the maxima in the amplitude roughly coincide with the
times of solar maximum. PC-2 varies more slowly compared to PC-1 and peaks
approximately every 3 years. EOF-2 shows patches of opposite sign that are
centred on the geographic equator around Central America and the western
Pacific Ocean. The location and appearance of these patches could indicate
that the second mode is related to changes in the core field since geomagnetic
impulses have been reported around these areas [43, 5, 64, 13]. The third PC,
PC-3, combines a slow variation that peaks around 2002 and 2014 with an annual
oscillation, which is much weaker in amplitude compared to PC-1 but likely
also of external origin. EOF-3 exhibits a large number of small-scale positive
and negative patches. The remaining PCs (PC-4 through PC-6) show an
oscillatory behaviour with a period close to one year and the corresponding
EOFs exhibit a wide range of patterns, which are difficult to interpret.

Based on the results of the PCA, we tried to implement a post-processing step
that removes the modes that we assume are the result of a leakage of the
external field into the internal field model as is the case, for example, for
the first mode, since it is a global annual oscillation with an apparent solar
cycle dependency. Successful removal of these modes would provide a smoother
timeseries of the coefficients of the time-dependent internal field model,
which we could use to analyse the SA. Unfortunately, this approach did not
work in practice since the PCA is not able to isolate the annual oscillations
in terms of a single mode. Instead, we found that these oscillations are
visible in most, if not all, of the modes and cannot be removed to a
satisfactory level, including by further relaxing the temporal regularisation,
which makes the annual oscillations clearer. Nonetheless, we find the PCA
provides clear insight into the signals entering internal field models as the
temporal regularisation is relaxed.

As a simpler alternative to removing PCA modes, we resorted to filtering out
the annual oscillations by computing centred annual differences of monthly
values of the internal field at the CMB as given by Model-B to find the SV
and, by repetition, the SA. Fig. 15 shows time-longitude plots of the obtained
radial SA up to degree 10 for Model-B on the left and, by the same
computation, for Model-A and CHAOS-7.9 in the middle, and the difference
between Model-B and CHAOS-7.9 on the right.

Figure 15: Time-longitude plot of the radial SA up to degree 10 along the
geographic equator at the CMB as given by repeated computation of annual
differences of the radial field from Model-B, Model-A, CHAOS-7.9, and the
radial field difference between Model-B and CHAOS-7.9.

We see that, in comparison to Model-A and also CHAOS-7.9, which is similar to
Model-A in Fig. 15, there is more power in the SA of Model-B. The patterns
have sharper edges and there is generally more structure, which indicates an
improved temporal resolution. This is especially apparent during the CHAMP
period, until 2010, and in the longitude interval between
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$ for the entire model time span. We
produced a similar plot for a weakly regularised version of Reference and
found that it looks very similar to Fig. 15. We interpret this to mean that
the increase in resolution is mostly due to the reduced temporal
regularisation and not due to the co-estimation of the ionospheric field
model. We acknowledge that computing annual differences has the caveat of
removing genuine internal field signals that have frequencies which are
integer multiples of one oscillation per year. Further work is clearly needed
on better methods to remove the annual signal.

To summarise, we find that reducing the temporal regularisation of the
internal field model causes signals which we suspect are of external origin to
leak into the estimated internal field despite the co-estimation of an AMPS-
type model of the ionospheric field. It is possible that our approach to
parametrize the ionospheric field lacks terms that can account for the signals
like those seen in the first and third mode identified by the PCA (see
discussion in Sect. 5). To ensure that the co-estimation of the ionospheric
field model is not, in fact, introducing artefacts into the internal field
model, we performed the same analysis of principal components but using the
reference model. We found that the first two modes are similar, which suggests
that these modes originate from genuine signals in the magnetic dataset used
in this study.

##  5 Discussion

Our results show that the co-estimation of a climatological ionospheric field
as part of the CHAOS modelling approach takes into account previously
unmodelled signals in the polar regions and produces geomagnetic field models
that fit the magnetic input data for geomagnetically quiet conditions well. It
enables the construction of high quality models of the core field while using
vector field data at all latitudes and local time. Similar to the AMPS model,
our approach provides estimates of the average polar ionospheric currents that
are realistic in structure and able to vary in response to changes of the
external driving. However, this approach is only capable to represent the long
term average of the currents and not individual highly dynamic events.

One aim of this study is to answer the question of whether the co-estimation
of an AMPS-type ionospheric model could allow the construction of internal
field models that are less contaminated by the ionospheric field in the polar
regions. By comparing the recovered SV of Model-A and Reference at the CMB in
the northern polar region (see Fig. 12), we find that the co-estimation
reduces the leakage of ionospheric signals into the time-dependent internal
field model. However, the improvement in the recovered SV in the polar
regions, which is most apparent in the zonal terms of the high spherical
harmonic degrees, is relatively small and difficult to interpret due to
significant ionospheric signals that remain even when co-estimating the
ionospheric field model and due to the strong time averaging effect of the
temporal regularisation applied on the internal field model. We find that the
regularisation of the poloidal ionospheric field is important. This mostly
affects the zonal or near-zonal parts of the poloidal ionospheric model as
observed within an Earth-fixed frame, which are the terms that show the
largest difference between Model-A and Reference in the internal field model.
We acknowledge that the Reference model is a somewhat extreme case since it
uses vector data at high latitudes without accounting for polar ionospheric
signals in any way. To test whether the ionospheric leakage into the internal
magnetic field model can be further reduced in the polar regions, we derived
additional test models where the internal field was estimated using the scalar
component of the magnetic observations instead of vector data at polar
latitudes. However, since we found these test models to be very similar to the
models shown here, we preferred to present the models derived from vector data
at all latitudes.

In non-polar regions we find that the estimated ionospheric field model
captures the basic Solar-quiet pattern but lacks power on the dayside. Hence,
there is a possibility of leakage into the time-dependent internal field due
to the use of dayside data and an imperfectly modelled ionospheric field at
mid and low latitudes. Comparing internal field coefficients $g_{1}^{0}$ and
$g_{3}^{0}$, which are known to be the coefficients most affected by the
ionospheric field at mid and low altitude, we find that the timeseries of
these coefficients for Model-A, and similarly for Reference, are slightly
shifted with respect to CM6, by $25\text{\,}\mathrm{nT}$, while the shift is
smaller with respect to CHAOS-7.9. A possible future remedy could be to
estimate the time-dependent internal field only from nightside data, while the
ionospheric model is fit using data from all local times.

Whether the reduced ionospheric field leakage also leads to internal field
models that are better resolved in time, can only be investigated by reducing
the temporal regularisation that smooths the time-dependence of the internal
field model, in particular affecting the high spherical harmonic degrees. For
this reason we derived Model-B, which exhibits more power in the high degrees
of the SA. However, we find that the estimated SA is now dominated at the
large length-scales, approximately up to degree 8, by distinct annual
oscillations. Although these annual oscillations can be filtered out of the SA
estimates by using annual differences, the question of where these
oscillations come from and why they are not captured by the ionospheric model
remains.

A possible explanation for the leakage of annual oscillations into the
internal field model is related to the parametrization of the ionospheric
field, which is most likely not sufficient to account for all types of
ionospheric currents and their time-dependence. In addition, it is assumed
that field-aligned currents are radial, which is reasonable at high latitudes
but fails at mid latitudes. Hence, it is conceivable that the annual
oscillations result from seasonally varying currents at mid latitudes such as
interhemispheric field-aligned currents that produce magnetic signals at
satellite altitude at mid latitudes on the dayside. Accounting for these
currents is challenging since it requires the estimation of poloidal and
toroidal potentials of the ionospheric magnetic field within the measurement
shell traced by the satellite orbits [[, e.g.]]Fillion2023,Olsen1997. Apart
from ionospheric sources, other processes cannot be ruled out as the origin of
the annual oscillations.

Another limitation of our models concerns the treatment of electromagnetically
induced currents in Earth’s interior and oceans associated with variations of
the ionospheric field. Since our models are only estimated from satellite
data, both the inducing and induced ionospheric fields are internal with
respect to the input data. Hence, our estimated ionospheric field model also
contains the induced response. In principle, through an a-posteriori analysis,
the estimated ionospheric field model could be separated into induced and
inducing parts using the Q-response functions [[, e.g.]]Grayver2021 for a
given model of the Earth’s conductivity. However, this approach would not
affect the quality of the core field model or resolve the ambiguity between
the sources in the ionosphere and those in the core and lithosphere.

Instead of a-posteriori separating the induced and inducing parts of our
estimated ionospheric field model, an alternative approach is to include the
induced response during the modelling. For example, one could derive a new set
of AMPS-type ionospheric field basis functions that take into account the
induced counterpart via Q-response transfer functions. This would have the
advantage of allowing ground-based observations to be used during the
estimation of both the core and ionospheric fields, for example, using hourly
mean values. For these observations, the inducing ionospheric field sources
are then external, which would aid the separation of the core and ionospheric
fields.

##  6 Conclusions

In this study we successfully combined the climatological approach of the AMPS
model, which is suitable for modelling the ionospheric magnetic field in the
polar regions, and the CHAOS modelling framework to derive models of the
geomagnetic field that take into account internal and magnetospheric fields as
well as the climatological aspects of the ionospheric field. We used this new
approach to estimate a geomagnetic field model from satellite magnetic vector
data under geomagnetic quiet conditions.

The derived model, called Model-A, shows a good fit to the input vector data
and successfully removes obvious systematic errors related to ionospheric
signals in the polar regions, which were previously unaccounted for in the
CHAOS modelling framework. By investigating the effect of co-estimating the
ionospheric field on the internal field and its time variations in the polar
regions, we find only small differences, which are most visible in the zonal
terms of the high-degree spherical harmonics of the estimated SV. Importantly,
high latitude non-axisymmetric SV flux features stay mostly unchanged, which
adds to the evidence that they are of internal origin and therefore relevant
for studies of the core flow [32].

The distinct annual oscillations in the internal field from Model-B, which was
weakly regularised in time, could indicate that there remain ionospheric or
related induced signals in the modelled internal field at low-to-mid latitudes
despite the co-estimation an AMPS-type ionospheric field. This suggests
shortcomings of our ionospheric field parametrization in non-polar regions,
and it indicates that the noise present in time-dependent internal field
models, which mostly affects the low-order spherical harmonic coefficients, is
not only due to the leakage of magnetic field signals produced by high-
latitude currents. Identifying the physical origin of these signals and taking
them into account will be important to increase the resolution of internal
field models in time.

The ambiguity between the internal field and poloidal ionospheric field models
was reduced through a practical approach by regularising the time-averaged
divergence-free part of the ionospheric currents, although this involves
regularisation parameters that must be carefully chosen during model
construction. In the future not only satellite magnetic data but also ground-
based magnetic observations could be used in order to better resolve this
ambiguity. But this requires treatment of the internally induced field due to
the poloidal ionospheric field on the ground-based data.

###### Acknowledgements.

We gratefully acknowledge ESA for providing access to the Swarm L1b data and
the fully calibrated CryoSat2 magnetometer data. We wish to thank the German
Aerospace Center and the Federal Ministry of Education and Research for
supporting the CHAMP mission. Furthermore, we would like to thank the staff of
the geomagnetic observatories and INTERMAGNET for supplying high-quality
observatory data. Susan Macmillan (British Geological Survey) is gratefully
acknowledged for collating checked and corrected observatory hourly mean
values in the AUX OBS database. We also thank the editor and two anonymous
reviewers for helpful comments and suggestions, which clarified and improved
the manuscript. CK and CCF were funded by the European Research Council under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 772561). The study has been partly supported as part of Swarm
Data Innovation and Science Cluster (Swarm DISC) activities, funded by the ESA
contract No. 4000109587/13/I-NB. KL was funded by the Trond Mohn Foundation,
and the AMPS model is funded by ESA through Swarm DISC within the reference
frame of ESA contract No. 000109587/13/I-NB. CK developed the computer
software used here for modelling the geomagnetic field, derived and analysed
the test field models, and prepared the manuscript. CCF helped with the
presentation and interpretation of the results, and assisted with the outline
of the manuscript. KL helped with the interpretation of the results and
provided guidance on the implementation of the part of the developed computer
software that involves the modelling of the ionospheric magnetic field. NO
helped with the interpretation of the results and provided the computer
scripts on which parts of the developed modelling software are based. All
authors read and accepted the manuscript.

The data underlying this article are available in the following repositories:

  1. 1.

Swarm and CryoSat2 data are available through ESA at https://swarm-
diss.eo.esa.int/#.

  2. 2.

CHAMP data are available from [50].

  3. 3.

Ground magnetic observatory data from INTERMAGNET are available at
ftp://ftp.nerc-murchison.ac.uk/geomag/Swarm/AUX_OBS/hour/ or via the virtual
research platform VirES https://vires.services/.

  4. 4.

$\mathit{RC}$ index is available at http://www.spacecenter.dk/files/magnetic-
models/RC/.

  5. 5.

$\mathit{SML}$ index is available at https://supermag.jhuapl.edu/info/.

  6. 6.

Solar wind speed, interplanetary magnetic field and $\mathit{Kp}$ index are
available through the GSFC/SPDF OMNIWeb interface at
https://omniweb.gsfc.nasa.gov/ow.html.

  7. 7.

CHAOS-7 model and its updates are available at
http://www.spacecenter.dk/files/magnetic-models/CHAOS-7/.

Files containing the estimated parameters of our models and the processed data
are available at https://doi.org/10.11583/DTU.24025596.

## References

  * [1] P. Alken et al.  “International Geomagnetic Reference Field: the thirteenth generation”  In _Earth, Planets and Space_ 73.1 Springer ScienceBusiness Media LLC, 2021  DOI: 10.1186/s40623-020-01288-x
  * [2] Jonathan Aurnou, Steven Andreadis, Lixin Zhu and Peter Olson  “Experiments on convection in Earth’s core tangent cylinder”  In _Earth and Planetary Science Letters_ 212.1-2 Elsevier BV, 2003, pp. 119–134  DOI: 10.1016/s0012-821x(03)00237-1
  * [3] Julien Baerenzung, Matthias Holschneider, Jan Saynisch-Wagner and Maik Thomas  “Kalmag: a high spatio-temporal model of the geomagnetic field”  In _Earth, Planets and Space_ 74.1 Springer ScienceBusiness Media LLC, 2022  DOI: 10.1186/s40623-022-01692-5
  * [4] Julien Baerenzung, Matthias Holschneider, Johannes Wicht, Vincent Lesur and Sabrina Sanchez  “The Kalmag model as a candidate for IGRF-13”  In _Earth, Planets and Space_ 72.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1186/s40623-020-01295-y
  * [5] A. Chulliat and S. Maus  “Geomagnetic secular acceleration, jerks, and a localized standing wave at the core surface from 2000 to 2010”  In _Journal of Geophysical Research: Solid Earth_ 119.3 Wiley-Blackwell, 2014, pp. 1531–1543  DOI: 10.1002/2013jb010604
  * [6] Arnaud Chulliat, Pierre Vigneron, Erwan Thébault, Olivier Sirol and Gauthier Hulot  “Swarm SCARF dedicated ionospheric field inversion chain”  In _Earth, Planets and Space_ 65.11 Springer ScienceBusiness Media LLC, 2013, pp. 1271–1283  DOI: 10.5047/eps.2013.08.006
  * [7] C.. Constable  “Parameter estimation in non-Gaussian noise”  In _Geophysical Journal International_ 94.1 Wiley Online Library, 1988, pp. 131–142
  * [8] João Domingos, Maria Alexandra Pais, Dominique Jault and Mioara Mandea  “Temporal resolution of internal magnetic field modes from satellite data”  In _Earth, Planets and Space_ 71.1 Springer ScienceBusiness Media LLC, 2019  DOI: 10.1186/s40623-018-0983-5
  * [9] James W. Dungey  “Interplanetary magnetic field and the auroral zones”  In _Physical Review Letters_ 6.2 American Physical Society (APS), 1961, pp. 47–48  DOI: 10.1103/physrevlett.6.47
  * [10] J.. Emmert, A.. Richmond and D.. Drob  “A computationally compact representation of Magnetic-Apex and Quasi-Dipole coordinates with smooth base vectors”  In _Journal of Geophysical Research: Space Physics_ 115.A8, 2010  DOI: 10.1029/2010JA015326
  * [11] M. Fillion, G. Hulot, P. Alken and A. Chulliat  “Modeling the climatology of low- and mid-latitude F-region ionospheric currents using the Swarm constellation”  In _Journal of Geophysical Research: Space Physics_ 128.5 American Geophysical Union (AGU), 2023  DOI: 10.1029/2023ja031344
  * [12] C.. Finlay et al.  “Challenges handling magnetospheric and ionospheric signals in internal geomagnetic field modelling”  In _Space Science Reviews_ 206.1-4 Springer Nature, 2016, pp. 157–189  DOI: 10.1007/s11214-016-0285-9
  * [13] C.. Finlay et al.  “The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic Anomaly”  In _Earth, Planets and Space_ , 2020  DOI: 10.1186/s40623-020-01252-9
  * [14] E. Friis-Christensen, C.. Finlay, M. Hesse and K.. Laundal  “Magnetic field perturbations from currents in the dark polar regions during quiet geomagnetic conditions”  In _Space Science Reviews_ 206.1-4 Springer Nature, 2017, pp. 281–297  DOI: 10.1007/s11214-017-0332-1
  * [15] E. Friis-Christensen, H. Lühr and G. Hulot  “Swarm: a constellation to study the Earth’s magnetic field”  In _Earth, Planets and Space_ 58.4 Springer ScienceBusiness Media LLC, 2006, pp. 351–358  DOI: 10.1186/bf03351933
  * [16] GFZ Section 2.3  “CHAMP: overview of final ME products and format description”  In _Scientific Technical Report - Data; 19/10; ISSN 2190-7110_ GFZ German Research Centre for Geosciences, 2019  DOI: 10.2312/GFZ.B103-19104
  * [17] J.. Gjerloev  “A Global Ground-Based Magnetometer Initiative”  In _Eos, Transactions American Geophysical Union_ 90.27 American Geophysical Union (AGU), 2009, pp. 230  DOI: 10.1029/2009eo270002
  * [18] J.. Gjerloev  “The SuperMAG data processing technique”  In _Journal of Geophysical Research: Space Physics_ 117.A9 American Geophysical Union (AGU), 2012  DOI: 10.1029/2012ja017683
  * [19] K.. Gorski et al.  “HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere”  In _The Astrophysical Journal_ 622.2 American Astronomical Society, 2005, pp. 759–771  DOI: 10.1086/427976
  * [20] Alexander V. Grayver, Alexey Kuvshinov and Dieter Werthmüller  “Time-domain modeling of three-dimensional Earth’s and planetary electromagnetic induction effect in ground and satellite observations”  In _Journal of Geophysical Research: Space Physics_ 126.3 American Geophysical Union (AGU), 2021  DOI: 10.1029/2020ja028672
  * [21] Richard Holme and Jeremy Bloxham  “The treatment of attitude errors in satellite geomagnetic data”  In _Physics of the Earth and Planetary Interiors_ 98.3-4 Elsevier, 1996, pp. 221–233
  * [22] P.. Huber  “Robust statistics”, Wiley Series in Probability and Statistics - Applied Probability and Statistics Section Series  Wiley, 2004  URL: https://books.google.dk/books?id=e62RhdqIdMkC
  * [23] Loïc Huder, Nicolas Gillet, Christopher C. Finlay, Magnus D. Hammer and Hervé Tchoungui  “COV-OBS.x2: 180 years of geomagnetic field evolution from ground-based and satellite observations”  In _Earth, Planets and Space_ 72.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1186/s40623-020-01194-2
  * [24] T. Iijima and T.. Potemra  “Large-scale characteristics of field-aligned currents associated with substorms”  In _Journal of Geophysical Research_ 83.A2 American Geophysical Union (AGU), 1978, pp. 599  DOI: 10.1029/ja083ia02p00599
  * [25] Takesi Iijima and T.. Potemra  “The amplitude distribution of field-aligned currents at northern high latitudes observed by Triad”  In _Journal of Geophysical Research_ 81.13 American Geophysical Union (AGU), 1976, pp. 2165–2174  DOI: 10.1029/ja081i013p02165
  * [26] Larry Kepko et al.  “Substorm current wedge revisited”  In _Space Science Reviews_ 190.1-4 Springer, 2015, pp. 1–46
  * [27] J.. King and N.. Papitashvili  “Solar wind spatial scales in and comparisons of hourly wind and ACE plasma and magnetic field data”  In _Journal of Geophysical Research: Space Physics_ 110.A2, 2005  DOI: https://doi.org/10.1029/2004JA010649
  * [28] K.. Laundal, C.. Finlay, N. Olsen and J.. Reistad  “Solar wind and seasonal influence on ionospheric currents from swarm and CHAMP measurements”  In _Journal of Geophysical Research: Space Physics_ 123.5 American Geophysical Union (AGU), 2018, pp. 4402–4429  DOI: 10.1029/2018ja025387
  * [29] Karl M. Laundal, Christopher C. Finlay and Nils Olsen  “Sunlight effects on the 3D polar current system determined from low Earth orbit measurements”  In _Earth, Planets and Space_ 68.1 SpringerOpen, 2016, pp. 142
  * [30] Karl Magnus Laundal and Arthur D. Richmond  “Magnetic coordinate systems”  In _Space Science Reviews_ 206.1-4 Springer, 2017, pp. 27–59
  * [31] Vincent Lesur, Ingo Wardinski, Martin Rother and Mioara Mandea  “GRIMM: the GFZ Reference Internal Magnetic Model based on vector satellite and observatory data”  In _Geophysical Journal International_ 173.2 Oxford University Press (OUP), 2008, pp. 382–394  DOI: 10.1111/j.1365-246x.2008.03724.x
  * [32] Philip W. Livermore, Rainer Hollerbach and Christopher C. Finlay  “An accelerating high-latitude jet in Earth’s core”  In _Nature Geoscience_ 10.1 Nature Publishing Group, 2017, pp. 62
  * [33] S. Macmillan and N. Olsen  “Observatory data and the Swarm mission”  In _Earth, Planets and Space_ 65.11 Springer ScienceBusiness Media LLC, 2013, pp. 1355–1362  DOI: 10.5047/eps.2013.07.011
  * [34] J. Matzka, C. Stolle, Y. Yamazaki, O. Bronkalla and A. Morschhauser  “The geomagnetic Kp index and derived indices of geomagnetic activity”  In _Space Weather_ American Geophysical Union (AGU), 2021  DOI: 10.1029/2020sw002641
  * [35] Jürgen Matzka, Oliver Bronkalla, Katrin Tornow, Kirsten Elger and Claudia Stolle  “Geomagnetic Kp index”  GFZ Data Services, 2021  DOI: 10.5880/KP.0001
  * [36] S. Maus et al.  “Third generation of the Potsdam Magnetic Model of the Earth (POMME)”  In _Geochemistry, Geophysics, Geosystems_ 7.7 American Geophysical Union (AGU), 2006  DOI: 10.1029/2006gc001269
  * [37] R.. McPherron, C.. Russell and M.. Aubry  “Satellite studies of magnetospheric substorms on August 15, 1968: 9. Phenomenological model for substorms”  In _Journal of Geophysical Research_ 78.16 American Geophysical Union (AGU), 1973, pp. 3131–3149  DOI: 10.1029/ja078i016p03131
  * [38] Christer Van Der Meeren et al.  “Apexpy version 1.1.0”  Zenodo, 2021  DOI: 10.5281/ZENODO.4585641
  * [39] P.. Newell and J.. Gjerloev  “Evaluation of SuperMAG auroral electrojet indices as indicators of substorms and auroral power”  In _Journal of Geophysical Research: Space Physics_ 116.A12 American Geophysical Union (AGU), 2011  DOI: 10.1029/2011ja016779
  * [40] P.. Newell, T. Sotirelis, K. Liou, C.. Meng and F.. Rich  “A nearly universal solar wind-magnetosphere coupling function inferred from 10 magnetospheric state variables”  In _Journal of Geophysical Research: Space Physics_ 112.A1 Wiley Online Library, 2007  DOI: 10.1029/2006JA012015
  * [41] Nils Olsen  “Ionospheric F region currents at middle and low latitudes estimated from Magsat data”  In _Journal of Geophysical Research: Space Physics_ 102.A3 Wiley Online Library, 1997, pp. 4563–4576
  * [42] Nils Olsen, Giuseppe Albini, Jerome Bouffard, Tommaso Parrinello and Lars Tøffner-Clausen  “Magnetic observations from CryoSat-2: calibration and processing of satellite platform magnetometer data”  In _Earth, Planets and Space_ 72 Springer, 2020, pp. 1–18
  * [43] Nils Olsen and Mioara Mandea  “Investigation of a secular variation impulse using satellite data: The 2003 geomagnetic jerk”  In _Earth and Planetary Science Letters_ 255.1-2 Elsevier BV, 2007, pp. 94–105  DOI: 10.1016/j.epsl.2006.12.008
  * [44] Nils Olsen, Mioara Mandea, Terence J. Sabaka and Lars Tøffner-Clausen  “CHAOS-2 — a geomagnetic field model derived from one decade of continuous satellite data”  In _Geophysical Journal International_ 179.3 Oxford University Press (OUP), 2009, pp. 1477–1487  DOI: 10.1111/j.1365-246x.2009.04386.x
  * [45] Nils Olsen, Mioara Mandea, Terence J. Sabaka and Lars Tøffner-Clausen  “The CHAOS-3 geomagnetic field model and candidates for the 11th generation IGRF”  In _Earth, Planets and Space_ 62.10 Springer ScienceBusiness Media LLC, 2010, pp. 719–727  DOI: 10.5047/eps.2010.07.003
  * [46] Nils Olsen et al.  “CHAOS—a model of the Earth’s magnetic field derived from CHAMP, Ørsted, and SAC-C magnetic satellite data”  In _Geophysical Journal International_ 166.1 Blackwell Publishing Ltd Oxford, UK, 2006, pp. 67–75
  * [47] Nils Olsen et al.  “The CHAOS-4 geomagnetic field model”  In _Geophysical Journal International_ 197.2 Oxford University Press, 2014, pp. 815–827
  * [48] A.. Richmond  “Ionospheric electrodynamics using magnetic apex coordinates”  In _Journal of geomagnetism and geoelectricity_ 47.2 Society of GeomagnetismEarth, PlanetarySpace Sciences, 1995, pp. 191–212  DOI: 10.5636/jgg.47.191
  * [49] Guillaume Ropp, Vincent Lesur, Julien Baerenzung and Matthias Holschneider  “Sequential modelling of the Earth’s core magnetic field”  In _Earth, Planets and Space_ 72.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1186/s40623-020-01230-1
  * [50] Martin Rother and Ingo Michaelis  “CH-ME-3-MAG - CHAMP 1 Hz combined magnetic field time series (level 3)”  GFZ Data Services, 2019  DOI: 10.5880/GFZ.2.3.2019.004
  * [51] Terence J. Sabaka, Nils Olsen and Robert A. Langel  “A comprehensive model of the quiet-time, near-Earth magnetic field: phase 3”  In _Geophysical Journal International_ 151.1 Oxford University Press (OUP), 2002, pp. 32–68  DOI: 10.1046/j.1365-246x.2002.01774.x
  * [52] Terence J. Sabaka, Nils Olsen and Michael E. Purucker  “Extending comprehensive models of the Earth’s magnetic field with Ørsted and CHAMP data”  In _Geophysical Journal International_ 159.2 Oxford University Press (OUP), 2004, pp. 521–547  DOI: 10.1111/j.1365-246x.2004.02421.x
  * [53] Terence J. Sabaka, Nils Olsen, Robert H. Tyler and Alexey Kuvshinov  “CM5, a pre-Swarm comprehensive geomagnetic field model derived from over 12 yr of CHAMP, Ørsted, SAC-C and observatory data”  In _Geophysical Journal International_ 200.3 Oxford University Press (OUP), 2015, pp. 1596–1626  DOI: 10.1093/gji/ggu493
  * [54] Terence J. Sabaka, Lars Tøffner-Clausen and Nils Olsen  “Use of the Comprehensive Inversion method for Swarm satellite data analysis”  In _Earth, Planets and Space_ 65.11 Springer, 2013, pp. 2
  * [55] Terence J. Sabaka, Lars Tøffner-Clausen, Nils Olsen and Christopher C. Finlay  “A comprehensive model of Earth’s magnetic field determined from 4 years of Swarm satellite observations”  In _Earth, Planets and Space_ 70.1 Springer ScienceBusiness Media LLC, 2018  DOI: 10.1186/s40623-018-0896-3
  * [56] Terence J. Sabaka, Lars Tøffner-Clausen, Nils Olsen and Christopher C. Finlay  “CM6: a comprehensive geomagnetic field model derived from both CHAMP and Swarm satellite observations”  In _Earth, Planets and Space_ 72.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1186/s40623-020-01210-5
  * [57] D. Saturnino, M.. Pais and J. Domingos  “The signature of geomagnetic field external drivers in virtual observatory 30-day means derived from Swarm data”  In _Journal of Geophysical Research: Space Physics_ American Geophysical Union (AGU), 2021  DOI: 10.1029/2021ja029579
  * [58] N. Schaeffer, D. Jault, H.-C. Nataf and A. Fournier  “Turbulent geodynamo simulations: a leap towards Earth’s core”  In _Geophysical Journal International_ 211.1 Oxford University Press (OUP), 2017, pp. 1–29  DOI: 10.1093/gji/ggx265
  * [59] P. Schwintzer, L. Grunwaldt and C. Förste  “CHAMP reference systems, transformations and standards”, 2002
  * [60] Andrey Sheyko, Christopher Finlay, Jean Favre and Andrew Jackson  “Scale separated low viscosity dynamos and dissipation within the Earth’s core”  In _Scientific Reports_ 8.1 Springer ScienceBusiness Media LLC, 2018  DOI: 10.1038/s41598-018-30864-1
  * [61] R.. Shore et al.  “Decadal period external magnetic field variations determined via eigenanalysis”  In _Journal of Geophysical Research: Space Physics_ 121.6 American Geophysical Union (AGU), 2016, pp. 5172–5184  DOI: 10.1002/2015ja022066
  * [62] M. Sugiura and T. Kamei  “Equatorial Dst index 1957–1986, IAGA Bull., 40”  In _by A. Berthelier and M. Menville (Int. Serv. Geomagn. Indices Publ. Off., Saint Maur, 1991)_ , 1991
  * [63] Lars Tøffner-Clausen and Jonas Nielsen  “Swarm - Level 1b processor algorithms”, 2019
  * [64] J. Torta, F. Pavón-Carrasco, Santiago Marsal and Christopher C. Finlay  “Evidence for a new geomagnetic jerk in 2014”  In _Geophysical Research Letters_ 42.19 American Geophysical Union (AGU), 2015, pp. 7933–7940  DOI: 10.1002/2015gl065501
  * [65] D.. Weimer  “An empirical model of ground-level geomagnetic perturbations”  In _Space Weather_ 11.3 American Geophysical Union (AGU), 2013, pp. 107–120  DOI: 10.1002/swe.20030
  * [66] Y. Yamazaki and A. Maute  “Sq and EEJ—A Review on the Daily Variation of the Geomagnetic Field Caused by Ionospheric Dynamo Currents”  In _Space Science Reviews_ 206.1-4 Springer ScienceBusiness Media LLC, 2016, pp. 299–405  DOI: 10.1007/s11214-016-0282-z