A Swarm lithospheric magnetic field model to SH degree 80

The Swarm constellation of satellites was launched in November 2013 and since then has delivered high-quality scalar and vector magnetic field measurements. A consortium of several research institutions was selected by the European Space Agency to provide a number of scientific products to be made available to the scientific community on a regular basis. In this study, we present the dedicated lithospheric field inversion model. It uses carefully selected magnetic field scalar and vector measurements from the three Swarm satellites between March 2014 and December 2015 and directly benefits from the explicit expression of the magnetic field gradients by the lower pair of Swarm satellites. The modeling scheme is a two-step one and relies first on a regional modeling approach that is very sensitive to small spatial scales and weak signals which we seek to describe. The final model is built from adjacent regional solutions and consists in a global spherical harmonics model expressed between degrees 16 and 80. The quality of the derived model is assessed through a comparison with independent models based on Swarm and the CHAMP satellites. This comparison emphasizes the high level of accuracy of the current model after only 2 years of measurements but also highlights the possible improvements which will be possible once the lowest two satellites reach lower altitudes.


Background
The magnetic field of the Earth's crust has spatial scales ranging from meters to 1000s of kilometers [e.g., Thébault et al. (2010b) for a review]. It is generated by the (induced or remanent) magnetized rocks and depends on parameters such as their mineralogy, temperature, chemical alteration, age, deformation as well as on the Earth's main field orientation (Dunlop and Özdemir 2007). The lithospheric magnetic field is thus of great interest for geophysical explorations and tectonic interpretations. After the early success of the Magsat and Oersted missions (e.g., Purucker et al. 2002) and that of the CHAMP mission (e.g., Maus 2010), its accurate global mapping from space is currently one of the main goals of the innovative Swarm mission launched by the European Space Agency (ESA) on November 22, 2013(Friis-Christensen et al. 2006).
The Earth's magnetic field is the sum of fields of internal and external origins that operate on a wide range of spatial and temporal scales . The three major sources are (1) the main field, which is generated by electrical currents in the liquid outer core, (2) the external fields produced by currents in the ionosphere and the magnetosphere, and (3) the lithospheric field. Since the core and external fields mask the contributions of the lithosphere, current practice is to map the lithospheric field in spherical harmonics (SH) either by a joint inversion with other source magnetic fields or only after the measurements are corrected for the core and external fields. The first approach, hereafter referred to as the comprehensive inversion (CI), relies on mathematical parameterizations that describe each source field on average . The second approach is sequential (sequential inversion, SI) and relies on designing, testing, and applying different corrections by trials and errors (e.g., Maus et al. 2008). The SI requires constant supervision of the operator. However, since it is purposefully developed for modeling the lithospheric field, it is also expected to provide more detailed spatial structures of this source field.

Open Access
*Correspondence: erwan.thebault@univ-nantes.fr 1 Laboratoire de Planétologie et de Géodynamique de Nantes, UMR CNRS 6112, University of Nantes, 2 chemin de la Houssinière, 44300 Nantes, France Full list of author information is available at the end of the article For the identification of the large scales of the lithospheric magnetic field, the data measured by the low Earth orbiting satellites are crucial because they provide, in a statistical sense, homogeneous measurements on a global scale. Unfortunately, the lithospheric field strongly attenuates with altitude and its mapping is challenging. The signal can furthermore be masked by time-varying external fields that create anisotropic errors which are aligned along the satellite orbits. Attempts to correct for these contaminations can be done by along-track analysis (Maus et al. 2006a;Thébault et al. 2012) or statistically ), but the transient nature of these disturbances makes these approaches imperfect. As a result, satellite lithospheric field models often show disagreements in their sectorial SH terms. At present, cross-comparisons between lithospheric field models suggest that the lithospheric field can be imaged robustly on a global scale up to SH degree 85 [see for instance Olsen et al. (2014)].
A major breakthrough is expected from the Swarm satellite mission thanks to the orbital configuration of its three polar orbiting satellites (Friis-Christensen et al. 2006). The east-west (EW) separation of the lowest Swarm A and C satellites can be exploited to approximate the EW magnetic field gradient, and differences between successive vector measurements along each individual satellite track can be used to estimate north-south gradient (NS) (e.g., Maus et al. 2006a;Kotsiaros and Olsen 2014). These EW and NS components provide 5 out of the 6 components required to fully characterize the gradient of the magnetic field. From a mathematical viewpoint, such gradient components do not provide more information than the vector measurements. However, they help avoiding undesired external field signals and amplify the signal-to-noise ratio of the lithospheric field over wavebands that are poorly constrained by just vector and scalar measurements (Kotsiaros and Olsen 2012). The proof of concept of Swarm was recently validated by the preliminary Swarm Initial Field Model (SIFM, Olsen et al. 2015).
In this paper, we use 22 months of Swarm measurements and derive a new lithospheric field model by applying a procedure that was developed during the preparation phase of the Swarm mission . This dedicated lithospheric field inversion (DLFI) is the one currently used to produce one of the official ESA level 2 candidate products in the framework of the Swarm Satellite Constellation Application and Research Facility (SCARF; Olsen et al. 2013). It relies on a selection of the measurements to minimize the contributions of transient external fields ("Background" section) and on corrections for the main and average external magnetic fields ("Swarm data selection" section). The lithospheric field is then modeled using regional functions ("Data modeling" section) that are finally transformed into a unique set of SH Gauss coefficients. The output model is discussed and compared to available lithospheric field models ("Results" section).

