A Spherical Harmonic Model of Earth's Lithospheric Magnetic Field up to Degree 1050

Detailed mapping of the Earth's magnetic field brings key constraints on the composition, dynamics, and history of the crust. Satellite and near‐surface measurements detect different length scales and are complementary. Here, we build a model, selecting and processing magnetic field measurements from the German CHAMP and ESA Swarm satellites, which we merge with near‐surface scalar anomaly data. We follow a regional approach for modeling the magnetic field measurements and next transform the series of regional models into a unique set of spherical harmonic (SH) Gauss coefficients. This produces the first global model to SH degree 1050 derived by inversion of all available measurements. The new model agrees with previous satellite‐based models at large wavelengths and fits the CHAMP and Swarm satellite data down to expected noise levels. Further assessment in the geographical and spectral domains shows the model to be stable when downward continued to the Earth's surface.

4 km altitude above the Earth's ellipsoid. Alternative global scalar anomaly grids, such as the Earth Magnetic Anomaly Grid at 2 arc minutes resolution (EMAG2, Maus et al., 2009) and its latest update, the EMAG2v3 grid, were also constructed along the same lines (Meyer et al., 2017).
Both WDMAM and EMAG2 scalar anomaly grids were next converted into sets of spherical harmonic (SH) Gauss coefficients, using linearization and regularization (Lesur et al., 2016;Maus, 2010) to recover models of the vector magnetic field. However, no global inversion combining all available near-surface gridded data and satellite direct measurements have yet been performed to best take advantage of the data and build a fully consistent joint model.
Here, we carry such a joint inversion, using CHAMP and Swarm measurements combined with the WD-MAM-2 anomaly grid to build a model of the Earth's lithospheric magnetic vector field with a 40-km horizontal spatial resolution. For the first time, all measurements and data types (satellite vector, scalar, and gradient data, as well as aeromagnetic and marine compilations) are simultaneously used. The method involves an iterative piecewise approach using the Revised Spherical Cap Harmonic functions (R-SCHA, Thébault, 2006) to allow independent regional data analyses, build regional models and combine them into a global SH model. This makes it possible to best take advantage of all available measurements at reduced and manageable numerical costs.
This study is organized as follows. Section 2 describes the CHAMP and Swarm satellite data, the selection and processing methods we use, and provides a summary of the WDMAM scalar grid. Section 3 describes the regional mathematical representation of the magnetic field and the inversion approach. Inversion statistics are next presented in Section 4, where we also assess and discuss the final SH model, which we compare to earlier models.

Swarm and CHAMP Satellite Data Selection
We use 1 Hz Level 1b vector and scalar measurements of the lowest Swarm A and C satellites between March 1, 2014 and March 31, 2020, and between 430 and 460 km altitude. Data are selected only when the instruments are in nominal mode and when the thrusters are off. They are next sub-sampled every 5 s, which corresponds to about 35 km along-track spacing at satellite altitude.
Only magnetic field measurements corresponding to magnetospheric quiet times are used. These are identified by the planetary Kp index being lower than 2° in value both at time of measurement and during the previous three hours. This selection is refined by further requesting the Disturbed Storm Time Index (Dst) to not exceed 5 nT in absolute value and not vary by more than 5 nT during the three previous hours. We additionally require the interplanetary magnetic field (IMF) to be less than 5 nT in absolute value for the y-component and positive for the z-component, for a better rejection of measurements corresponding to precursors of strong events such as geomagnetic storms and substorms (Richardson & Cane, 2011).
The statistical properties of the remaining perturbation fields after selection are different at low and high latitudes. We separate the selected data set into midlatitudes (between −60° and 60° magnetic latitudes MLAT, as defined by the Quasi-Dipole coordinates at satellite altitude, Richmond, 1995) and polar-latitude regions (between 50° and 90° and −90° and −50° MLAT for the north and south caps, respectively). The polar and midlatitude subsets overlap by 10° in latitude. This precaution allows more robust further dedicated data corrections, mitigates edge effects, and ensures a better continuity between mid and polar regions (see also Thébault et al., 2013 for the algorithm description). Nighttime measurements between 21:00 and 5:00 local time (LT) are kept at midlatitudes to minimize contamination from the pseudo-periodic diurnal ionospheric Sq field (see Chulliat et al., 2013, for instance). In the polar regions, the measurements are taken in the dark and when the Sun is at least 10° below the horizon.
We also include a selection of the latest calibrated CHAMP measurements (GFZ Section 2.3, 2019) from January 2001 to September 2010 between 250 and 480 km altitude. As is done for the Swarm data, the CHAMP data are split into mid and high magnetic latitude datasets with a 10° overlap between regions. At midlatitudes, the measurements are selected between 22:00 and5:00 LT and all local times are considered in the polar areas with the condition that the satellite was in the dark with the Sun located at least 5° below the horizon. A final data selection is made against the Kp and the Dst magnetic indices and the IMF values following the same criteria as those considered for the Swarm data.

