Magnetometer data from the GRACE satellite duo

This paper describes and discusses the preprocessing and calibration of the magnetic data taken by the navigational magnetometers onboard the two GRACE satellites, with focus on the almost 10 years period from January 2008 to the end of the GRACE mission in October 2017 for which 1-Hz magnetic data are available. A calibration of the magnetic data is performed by comparing the raw magnetometer sensor readings with model magnetic vector values as provided by the CHAOS-7 geomagnetic field model for the time and position of the GRACE data. The presented approach also accounts for magnetic disturbances produced by the satellite’s magnetorquer and for temperature effects, which are parametrized by the Sun incident angle. The root-mean-squared error of the difference between the calibrated data and CHAOS-7 model values is about 10 nT, which makes the GRACE magnetometer data relevant for geophysical investigations.


Introduction
Most satellites are equipped with magnetometers as part of their attitude control system. The measurements taken by these platform magnetometers are not only useful for navigational purposes, but may also be used for scientific investigations, provided they have been properly calibrated. In recent years there has been an increased focus on the utilization of platform magnetometer data, starting from the pioneering work done by Brian Anderson and his team in calibrating and interpreting magnetometer data from the IRIDIUM satellites (e.g. Anderson et al. 2000) to recent work done with CryoSat-2 (e.g. Olsen et al. 2020) and DMSP (e.g. Alken et al. 2020). A study on "Feasibility of using platform magnetometers to observe and detect Space Weather events" (Doornbos et al. 2018) is another example of recent activities regarding scientific utilization of data from non-scientific satellite magnetometers.
Proper calibration of platform magnetometer data is an endeavour that needs to be adapted to the specific satellite mission. There is no "one-size-fits-all" solution. A procedure that is successful for one mission is not necessarily the optimal approach for other satellites. In particular, the correction for magnetic disturbances due to satellite and payload operation-for instance how to account for magnetic signatures from magnetorquer activity and solar array currents-usually needs to be adapted for the specific satellite in consideration.
This paper describes the preprocessing and calibration of platform magnetometer data collected by the GRACE (Gravity Recovery and Climate Experiment) satellite duo (Tapley et al. 2004). The two GRACE satellites were launched on 17 March 2002 into a near-polar orbit with altitudes of about 470 km, decreasing in particular after 2013 to about 330 km and lower towards the end of the mission in October 2017 (cf. Fig. 1). Orbital inclination of 89.0 • implies a Local Time (LT) drift rate of the ascending node (LTAN) of 2.25 h/month; all local times are thus covered every 5.3 months. The two satellites, denoted as GRACE-A and GRACE-B, follow each other with an along-track separation of approximately 220 km.
The main scientific objective of the GRACE tandem is to measure gravity anomalies with very high precision (Tapley et al. 2004), enabling to determine changes in distribution of masses around the planet. Data from the GRACE satellites are therefore useful for studying Earth's geology, hydrology, ocean, and climate as well as processes in Earth's thermosphere and ionosphere. In addition, each of the two GRACE satellites is equipped with a magnetometer used for navigational purposes. In this paper, we present our effort to calibrate the magnetic data collected by the two GRACE satellites and to provide calibrated magnetic data as daily files of 1 Hz temporal resolution.
Similar to GRACE, the GRACE Follow On (GRACE-FO) mission also consists of two identical satellites, launched in May 2018. But these satellites carry improved magnetometers compared to those flying on GRACE. In particular, their 16-bit digitization results in less noisy magnetometer data, as demonstrated by Stolle et al. (2021). Since GRACE-FO flies simultaneously with the high-precision Swarm satellites, magnetic data from GRACE-FO are less interesting for investigations of Earth's internal magnetic field but of more value for studying the dynamics of external field variations like field-aligned currents (FACs). However, the unique value of GRACE magnetic measurements, despite larger instrument noise, is the fact that they bridge the data gap between CHAMP and Swarm, with a solid margin of more than 2 years on each side of this gap. This makes GRACE magnetic data along with similar platform magnetometer data from the CryoSat-2 satellite interesting for modelling Earth's core field, as shown, e.g. by Kloss et al. (2021).