Swarm data selection
The algorithm we employ is sketched in Fig. 1. It has been constantly modified since the beginning of the launch of Swarm and is applied to its measurements for the first time. The numerical protocol begins with the Swarm data selection (process 1). We select the 1-Hz level 1b measurements of Swarm A and C between March 1, 2014, corresponding to the satellites insertion into their final orbits, and December 31, 2015. This dataset is downsampled by taking only one vector and scalar measurement every 5 s. We use the magnetic data with baseline 0408 for Swarm A and baseline 0409 for Swarm C. We select vector measurements of Swarm A and C measured by the vector field magnetometers (VFM) and the scalar measurements acquired by the absolute scalar magnetometer (ASM) in nominal mode only. However, due to an anomaly of the ASM on Swarm C [see Fratter et al. (2016)], we select the ASM scalar data in nominal mode only up to November 2014 and thereafter use the field strength measured in nominal mode by the calibrated VFM magnetometer. Measurements of satellite B, taken at an average altitude of 510 km, are not used as they are not as sensitive to the lithospheric magnetic signal as the Swarm A and C satellites.
Middle-and high-latitude measurements are processed separately in order to account for the different statistical properties of the external source fields. High-latitude measurements, defined as those covering regions of poleward magnetic latitudes larger than |±52°|, are selected when the Sun is at least 10° below the horizon, regardless of the local time (LT). In contrast, low-latitude measurements in the magnetic latitude range of −52° to 52° are selected for LT between 23:00 and 5:00. For the purpose of the along-track pre-processing (detailed hereafter), we construct a provisional datasets where the latitude limit of the selected measurements extends up to |±62°| for low-latitude measurements and down to |±42°| for the polar measurements. These extra measurements are rejected after the dedicated corrections.
We select vector and scalar measurements only when the hourly Dst index is lower than ±5 nT, with variations lower than 5 nT over the three previous hours, in order to minimize the external fields generated by the ring current. Similarly, measurements associated with Kp index lower than 2 − during the considered time interval and the previous 3-h one are selected. The measurements corresponding to reduced magnetic field perturbation are selected using the interplanetary magnetic field (IMF) values of |IMFBy| <2 nT and IMFBz >0. Finally, we reject all selected segments of satellite track smaller than 10° angular length.

