Co-estimating geomagnetic field and calibration parameters: modeling Earth's magnetic field with platform magnetometer data

Models of the geomagnetic field rely on magnetic data of high spatial and temporal resolution. The magnetic data from low-Earth orbit satellites of dedicated magnetic survey missions such as CHAMP and Swarm play a key role in the construction of such models. Unfortunately, there are no magnetic data from such satellites after the end of CHAMP in 2010 and before the launch of Swarm in late 2013. This limits our ability to recover signals on timescales of 3 years and less during this gap period. The magnetic data from platform magnetometers carried by satellites for navigational purposes may help address this data gap provided that they are carefully calibrated. Earlier studies have demonstrated that platform magnetometer data can be calibrated using a fixed reference field model. However, this approach can lead to biased calibration parameters. An alternative has been developed in the form of a co-estimation scheme which consists of simultaneously estimating both the calibration parameters and a model of the internal geomagnetic field. Here, we develop a scheme, based on the CHAOS framework, that involves the co-estimation of a geomagnetic field model along with calibration parameters of platform magnetometers. Using our implementation, we are able to derive a geomagnetic field model from 2008 to 2018 with satellite magnetic data from CHAMP, Swarm, secular variation data from ground observatories, and platform magnetometer data from CryoSat-2 and GRACE. Through experiments, we explore correlations between the estimates of the geomagnetic field and the calibration parameters, and suggest how these may be avoided. We find that platform magnetometer data provide additional information on the secular acceleration, especially in the Pacific during the gap period. This study adds to the evidence that it is beneficial to use platform magnetometer data in geomagnetic field modeling.