The GRACE magnetometer measurements
Each of the two GRACE satellites carries a 3-axis fluxgate magnetometer, mounted at the tip of a short (about 0.5 m long) boom ( Fig. 2) located below the satellite pointing in nadir direction. These magnetometers are of the same type as that flown on the ROSAT satellite (Trümper 1982) and sample the magnetic field with 12-bit resolution. This coarse solution (most other satellite are equipped with 16 bit magnetometers) implies that the vertical magnetic field (roughly aligned with sensor axis 3) measures the magnetic field in 2 12 steps separated by d = 25.88 nT between − 53, 000 nT +d = − 52, 974 nT and + 53, 000 nT. In Antarctic regions close to the magnetic pole where the vertical magnetic field exceeds − 52, 974 nT, the magnetic data are saturated and have been flagged and discarded from the further analysis. The 12-bit digitization of the GRACE magnetometer data result in increased noise when analysing the data: the quantization noise produced by a N = 12-bit digitizer sampling in steps of size d results in a root-mean-squared (rms) noise (Schwartz 1980) of rms = d/ √ (N ) = 25.88 nT/ √ (12) ≈ 7.5 nT. From the launch of the satellites in 2002 until the end of 2007 the GRACE magnetometer data are available at 10-s sampling rate (1-s data are available for shorter time periods on selected days), while from 1 January 2008 to end of mission in October 2017 the instruments provided data at 1-s sampling rate. In this article, we focus on the 10-year period 2008 to 2017 for which 1-s data are available (cf. Fig. 1).
To illustrate some characteristics of the GRACE magnetometer data, including challenges and opportunities, we show in Fig. 3 data collected by GRACE-A on 2 February 2014. The top panel of the figure shows the raw (uncalibrated) readings B 3 of the 3rd magnetometer sensor for the whole day as given in the Level-1b data files. Saturation of the magnetometer at values below − 52, 974 nT (highlighted in red) is obvious between 02:00 and 05:00 UTC, and again between 14:00 and 17:00 UTC. Also seen are magnetic field variations of a few thousands of nanotesla, for instance during the few minutes before 12:00 UTC.
A 20-min time segment from 11:30 to 11:50 UTC (grey region in upper panel) is presented in the second panel of the figure. The observed magnetic field B 3 (blue curve) shows fluctuations of several thousands of nanotesla. Since the third magnetometer axis is roughly aligned with nadir, B 3 should approximately sample the radial magnetic field component B r . Model values of B r as given by the CHAOS-7 magnetic field model (Finlay et al. 2020) are shown for comparison. The observed field B 3 (blue curve) follows the general trend of the model values (red curve) but is shifted by about 2000 nT, which hints at a magnetometer bias of roughly 2000 nT in that component. But the scatter of the B 3 magnetometer readings of several thousands of nanotesla is likely of nongeophysical origin. A more detailed analysis revealed that GRACE magnetometer data are heavily disturbed by magnetorquer action.
Each GRACE satellite is equipped with three magnetorquers (MTQs) for attitude control. See Fig. 2 for the location of the MTQs. These magnetorquers produce magnetic fields of up to several thousands of nanotesla at the location of the magnetometers, and correcting for these magnetic disturbances is essential before GRACE magnetometer data can be used for geophysical investigations. Fortunately, MTQ current values are provided in the Level-1b data for the time instant of each magnetometer reading, and thus these disturbances can be accurately modelled and removed from the measurements, despite their magnitude.
Following Olsen et al. (2020), we model the magnetic disturbance vector B MTQ at the magnetometer produced by magnetorquer currents I MTQ (separate currents for each of the three magnetorquer coils) as: where M is a 3 × 3 matrix that describes the magnetic disturbance at the magnetometer location caused by I MTQ . As discussed below, we co-estimate the 9 elements of M as part of the calibration effect and obtain (1a)  Table 2 of Olsen et al. (2020) is incorrectly given as nT/A. The correct unit is nT/mA.) Typical GRACE magnetorquer current strengths have values of up to 100 mA which results in magnetic field disturbances of several thousand nanotesla. Accounting for magnetorquer disturbance is therefore crucial for GRACE. The third panel of Fig. 3 shows B MTQ,3 for the same 20-min time segment as in the second panel; the excellent agreement with the observed scatter of B 3 (blue curve of second panel) is obvious. Removing the modelled MTQ magnetic disturbances from the magnetometer observations leads to the smooth orange curve shown in the second panel of the figure. The good agreement with the CHAOS-7 model value of B r -apart from the offset of about 2000 nT-is encouraging for using GRACE magnetometer data for geophysical purposes.
The bottom figure panel presents the fully calibrated and corrected downward magnetic field component B C after removal of the CHAOS-7 model values. Most of the remaining scatter (of a few tens of nanotesla peakto-peak) in these 1-Hz data is probably due to the 26-nT digitization steps of the magnetometer, while some of the spikes can be attributed to strong magnetorquer action: current strengths larger than 40 mA (highlighted in orange) occur, on average, during less than 10% of the time but are more frequent in specific geographic regions. Some of the spikes can be attributed to magnetorquer current saturation (currents > 119 mA, highlighted in red). To reduce the scatter further and remove remaining spikes, we applied a 11-s median filter; the resulting data, shown by the green curve, are those provided in the final calibrated data product (see Table 3 below).
One of the challenges of modelling Earth's magnetic core field is the 3-year gap in satellite data coverage between the atmospheric re-entry of the CHAMP satellite in September 2010 and the launch of the Swarm satellite trio in November 2013. Although platform magnetometer data from CryoSat-2 (available at 4 s sampling rate since August 2010, Olsen et al. 2020) can close this gap, Cryo-Sat-2 magnetic data are of reduced quality during the first years of the mission before an in-flight software update allows for more stable-and thus less disturbed-satellite operation. GRACE platform magnetometer data provide an interesting alternative due to several reasons: firstly, the 3-year gap in high-precision satellite magnetometer coverage happened during the second half of the GRACE mission and thus well within its operational phase; secondly, GRACE provides 1-s magnetometer data (and thus at higher sampling rate than available for CryoSat-2); and thirdly, the availability of data from two spacecraft following each other with a distance of 220 km (and thus performing magnetic field measurements at about the same location but separated by 30 s) gives additional opportunities for geophysical studies.
This paper describes details of the preprocessing and calibration steps applied to the GRACE magnetometer data. Information on the GRACE data used in this study and how they have been preprocessed is provided in "Data and data preprocessing" section. "Calibration of magnetometer data" section describes the mathematical model used to calibrate the magnetometer data, including how to account for spacecraft disturbance fields and non-linear magnetometer effects. Results are discussed in "Results and discussion" section. The paper concludes with a summary and outlook in "Summary and outlook" section.
The calibration effort described here is performed using an a priori model of Earth's magnetic field-for that we used the CHAOS-7 model of Finlay et al. (2020). Alternatively, a geomagnetic field model can be co-estimated with (at least some of ) the magnetometer calibration instead of using an a priori field model. Alken et al. (2020) and Kloss et al. (2021) describe such an approach, and Kloss et al. (2021) applied it to GRACE magnetic data after correcting for spacecraft disturbance fields and magnetometer non-linearities using (a preliminary version of ) the "common" calibration parameters determined in the present paper.