Correction with SH models
In a second stage of the processing, we correct the measurements for the main and external magnetic fields (process 2a, 2b, and/or process 3 in Fig. 1). Explicit models relying on the parameterization of the main internal and external fields in SH are most appropriate for such data correction. The Swarm processing chains developed within the SCARF consortium ) provide various models for these source fields (Hamilton 2013;Rother et al. 2013;Sabaka et al. 2013). However, data corrections based on models derived from different dataset Fig. 1 Flowchart of the processes runs to generate a lithospheric field model from the level 1b Swarm magnetic field measurements may lead to bias and errors (e.g., Sabaka and Olsen 2006). For example, the CI core field model relies on a selection that uses both night and day data, which is different from the strict LT and Sun elevation selection criteria we use here. This model also attempts to separate the core and the internal fields induced in the Earth's mantle by the external fields . As a consequence, it is preferable to use either the full CI model or a dedicated core field model together with an additional model for the ionospheric field , to ensure that the contribution of the ionosphere induced in the Earth's mantle is also corrected for. Additional magnetospheric field corrections can then also be implemented based on the so-called Model of the MAgnetosphere (MMA) SH model (Hamilton 2013). As an attractive alternative to the CI model the Dedicated COre field model (DCO) can, for instance, be used. It is derived from night-side data and does not separate core and internally induced fields . It thus provides the correction sought with a minimum risk of bias. For the current release of L2 models to ESA, the work flow was designed to ensure that DLFI model is only computed once the DCO and CI models are available. For the present study, we use the DCO model to correct the measurements for the Earth's core field over the time span of interest.
In order to complement this first-order correction and to allow the DLFI chain to also operate independently of the SCARF sequences of model production, an alternative processing block (process 3) is used. We derive an additional built-in model of the Earth's main and lithospheric internal field up to SH degree 60, the main field secular variation and acceleration up to SH degrees 10 and 8, respectively, and a simple model for the magnetospheric field to SH degree 1 scaled with the hourly Dst index [details about this relatively standard parameterization may be found in Thébault et al. (2010a)]. The internal part of the model up to SH degree 15 and its magnetospheric part can then be used to correct the measurements for both these main and magnetospheric field contributions.
For the present study, since the DCO model had been made available to us, a combined approach is in fact used. The data are first corrected for the DCO model over the time span of interest (process 2a), and the residuals are next used to compute a model itself used to further correct the data in the way described above (process 3). The root-mean-square (RMS) of the final residuals is then plotted within 1° longitudinal bands over all latitudes to identify satellite tracks showing abnormal RMS values.