Introduction
The Earth's magnetic field is a superposition of many sources. By far, the largest contribution comes from within the Earth at a depth of more than 3000 km. There, in the outer core, a liquid iron alloy is rapidly moving and thus advecting, stretching, arXiv:2103.09078v1 [physics.geo-ph] 16 Mar 2021 and maintaining the ambient magnetic field against dissipation in a process called the Geodynamo. Earth's core dynamics are not fully understood, but can be studied using time-dependent geomagnetic field models. Such models are constructed using measurements of the magnetic field taken at and above Earth's surface.
The study of core processes on decadal or longer timescales requires long timeseries of magnetic vector data with high spatial and temporal resolution. Along with ground-based magnetic observatories, low-earth orbit satellites from dedicated magnetic survey missions such as CHAllenging Minisatellite Payload (CHAMP, 2000(CHAMP, -2010 and the Swarm trio (since 2013) provide such data. However, other than scalar data from Ørsted, no high-quality calibrated magnetic vector data from satellites are available between the end of the CHAMP mission in September 2010 and the launch of the Swarm satellites in November 2013. This data gap not only cuts in two an otherwise uninterrupted time-series of high-quality magnetic satellite data since the year 2000, but also limits our ability to derive accurate core field models that resolve temporal changes of the magnetic field on timescales of a few years and less in the gap period. To address the issue, one can utilize the crude magnetometers that are carried by most satellites for navigational purposes, the so-called platform magnetometers. Although not a substitute for dedicated highquality magnetic survey satellites, platform magnetometers can supplement ground observatory data in gaps between dedicated missions and help improve the local time data coverage of simultaneously flying high-quality magnetic survey satellites.
Satellite-based magnetic vector data need to be calibrated to remove magnetometer biases, scale factors, and non-orthogonalities between the three vector component axes (Olsen et al., 2003). Comparing the vector magnetometer output with a magnetic reference field allows the estimation of these calibration parameters. On dedicated survey mission satellites, the reference is a second, absolute scalar, magnetometer mounted in close proximity to the vector magnetometer and measuring the magnetic field intensity. However, non-dedicated satellites carrying platform magnetometers are typically not equipped with such scalar reference magnetometers. In this case, it is possible to use a-priori geomagnetic field models like CHAOS Olsen et al., 2006) or the IGRF (Thébault et al., 2015) as reference. Such an approach has been successfully used, e.g., by Olsen et al. (2020) for calibrating data from the CryoSat-2 magnetometer, but use of a fixed reference field model is not without risks and could lead to biased calibration parameters.
An alternative venue has been explored by Alken et al. (2020), who combined high-quality magnetic data from CHAMP and Swarm with platform magnetometer data from CryoSat-2 and several satellites of the Defense Meteorological Satellite Program (DMSP) to estimate a model of the internal field and the required calibration parameters for each satellite simultaneously. Ideally, such a co-estimation scheme eliminates the need for a-priori geomagnetic field models, but Alken et al. (2020) fall short by co-estimating only the internal field while still relying on a fixed model of the external field. Nevertheless, their study convincingly demonstrated that platform magnetometer data provide valuable information about the time-dependence of Earth's magnetic field.
In this study, we followed Alken et al. (2020) and developed a co-estimation strategy but within the framework of the CHAOS field model series. Our implementation differs in three important aspects. First, we estimated both the internal (core and crust) and external (magnetospheric) geomagnetic field contributions in contrast to only the internal field. This way, we avoided having to remove a fixed external field model from the satellite data prior to the model parameter estimation. Following the methodology of the CHAOS model, we did use a prior external field model for processing the ground observatory data which we used in addition to the satellite data. Second, we used the platform magnetometer data from CryoSat-2 and, instead of DMSP, data from the Gravity Recovery and Climate Experiment (GRACE) satellite pair. Finally, to reduce the significant correlation between the internal axial dipole and the calibration parameters during periods of poor coverage of high-quality magnetic data, we excluded platform magnetometer data from determining the internal axial dipole (its time variation is well resolved with ground observatory data during the gap period, while its absolute value is constrained by Swarm and CHAMP data on both sides of the gap) rather than controlling the temporal variability of the internal axial dipole through an additional regularization as done by Alken et al. (2020).
The paper is organized as follows. In the first part, we present the datasets and the data processing. Next, we describe the model parameterization and define the calibration parameters, which are similar to those used for the Ørsted satellite (Olsen et al., 2003). We go on by presenting a geomagnetic field model derived from high-quality calibrated data from the CHAMP and the Swarm satellites as well as ground observatory secular variation data and supplemented this with previously uncalibrated platform magnetometer data from CryoSat-2 and GRACE, spanning a 10 year period from 2008 to 2018. Finally, we explore in a series of experiments the effect of co-estimating an external field, the trade-off between the internal dipole and the calibration parameters, and the importance of including dayside platform magnetometer data when estimating calibration parameters. We conclude the paper by looking at the secular acceleration of our model, paying particular attention to the data gap between 2010 and 2013.

Data and data processing
We used calibrated magnetic data from the Swarm satellites Alpha (Swarm-A) and Bravo (Swarm-B), and from the CHAMP satellite from January 2008 to the end of December 2017, supplemented with five datasets of uncalibrated magnetic data from the three platform fluxgate magnetometers (FGM) on-board the CryoSat-2 satellite (CryoSat-2 FGM1, CryoSat-2 FGM2 and CryoSat-2 FGM3), the one on-board the first GRACE satellite (GRACE-A), and the other one on-board the second GRACE satellite (GRACE-B). In addition to the satellite data, we included revised monthly mean values of the SV from ground observatories to contribute to the Earth's internal time-dependent field. Details of the datasets are given in the following.

Absolute satellite data from scientific magnetometers
The satellite data from scientific magnetometers are in general of high quality in terms of accuracy, precision and magnetic cleanliness. The high standard of the data is achieved by low noise instruments that are mounted together with star cameras on an optical bench further away from the spacecraft body at the center of a several meter long boom. The data are regularly calibrated in-flight with a second absolute scalar magnetometer placed at the end of the boom and carefully cleaned from magnetic disturbance fields originating from the spacecraft body.
From the CHAMP mission, we used the Level 3 1 Hz magnetic data, version CH-ME-3-MAG (Rother and Michaelis, 2019), between January 2008 and August 2010, downsampled to 15 s, and only when attitude information from both star cameras was available. From the Swarm mission, we used the Level 1b 1 Hz magnetic data product, baseline 0505/0506, from the Swarm-A and Swarm-B satellites between November 2013 and December 2018, also downsampled to 15 s. Here, we worked with vector data from CHAMP and Swarm in the magnetometer frame.

Relative satellite data from platform magnetometers
Relative satellite data refer to the raw sensor output from platform magnetometers. The data have to be corrected and calibrated before they can be used in geomagnetic field modeling. The correction of the data accounts for temperature effects, magnetic disturbances due to solar array and battery currents, magnetorquer activity, as well as non-linear sensor effects, whereas the calibration removes magnetometer biases, scale differences, and non-orthogonalities between the three vector component axes.
From CryoSat-2, we took magnetic data, baseline 0103, from the three platform magnetometers as described in Olsen et al. (2020) from August 2010 to December 2018 and only when the attitude uncertainty q error was below 40 . Since the purpose of this paper is the co-estimation of calibration parameters for the platform magnetometers, we processed the dataset using the original calibration parameters to undo the calibration step that has been performed by Olsen et al. (2020) but keeping the applied correction for magnetic disturbances from the spacecraft and its payload. This way, we obtained essentially uncalibrated data while still retaining the corrections for magnetic disturbances, temperature effects and nonlinearities. In a pre-whitening and data reduction step, we computed residuals to the CHAOS-6-x9 model in the uncalibrated magnetometer frame, removed those larger than 1000 eu (quasi nanoTesla, in the following referred to as engineering units) in absolute value to discard gross outliers, computed component-wise robust mean values of the residuals in 1 min bins to reduce the original 4 s sampled data to 1 min values, and added the CHAOS-6-x9 model values back. Fig. 1 shows an example of the raw vector residuals ∆E of CryoSat-2 FGM1 in the uncalibrated magnetometer frame over 3 h on March 24, 2016. In a similar way, we processed the 1 Hz data from the GRACE satellites, baseline 0101, to obtain 1 min uncalibrated but corrected vector data between January 2008 and October 2017 (GRACE-A) and August 2017 (GRACE-B) (Olsen, 2020).
The computation of 1 min values served two purposes. First, to reduce the random noise of the magnetometers by taking the average of successive values and, second, to decrease the number of platform magnetometer data, so that a fair amount of absolute satellite data was able to guide the co-estimation of the calibration parameters.

Ground observatory data
In addition to satellite data, we added annual differences of monthly mean values from 162 ground observatories to help determine the time changes of the core field  Figure 1 Residuals of raw vector data from CryoSat-2 FGM1 with respect to the CHAOS-6-x9 model values in the uncalibrated magnetometer frame for an example period of 3 h on March 24, 2016. The gap in the raw data between 10:30 and 10:40 is due to the rejection of data with poor attitude information (qerror > 40 ).
(secular variation). Following Olsen et al. (2014), we computed revised monthly means as Huber-weighted averages of the hourly observatory mean values from the AUX OBS database (Macmillan and Olsen, 2013) at all local times after removing estimates of the ionospheric field of the CM4 model (Sabaka et al., 2004) and the large-scale magnetospheric field of CHAOS-6-x9, including their internally induced parts.

Satellite data selection
We organized the satellite data according to quasi-dipole (QD) latitude (Richmond , 1995) into a non-polar (equal to and equatorward of ±55°) and a polar (poleward of ±55°) data subset. From each subset, we selected data under quiet geomagnetic conditions. Specifically, we selected data from the non-polar subset that satisfied the following criteria: • Low geomagnetic activity as indicated by the planetary activity index Kp smaller than or equal to 2 o ; • Dark condition as indicated by a solar zenith angle greater than 100°for the Swarm and CHAMP satellites (i.e., sun at least 10°below the horizon). From CryoSat-2 and GRACE, we used data from dark and sunlit regions, since we found that this leads to better determined calibration parameters; • Slow change of the magnetospheric ring current as indicated by the RCindex (Olsen et al., 2014) rate of change in absolute terms being smaller than 2 nT h −1 . From the polar subset, we kept data according to the following criteria: • Dark condition except in the case of platform magnetometers on-board CryoSat-2 and GRACE, where we also used sunlit data; • RC-index rate of change in absolute terms smaller than or equal to 2 nT h −1 ; • The merging electric field at the magnetopause E m = v 4/3 B 2/3 T sin|Θ|/2, where v is the solar wind speed, B T = B 2 y + B 2 z is the interplanetary magnetic field in the y-z-plane of the Geocentric Solar Magnetic (GSM) coordinates, and Θ = arctan(B y /B z ), was on average smaller than 2.4 mV m −1 over the previous 2 h; • The interplanetary magnetic field component B z in GSM coordinates was on average positive over the previous 2 h. Fig. 2 shows a stacked histogram of the number of data for each satellite after the data selection. It can be clearly seen that platform magnetometer data are the main contributor to the number of data in the gap period, whereas it is comparable to the number of data from CHAMP and the Swarm satellites in the time before and after the gap. The ground observatories contribute approximately 130 monthly mean values of the SV each month throughout the entire model time span, which is much less than the monthly average number of satellite data.

Model parameterization and estimation
We are interested in the magnetic field vector B on length scales smaller than Earth's circumference and time scales that are much longer than the time it takes light to traverse these distances (Backus et al., 1996;Sabaka et al., 2010). On these scales, the displacement current can be neglected and the magnetic field is governed by Ampere's law. We assume that the measurements of Earth's magnetic field are taken in a region free of electrical currents and magnetized material, such that the field is irrotational, which allows us to introduce a scalar potential V to represent the magnetic field as the gradient of the potential B = −∇V . The potential consists of two terms V = V int + V ext that describe internal sources such as the time-dependent core-generated field and the assumed static lithospheric field, and external sources that we assume are mainly magnetospheric in origin for our chosen data selection criteria and have an internally induced counterpart associated with them (by selecting data from dark regions, we minimize ionospheric field contributions).
To describe the geomagnetic field, we use an Earth-fixed frame of reference whose point of origin coincides with the Earth's center and in which the position vector r is given in spherical polar coordinates by the radial distance r as measured from the origin (radius), the angular distance θ (co-latitude) as measured from the north polar axis, and the azimuthal angular distance φ (longitude) as measured from the Greenwich meridian. In the following, we refer to that system as the Radius-Theta-Phi (RTP) reference frame.
In spherical coordinates the scalar potential can be expressed as a weighted sum of solid harmonics, which are harmonic functions of the spatial coordinates. Our modeling approach follows that of earlier models of the CHAOS model series (Finlay et al., 2016Olsen et al., 2006Olsen et al., , 2014 and consists of describing the geomagnetic field with the help of a scalar potential whose exact form depends on a set of coefficients that multiply the solid harmonics. The coefficients are estimated by minimizing a quadratic cost function in the residuals, the difference between the magnetic observations and the magnetic data calculated with the model. We used two kinds of residuals: the components of vector differences in the RTP frame (vector residuals) and the difference of vector magnitudes (scalar residuals). More specifically, we computed vector residuals of the non-polar satellite data, scalar residuals of the polar satellite data, and vector residuals of the ground observatory SV data at all QD latitudes.

Internal field parameters
The scalar potential of the internal sources is given by where a = 6371.2 km is the chosen spherical reference radius of the Earth, n and m are, respectively, the spherical harmonic degree and order, N int is the truncation degree, g m n (t) and h m n (t) are the Gauss coefficients in nanoTesla (nT) for a given n and m, and P m n (cos θ) are the Schmidt quasi-normalized associated Legendre functions. We truncated the formally infinite sum of solid harmonics at N int = 50 and expanded the Gauss coefficients of degree n ≤ 15 in time using sixth-order B-splines (De Boor , 1978), while we kept the higher degree coefficients (n > 15) constant in time where g m n,j (similarly for h m n,j ) is the coefficient of B 6,j (t)-the jth function of the B-spline basis that has knots at 6-month intervals and six-fold multiplicity at the model endpoints in t s = 2008.0 and t e = 2018.0 in years. For the purposes of testing the co-estimation of calibration parameters here, a truncation of the time-dependent internal field at degree N int = 15 was deemed sufficient.

External field parameters
The scalar potential of the external sources consists of two terms V ext = V SM +V GSM that are designed to account for near and remote magnetospheric sources. We use the Solar Magnetic (SM) coordinate system to parameterize near magnetospheric sources where θ SM and φ SM are, respectively, the SM co-latitude and longitude, q m n,SM and s m n,SM are the Gauss coefficients with respect to the SM coordinate system, ∆q m 1,SM (t) and ∆s m 1,SM (t) are the RC-baseline corrections, and R m,s n,SM and R m,c n,SM are modification of the solid harmonics that account for the time-dependent transformation from the SM to the geographic coordinate system and include internally induced contributions based on the diagonal part of the Q-response matrix that has been derived from a 3D conductivity model of Earth . The external Gauss coefficients with n = 1 have a specific time-dependence in the form of where (t) and ι(t) are the respective internal and external part of the RC-index linearly interpolated from hourly values. The RC-baseline corrections were estimated in bins of 30 days except in the gap period, where we used a single bin from August 2010 to January 2014 to reduce the strong co-linearity between the calibration parameters and the baseline corrections that earlier tests had revealed. The remote magnetospheric sources, the currents at the magnetopause and in the magnetotail, are taken into account by a purely zonal potential in the GSM coordinate system up to degree 2 where q m n,GSM and s m n,GSM are Gauss coefficients that are constant in time with respect to the GSM coordinate system, and R m,c n,GSM are modifications of the solid harmonics similar to corresponding terms in Eq. (3) but for the GSM coordinates.

Alignment parameters
Using satellite data in the vector field magnetometer frame (VFM) requires an additional step, called data alignment, which involves determining alignment parameters that describe the rotation of the magnetic field vector B VFM in the VFM frame to B CRF in the common reference frame (CRF) of the satellite. Once in the CRF, the vector components can be combined with the attitude information from the star camera and rotated into the RTP frame for computing the vector residuals. We performed the data alignment for CHAMP, Swarm, CryoSat-2, and GRACE.
The alignment parameters are usually parameterized in the form of Euler angles α, β, and γ. We adopted the 1-2-3 convention of the Euler angles to align the magnetic field where the rotation matrix is a combination of the three rotations Following the alignment, we applied another rotation matrix R CRF RTP to rotate the field components from the CRF to the RTP reference frame which depends on position and time. That rotation matrix was computed by combining the quaternions that express the rotation from the CRF to the Earth-fixed Earth-centered North-East-Center (NEC) frame with quaternions that describe the change from the NEC to the RTP reference frame. For each satellite dataset, we parameterized the Euler angles in time as a piecewise constant function using a sequence of 30 day bins.

Calibration parameters
The calibration can be viewed as an extension of the data alignment which makes it possible to use platform magnetometer data in geomagnetic field modeling. We performed the calibration for CryoSat-2 and the GRACE satellites.
We assume that the platform magnetometer is a linear vector field magnetometer, which provides information about the desired local magnetic field vector B VFM (units of nT) in the form of the sensor output E = (E 1 , E 2 , E 3 ) T (units of eu), which typically consists of components that are measured relative to three biased and non-orthogonal axes employing different scale factors (Olsen et al., 2003). More specifically, the sensor output in the magnetometer frame is related to the local magnetic field through where is the diagonal matrix of sensitivities or scale factors s = (s 1 , s 2 , s 3 ) T (units of eu nT −1 ), is the matrix that projects the orthogonal components of magnetic field vector B VFM onto three non-orthogonal directions defined by the non-orthogonality angles is the offset or bias vector (units of eu). Combining the calibration step in Eq. (9), the alignment step involving the Euler angles in Eq. (6) and the change of frame in Eq. (8), yields an equation that transforms the uncalibrated sensor output E into calibrated, aligned field components in the RTP frame We estimated the nine basic calibration parameters and the three Euler angles in bins of 30 days. For data equatorward of ±55°QD latitude, we performed a vector calibration using the component residuals of B RTP for estimating the model parameters (see Sec. Model parameter estimation). In contrast, for data poleward of ±55°QD latitude, we performed a scalar calibration by using the residuals of the vector magnitude, in which case the rotation matrices from the VFM to the RTP frame including the Euler angles disappear at the expense of loosing the ability to estimate the Euler angles. Tab. 1 summarizes the different parts of the model and the corresponding number of parameters.

Model parameter estimation
The geomagnetic field model parameters p, the Euler angles q, and the calibration parameters e were derived by solving the least-squares problem where m = (p T , q T , e T ) T is the entire model parameter vector, and Φ is the cost function which penalizes a quadratic form in the residuals-the difference between the computed geomagnetic field model values g(p) and the calibrated, aligned magnetic data d(q, e)-using the inverse of the data covariance matrix C d , and a quadratic form in the model parameter vector using the regularization matrix Λ. For the definition of the matrices C d and Λ, see, respectively, Secs. Data weighting and Model regularization.
The least-squares solution m * in Eq. (15) is found through an iterative quasi-Newton method, which consists of updating the model parameter vector m k at iteration k using m k+1 = m k + ∆m together with where d k = d(q k , e k ), g k = g(p k ), and G k is a matrix with entries corresponding to the partial derivative of the ith residual with respect to the jth model parameter 1 Entries related to g 0 1 B-spline coefficients and platform magnetometer data are zero.
evaluated at iteration k (Tarantola, 2005, p. 69). Some entries of G k are zero owing to data subsets that do not provide information on parts of the model. For example, the scalar data do not constrain the Euler angles and the vector data from one magnetometer do not constrain the Euler angles associated with another magnetometer. With the same idea in mind, we modified entries of G k to prevent some data subsets from constraining certain parts of the internal field model. In particular, we set entries to zero for the following criteria: 1 The row index of the matrix entry corresponded to dayside data from a platform magnetometer, on-board CryoSat-2 or GRACE, and the column index corresponded to model parameters that describe the internal and external magnetic field. Therefore, the dayside data were only used to constrain the Euler angles and calibration parameters of the respective platform magnetometer. 2 The row index of the matrix entry corresponded to data from a platform magnetometer, on-board CryoSat-2 or GRACE, and the column index corresponded to the B-spline parameters that parameterize the g 0 1 Gauss coefficient of the internal field in time. Therefore, no platform magnetometer data were used to constrain the B-spline coefficients of the axial dipole which we believe are well determined using ground observatory data. Tab. 2 gives an overview of whether or not certain datasets constrained specific parts of the model. Nevertheless, we used the full model description in the forward evaluation to compute the residuals.
The iterative procedure described in Eq. (17) requires a starting model m 0 to initialize the model parameter estimation. We initialized the internal field model parameters using the corresponding part of CHAOS-6-x9, while we set the external field model parameters to zero. To initialize the Euler angles, we used the values from CHAOS-6-x9 in case of Swarm and CHAMP satellites, or set the angles to zero in case of CryoSat-2 and the GRACE satellite duo. For the calibration parameters, we simply set the offsets and non-orthogonalities to zero and the sensitivities to one over the whole time span. The parameter estimation usually converged after 10-15 iterations. We also tested other starting models, e.g. random calibration parameters, but found that our choice had little impact on the converged model parameters other than increasing the number of necessary iterations.

Data weighting
For the vector components of the non-polar satellite data, we used a covariance matrix that accounts for the attitude uncertainty of the star cameras with respect to the B23 reference frame defined by unit vectors in the direction of B, n × B, and n × (n × B), where n is an arbitrary unit vector not parallel to B that we chose to be the third CRF base vector, σ 2 is the variance of an isotropic instrument error and ψ 2 is the variance associated with random rotations around the three reference axes (Holme and Bloxham, 1996). Tab. 3 summarizes the values of σ and ψ for the different satellite datasets. We scaled the diagonal entries of the covariance matrix with Huber weights (Constable, 1988;Sabaka et al., 2004) that we calculated for each component in the B23 reference frame to downweight data points that greatly deviated from the model evaluated at the previous iteration. After inverting and rotating the Huber-weighted covariance matrix of the individual data point into the RTP frame, we arranged them into a block-diagonal matrix completing the desired inverse data covariance matrix C d −1 . In case of the vector magnitude of the polar satellite data, we simply used σ 2 scaled with Huber weights as variance. The covariance of the ground observatory SV vector data was derived from detrended residuals to the CHAOS-6-x9 model, including the covariance between vector components at a given location.

Model regularization
The regularization in the form of the matrix Λ in Eq. (15) is designed to ensure the convergence of the model parameter estimation by limiting the flexibility of the model. The regularization matrix is block diagonal and consists of the blocks Λ int , Λ ext , and Λ cal , which regularized the internal, external, and the calibration parameters, respectively. We did not regularize the Euler angles, such that corresponding blocks in the regularization matrix are zero.
Turning to the internal part of the model, following the example of earlier models in the CHAOS series, we designed a regularization based on the square of the third time-derivative of the radial field component B r integrated over the core mantle boundary (CMB) and averaged over the entire model time span where c = 3485.0 km is the chosen spherical reference radius of the CMB, Ω(c) denotes the CMB given as the spherical surface of radius c, and dΩ = sin θdθdφ is the surface element for the integration. Furthermore, we set up a regularization of the internal field based on the square of the second time-derivative of the radial component integrated over the CMB at the model start time t s and similarly for the end time by replacing t s with t e . Returning to Eq. (20), thanks to the orthogonality of spherical harmonics on the surface of the sphere, carrying out the spatial integration leads to where w Ω = (n+1) 2 2n+1 a c 2n+4 is a spatial factor that follows from the surface integration and · t = 1 te−ts te ts dt denotes the time average. Utilizing the fact that the time-dependence of the Gauss coefficients is given by sixth-order B-splines, terms such as can be written as a quadratic form in g m n = (g m n,1 , g m n,2 , . . . ) T , the vector of the spline coefficients of g m n , using the matrix A t that has entries corresponding to the time averages of products of the third time-derivative of the B-splines. While the time-derivatives of the B-splines are known analytically, we approximated the time average numerically by a Riemann sum of rectangles. A similar computation of Eq. (21), now evaluating the derivatives only at the endpoints instead of averaging in time, yields matrices A ts jj =B 6,j (t s )B 6,j (t s ) and A te jj =B 6,j (t e )B 6,j (t e ). Finally, based on the physical quantities in Eqs. (20) and (21), we devised a blockdiagonal regularization matrix for the internal magnetic field model where n and m run over the degree and order in the spherical harmonic expansion of the internal field in Eq. (1); w m (m) and w tp (n) are functions which control the regularization strength based on the degree and order of the internal Gauss coefficients; λ t , λ ts , and λ te are parameters that, respectively, set the regularization strength over the entire model time span, at the model start time and end time. Following Finlay et al. (2020), in order to relax the regularization at higher spherical harmonic degree, we defined w tp (n) as a tapered window which gradually reduces from one to 0.005 where n min = 3 and n max = 6 are the chosen limits of a half-cosine taper τ (n) = 0.995 2 1 + cos π n − n min n max − n min + 0.005.
In contrast to Finlay et al. (2020), who used n max = 11 to achieve stable power spectra with more power in the time-dependence of the high-degree coefficients without causing instabilities, we were able to further decrease the upper limit of the taper. The magnetospheric and ionospheric field and their induced counterparts may also cause the estimation of the internal field parameters to become unstable. Our experience shows that it is typically the zonal harmonics that become unstable first if the regularization is not sufficiently strong. Therefore, in addition to the degree-dependent temporal regularization, there is a special treatment of zonal and non-zonal spherical harmonics based on Note that the regularization of the internal field model only constrains the timederivatives of the field but not the field itself.
Turning to the external part of the model, we regularized only the bin-to-bin variability of the three RC baseline corrections ∆q 0 1,SM , ∆q 1 1,SM , and ∆s 1 1,SM in Eq. (3) using a quadratic form in the first forward difference of neighboring bins. The forward difference was calculated with the matrix whose number of columns is equal to the number of bins that comprise each RCbaseline correction. Taken together, the regularization matrix for all parameters related to the external field model reads where ⊗ is the Kronecker product, I 3 is the unit matrix of size three corresponding to the three RC-baseline corrections, D 2 = D T D is the coefficients matrix that determines the quadratic form, additional zeros on the diagonal indicate the other unregularized model parameters of the external field, and λ ext is the chosen regularization parameter.
Turning to the calibration parameters, we regularized a quadratic form in the bin-to-bin variability of each calibration parameter for the five platform magnetometers (three on CryoSat-2 and one on each of the two GRACE satellites). The regularization matrix Λ cal is block-diagonal with each block Λ cal,i , i = 1, . . . , 5, corresponding to the calibration parameters for each of the five platform magnetometers. The regularization matrix can be written as where we define the regularization parameters λ b,i , λ s,i and λ u,i to control the temporal smoothness of the offsets, sensitivities, and non-orthogonalities, respectively.

Results and discussion
We built two geomagnetic field models which span 10 years from the 1st of January 2008 to the 31st of December 2018, but differ in the use of platform magnetometer data to constrain the field model parameters.
The first model, Model-A, was derived with data from the Swarm-A, Swarm-B, and CHAMP satellites, and the monthly SV data from ground observatories. It served as a reference model, which allowed us to identify differences to models which were derived using platform magnetometer data in addition. Considering the model parameterization, regularization, and estimation, Model-A is very similar to the CHAOS model series. In fact, the parameterization of the geomagnetic field and the alignment parameters of the satellite data are identical, except for the lower truncation degree of the internal field and the longer bins of the alignment parameters and RC-baseline corrections in Model-A. A notable difference is the use of gradient data in the CHAOS model. The strong temporal regularization of the high-degree Gauss coefficients of the time-dependent internal field has been relaxed in the newly released CHAOS-7 model through a taper, which we also used here. For Model-A, we tuned the regularization, such that the model parameters matched the ones of the CHAOS-6-x9 model as close as possible. Tab. 4 shows the numerical values of the regularization parameters.
The second model, Model-B, is our preferred model and was derived with data from Swarm-A, Swarm-B, CHAMP, monthly ground observatory SV data, and, as opposed to Model-A, platform magnetometer data from CryoSat-2 FGM1, CryoSat-2 FGM2, CryoSat-2 FGM3, GRACE-A, and GRACE-B. In addition to Model-A and Model-B, we built test models in a series of experiments to investigate the effect of platform magnetometer data on the estimation of the geomagnetic field model. Details of the test models are given below. The regularization parameters are the same for all the presented models, i.e., Model-A, Model-B, and the test models. Table 4 Chosen numerical values of the regularization parameters. The values are valid for all the models built in this paper insofar as the regularization terms are applicable to the specific model. 1 Not applicable to Model-A, which was not derived from platform magnetometer data.

Fit to satellite data and ground observatory SV data
We begin with reporting on the fit of Model-B to the satellite data and ground observatory SV data. The histograms of the scalar and vector residuals for each dataset are shown in Fig. 3. The residuals of Swarm-A, Swarm-B, CHAMP and the ground observatories show narrow and near-zero centered peaks, which demonstrate the high-quality and low-noise level of these datasets. In contrast, the peaks are broader for CryoSat-2 and even more in the case of GRACE, which is, as expected, due to the higher data noise level. By separating the residuals poleward of ±55°Q D latitude from the ones equatorward, we find that peaks are broader at polar QD latitudes for all datasets, which is a result of unmodeled magnetic signal of the polar ionospheric current system. Also, the histograms of the GRACE residuals are biased toward negative values. Upon further investigation, we found a local time-dependence especially visible in the scalar residuals, which could indicate that signals from solar array and battery currents have not been fully removed from the GRACE datasets used here. The residual statistics are summarized in Tab. 5 for the satellite data and Tab. 6 for the ground observatory SV data. Fig. 4 shows the time-series of the SV components at six chosen ground observatories together with the computed values from Model-A and Model-B. Overall, the fit of Model-A and Model-B to the ground observatory SV data is good, as expected, for the first five observatory SV shown since these data were used in the model parameter estimation. The computed values of Model-A and Model-B differ especially during the gap from 2010 to 2014, where Model-B can make use of platform magnetometer data in addition to the ground observatory SV data, while Model-A only relies on the ground observatories. That shows that platform magnetometer data contribute to the internal field model especially when there is a lack of calibrated satellite data from CHAMP and Swarm. Perhaps even more convincing is the performance of both models when compared to a dataset not used in the inversion. With the SV data from Saint Helena, we show such an independent dataset in the last row of Fig. 4. Although both models fit Saint Helena well, Model-B performs slightly better in the radial SV in 2013 and the azimuthal SV at least in the first half of the gap period, until 2012.   Figure 3 Histograms of the residuals of each satellite and ground observatory SV data using Model-B. The histograms have been normalized to have unit area. Computed statistics are shown in Tab. 5 for the satellite data and Tab. 6 for the ground observatory SV data. Table 7 Mean values of the calibration parameters for each platform magnetometer dataset. The time-series are shown in Fig. 5. To summarize, with Model-B we built a model that fits both the satellite and ground observatory SV data to a satisfactory level, which shows that platform magnetometer data can be successfully used in geomagnetic field modeling.

Calibration parameters
We document the estimated calibration parameters of each platform magnetometer dataset by showing the time-series in Fig. 5 and the respective mean values in Tab. 7. In Fig. 5, the rows of panels correspond to the CryoSat-2 (top three) and GRACE (bottom two) platform magnetometer datasets, and the columns of panels show the offsets (left), sensitivities (middle), and non-orthogonality angles (right). Since Alken et al. (2020) also used magnetic data from the three platform magnetometers on-board CryoSat-2, it is possible to compare the estimated calibration parameters. First, comparing the time-averaged values of the calibration parameters (Tab. 7 here and Tab. 4 in Alken et al. (2020)), we find that the nonorthogonalities are equal to within 0.01°and the offsets to within 1 eu. The averaged values of sensitivities are equal to within 1 × 10 −4 eu/nT (notice that Alken et al. (2020) use the reciprocal of the sensitivity). In terms of the temporal variability, we find that our estimated calibration parameters have amplitudes that are smaller, or equal in case of the offsets, which is likely due to a difference in the regularization strength. In

Results of the experiments
We conducted a series of experiments in which we changed the model estimation, parameterization, and data selection with the goal to investigate and document difficulties when dealing with platform magnetometer data in a co-estimation scheme. This section also justifies the modeling strategies that went into the construction of our preferred geomagnetic field model, Model-B. In a first experiment, we allowed the nightside platform magnetometer data to participate in the estimation of the axial dipole coefficient of the time-dependent internal field. That is, we derived a test model, Model-C, identical to Model-B but left the matrix of partial derivatives G unchanged so that the entries corresponding to the B-spline coefficients g 0 1,j were non-zero and, thus, the satellite data contributed to the estimation of the internal dipole coefficients. On the left of Fig. 6, we show the time-derivative of g 0 1 as a function of time computed with Model-B and Model-C, while, on the right, we show s 1 of GRACE-A as an example of the calibration parameters. In contrast to Model-B, Model-C features a conspicuous detour of the time-derivative of the g 0 1 coefficient in the gap between CHAMP and Swarm data (blue-shaded region). Although we only show s 1 of GRACE-A in Fig. 6, we find that all three sensitivities of each platform magnetometer differ in the gap period between Model-C and Model-B. The other internal Gauss coefficients also deviate but to a lesser extent. Interestingly, other model parameters such as the offsets, nonorthogonality angles, Euler angles and external field parameters seem qualitatively unaffected. The same correlation between the internal axial dipole coefficient and the sensitivities has been reported by Alken et al. (2020) who show that this effect can be mitigated either by including large amounts of previously calibrated data or through the use of a regularization that favors a linear time-dependence of the internal dipole during the gap period. Due to the lack of additional calibrated data and our interest in the high-degree SA during the gap that such a regularization affects by redistributing power to higher degrees, we chose to set the dependence of g 0 1 , the most affected internal Gauss coefficient, on the satellite platform magnetometer data to zero. In other words, we completely relied on the ground observatory SV data and the temporal regularization to estimate the time-dependence of g 0 1 in the gap period.
In a second experiment, we built a test model, Model-D, which uses 30 day bins of the RC-baseline corrections consistently over the whole model time span in contrast to Model-A and Model-B, which use a single bin spanning the entire gap period. As an example, Fig. 7 shows the RC-baseline correction ∆q 0 1 on the left and the calibration parameter s 1 of GRACE-A on the right, computed with Model-D and Model-B. In Model-D, ∆q 0 1 has a noticeable peak during the gap period that is much larger in value than the variation during CHAMP or Swarm times while the sensitivity is slightly offset to higher values. We find the same behavior for all RC-baseline corrections and calibration parameters, although most prominently for the sensitivities. Again, other model parameters seem unchanged, which indicates that there is a significant correlation between the RC-baseline corrections and the calibration parameters of the platform magnetometers. Using a single bin for the RC-baseline corrections in the gap period helps to reduce that effect. As a final comment regarding Model-C and Model-D, we performed a simulation combining both experiments; that is, we determined g 0 1 with the platform magnetomter data and estimated the RC-baseline corrections in 30 day bin over the entire model time span. In this case, we observed deviations from Model-B which were identical to those shown in Figs. 6 and 7 but, now, affected the internal axial dipole, the RCbaseline corrections, and the sensitivities all at the same time. In an effort to analyze the relationship between the calibration and the other model parameters in a quantitative manner, we also investigated the model correlations ρ ij = C ij / C ii C jj based on the entries of the model covariance matrix evaluated with the converged model parameters (Tarantola, 2005, p. 71). Unfortunately, the analysis revealed a large number of small correlations, which are difficult to interpret. Therefore, we did not make significant use of it in the modeling and preferred to rely on experiments to guide our modeling strategy. In a final experiment, we derived a test model, Model-E, by only using nightside platform magnetometer data as opposed to Model-B, where the calibration parameters were determined from dayside and nightside platform magnetometer data. Fig. 8 shows the calibration parameters for GRACE-A computed with Model-B (thick lines) and Model-E (thin lines). In the case of GRACE-A, using dayside data to determine the calibration parameters considerably changes the sensitivities and non-orthogonalities as can be seen, for example, when looking at s 1 , s 2 or u 3 . In particular for s 2 , there is a vertical shift of approximately 200 × 10 −6 eu/nT, which translates to 10 nT in a magnetic field of 50,000 nT. Irrespective of the platform magnetometer, the experiment shows that the local time coverage of the data plays an important role in determining the calibration parameters. The importance of using both day and nightside data becomes clear when appreciating that the orbital plane of the satellites is slowly drifting in local time. Under a possible nightside data selection criteria, the drift leads to the selection of data from either the ascending or descending parts of the orbit at a time. For example, if the ascending node of the orbit is on the nightside, then the platform magnetometer collects data of the magnetic field that mostly points along the direction of flight, in agreement with the predominant dipolar field configuration, until the ascending node crosses over to the dayside placing the descending part of the orbit on the nightside. Now, the observed magnetic field mostly points against the direction of flight. In the case of CryoSat-2, it takes the ascending node 8 months and GRACE around 11 months to traverse the nightside, which is longer than the monthly bins used for estimating the calibration parameters. Hence, the data of each bin will be collected either from the ascending or descending nodes with the respective bias of the field direction. Instead, by using both nightside and dayside, we ensured that the data within each bin covered a broad range of local times to excite the platform magnetometer from various directions, which we believe improves the estimation of the calibration parameters. Nevertheless, we did not use any dayside data to constrain the geomagnetic field model since we do not account for the strong ionospheric sources on the dayside. Those ionospheric sources, however, may contaminate the calibration parameters.

Secular acceleration
One motivation for using platform magnetometer data has been the growing interest in SA pulses, enhancements of the SA that occur on sub-decadal time scales and are seen most prominently at low latitudes. These pulses have been reported by several studies (Chulliat and Maus, 2014;Olsen and Mandea, 2007) and are thought to reflect the dynamical processes in the Earth's outer core.
To further study SA pulses and the SA in general, accurate internal field models are needed, which rely on long and continuous time-series of satellite data to give a global picture. When supplemented with high quality satellite data, platform magnetometer data may play an important role in providing those models.
To investigate the effect of platform magnetometer data on the recovered SA, we show in Fig. 9 time-longitude maps of the radial SA on the Equator at the CMB computed with Model-B (left) and Model-A (center) alongside the difference map (right). Recall that Model-B is partly based on platform magnetometer data in contrast to Model-A, so that the difference of the two reflects the use of these data. Both models show the SA pulses in 2009, 2013 and most recently in 2017 as enhancement of the radial SA on the Equator. Of special interest is the pulse in 2013, right in between periods of high-quality magnetic data from the CHAMP and Swarm missions. In the difference map, the SA during CHAMP and Swarm period is largely unchanged, which suggests that the effect of the CryoSat-2 and GRACE data is rather minimal during these times. In contrast, the SA in the gap period is distinctly different for the two models. Differences that are large in absolute value seem to be concentrated around 0°and 180°longitude on the Equator which coincides with the Pacific and the region in the South Atlantic close to Central Africa. The geographical location of the differences is more clearly seen in Fig. 10, which shows global maps of the radial SA at the CMB during the SA pulses in 2009, 2013 and 2017. Again, the difference between Model-B and Model-A is small in 2009 and 2017, i.e. during CHAMP and Swarm times, but large in 2013 in the middle of the gap period. The regions with the largest differences are located in the Southern hemisphere and the Equatorial region with prominent examples in the West and South Pacific Ocean, and Central Africa. Our findings seem to indicate that the platform magnetometers have the desired effect of balancing the uneven spatial distribution of the ground observatory network in the gap period.

Conclusions
In this study, we present a co-estimation scheme within the framework of the CHAOS field model series that is capable of estimating both a geomagnetic field model and, at the same time, calibration parameters for platform magnetometers. This approach enables us to use platform magnetometer data to supplement highquality magnetic data from magnetic survey satellites and removes the requirement for utilizing a-priori geomagnetic field models to calibrate platform magnetometer data.
We followed Alken et al. (2020) but went further in that we co-estimated a model of not only the internal field but also the external field. The co-estimation scheme relies on absolute magnetic data which we took from CHAMP, Swarm-A, Swarm-B and the monthly SV data from ground observatories between 2008 and 2018. Magnetic data from five platform magnetometers were used: three on-board CryoSat-2 and one on-board each of the GRACE satellite pair. This allowed us to considerably improve the geographical and temporal coverage of satellite data after CHAMP and before the launch of the Swarm satellites.
We successfully co-estimated a geomagnetic field model along with calibration parameters of the five platform magnetometers. The misfit to the high-quality satellite data and ground observatory SV data was similar to that for models derived without including platform magnetometer data, and the good fit to an independent ground observatory dataset from Saint Helena provide evidence that our modeling approach performs well.
In a series of experiments we investigated the trade-offs when co-estimating calibration and geomagnetic field model parameters. We found that the calibration parameters strongly correlate with the internal axial dipole and the RC-baseline corrections of the external field during the gap period, when there is less high-quality data available. By preventing platform magnetometer data from contributing to the internal axial dipole and using constant RC-baseline corrections throughout the entire gap period, we successfully avoided those complications.
Our experiments showed that including platform magnetometer data leaves the SA signal practically unchanged during the CHAMP and Swarm period but leads to differences in the gap period. The difference in the recovered SA signal is stronger in the West and South Pacific, where only a few observatories are located, which suggests that platform magnetometer data help to improve the global picture of the SA. Based on our investigations, we find that it is worthwhile to include platform magnetometer data in internal field modeling, in particular from CryoSat-2 given the relative low noise level.

Abbreviations
CHAMP CHAllenging Minisatellite Payload CMB Core mantle boundary CRF Common reference frame DMSP Defense Meteorological Satellite Program FGM Fluxgate magnetometer GRACE Gravity Recovery and Climate Experiment GSM Geocentric solar magnetic NEC North-east-center nT NanoTesla QD Quasi-dipole RTP Radius-theta-phi SM Solar magnetic VFM Vector field magnetometer Declarations Availability of data and materials The datasets supporting the conclusions of this article are available in the following repositories: Swarm and CryoSat-2 data are available from https://earth.esa.int/web/guest/swarm/data-access; The GRACE data are available from ftp://ftp.spacecenter.dk/data/magnetic-satellites/GRACE/ CHAMP data are available from https://isdc.gfz-potsdam.de/champ-isdc; Ground observatory data are available from ftp://ftp.nerc-murchison.ac.uk/geomag/Swarm/AUX_OBS/hour/;