Data and data preprocessing
GRACE data spanning the 118 months between January 2008 and October 2017 (for GRACE-A), and September 2017 (for GRACE-B), respectively, have been analysed in this study. We use raw magnetometer sensor readings E and magnetorquer currents I MTQ as given in the GRACE MAG1B L1b data files, sampled at a frequency of 1 Hz. Satellite positions are taken from the GNV1B files, sampled at 5 s. Spacecraft attitude is available also every 5 s in the SCA1B files. All these data files are publicly available at ftp://isdcftp.gfz-potsdam.de/grace/ Level-1B/JPL/. Position and attitude data have been spline-interpolated to the 1-s time instants of the magnetometer readings. If one of the 3 magnetometer sensor readings was outside of ± 52, 974 nT, indicating saturation of the sensor, the corresponding values of all 3 components have been flagged and excluded from the further analysis. Time series of electric currents of the two solar arrays ( I SA1 , I SA2 ) and of battery ( I Batt ), sampled at 3 s, have been kindly provided by Krzysztof Snopek (GFZ Potsdam).
For the estimation of the magnetometer calibration parameters, we use measurements downsampled to 60 s. Data from all local times and geomagnetic activity conditions were included. However, since the CHAOS-7 model that we use as reference model field does not attempt to describe the strong ionospheric magnetic field signatures at polar latitudes, we only use data at quasi-dipole (QD) latitudes (Emmert et al. 2010) equatorward of 60 • for estimating the calibration parameters (which, however, subsequently are applied to data from all latitudes and with full 1 Hz time resolution). For GRACE-A this yields 10.5 mio data points (corresponding to 3 × 2.62 mio vector triplets plus 2.62 mio intensity data) for estimating the magnetometer calibration model; for GRACE-B the corresponding numbers are 9.56 mio data points (2.39 mio vector triplet splus 2.39 mio intensity data). Following Olsen et al. (2003) and Olsen et al. (2020), the relationship between the raw magnetometer sensor readings E and the calibrated and aligned magnetic field vector B NEC in the NEC (North-East-Center) frame is given by where b is the vector of 3 magnetometer offsets, S −1 is the inverse of the diagonal matrix S of 3 scale values; P −1 is the inverse of the lower triangular matrix P of 3 (2a)