Satellite Data Correction
We correct the Swarm and the CHAMP measurements for the internal core field and its temporal variations up to SH degree 15 and for the external primary and internal induced magnetospheric fields using the CHAOS-7.5 model (Finlay et al., 2020). This model covers the full CHAMP and Swarm data epochs and was derived with the same calibration versions of the CHAMP and Swarm data. After these corrections, vector and scalar gradient data sets are also built. For Swarm, across track gradients are constructed from satellites A and C by selecting pairs of data measured at the same Universal Times at all latitudes. For CHAMP and Swarm, along-track gradients are constructed considering 15 s time spacing between measurements along the satellite orbits, following Finlay et al. (2020).
The corrected data show large-scale remaining structures at mid-latitudes and significant transient perturbations in the polar areas blurring the lithospheric field signal. At mid-latitudes, these signals are caused by rapid magnetospheric and ionospheric fields not accounted for by the CHAOS-7.5 model correction. In polar areas, perturbations along satellite orbits are mostly due to the magnetic field of Polar Electrojet and Field Aligned Currents (FAC). This calls for dedicated corrections using additional filters along the Swarm and CHAMP satellite trajectories.
Such dedicated corrections, only applied to portions of the sphere, require some care to minimize spectral leakage and avoid filtering out some of the signal carried by the large-scale lithospheric field structures (Thébault et al., 2012). To minimize such effects, we first correct all measurements with a lithospheric field model expanded to SH 80 derived from 5 years of Swarm measurements following a procedure described in Thébault et al. (2016). This model, which we later reinstate, is available on the ESA Swarm Payload Data Ground Segment (PDGS) under the file name SW_OPER_MLI_SHAi2D_00000000T000000_99999999T99 9999_0501.shc.
At midlatitudes, the perturbation field contaminates mostly the zonal terms (Maus et al., 2008). To deal with this, we next correct the data along each satellite track in the geomagnetic reference frame for an external SH degree 1 field and its internal induced part, and an additional SH degree two external zonal term. In polar regions, the transient external fields are more effectively reduced with a non-parametric filter such as the Singular Spectral Analysis (SSA; see Thébault et al., 2017). For the SSA, the perturbation field is mostly carried by the first two singular values of the SSA. A final screening is performed to manually reject satellite track measurements with obvious first-order issues such as spikes or outlier values.
After all corrections, only measurements between 55° and 90°, and between −90° and −55° MLAT for the north and south caps, respectively, and between −55° and 55° MLAT at midlatitudes are kept, thus removing the overlap areas between mid and high magnetic latitudes regions that were considered to minimize edge effects between regions. The contributions of the lithospheric field magnetic field to SH degree 80 are then added back to all data.

Near-Surface Data
Several near-surface magnetic anomaly compilations are publicly available. The EMAG2 scalar anomaly grid (Maus et al., 2009) is not considered in this work as it misses some recent regional updated compilations. Its latest version, the EMAG2v3 (Meyer et al., 2017) is not considered either, since parts of oceanic areas contain wide regional data voids and sharp transitions. To minimize numerical instabilities and spatial ringing in our model, we select the WDMAM-2 grid that offers the most extensive geographical coverage (Lesur et al., 2016). About 3.5% of the data sets are synthetic but built from reliable and carefully calibrated geophysical priors (Dyment et al., 2015). As is the case for all global scalar anomaly grids, the large wavelengths are imposed by a satellite-based model. In the WDMAM-2, the structures corresponding to SH degrees between 16 and 100 originate from the GRIMM_L120 lithospheric field model downward continued to 5 km altitude (Lesur et al., 2013). These large-scale structures might not be fully compatible with our selected and processed measurements. We describe in the next section an approach designed to minimize this issue.
The WDMAM-2 grid is provided as a grid with a 5 × 5 km cell size at 5 km altitude above the WGS84 geoid (Kovalevsky et al., 2012). We reduce the 25,927,200 total intensity anomaly data of the initial grid to about a 0.15 × 0.15° cell size equal area grid using a cubic B-spline interpolation (about 15 km spacing). We are therefore left with 3,240,000 intensity anomaly values converted from the WGS84 geodetic to the geographic Earth-Centered Earth-Fixed (ECEF) coordinate system. This reduction procedure is necessary to avoid deriving a lithospheric magnetic field model too biased toward near-surface measurements.