Empirical corrections
The magnetospheric field signal from the ring current, which is at least of the same order as the magnitude of the lithospheric field strength, is dynamic and varies from one orbit to the next. The Dst is an hourly index, and its parameterization for the modeling of the magnetospheric field in SH to degree 1 (in process 2b or/and 3) is not sufficient to account for all its time and space variations. This is especially true for time variations shorter than 1 h which are responsible for systematic offsets between the magnetic field measurements on nearby satellite orbits. To correct for this, we introduce an additional processing block (process 4). For each low-latitude satellite track, we fit an SH external dipole (three parameters) with its internally induced zonal dipole part (one parameter) in the inclined internal dipole reference frame and remove this model from the measurements. Such along-track filtering must be applied with care as it also tends to filter out genuine lithospheric magnetic field structures (e.g., Maus et al. 2006a, b) and to introduce artifacts (Thébault et al. 2012). To minimize its adverse effect, we first subtract the lithospheric field contribution from SH degrees 16-45 of the model derived in process 2 from the scalar and vector measurements before estimating the alongtrack correction. The along-track estimation of the external field parameters is performed using jointly Swarm A and Swarm C measurements whenever possible. Since Swarm A and Swarm C travel along slightly different orbits, this joint inversion introduces a longitudinal constrain that minimizes the error of correction. The alongtrack primary and induced external field correction is the one we apply to the selected and corrected data.
At mid-latitudes, this processing reduces considerably the offsets between the adjacent satellite tracks but not perfectly particularly along the magnetic equator and South East Asia (see Fig. 2). However, initial tests with this Swarm dataset demonstrated that more complex data processing along the satellite tracks would distort the lithospheric field structures and lead to models with less power than expected.
At polar latitudes, this SH along-track analysis is also imperfect. The data selection process minimizes the ionosphere diurnal variations but cannot avoid all ionospheric field contributions. In particular, the field produced by the polar electrojets persists at night times. It is confined within the auroral oval (between 55° and 75° absolute latitude) during magnetically quiet times, and along-track correction with low SH degree models is not as efficient. Maus et al. (2008) propose an additional line leveling procedure that consists in minimizing the differences at crossover points between satellite tracks. The leveling is effective for mapping the lithospheric field. However, the corresponding inverse problem is illposed because one can always add a constant to the satellite tracks without affecting the difference between the tracks (Menke 2012). Setting this constant value a priori is challenging, and an erroneous value can lead to a lithospheric field model with a biased offset. For this reason, we do not apply the line leveling procedure for the present model.
Rather, at polar latitudes, we stick to the SH along-track correction and implement an additional empirical correction based on singular spectral analysis (SSA). SSA is a numerical method that decomposes the signal into its principal components (Golyandina and Zhigljavsky 2013). The vector and scalar measurements of each Swarm A and C satellite are analyzed independently track by track. For each polar portion of a satellite half orbit, the measurements are treated as a series B i (r), where r represents the location of the ith magnetic component (vertical, north, east, or scalar). In a first step, the covariance matrix of the series B i (r) is computed and decomposed by singular value decomposition (SVD). This SVD provides a set of eigenvalues sorted by order of magnitude together with a matrix of empirical orthogonal functions (EOF). The projection of the signal B i (r) onto each EOF that is multiplied by the corresponding set of eigenvalues is then used to identify and filter out unwanted features. The filtering is done without a priori information about the nature or the shape of magnetic field structures; the procedure is model independent. This flexibility is both the main advantage and the main disadvantage of SSA. Compared to SH filtering, no information is required with respect to the typical length scale of the noise along the satellite orbit. SSA can handle automatically changes in the spatial scales of the external field signal due to changes in the magnetic activity. Contrary to techniques based on Fourier analyses [see, for instance, Langel and Hinze (1998)], SSA is also less prone to severe ringing and aliasing that often produce artificial NS oscillations in the corrected measurements. However, SSA first requires the identification of the set of EOFs that carry the external fields. This is a crucial yet fairly arbitrary step that is left to the operator who must decide which sets of EOFs best represent the signal to be filtered out.
The SSA filtering is applied on polar data after correction for the main, magnetospheric, and lithospheric fields up to SH degree 45 and after implementing the dedicated along-track SH correction. Synthetic tests and inspection of the residual plots (see Fig. 3) reveal that the contaminating external fields are mostly represented by the first two EOFs in the polar regions. This, we note, is often violated for satellite track segments shorter than about 8° (typically corresponding to the wavelength at SH degree 45) if one is not careful to first remove a lithospheric field model to SH degree 45 in the way we did. This additional SSA along-track correction is applied to all selected data in the high-latitude regions (but not at the mid-latitude regions). Figure 3 shows the measurements after the SH and SSA corrections above the polar regions. Once this series of dedicated correction are made, the lithospheric field model from SH degree 16 to 45 that was first subtracted is added back and the measurements temporarily selected above magnetic latitude |±52°| in the mid-latitude dataset and below magnetic latitude |±52°| in the polar dataset are removed.
Across-track vector and scalar differences are next computed by selecting Swarm A and C data measured at the same universal time (UT) with the additional condition that the distance between the measurements does not exceed 1.4° in great circle distance (this corresponds to the angular distance between Swarm A and C at the equator crossings). These across-track differences approximate the gradient of the measurements. However, we do not divide the difference by the actual distance between the measurements so that vector, scalar, and "gradient data" are all in units of nanotesla. This leads to an inverse problem statistically easier to handle (e.g., Thébault et al. 2013). The across-track difference computed from the Swarm A and C synchronous measurements is not exactly EW oriented and contains a significant amount of NS contributions near the geographic poles where Swarm A and C no longer fly side by side. This construction of the gradients aims at canceling out all large-scale contributions but also very rapid and transient field fluctuations related to external field sources remaining in the data. This differs from the approach chosen by Olsen et al. (2015) where a more exact separation of EW and NS gradients is sought numerically. Their approach does not preclude the contamination of gradients by rapid external field variations (under a few seconds) although it opens the interesting, but not yet fully exploited, possibility to build an optimal combination from all gradient components depending on their sensitivity to the lithospheric field (Kotsiaros and Olsen 2012).
Differences along each satellite A and C tracks that approximate the NS gradients at equator crossing are computed by selecting data measured by the same satellite with a time stamp difference of 20 s. This time difference corresponds to about 1.2° or 140 km (at 460 km altitude) spacing between the measurements. This distance is chosen smaller but close to the angular EW separation of Swarm A and C of about 1.4° at the equator. Kotsiaros et al. (2015) performed a retrospective analysis of the last 2 years of the CHAMP satellite measurements and concluded that computing NS gradients from measurements separated by not more than 30 s is close to optimum. Contrary to the across-track gradients, we found that uncorrected external fields were enhanced at high latitudes when including the satellite A and C alongtrack gradients, particularly in the auroral ovals. For this reason, along-track differences were considered only in mid-latitude regions and not in polar regions.