Calibration of magnetometer data
non-orthogonalities; B CRF is the magnetic field vector in the "Common Reference Frame" (CRF) of the satellite's star tracker (STR); the rotation matrix R A describes the rotation between the magnetometer (FGM) frame and CRF; the rotation matrix R J2000 contains the attitude information provided by the star tracker; and R NEC describes the rotation from the J2000 to the NEC frame. The matrices R NEC and R J2000 are assumed to be known ( R NEC is determined using IERS algorithms (Vallado et al. 2006) and R J2000 is derived from the STR attitude) and thus only the 9 elements of the 3 × 3 matrix A = R A P −1 S −1 (which involves the 3 Euler angles describing R A , the 3 non-orthogonality parameters describing matrix P −1 and the 3 scale values describing matrix S −1 ) plus the 3 elements of b = −Ab , that means in total 12 parameters, have to be determined in order to calibrate the magnetometer data for a given time interval. The 12 "classical" calibration parameters (3 offsets, 3 scale values, 3 non-orthogonalities and 3 Euler angles) can be determined from these 12 model parameters via QL decomposition of A . See Olsen et al. (2020) for details of this parametrization.
Knowledge of B CRF (e.g. as given by a model of Earth's magnetic field) plus observations of the magnetometer sensor readings E allows for estimating the 9 elements of A and the 3 elements of b as a linear inverse problem. This approach has been used for calibrating the Cryo-Sat-2 magnetic data ). However, for calibrating GRACE magnetometer data we instead chose to directly estimate the 3 offset values, 3 scale values, 3 nonorthogonalities and 3 Euler angles by solving a non-linear inverse problem. Although slightly more complicated, this approach allows us to apply separate temporal regularizations of the month-to-month variability of the Euler angles and of the other 9 parameters, enabling us to better account for boom bending, i.e. a mechanical variation of the coordinate frames of the magnetometer (at the tip of the boom) and of the star imager (fixed to the spacecraft body).

Accounting for time dependence, non-linearities, sensor cross-talk and s/c disturbance fields
To account for the magnetic disturbance of magnetorquer activity, solar arrays, batteries as well as of nonlinearities in sensor output and cross-talk between the three sensors we augment the linear equation of Eq. 2c as where the 3 × 3 matrix M (i.e. 9 parameters) describes the magnetic disturbance caused by the magnetorquer coil currents I MTQ (separate currents for each of the three magnetorquer coils); b SA1 , b SA2 and b Batt (i.e. 3 × 3 = 9 parameters) describe the magnetic disturbances due to the electric currents I SA1 , I SA2 and I Batt of the two solar arrays and the batteries, respectively; the vectors ξ and η describe quadratic, and, respectively, cubic, effects in the sensor output E , and b ADC is a 3 × 3 diagonal matrix with elements b k,ADC , k = 1 − 3 , to account for zero-offsets in the analog-digital converter (ADC) of the magnetometer.
The ith components ( i = 1 − 3 ) of ξ and η are parametrised according to where we have scaled the sensor output E k according to Ê k = E k /E 0 with E 0 = 10 4 nT in order to achieve nondimensional quantities Ê k of order O(1), thereby reducing numerical instabilities in estimating these parameters.
To account for a possible zero-offset in the Analog-Digital-Converter (ADC) of the magnetometer we co-estimate a bias correction b k,ADC for all three magnetometer sensors, resulting in additional 3 model parameters.

Accounting for sun-related magnetic disturbance fields
The calibration parameters of fluxgate magnetometers are known to depend on the temperature of sensor and electronics, an effect that is accounted for, e.g. in the calibration of Swarm (Tøffner-Clausen et al. 2016) and Cryo-Sat-2 ) magnetometer data for which measurements of magnetometer temperatures are available. However, such temperature data are unfortunately not available for GRACE. Not accounting for the temperature dependence of the magnetometer readings results in magnetic data that show a pronounced variation with the Local Time of the satellite orbit, likely due to the heating caused by Sun illumination. To mitigate the lack  Fig. 4, as proxy for magnetometer temperature. The "azimuth" angle α varies rapidly within one orbit (i.e. within about 1.5 h), while the "elevation" β varies slowly with orbital LTAN. Figure 5 shows the time evolution of the slowly changing Sun incident angle β for both GRACE satellites. Their nominal flight configuration (since December 2005) as "tandem" means that the lead satellite GRACE-B is flying in +x-direction, while the trailing GRACE-B satellite flies in −x-direction. The β angle of the two satellites therefore follow each other closely, but have opposite sign: if GRACE-A is illuminated by the Sun "from left" ( −y-direction, i.e. β ≈ 90 • ) the corresponding β angle for GRACE-B is β ≈ −90 • corresponding to illumination "from right" ( +y-direction). After 2014, however, the two satellites "swap order", and correspondingly flight direction, during four periods of a few months duration each. These periods are indicated by grey-shaded areas in Fig. 5. As a consequence, GRACE-A is always illuminated "from left" ( β > 0 ) after mid-2014, while GRACE-B is always illuminated "from right" ( β < 0).
In case of the Swarm satellite trio, a magnetic disturbance field depending on the Sun incident angles α and β was found, and thermoelectric currents in the groundings of the blankets covering the magnetometers have been identified as the likely cause of these disturbances. These disturbances have been successfully modelled by co-estimating offset values depending on α and β together with the other magnetometer calibration parameters-see Tøffner-Clausen et al. (2016) for details.
In case of the GRACE satellites the situation is more complicated since α and β are used as proxy for instrument temperature, which not only affects sensor offsets, but also scale values (the impact of temperature changes on non-orthogonalities is found to be small). In addition, a temperature-related bending of the magnetometer boom will result in Euler angles that depend on the Sun incident angles α and β.
We therefore include an explicit dependence of the three offsets b = b 1 , b 2 , b 3 , of the three scale values S 1 , S 2 , S 3 , and of the three Euler angles ǫ 1 , ǫ 2 , ǫ 3 on Sun incident angles α, β by expanding each of these 9 parameters in terms of spherical harmonics: where x is any of the 9 calibration parameters offsets, scale values and Euler angles, P m n are associated Schmidtnormalized Legendre functions with N = 8 as the maximum degree and M = 2 as the maximum order of the spherical harmonic expansion. Since co-estimating terms of high spherical harmonic order m may absorb small-scale along-track geophysical magnetic field variations (like the Equatorial Electrojet on the dayside at the magnetic equator) we limit our expansion to rather low degrees of m ≤ M = 2 , and x m,c n and x m,s n are expansion coefficients corresponding to cos(mφ)P m n (sin β) and sin(mφ)P m n (sin β) , respectively. This expansion yields (5) x m,c n cos mα + x m,s n sin mα P m n (sin β), 38 parameters for each of the 9 calibration parameters, resulting in 9 × 38 = 342 additional model parameters in total. Preliminary calibration efforts revealed rather large changes, in particular of the Euler angles for months when the GRACE satellites fly in non-nominal, "reversed", configuration (indicated by the grey shadow time periods in Fig. 5 and following), compared to nominal flight direction. The reason for this is presently unknown. To account for this effect we solve for 12 "reverse flying" correction parameters, describing corrections of the offsets, scale values, non-orthogonalities and Euler angles for time periods of "reverse flying satellites", assumed to be the same for those four periods.
These 12 "reverse flying" parameters together with the 9 parameters describing magnetorquer disturbances (matrix M ), the 3 × 3 = 9 model parameters b SA1 , b SA2 and b Batt , the 18 + 30 = 48 parameters describing quadratic and cubic magnetometer sensor effects and cross-talks, the 3 ADC offset corrections, and the 342 parameters describing dependency on α, β result in a total of additional 9 + 9 + 48 + 3 + 342 + 12 = 423 model parameters. We refer to them as "common parameters" since they are assumed to be constant for the whole data period.