Data Inversion and Statistics
The scalar WDMAM-2 anomaly data and the CHAMP and Swarm vector, scalar, and gradient data (see Table 1) are regionally modeled with the R-SCHA mathematical functions. The Earth is tiled with 700 caps placed on the nodes of an equal area grid. A description of the way this procedure was applied globally on the CHAMP measurements (third calibration) is provided and illustrated in Thébault (2006, their Figure 1) and the way the algorithm has been extended for taking advantage of gradient data is detailed in Thébault et al. (2013). The same general procedure is used here.
We solve a least-squares inverse problem within each spherical cap. In the spherical cap's reference frame where α = 6,371.2 km is the Earth's reference radius, r is the radius, θ is the colatitude, and ϕ is the longitude, the R-SCHA mathematical basis functions express the magnetic potential V as a sum of the Legendre (V 1 ) and Mehler (V 2 ) potentials where the Legendre potential writes    Note that when considering such regional modeling, both sets of parameters E m p G and , E e m n k G are needed to account for magnetic field sources creating a field not fulfilling the null condition at the edges of the modeled region. This generally happens because magnetic field sources outside the modeled region contribute to the magnetic field measured inside this region (see Thébault et al., 2006).
The magnetic field (by convention the negative gradient of potential V) is solved within caps of 8° half aperture (about 1,800 km diameter at Earth's mean radius) and between the minimum and maximum available data altitudes in each region. The maximum expansion index k in the potential V 1 depends on the data availability on the cap's horizontal surface and the p index in the potential V 2 depends on the data availability with the altitude. We therefore restrict the Mehler's solution to a maximum index p = 4 and the Legendre's solution to a maximum index k = 150, for consistency with the near-surface and satellite data coverage on the sphere and with their altitude. We also select only degrees For the construction of the R-SCHA basis functions, we apply a procedure designed to both optimize the overall computation time and minimize the a priori large-scale satellite-based information contained in the WDMAM-2 anomaly field data. We rely on the statistical equivalence between the real-valued R-SCHA degrees n k and the integer SH degree n. Since the CHAMP and Swarm satellite measurements are available between 250 and 460 km altitude, the signal-to-noise ratio at length scales smaller than these altitudes is low. To account for this, R-SCHA basis functions corresponding to  300 E k n are not computed at satellite altitudes (typically corresponding to spatial structures smaller than about 125 km; that is, half of the minimum altitude found to be 250 km in the CHAMP satellite measurements). This approximation saves a significant amount of computational time. It was fully tested on synthetic simulations and is also justified considering the decrease of the Legendre potential values with increasing radius r.
The analogy between the standard SH and the R-SCHA degrees is further used to downweight the largescale GRIMM_L120 lithospheric field model contribution present by construction in the WDMAM-2 scalar anomaly data. This is done to avoid unnecessarily constraining our model by the information provided by GRIMM_L120. We define a weight function for the Legendre basis functions with a cutoff wavelength E c λ = 800 km (i.e., SH degree 100) using a high-pass Butterworth filter of order k = 6 At R-SCHA degrees n k ∼ 15 (i.e.,  2670km E k λ ) this weight is equal to 0.5 and at degrees n k ∼ 100 (i.e.,  400km E k λ ) it reaches 1. We multiply the Legendre basis functions computed at the near-surface scalar altitudes by the weight values   E k w λ for each degree n k . In contrast, the basis functions computed at satellite altitudes are not downweighted, to allow CHAMP and Swarm data to constrain the large lithospheric field wavelengths. The WDMAM-2 scalar anomaly grid is therefore downweighted in the GRIMM_L120 spectral band and smoothly weighted toward unity at the equivalent SH degree 100, beyond which the grid no longer contains the a priori GRIMM_L120 truncated model. Completely canceling out the WDMAM-2 contributions in the GRIMM_L120 spectral band (i.e., setting    0 ) was not possible as this was found to lead to instabilities in the inverse problem.
For the data inversion, scalar gradients and scalar anomaly data are linearized along the Earth's main field. We follow Lesur et al. (2016) and linearize the WDMAM-2 scalar anomaly data at each location along the CM4 vector main field model (Sabaka et al., 2004) for the 1990 reference epoch of WDMAM-2. The CHAMP and the Swarm magnetic field scalar and scalar gradient data are linearized along the CHAOS-7.5 time-varying core field model (Finlay et al., 2020) up to SH degree 15 at each data location and corresponding epoch between 2000 and 2020.
A regional regularization is applied in the inverse problem. Indeed, Backus (1970) showed that uniqueness is not guaranteed when trying to recover the internal geomagnetic vector field from just the knowledge of its intensity. When producing models, this leads to a "Backus effect" responsible for noise amplification in the direction perpendicular to the main field, particularly near the magnetic equator (e.g., Khokhlov et al., 1997). Global models involving SH representation usually rely on a global regularization (e.g., see Equation 39 in Maus, 2010). Here, however, the regional approach used allows us to account for the Backus effect at a regional scale and to minimize the crustal field model components that are perpendicular to the main ambient field only for spherical caps whose centers are located within a magnetic latitude band ±40° around the magnetic equator.
For the numerical least-squares inversion, we apply a reweighting scheme using Huber weights to minimize the leverage of discrepant magnetic field measurements (Constable, 1988). However, the usual reweighting scheme involving many iterations was not found to be optimal when merging the different data types (and data altitudes). The reason for this is that the WDMAM-2 data contains information at horizontal wavelengths smaller than 40 km that does not comply with the Huber statistics (based on a mixture of Gaussian and Laplacian error distributions). In the full iterative reweighting scheme, this leads to regional models with artificially low power in regions where the near-surface magnetic anomalies have complex and strong signal at length scales smaller than 40 km. In practice, a maximum number of two iterations was found to be more appropriate for preserving the signal amplitude and resolution at Earth's mean radius.
Each of the 700 R-SCHA models is described by 5,379 regional parameters. They can be derived independently, allowing parallelization, independent assessment, and testing of each regional solution. This also saves computational and memory resources otherwise hardly manageable in a single global SH inversion. The series of regional R-SCHA models expressed at Earth's mean radius was finally transformed into a consistent set of SH Gauss coefficients to degree and order 1100 (using Thébault et al., 2016). However, the model is released only up to SH degree 1050, considering that the last harmonics may be contaminated by spectral leakage of structures smaller than the maximum model resolution.

Results and Discussion
The satellite data are successfully fitted by the regional models to a root mean square residuals lower than 4 nT in polar areas (Table 1). These residuals reflect the presence of remaining transient external perturbation field in the preprocessed data, despite the carefully applied dedicated correction. A few large-scale residual structures also remain at mid-latitudes, following the magnetic equator, as a result of the more conventional and less stringent correction in terms of low SH functions. The mean residuals for the CHAMP intensity data are slightly biased toward negative values while Swarm intensity data residuals are close to zero. This peculiarity has already been observed when combining data from different satellite missions (see, e.g., Finlay et al., 2020, their Tables 6 and 7) and could result from the way their instruments are calibrated in space. Figures 1a and 1b show example scalar anomaly fields predicted by our final global model over Australia and the North Atlantic Ocean, as well as residuals between these model predictions and the original WD-MAM-2 grid. Residuals are mostly located over complex crustal magnetic field structures containing spatial wavelengths smaller than 40 km. These near-surface residuals also display large-scale structures testifying for disagreements between our model and the GRIMM_L120 lithospheric field model used to complement ground data in the WDMAM-2 grid. Some of these features are clearly apparent over Australia but also offshore of North Africa and Spain, in the form of North-South stripes with an East-West change of polarity. Such structures are known to be artifacts related to external field contamination when using satellite measurements (Thébault et al., 2012(Thébault et al., , 2017. The fact that they show up in our residual maps confirms that these GRIMM_L120 artifactual structures do not comply with the CHAMP and additional Swarm satellite measurements we used and that they are properly rejected by our model. This a posteriori validates our procedure to downweight contributions from GRIMM_L120 when using the WDMAM-2 grid. Over oceans, residuals are locally larger and structured along oceanic chrons. Figure 1b illustrates this point. It displays magnetic field structures and residuals with sharp transitions over the North Atlantic region. Some of these transitions are well covered and documented by marine and airborne measurements but not captured by our model. It is for instance the case over the Reykjanes ridge, south of Iceland, and down to the Bight oceanic transform fault. Accounting for these would require a higher resolution SH model. In other regions, such as those highlighted using light gray shading in Figure 1b, residuals also reveal small-scale patterns. They correspond to additional information added to the WDMAM-2 grid to fill oceanic data gaps using the magnetization model of Dyment et al. (2015). Significant residuals in oceans are globally well correlated with the location of these synthetic data because such discontinuous and geometrical patterns are difficult to resolve with continuous mathematical functions.
A map of the radial component predicted at Earth's mean radius by our SH vector lithospheric field is finally shown in Figure 1c. Magnitude is generally weaker near the equator than near the pole, as expected and as had been observed in the WDMAM-2 scalar anomaly grid. This is because the lithospheric field induced by the Earth's main field dominates in the continental crust and increases with absolute latitude (Maus & Haak, 2002). Figure 2a shows the power spectrum of the SH model plotted at Earth's surface between SH degrees 16 and 185. This range allows visual comparisons with the spectra of the MF7 model derived from the final 2 years of the CHAMP data (Maus et al., 2008), the LCS-1 model derived from the CHAMP and Swarm satellite observations (Olsen et al., 2017), and the WDMAM-2 model (Lesur et al., 2016). Between SH degrees 16 and 100, the present model closely follows the MF7 and LCS-1 model and deviates from the WDMAM-2 model. From SH degrees 100 onwards, the model then closely follows the WDMAM-2 model. This behavior is what was expected and further suggests that the model properly accounts for both the satellite and near-surface data we used. Figure 2b displays the power spectrum of our model over its entire range, which we now compare to the spectra of the NGDC-720 model derived from the EMAG2 grid (Maus, 2010), of the Enhanced Magnetic Model (EMM-2017) computed from the EMAG2v3 grid (Meyer et al., 2017) and of the WDMAM-2 model (Lesur et al., 2016). All spectra show some varying amount of power dip between SH degree 120 and about 200. This results from the lack of constraints provided by the satellite and near-surface data in this degree range. For larger degrees, all spectra have comparable levels up to SH degree 300, beyond which their slopes significantly differ. Though derived from different datasets and ocean models, the NGDC-720 and the WD-MAM-2 display similar spectra up to about SH degree 600, before becoming flat. In contrast, the spectrum of our model exhibits a quasi-constant negative slope up to SH degree 830 beyond which the slope becomes steeper. A decreasing, rather than flat, power spectrum is consistent with remaining data gaps and testifies for the stability of the model.
We also tested our model against a statistical representation of its spectrum between degrees 16 and 1050 (shown in black with its 95% Confidence Interval [CI] in Figure 2b). The statistical model is parameterized in terms of a mean apparent crustal magnetization, a mean magnetic crustal thickness, and a power-law related to the magnetic sources, following Thébault and Vervelidou (2015). The least-squares fit between the observed and the statistical spectra leads to a mean magnetization estimate of 0.4 ± 0.3 A m −1 (95% CI), a mean magnetic thickness equal to 36 ± 6.5 km (95% CI) and a power-law equal to −1.36 ± 0.32 (95% CI). These estimates are realistic and in agreement with independent geophysical values (see Thébault & Vervelidou, 2015). The observed power spectrum lies outside the 95% CI of the statistical fit only for SH degrees beyond 830, where the observed power spectrum has a steeper slope because of data gaps.
The model presented in this work is the first successful lithospheric magnetic field model derived up to SH degree 1050 and based on a joint inversion of publicly available near-surface global compilations and CHAMP and Swarm satellite data. Using the same approach, additional models with increasing resolution could next be built as more satellite and near-surface measurements become available. An important target will be to constrain the spectral gap between SH degrees 120 and 180. Different scenarios are currently explored for the evolution of the ongoing Swarm mission. Measuring the magnetic field scalar and vector gradients at steadily decreasing altitudes during the Swarm A and C satellites re-entry will provide invaluable low altitude measurements for filling this gap. In parallel, additional efforts to compile and add more marine and airborne measurements in global compilations such as the WDMAM are also much needed to improve the robustness of the still poorly constrained intermediate lithospheric field wavelengths in these compilations.

Data Availability Statement
The Swarm data are available at ftp://swarm-diss.eo.esa.int/Level1b/Entire_mission_data/MAGx_LR/ or by any web browser at https://swarm-diss.eo.esa.int/#swarm%2FLevel1b%2FEntire_mission_data%2FMAGx_ LR. The World Digital Magnetic Anomaly Map data are available at www.wdmam.org. The CHAOS-7.5 core  Maus et al., 2008), of the LCS-1 model (blue, Olsen et al., 2017), and of the WDMAM-2 model (green, Lesur et al., 2016); (b) Full spectra of the new SH model (up to SH degree 1050, red), of the NGDC-720 model (up to SH degree 720, blue, Maus, 2010), of the WDMAM-2 model (up to SH degree 800, green, Lesur et al., 2016) and of the EMAG2v3/EMM-2017 model (up to SH degree 790, pink, Meyer et al., 2017). The black curves are the maximum likelihood of a statistical model and its 95% CI best fitting the new SH model using the statistical expression of Thébault and Vervelidou (2015). field model is available at https://www.spacecenter.dk/files/magnetic-models/CHAOS-7/. The lithospheric field model to degree 1050 presented in this study is available in spherical harmonics at the Zenodo CERN's Data Center https://zenodo.org/record/5546528.