Polar ionospheric currents and high temporal resolution geomagnetic field models

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 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. 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.


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 km to 1000 km altitude. The ionospheric field undergoes daily, seasonal and solar cycle variations, which depend on solar activity and illumination (Yamazaki et al. 2016). 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 . 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 (Livermore et al. 2017), while highly time-dependent turbulent polar vortices are expected within the tangent cylinder (Aurnou et al. 2003;Schaeffer et al. 2017;Sheyko et al. 2018). 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 Olsen et al. 2009Olsen et al. , 2010Olsen et al. , 2006Olsen et al. , 2014, 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 Ropp et al. 2020 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. Friis-Christensen et al. 2017).
In the comprehensive modelling approach (Sabaka et al. 2002(Sabaka et al. , 2004(Sabaka et al. , 2015(Sabaka et al. , 2018(Sabaka et al. , 2020, 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 (Sabaka et al. 2020), 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 (Richmond 1995), 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 h, 12 h, 8 h, and 6 h, which are further modulated with annual and semiannual harmonics and scaled by a three-monthly average of F10.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 (Chulliat et al. 2013) 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 (Friis-Christensen et al. 2006). 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 (Baerenzung et al. 2022 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. Lesur et al. 2008 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 (Laundal et al. 2018) 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 (Richmond 1995) 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.

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 Hz magnetic data, product CH-ME-3-MAG (Rother et al. 2019), between January 2001 and August 2010, which we downsampled to 1 min values. We selected data according to the recommended quality flags that are provided in the distributed CHAMP data product files (GFZ Section 2.3 2019). 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 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 Olsen et al. 2020. We reduced the dataset to 1 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 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 min resolution. By using 20 s intervals for the linear fit, we followed Olsen et al. 2020, who recommends averaging over five successive values to reduce the intrinsic noise. In addition, we removed data if the attitude uncertainty qerror, which is provided in the CryoSat2 data product files, was larger than 40 arcseconds.
From the Swarm mission, we made use of the Level 1b 1 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 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 nT. We note that this approach also removed magnetic signals in the data associated with field-aligned currents, which can reach several thousands of 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 Kp index (Matzka et al. 2021a;Matzka et al. 2021b) was below 2o and the absolute rate of change of the RC index (Olsen et al. 2014), a quantitative measure of the magnetic disturbance at equatorial and mid-latitudes similar to the Dst index (Sugiura et al. 1991), was below 2 nT h −1 . Furthermore, we selected data if, on average over 2 h prior to the time of measurement, the Newell coupling function [Newell et al. 2007, for the exact definition used in this study, see Eq. (17)], 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,472,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. 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.

MODEL PARAMETRIZATION AND ESTIMATION
In this paper, we largely follow the modelling approach of the CHAOS geomagnetic field model series (Olsen et al. 2009(Olsen et al. , 2010(Olsen et al. , 2006(Olsen et al. , 2014, version CHAOS-7.9 ). 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 (Laundal et al. 2018). 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.

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 int , such that B int = −∇V int . In spherical coordinates V int is given by where a = 6371.2 km is the mean surface radius of the Earth, g m n are the spherical harmonic coefficients of degree n and order m, Y m n are the spherical harmonic functions, and N int = 55 is the chosen truncation degree to limit the spatial resolution of the model. The spherical harmonic functions are defined as where θ and ϕ are respectively the geocentric colatitude and longitude, and P m n are the associated Legendre functions using the Schmidt semi-normalization. We allow the spherical harmonic coefficients for n ≤ 20 to be time-dependent using a basis of 6th-order B-splines to account for the slow time changes of the internal field where B 6,k (k = 1, . . . , K) are the B-spline basis functions defined on the model interval using a sequence of knots with a 0.5 yr knot spacing and a 6-fold knot multiplicity at the model endpoints. The coefficients for 21 ≤ n ≤ N int are kept constant to represent the high-degree part of the assumed static lithospheric field.

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.

Ionospheric field
Following the approach of Laundal et al. 2016 we write the ionospheric magnetic field as where the poloidal magnetic field B pol , written in terms of the scalar potential V 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 B tor , written in terms of the potential T 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 Richmond 1995. There are two of these systems: Quasi-Dipole (QD) and Modified-Apex (MA). In QD coordinates the latitude is defined as where positive (negative) values refer to the northern (southern) magnetic hemisphere, h is the geodetic height of the point of interest, and hA 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 . In Modified-Apex (MA) coordinates the latitude is defined as where hR is a chosen reference height for the mapping, which we set to hR = 110 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 (Meeren et al. 2021), which is a wrapper of the Fortran library by Emmert et al. 2010. As the reference model for the field line tracing, we used the 13th generation of the International Geomagnetic Reference Field where θQD = π 2 − λQD and θMA = π 2 − λMA are the QD and MA colatitudes, respectively. We chose to truncate the spherical harmonic representations at N ion = 45 and N tor = 65. In addition, we used a maximum spherical harmonic order of M = 3 for both potentials in agreement with Laundal et al. 2018. Instead of the QD and MA longitudes, we used the Magnetic Local Time (MLT) where ϕnoon is the QD longitude of the subsolar point, computed on a sphere with radius r ≫ 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 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 −∇V 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 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 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 where {d1, d2, f1, f2} are the non-orthogonal base vectors for the magnetic apex coordinate systems ), andk is a unit vector in the geodetic upward direction. Here, we used thatr ≈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 |f1 × f2| based on the assumption that the vertical component of the ionospheric field scales with the linear dimension of the horizontal current system (Richmond 1995). The ionospheric magnetic field can be related to an electric sheet current density (in units of A m −1 ) that flows at a fixed height, chosen to be hR, written in the form of where J df is the divergence-free part of the sheet current density associated with B pol and J cf is the curl-free part associated with B tor . The potentials of the sheet current density parts are which were derived by treating the apex coordinates as if they were orthogonal, following the approach of Laundal et al. 2018. The curl-free part of the sheet current density can be furthermore related to an upward current density Ju (in units of A m −2 ) through Ju = −∇ · J cf , a statement of current continuity, which yields at the reference height At polar latitudes, where the magnetic field lines are close to vertical, Ju can be interpreted as field-aligned currents and J cf as the horizontal closure of these currents in the form of a sheet current. Instead of parametrizing the expansion coefficients g m,ion n and T m,ion n in time explicitly, we followed the climatological approach of the AMPS model (Laundal et al. 2018) and wrote these coefficients as linear combinations of external driving parameters Xi (i = 1, . . . , 19) so that g m,ion n (t) = g m,ion n,0 and similarly for T m,ion n . The Xi are combinations of solar wind parameters and IMF components that have been found suitable for characterising the external driving of the ionospheric current system (Laundal et al. 2018;Weimer 2013 which are all functions of time. The terms in Eq. (14) involve the clock angle where the components of the IMF, BIMF,y and BIMF,z, are with respect to the GSM frame; the dipole tilt angle whereŝ is a unit vector in the direction of the sun andm dip is the dipole moment of the IGRF magnetic field, parametrizes seasonal effects; the solar wind-magnetospheric coupling function (Newell et al. 2007) where Bt = B 2 IMF,y + B 2 IMF,z (given in nT) and vsw (given in km 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 τ = 10 −3 |vsw| 4/3 Bt 2/3 cos 8/3 θc 2 solar flux, sfu ≡ 10 −22 W m −2 Hz −1 , parametrizes solar cycle variations. Expanding on the original parametrization of AMPS, we included as X19 the SML index (Newell et al. 2011), developed by the SuperMAG initiative (Gjerloev 2009(Gjerloev , 2012, to parametrize indirectly driven currents in the polar ionosphere. With typically more than 100 contributing ground-based magnetometer stations, the SML index can be considered an extension of the traditionally used 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.

Magnetospheric field
The magnetic field B mag produced by electric currents in the magnetosphere can be separated into contributions due to the ring current in the near-Earth magnetosphere, B near , and the currents in the remote magnetosphere, B far , so that B mag = B near + B 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 B near = −∇V 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 . The scalar potential is given by where N near = 2 is the chosen truncation degree,q m,SM 1 are constant regression parameters multiplying the RC index, which consists of an internal part, RC i, and an external part, RC e, so that RC = RC i + RC e. We estimated the spherical harmonic coefficients q m,SM n 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 RC index into internal and external parts. The reader is referred to Finlay et al. 2020 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 B far = −∇V 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 ). The potential is given by where N mag = 2 is the chosen truncation degree. The treatment of the internally induced part is similar to the approach used for B near .

Alignment parameters
We estimate alignment parameters in the form of three Euler angles α, β and γ 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 where BCRF and BVFM are column vectors that contain the magnetic field components with respect to the CRF and the VFM frame, respectively, and R1, R2 and R3 are rotation matrices given by 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.

Model estimation
The model parameters are arranged into a column vector m = [p T , q T ] T , where the column vector p contains the parameters of the geomagnetic field model and the column vector 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 where g(p) is a column vector containing the model estimates of the magnetic vector components, d(q) is a column vector containing the aligned magnetic vector observations expressed in terms of spherical geocentric components, C −1 d is the inverse of the data error covariance matrix, and Λ is the model regularisation matrix. At each iteration k the model parameter vector is updated through where g k = g(p k ), d k = d(q k ) and G k is the matrix of partial derivatives of the residuals with respect to the model parameter vector

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 B,n × B and B × (n × B), wheren is the star camera bore sight assumed not parallel to B. In this reference frame the data error covariance matrix is diagonal and given by Holme et al. 1996 where σ (in nT) is an isotropic instrument error in the vector component magnitudes, χ (in radians) is an error in the attitude aboutn, ψ (in radians) is an error in the attitude about the two axes perpendicular ton. In the B23 frame, we multiplied the inverse of the data error covariance matrix, which is diagonal, by the Huber weights (Constable 1988;Huber 2004), 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 θ 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. The star camera bore sightn, 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 ofn is in fact arbitrary since we assume ψ and χ to be equal. However,n is important in the case of CHAMP when only one of the two head units was active at a time.

Model regularisation
The model regularisation matrix Λ 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 timederivative of the radial internal field at the Core-Mantle Boundary (CMB) (r = 3485 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, ts = 2001.0 and te = 2022.0 in units of decimal years, and averaged over the CMB. The corresponding regularisation parameters are λt, λt s and λt e , respectively. The time variation of each RC -baseline correction, {q −1,SM 1 , q 0,SM 1 , q 1,SM 1 }, was minimised using a quadratic norm of the bin-to-bin differences, which is scaled by the regularisation parameter λ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 B 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 where S(r0) is the spherical surface of radius r = r0 ≡ a + hR, J df s with s ∈ {r, θ, ϕ} 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 ti 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 λ 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 θMA. The associated regularisation matrix Λ tor is diagonal with entries n(n+1) 2n+1 , which depend on the degree of the expansion coefficients T m,ion n . The associated regularisation parameter is λ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.
The iterative minimisation for the model estimation was initialised with a starting model m0, 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 (Schwintzer et al. 2002) and Swarm (Tøffner-Clausen et al. 2019), and during in-flight calibration for CryoSat2-1 . Convergence was typically reached after 15 iterations. 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.

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 ∆Br, ∆B θ and ∆B ϕ , and the difference in scalar magnitude, ∆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 nT width. 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 θ and ϕ 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 (|λQD| > 55 • ) and non-polar (|λQD| ≤ 55 • ) latitudes. 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 ∆B θ , ∆B ϕ , and ∆F for CHAMP in the polar regions are 17.35 nT, 18.89 nT, and 6.20 nT for Model-A and 21.07 nT, 23.14 nT, and 8.52 nT for the reference model, respectively, which corresponds to a reduction of approximately 20 %. 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 (Gorski et al. 2005) 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. 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. Dungey 1961). 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. Yamazaki et al. 2016). 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 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 |λ QD | > 55 • and non-polar QD latitude to the opposite.

Model-A Reference
Mean ( Network (  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 (Sabaka et al. 2020), the sum of primary and secondary parts, at satellite altitude at 450 km during noon in Greenwich on March 21, 2018. 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.

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. 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 (Iijima et al. 1978;Iijima et al. 1976). For northward IMF (θc = 0 • ) 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 θc = 90 • leads to stronger currents than a counter-clockwise rotation to θc = −90 • , 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 (θc = 180 • ) 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 SML index and only for purely northward IMF. With decreasing 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 (SML = −40 nT) to the midnight sector (SML = −160 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 (Kepko et al. 2015;McPherron et al. 1973). Again, however, the case for SML = −160 nT should be considered an extrapolation given that there was a relatively small amount of data with large negative values of the SML index available during the modelling thanks to the quiet-time data selection (for the distribution of the 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 SML index and the ϵ coupling function probably compete to some extent for the same signal. So a part of the structure which may go into the ϵ terms, if the SML index was not included in the parametrization, may now be contained in the 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 SML index, because some of the signal for southward IMF conditions is now contained in the 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 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 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.

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. 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 highdegree 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. 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. Sabaka et al. 2013). 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 ≈ 0) and near-sectorial (m ≈ 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  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 (Livermore et al. 2017), 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.

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 λt, λt s and λt e by a factor of 80 to 0.0125, 1.25 × 10 −4 and 1.25 × 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. 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ġ 0 1 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ġ 0 1 , thė h 1 12 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 , GRIMM (Lesur et al. 2008), or CM (Sabaka et al. 2020], low resolution basis functions [e.g., POMME (Maus et al. 2006), CovObs (Huder et al. 2020)], or by including other internal sources [e.g., Kalmag models (Baerenzung et al. 2022), or sequential models of Ropp et al. 2020].
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. Shore et al. 2016) and magnetic observations made by satellites (e.g. Domingos et al. 2019;Saturnino et al. 2021). 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 (Gorski et al. 2005), covering the entire Earth's surface at approximately 2 • 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 X of size 201 × 32400. Finally, we centred X by removing the column-wise mean value,X, such that For a data matrix such asX, 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 ith mode obtained through the PCA consists of the Empirical Orthogonal Function (EOF), vi, which represents the spatial pattern of the time variation, and the principal component (PC) yi =Xvi, which is the timeseries of variance σ 2 i . These modes are typically sorted in decreasing order of variance and can be used to reconstruct the data matrix through or, if a subset of modes is chosen, to perform a partial reconstruction. We applied the PCA onX 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 % of the variance (13 modes account for 90 %). 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 0 2 except that the zero lines in latitude are shifted slightly northward. Given the spatiotemporal 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ġ 0 1 in Fig. 13 (phase shift betweenġ 0 1 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 (Chulliat et al. 2014;Finlay et al. 2020;Olsen et al. 2007;Torta et al. 2015). 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. 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 • and 90 • 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.

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 0 1 and g 0 3 , 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 2-5 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. Fillion et al. 2023;Olsen 1997). 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 Qresponse functions (e.g. Grayver et al. 2021) 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.

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 nonaxisymmetric 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 (Livermore et al. 2017).
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.

ACKNOWLEDGMENTS
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.

DATA AVAILABILITY
The data underlying this article are available in the following repositories: (i) Swarm and CryoSat2 data are available through ESA at https://swarm-diss.eo.esa.int/#. (ii) CHAMP data are available from Rother et al. 2019. (iii) 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/.
(vi) Solar wind speed, interplanetary magnetic field and Kp index are available through the GSFC/SPDF OMNIWeb interface at https: //omniweb.gsfc.nasa.gov/ow.html.
Files containing the estimated parameters of our models and the processed data are available at https://doi.org/10.11583/DTU. 24025596.