Model estimation procedure
Temporal variation of the 12 "classical" calibration parameters is accounted for by estimating them in monthly bins. Since we use GRACE data spanning 118 months (January 2008 to October 2017) for GRACE-A (and 1 month less for GRACE-B), i.e. for time instants t k , k = 1 − K with K = 118 , this yields 118 × 12 = 1416 model parameters (1404 for GRACE-B).
Adding the 423 "common parameters" which are assumed to be constant for the entire time period results in a total of 1839 model parameters for GRACE-A (1827 for GRACE-B).
These model parameters are estimated by performing a combined vector and scalar calibration, minimizing the squared difference | B| 2 + w F | F | 2 averaged over all data points, with where the calibrated magnetic field vector B CRF depends on the sensor output E and on the model parameters as indicated by Eq. 3 and w F is the relative weight given to the scalar residuals | F | 2 compared to the vector residuals | B| 2 . The reference (model) values B mod,CRF are given for each satellite position and time instant by the CHAOS-7 field model (Finlay et al. 2020); they consist of contributions from the core, lithosphere and large-scale magnetosphere (plus Earth-induced counterparts) and are rotated from NEC to the CRF frame: Note that a scalar calibration (i.e. minimizing the average of Eq. 6b) constrains all model parameters but the Euler angles, while a vector calibration (i.e. minimizing the average of Eq. 6a) provides information regarding all model parameters. Although this at first glance seems to be advantageous, mechanical bending of the GRACE magnetometer boom on time-scales not parametrized in the model (either explicitly due to the monthly bins or indirectly due to dependence of Euler angles on Sun incident angles α and β ) might degrade the estimation of all the other calibration parameters when performing a pure vector calibration. A combined vector and scalar calibration, with an empirically determined scalar weight of w F = 5 , mitigates this problem and has been found to be suitable for calibrating the GRACE magnetometer data. We use an Iteratively Reweighted Least-Squares approach with Huber weights (Huber 1981), minimizing the Chisquared misfit where the model vector m contains the 1839 model parameters and the data residual vector e = d obs − d mod is the difference between the data vector d obs (containing all elements of B CRF as well as F = |B CRF | ) and the model prediction vector d mod (containing all elements of B mod,CRF and F mod ). The diagonal weight matrix W contains the Huber weights (to account for data outliers) and is a block diagonal matrix which regularizes the month-to-month variation (first differences) of the 12 "basic" model parameters offsets, scale values, non-orthogonalities and Euler angles. We use different regularization parameters to constrain the temporal variation of the offsets ( b = 1 · 10 5 (month/nT) 2 ), scale values ( S = 1 · 10 14 month 2 ), non-orthogonalities ( u = 1 · 10 14 (month/rad) 2 ) and Euler angles ( ǫ = 1 · 10 13 (month/rad) 2 ) due to their different magnitude and expected temporal variability; for instance we allow for larger month-to-month variations of the Euler angles to account for a likely boom bending, compared to the non-orthogonalities which are expected to have less temporal variability. The values of these regularization (8) χ 2 = e T W e + m T �m, parameters were found by visual inspection of the month-to-month variations of the estimated calibration parameters and the achieved rms data misfit.