Data weighting
The scalar, vector, and gradient data we use are statistically different and undergo various types of processing in the algorithm. Defining which weight should be affected to each dataset is therefore not trivial. For example, it is known that a Sun-related perturbation affects the VFM vector measurements on all three Swarm satellites and that an empirical model is currently being used to correct and calibrate the official 0408 and 0409 vector data, using an approach similar to the one initially proposed by Lesur et al. (2015) (see https://earth.esa.int/web/ guest/swarm/data-access/dataset-history and related documents for more details). In addition, because of an anomaly affecting the ASM instrument on Swarm C, no more absolute scalar data are being acquired on that satellite since November 5, 2014, and calibration of Swarm C data must also rely on absolute scalar ASM measurements made on the nearby Swarm A satellite. Scalar measurements are also less sensitive to external field variations and errors due to the satellite attitude uncertainties (Holme and Bloxham 1996). Overall, it is thus fair to state that the error budget of vector measurements must a priori be considered as being larger than that of scalar measurements. On the other hand, it is important to avoid putting too much weight on scalar measurements in order to limit possible Backus effects near the equatorial regions (Backus 1970). More generally, we note that regional geomagnetic field models also can suffer from similar non-uniqueness issues when derived from scalar measurements only. To guarantee uniqueness of the derived models, regional modeling thus also requires vector measurements in all spherical caps, even in polar regions where they are comparatively noisier than in the mid-latitude ones.
Across-track gradient data are less sensitive to largescale external fields, in particular to those responsible for leveling issues between nearby measurements. Therefore, putting more weights on the across-track gradient data seems an attractive option. However, handling the gradient data is not trivial because NS and EW differences of Swarm measurements only provide information about the horizontal and not on the vertical gradient. In addition, since gradient data are mostly sensitive to the small spatial scales of the lithospheric fields, vector and scalar measurements are much needed for constraining the larger-scale lithospheric field structures (Friis- Christensen et al. 2006). For these two reasons, it is not certain that the incomplete gradient data carry enough information to completely constrain the vector magnetic field at every data location. The different experimental and processing errors thus justify weighting each data type differently. We do not try to set realistic variances to each data type, but we rather define weights with the priority order w δF ≥ w F ≥ w δV ≥ w V , with w δF the weight on the scalar gradients, w F the weight on the scalar measurements, w δV the weight on the vector gradients, and w V the weight on the vector measurements, respectively. Setting the numerical values is arbitrary, but taking into account the data quality and theoretical limitations, we choose w δF = 10, w F = 5, w δV = 5, w V = 1. These a priori weights are used in the inverse problem.

Regional inversions
The selected and corrected scalar and vector measurements are iteratively modeled (process 5 in Fig. 1) on regional scales using the revised spherical cap functions (Thébault 2006). Full vector information is not necessary to recover the poloidal field of internal origin, provided the external field signal has been removed correctly (Backus et al. 1996). In principle, it is then sufficient to use either the vertical or the horizontal vector field components. The vertical component, however, is less disturbed by the magnetic field generated by the fieldaligned currents in the polar regions. We therefore avoid using the horizontal component in those regions. In practice, we simply use the full vector and its gradient in combination with the total field intensity in the equatorial band limited by latitudes ±30° and discard the vector horizontal components at poleward latitudes.
The inverse problem is solved independently for 600 spherical caps that are placed on the nodes of an equal area grid on the sphere at 460 km altitude (median altitude of the selected Swarm measurements). All caps overlap in space so that adjacent regional models share a common subset of data. The half-aperture of each cap is 8°, and the maximum horizontal spatial resolution is set to about 350 km (corresponding roughly to SH degree 120). The inverse problem is solved with an iteratively reweighted least-squares (IRLS) method with the a priori weights w δF , w F , w δV , and w V , set to each data type. The weights are then iteratively and numerically reallocated following a hypothesis on the error distribution known as the Huber distribution (Huber 1981). The Huber distribution is chosen so that the error is assumed to follow a Gaussian distribution within one standard error but a Laplace distribution outside. This longer-tailed Laplace distribution minimizes the leverage of aberrant values in the measurements and provides more robust solutions for the local models. For each spherical cap, the algorithm computes the statistical misfit information between the model and each data type (Fig. 4).

Transformation in spherical harmonics
The regional models are finally converted into a unique set of SH Gauss coefficients. We designed two transformations, one using a fast spherical transform (FST in process 6a in Fig. 1) and one relying on a least-squares inversion with regularization (process 6b). In both cases, once the regional inversions are completed for all caps, the regional parameters are used to compute the vertical component of the lithospheric field on the nodes of a grid. The arithmetic mean is considered in the overlapping areas between two adjacent models. This set of values is the one used as input for the SH model computation. For the FST algorithm, the forward problem is computed on the 180,901 nodes of a Gauss-Legendre grid at the median altitude (about 460 km) of the measurements. In principle, and from the sampling theorem on the sphere, this allows the lithospheric field to be estimated up to SH degree 300. This fast track algorithm is the one we used for many dry runs during which we tested the efficiency of the processing chain and adjusted the different corrections and parameters. For the production of the final model, however, we rely on the least-squares inversion. The forward problem using the regional models is computed on the nodes of an equal area grid (also at the median altitude of the measurements, 460 km). The number of points is defined so that a model can be estimated without aliasing at least up to SH degree 150, but for the present study we restrict ourselves to model estimation up to SH degree 90. The leastsquares inversion is regularized starting from SH degree 70. Maus et al. (2006b) proposed an empirical power model for each Gauss coefficient that is proportional to (1 − m 2 /l 2 ), where l is the degree and m the order of the Gauss coefficient. Their regularization scheme leads to low SH orders having more weights than high SH orders. This energy distribution among orders is arbitrary but consistent with the fact that the noise is correlated along the satellite polar orbits, thus involving the sectorial Gauss coefficients being in general noisier than the tesseral and/or zonal terms. We follow the same empirical distribution of power among individual coefficients. However, we also impose that the Lowes-Mauersberger power spectrum of the individual Gauss coefficients should be consistent with the level of the statistical power spectrum proposed by Thébault and Vervelidou (2015). This statistical expression, plotted in Fig. 5, depends on the mean magnetization of the lithosphere, its thickness, and a power law with parameters numerically estimated using a best fit to the power spectrum of the observational NGDC-720 model (Maus 2010).
In practice, the damping of each Gauss coefficient is defined by computing the ratio between its power as estimated without regularization and its power as expected following the statistical form. When this ratio exceeds 5, the coefficient is considered as entirely contaminated by noise and the corresponding diagonal element of the SH design matrix is divided by the corresponding ratio to bring its power down to its expected statistical level (Maus et al. 2006b). The overall misfit between the final regularized DLFI model obtained in this way and the predicted values by the regional models is 0.18 nT.

Results
In the DLFI algorithm, the modeling of the selected and corrected Swarm measurements is performed in two steps. The first step involves the iterative regional modeling over the Earth at the median altitude of the measurements. Table 1 shows the Huber weighted misfit and mean residuals between the data and the regional models for the different data types (at the end of process 5 in Fig. 1). The data misfit is below 1 nT at mid-latitudes and does not exceed 3.1 nT at high latitudes. As expected, the noise level in polar regions is larger for the scalar and vertical field components than for the gradient data. The scalar gradient misfit is particularly low at mid-latitude, thus confirming the good quality of the measurements. The mean residual values are close to zero for all components. The small offset from zero for the gradient of the scalar and the vertical components likely indicates a small difference of baseline between the VFM and ASM onboard the Swarm A and C satellites . The geographic distribution of the RMS in Table 1 is detailed in Fig. 4. The misfits for the four data types are significantly larger over North America and East Antarctica. Interestingly, the misfits of the scalar and vector measurements are similar, but the scalar gradient data have comparatively lower noise in polar regions than the scalar, the vector, and the vector gradient data. This is a clear indication that the general misfit is primarily caused by errors in the correction for the external fields and not by the raw data uncertainty.
The regional models are then converted into a unique set of Gauss coefficients. Assessing the quality and robustness of a SH model or comparing it to previous models is no trivial matter. Here, we rely on statistical criteria defined in the spectral domain such as the shape of the Lowes-Mauersberger power spectrum, the degree correlation, and the azimuthal power spectrum [see Finlay et al. (2010), for definitions]. These quantities are also computed for two analogous models (SIFM, Olsen et al. 2015) and MF5 (Maus et al. 2007) and compared to the MF7 model (Maus 2010) used here as a reference. These analyses are convenient as they make it possible to identify coefficients that most differ among models. The SIFM was derived from 14 months of Swarm measurements Fig. 4 Misfit between each regional model and the measurements for the vector (a), the scalar (b), the gradient of the vertical (c), and the gradient of the scalar components (d) and is to date the only published Swarm model relevant for such comparisons. The MF5 model was derived from single CHAMP satellite data, using a similar number of measurements and at a median altitude comparable to the one of the Swarm measurements considered here. In this respect, it is better suited for comparison purposes than more recent models, such as MF6 (Maus et al. 2008). MF7, finally, is used as a reference because it was derived from data acquired toward the end of the CHAMP mission, at very low altitudes (about 300 km on average), and displays the best lithospheric signal-to-noise ratio (thanks also to geometric amplification of the lithospheric signal with downward altitudes).
The Lowes-Mauersberger power spectra of the DLFI (blue curve in Fig. 5, left) and MF7 (black curve) models are similar up to SH degree 70 (correlation larger than 0.85) and compare visually well to that of the SIFM (red), even though we observe a small offset of the DLFI spectrum with respect to the two other spectra. The degree correlation (Fig. 5, center) shows that the DLFI model better correlates with the MF7 model than the SIFM. This indicates that the difference in amplitude, which could be due to the filtering of the measurements along the satellite tracks, is a systematic feature. The azimuthal power spectra of the three models are also in good visual agreement. Nonetheless, we note that the DLFI model has less power than the SIFM and MF7 models in the sectorial harmonics (m/n = ±1). This is because the damping process we applied to the Gauss coefficients mostly affects these coefficients. Recall, however, that modeling errors are mostly concentrated in these sectorial harmonics and that the MF7 and SIFMs could thus both have too much power in these terms. Based on these three criteria, the SIFM and DLFI models are also in much better agreement with the MF7 model than the MF5 model. This illustrates well the Swarm constellation advantage of approximating the EW and NS magnetic field gradients using the combined Swarm A and C measurements.   Fig. 6. This comparison reveals that at mid-latitudes, the DLFI model misses the smallest scales of the two strongest Kursk and Bangui magnetic anomalies visible from space. The reason is not clear, but this might be an effect either of the along-track correction or of the iterative Huber weighting scheme which penalizes extreme but significant magnetic field values too strongly. Otherwise, the residual map does not reveal any systematic strong or large-scale structures. Most differences are found in polar regions, near the polar electrojets, as expected. This shows that in those regions the DLFI model contains small-scale structures that probably have too much power when the model is downward continued to the Earth's mean surface.

Conclusion
The DLFI model, which uses 22 months of Swarm satellite measurements, suggests that the data collected so far are (globally) sensitive to crustal field variations at least up to degree 80 with some geographic disparities. The lithospheric map we produced compares well with those derived from previous CHAMP and Swarm models. Differences are strongest in the northern and southern polar regions where rapid and significant non-lithospheric contributions in the measurements lead to arbitrary offsets between adjacent tracks that persist despite dedicated corrections and filtering. The perturbation of the small scales is particularly problematic when using regional modeling functions that were purposefully designed to image these structures. At present, a cleaner separation of the lithospheric signal from ionospheric and magnetospheric noise sources cannot be done simply. On the one hand, increasing the data volume is contradictory with the strategy of the DLFI which is to identify and select the least disturbed satellite tracks by relying on drastic selection criteria. On the other hand, applying heavier processing or damping increases the risk of systematic errors. The explicit use of the Swarm gradient configuration substantially improves the situation by allowing lighter corrections to be applied on the vector and scalar data along the satellite track than would have been needed with a single satellite mission. Figures 4  and 6, however, show that the DLFI model should still be regarded with some caution at the Earth's mean radius and auroral latitudes. Future significant improvements can nevertheless be expected from the progressively decaying altitudes of the Swarm satellites, which will increase the signal-to-noise ratio, and from data accumulation at all LTs, which will better average out the contributions of external magnetic field perturbations in polar regions (Maus et al. 2006b). From a methodological viewpoint, we find promising the application of the singular spectrum analysis (SSA) to correct the measurements in the polar regions. In this study, SSA was applied on a track-by-track basis. But patterns common to adjacent satellite paths could be identified by using instead multitaper SSA. However, improvements of the lithospheric field model with datasets similar to the one considered in this paper will remain challenging because the RMS of the crustal field beyond SH degree 80 is <0.1 nT at 460 km altitude (the median altitude of the Swarm measurements after 2 years). This value is comparable to the instruments error budget and is smaller than the data calibration uncertainties . Therefore, any systematic error larger than this threshold is likely to produce a map at the Earth's mean radius with biased structures at short length scales.