Results and discussion
The left part of Table 1 lists Huber-weighted rms data misfits (considering all three magnetic vector components and magnetic field intensity) applied to the full data set spanning January 2008 to October 2017, while its middle and right part shows statistics separately for dayside and nightside part of the data.
Achieved rms misfit values when considering all data are generally slightly larger than 10 nT. Consistently for both GRACE satellites, the misfits of the scalar residuals are lowest while those of the East components are highest. Looking at statistics separately for dayside and nightside data (middle and right part of the table) reveals consistently lower rms data misfit during night. The cause for this could be contributions from unmodelled geophysical sources (e.g. from dayside ionospheric magnetic contributions which are not described by the CHAOS-7 model), larger magnetic disturbances when the spacecraft is illuminated by the Sun, and/or magnetometer temperature effects that are not properly described by the Sun incident angles α, β.
The achieved rms misfit of about 10 nT is higher than corresponding values found for magnetic data of the CryoSat-2 (of about 7 nT, see Olsen et al. 2020) and GRACE-FO (5 to 8 nT, see Stolle et al. 2021) satellites. A possible reason for this is lack of magnetometer temperature measurements for GRACE and the coarser digitization (12 bits compared to 16 bits for CryoSat-2 and GRACE-FO). For comparison, the CHAOS-7 model rms misfit statistics to data that were used in deriving that model are about 5 nT for CryoSat-2 and about 2 nT for Swarm (see Tables 6 and 7 of Finlay et al. 2020). Figure 6 shows the temporal variation of the 12 "basic" calibration parameters (3 offsets, 3 scale values, 3 nonorthogonalities and 3 Euler angles) in terms of the deviation from their mean values, for GRACE-A (top) and GRACE-B (bottom), respectively. The thin curves, presenting results for a calibration model that does not include the 12 correction parameters for "reverse flying satellites", reveal larger deviations during the periods (grey areas) of reverse flying satellites, in particular in the second Euler angle ǫ 2 and the non-orthogonality parameter u 2 . These fluctuations are absent (for GRACE-A), or at least considerably reduced (for GRACE-B) in a calibration model that includes these 12 additional correction Table 1 Huber-weighted mean and rms misfit (in nanotesla) for the various vector components, in the magnetometer frame (B 1 , B 2 , B 3 ) and in the NEC frame (B N , B E   There is a clear secular trend in some calibration parameters, for instance in the offsets b 3 for GRACE-A, b 1 for GRACE-B, and in the second Euler angle for both satellites. Apart from this trend there are variations of a few nanotesla in the offsets, of up to ± 100 ppm in the scale values, and of up to ± 0.01 • (corresponding to about ±36 arcsecs) in the non-orthogonalities and Euler angles. These variations are comparable in magnitude to those found when analysing data from other platform magnetometers like CryoSat-2  and GRACE-FO (Stolle et al. 2021). Unfortunately the variations, e.g. of the Euler angles are quasi-periodic and do neither follow a seasonal nor local time periodicity; attempts to fit annual and local time periods (up to 3rd sub-harmonics) to these variations were not successful. Figure 7 shows maps of the dependencies of the offsets, scale values and Euler angles on Sun inclination angle α, β . The dashed ellipses indicate when the satellite is in eclipse (defined as the Sun being more than 15 • below the x − y-plane, corresponding to zero solar array currents). Despite the differences between the two satellites there are remarkable similarities: most striking is probably the strong increase of scale values for |β| 60 • , corresponding to periods when the magnetometer (at the tip of the boom pointing nadir) is in full solar illumination, either from the left (when β ≈ + 90 • ) or from the right ( β ≈ − 90 • ). Scale values of fluxgate magnetometers are known to depend on temperature; for the highprecision science magnetometers onboard the Ørsted and Swarm satellites they change by up to about 10 ppm per degree change in temperature (for Ørsted, see Olsen et al. 2003) and by up to 30 ppm/ • C (for Swarm, Tøffner-Clausen, personal communication) while values of up to 100 ppm/°C have been found for the navigational magnetometers on CryoSat-2 . CryoSat-2 magnetometer temperatures vary by about °C, resulting in variations of scale values of up to 600 ppm. The GRACE magnetometers are probably more exposed to temperature variations due to their location at the tip of a boom, compared to the CryoSat-2 magnetometers which are located in the satellite body, and thus GRACE magnetometer temperature variations may exceed the 6 °C variations measured on CryoSat-2. This could explain the variations of the GRACE magnetometer scale values of up to 1000 ppm (middle panel of Fig. 7) which are larger than those found on CryoSat-2 but still of same order of magnitude, supporting their plausible origin in magnetometer temperature variations.
In order to investigate the robustness of the estimated common model parameters we derived separate models using data from 2008 to 2012 and 2013 to 2017, respectively. Maps similar to Fig. 7 but determined separately from these two data subsets are included in the Additional file 1: Figs. S12 and S13), and show that the dependence on α, β is robustly determined for most of the parameters-e.g. for the scale values and Euler angles-while other parameters, for instance the offset b 2 , are rather different for the two data subsets, indicating less robust estimates. This behaviour is seen for both GRACE satellites. Table 2 lists values of the other common parameters. Its left part, labelled "2008-2017", shows values estimated from the entire data set while its middle and right panel presents results of the independent estimations based on the two data subsets. In general the results obtained from the two data subsets are rather similar also for these common calibration parameters, which gives confidence that inclusion of these parameters is justified. However, co-estimation of a set of common calibration parameters should be considered as a tool to obtain an improved calibration of the magnetic data-similar to co-estimating of "nuisance parameters" in statistics-and interpretation of individual parameters should be made with care. Nevertheless, the excellent agreement of the estimated values of the matrix M used for correction of magnetorquer disturbances is striking. The differences between their values derived from the two data sets are below 1%, and also the values determined for the two GRACE satellites are rather similar. All this gives confidence in an accurate correction of the rather large magnetorquer disturbances on GRACE.
Also the obtained estimates of the ADC offset corrections b ADC are fairly consistent: that of the 1st sensor axis is about + 3.5 nT for GRACE-A and slightly lower for GRACE-B, while those for the other two axes are approximately zero.
Although we formally solve for 1839 model parameters (for GRACE-A), the number of degrees of freedom is considerably smaller due to the strong temporal regularization of the 118 × 12 = 1416 "basic" model parameters. The trace of the model residual matrix is 819, which means that actually only 45 % of the 1839 model parameters are determined by the data while the remaining 55 % are constrained by model regularization (i.e. the constraints put on the month-to-month variation of the 12 "basic" calibration coefficients). See Tarantola (1987); Sabaka et al. (2004); Alken et al. (2020) for details on the calculation of the model resolution matrix for data subsets. Since model regularization is only applied on the 1416 "basic" model parameters (the common parameters are not regularized) this results mean that effectively only 1416 − 819 = 597 (corresponding to about 42 %) of these "basic" calibration parameters are determined by the data   while the remaining fraction is constrained by model regularization.
While the rms data misfit values listed in the left column of Table 1 are mean values for the entire mission, the upper part of Fig. 8 shows the time evolution of daily mean values of the residuals at non-polar latitudes (QDlatitude equatorward of ± 60 • ), for the three vector components (B N , B E , B C ) and for scalar intensity F. There is no obvious systematic behaviour in the mean values, although there seems to be a tendency for larger deviations from zero during dawn/dusk, indicated by the vertical dashed red and blue lines.
Also daily misfit rms values, presented in the lower part of the figure, do not show obvious systematic temporal variations, perhaps apart from increased rms of B N during the transition between "nominal" and "reverse" flying satellites in 2014 and 2015.
However, a small but systematic Local Time variation in the rms misfit becomes visible when plotted in dependence on LT and separately for ascending and descending orbit parts, as presented for GRACE-A in the right part of Fig. 9. This is consistent with the statistics presented in Table 1 for nightside and dayside data, respectively. Local Time variations are less obvious for the daily mean values, shown in the left part of the figure.

Summary and outlook
We have described and discussed the preprocessing and calibration steps applied to the magnetic data measured by the platform magnetometers onboard the two GRACE satellites. This calibration is performed by comparing the magnetometer sensor readings with magnetic field values as given by the CHAOS-7 geomagnetic field model (Finlay et al. 2020) for the time and position of the satellites. Fully calibrated magnetic vector data, together with time and position, are provided as daily files in the Common Data Format (CDF, see Goucher and Mathews (1994) and http://cdf.gsfc.nasa.gov/) in a format similar to that used for distribution of the Swarm magnetic Level 1b data (e.g. Olsen et al. (2013) and http://www.earth .esa.int/swarm /). The product name follows that used for the Swarm Level 1b data: file GRACE_A_MAG_20140817T00000 0_20140817T235959_0102.cdf contains the magnetic data taken by the GRACE-A satellite on 17. August 2014. The last four digits of the filename-here 0102indicate product version 0102. Table 3 lists the content of this GRACE magnetic data product. In addition to time and position it contains the three components of the calibrated magnetic vector in the Common Reference Frame (CRF) of the satellite as well as in the NEC frame. To reduce the intrinsic noise of the magnetometers we applied a 11-s robust average of the magnetic field vectors; this is provided as variable B NEC . Scalar intensity F = |B NEC | is also included in the data file.
The calibrated GRACE magnetic data product files are freely available at ftp.spacecenter.dk/data/magnetic-satellites/GRACE/. Since the data format follows that used for the Swarm Level 1b data, a visualization of the GRACE magnetic data is possible by uploading the CDF files to the VirES for Swarm web service available at http://www. vires .servi ces.
To illustrate the quality of the calibrated magnetic data, we present in Fig. 10 the difference B = B obs − B mod between the observed magnetic field and CHAOS-7 model predictions in dependence on QD magnetic latitude, separately for dayside and nightside parts of the orbits. In this example, we show GRACE-A data from the 10-day period 10 to 19 August 2010. The magnetic signature of polar currents, in particular in the horizontal components at QD latitudes of ±(65 − 80) • , is obvious. Also the signature of the mid-latitude Sq current system in the dayside data is clearly visible. All this confirms that GRACE magnetic data, after calibration, can be used for geophysical investigations.
As a final example, we present in Fig. 11 the timeevolution of large-scale magnetospheric magnetic field contributions in March and April 2015 as determined with magnetic data of various satellites. We assume that the magnetospheric magnetic field disturbance δB(t) = −∇V ext (t) can be derived from a scalar potential of external origin which is expanded in terms of spherical harmonics: where a = 6371.2 km is Earth's mean radius, θ d is dipole co-latitude, T d is dipole Local Time (see e.g. Hulot et al. 2015), and q m n (t) and s m n (t) are the Gauss coefficients describing external sources. They are estimated in bins of 1.5 h length (which is the approximate orbital period of LEO satellites) from the observed horizontal magnetic field components δB θ d , δB φ d in the dipole coordinate frame, after removal of core and crustal field contributions as given by the CHAOS-7 model. By using only the magnetic horizontal components we do not attempt to separate external (magnetospheric) and internal (electromagnetic induced) contributions; the presented results thus include both.
Magnetic measurements from six satellites downsampled at 15 s are used for this: In addition to GRACE-A and GRACE-B (the orbital planes of which were at 05:45/17:45 LT in March 2015) we used data from CryoSat-2 (00:40/12:40 LT), Swarm Alpha and Charlie (07:50/19:50 LT) and Swarm Bravo (09:30/21:30 LT). The upper part of Fig. 11 shows individual estimates of q 0 1 (t) alone, determined separately using data from each of the six satellites. q 0 1 (t) monitors the strength of a symmetric magnetospheric ring-current in the dipole equatorial plane. One of the largest geomagnetic storms of the recent years occurred on 17 March 2015, denoted as "Saint Patrick's Day Storm" and is well captured by all Fig. 9 Local Time dependence of GRACE-A daily mean (left) and misfit rms (right) values at non-polar latitudes for the three vector components in the NEC frame and for scalar intensity, respectively satellites and by the Dst-index that is derived from geomagnetic ground observatory data. Although there is rather good agreement between these separate estimates of q 0 1 , there are also subtle differences, likely caused by deviations from a symmetric ring current. Since the satellites sample the magnetic field simultaneously in different Local Time sectors, their combined analysis allows for a determination of such deviations from a symmetric ring current. The lower panel of Fig. 11 shows estimated time series of the 8 coefficients of Eq. 9 with N = 2 . The symmetric ring-current, modelled by q 0 1 , is by far the largest source, but non-symmetric contributions are particularly enhanced during geomagnetic storms, as seen, e.g. in the coefficients q 1 2 and s 1 2 . Kuvshinov et al. (2021) performed a combined analysis of magnetic data from Swarm and CryoSat-2 and estimated external and induced spherical harmonic expansion coefficients in their attempt to determine the electrical conductivity of Earth's mantle. The results of the above analysis demonstrates that GRACE data are also capable of monitoring magnetospheric field variations. We hope that the scientific community will find GRACE magnetic data useful for scientific studies, both regarding sources in Earth's interior and space physics/ space weather applications. One example of research that used (a preliminary version of the) GRACE magnetic data is the geomagnetic field modelling study by Kloss et al. (2021), who jointly analysed magnetic data from the satellites Ørsted, CHAMP, CryoSat-2, GRACE and Swarm. Top: time evolution of q 0 1 (t) describing the strength of the symmetric ring-current as determined separately from magnetic observations of six satellites. Bottom: time evolution of spherical harmonic expansion coefficients up to degree N = 2 by joint analysis of data from these 6 satellites