Swarm SCARF Dedicated Lithospheric Field Inversion chain

Abstract

The forthcoming Swarm satellite mission is a constellation of three satellites dedicated to the study of the geomagnetic field. The orbital characteristics of the mission, which includes a pair of satellites flying side by side, has prompted new efforts in data processing and modeling. A consortium of several research institutions has been selected by the European Space Agency (ESA) to provide a number of Level-2 data products which will be made available to the scientific community. Within this framework, specific tools have been tailor-made to better recover the lithospheric magnetic field contribution. These tools take advantage of gradient properties measured by the lower pair of Swarm satellites and rely on a regional modeling scheme designed to better detect signatures of small spatial scales. We report on a processing chain specifically designed for the Swarm mission. Using an End-to-End simulation, we show that the tools developed are operational. The chain generates a model that meets the primary scientific objectives of the Swarm mission. We also discuss refinements that could also be implemented during the Swarm operational phase to further improve lithospheric field models and reach unprecedented spatial resolution.

Earth Planets Space, 65, 12571270, 2013

Swarm SCARF Dedicated Lithospheric Field Inversion chain

E. Thbault1, P. Vigneron1, S. Maus2, A. Chulliat1, O. Sirol1, and G. Hulot1

1Equipe de Gomagntisme, Institut de Physique du Globe de Paris, Sorbonne Paris Cit, Universit Paris Diderot, UMR CNRS 7154, 1 rue Jussieu F-75005 Paris, France

2Geophysical Data Center, NOAA, Boulder, Colorado, U.S.A.

(Received March 31, 2013; Revised July 10, 2013; Accepted July 22, 2013; Online published November 22, 2013)

The forthcoming Swarm satellite mission is a constellation of three satellites dedicated to the study of the geomagnetic eld. The orbital characteristics of the mission, which includes a pair of satellites ying side by side, has prompted new efforts in data processing and modeling. A consortium of several research institutions has been selected by the European Space Agency (ESA) to provide a number of Level-2 data products which will be made available to the scientic community. Within this framework, specic tools have been tailor-made to better recover the lithospheric magnetic eld contribution. These tools take advantage of gradient properties measured by the lower pair of Swarm satellites and rely on a regional modeling scheme designed to better detect signatures of small spatial scales. We report on a processing chain specically designed for the Swarm mission. Using an End-to-End simulation, we show that the tools developed are operational. The chain generates a model that meets the primary scientic objectives of the Swarm mission. We also discuss renements that could also be implemented during the Swarm operational phase to further improve lithospheric eld models and reach unprecedented spatial resolution.

Key words: Lithospheric magnetic eld, space magnetometry, Swarm.

1. Introduction

Large magnetic anomaly eld features with spatial dimensions of more than a few hundred kilometers can be inferred from Low Earth Orbiting (LEO) satellites. These length scales are typically dominated by the magnetic eld created by the spatial variation in structure and composition of the Earths crust and upper mantle (Thbault et al., 2010 for a review). The previous CHAMP (Reigber et al., 2002) satellite mission allowed impressive progress towards mapping the Earths lithospheric magnetic eld (e.g., Maus et al., 2008), understanding their geological sources (e.g., Hemant and Maus, 2005), and inferring other physical quantities such as the equivalent global Curie depth (Purucker et al., 2002) or the heat ow anomalies of the Earths crust in remote regions (e.g., Fox-Maule et al., 2005). The Earths lithospheric eld at large wavelengths (>400 km) is conceptually simple because mostly static in time (e.g., Thbault et al., 2009). Unfortunately, it is also one of the weakest detectable sources at satellite altitudes and is masked by all other constituents of the total Earths magnetic eld (e.g., Hulot et al., 2007). For instance, it represents about 0.01% of the total Earths magnetic eld strength at 400 km height (with a root mean square of about3.8 nT using the NGDC-720 model; Maus, 2010) and its magnitude rarely exceeds 20 nT.

Planetary magnetic elds are often expressed in terms of a nite series of spherical harmonics (SH) to a maximum

Copyright c

[circlecopyrt] The Society of Geomagnetism and Earth, Planetary and Space Sci

ences (SGEPSS); The Seismological Society of Japan; The Volcanological Society of Japan; The Geodetic Society of Japan; The Japanese Society for Planetary Sciences; TERRAPUB.

doi:10.5047/eps.2013.07.008

expansion degree. After ten years of the CHAMP mission, the most recent models describing the magnetic eld of the lithosphere, such as the MF7 model (an upgrade of Maus et al., 2008; http://geomag.org/models/MF7.html) expanded to SH degree 133, and the CHAOS-4 model (an upgrade of Olsen et al., 2010a) expanded to SH degree 100, agree reasonably well up to degree 85, thus providing a robust bound of about 470 km on the maximum horizontal spatial resolution. The geometrical attenuation with altitude and the prominence of external elds at satellite level are the two main factors limiting our ability to reach higher SH degrees and thus to detect smaller scale structures. However, for equivalent altitude datasets, the differences between MF7 and CHAOS-4 mainly lie in the way the different models deal with source eld separation and, more specically, with the rapidly time-varying external elds.

Two approaches can indeed be followed for the analysis of magnetic eld measurements from space. The rst one, referred to as the comprehensive inversion (CI) approach, attempts to solve for all known magnetic eld contributions through a grand inversion (e.g., Sabaka and Olsen, 2006; Sabaka et al., 2013; Olsen et al., 2010a). This approach requires some assumptions about the geometry of the dominant source elds to guide the model parameterization. Because of the time needed to properly sample the Earths surface and because of the satellite orbital precession through local time (LT), the CI approach relies on mathematical models describing the sources only on average in space and time. For the recent CHAMP mission, the complete coverage in space and LT was achieved after about 130 days (Olsen et al., 2010b, gure 8). As a result, in such models the external magnetic elds with higher time frequen-

1257

1258 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

cies may overlap with the small scale and weak lithospheric elds so that both contributions can be difcult to identify without ambiguity (e.g., Olsen et al., 2010c).

One may follow a second approach, referred to as the sequential approach, and focus on one specic constituent such as the lithospheric magnetic eld (e.g., Maus et al., 2008). The ambiguity with respect to internal and external eld separation can then be partly removed by adding prior information or by selecting measurements in a way that is more optimum. This is a complex and subjective procedure where the investigators decisions are needed all along the way to isolate the corresponding contributions. This approach successfully led to a suite of lithospheric eld models with unprecedented spatial resolution (SH degree 133 for the MF7 model). Even though specic problems related to spectral leakage can arise (Sabaka and Olsen, 2006; Thbault et al., 2012), implying that all SH coefcients are not necessarily well resolved, these models clearly revealed striking regional correlations between magnetic anomalies and continental geological bodies at small spatial scales.

One of the primary goals of the European Swarm satellite mission (Friis-Christensen et al., 2006) is to unveil the lithospheric eld signal to the highest possible horizontal spatial resolution. We therefore decided to develop a processing chain dedicated to the production of a lithospheric eld model. This so-called Dedicated Lithospheric Field Inversion Chain was proposed within the Swarm SCARF (Satellite Constellation Application and Research Facility, see Olsen et al., 2013) set up by a consortium of several research institutions and selected for implementation by ESA. It will be operated during the mission as one of more than a dozen Level 2 processing chains, which will also include a Comprehensive Inversion chain (producing an independent lithospheric eld model, see Sabaka et al., 2013) and two other chains: the Swarm SCARF Dedicated Ionospheric Field Inversion chain (Chulliat et al., 2013) and the Swarm SCARF equatorial electric eld inversion chain (Alken et al., 2013) developed by IPGP or in collaboration with NOAA.

Whichever the chosen modeling approach, the general idea is to carefully represent the selected and corrected data by a single set of parameters that allows computing the lithospheric magnetic eld values outside the dominion of data. The common approach relies on the SH series expansion that is the most natural solution of the Laplaces equation in the geocentric reference frame. Here we favor a regional modeling scheme which has several advantages. First, it makes it possible to bypass the problem of truncated spherical harmonics that can arguably be suboptimal for the purpose of highlighting small magnetic eld wavelengths. Indeed, splitting the global problem into subregions allows high regional spatial resolution to be reached with a manageable number of parameters. Second, correction errors result mostly from transient external magnetic elds that contaminate the data in a different way depending on the location. Regional modeling can adjust to this by treating the noise regionally; it also offers the possibility of constraining the inverse problem depending on how well the regional model parameters are resolved (e.g., Lesur and Maus, 2006). Finally, computing regionally the covariance

matrices of the errors for a better understanding of their correlation in space and in time is also possible, something that is demanding in terms of computer resources when building high degree SH models.

We report below on the Dedicated Lithospheric Field Chain as it was implemented during the preparation phase of the Swarm satellite mission. The block processes were duly tested in the course of the development and we do not discuss the details of the closed-loop simulations that were successfully completed. Rather, we focus on the end-toend simulation and describe in some details the algorithm, its advantages but also its shortcomings. Because the dedicated techniques rely on the operators judgment rather than unequivocal scientic observations (which often raises some contentious issues), we describe the difculties we met and the compromises we decided to make. We nally show that the chain is able to produce a lithospheric eld model that reaches the Swarm performance criteria (Swarm Level 2 Processing System Consortium, 2013).

2. Theoretical Background

Different regional modelling techniques are currently used or under development in the eld of geomagnetism (see Schott and Thbault, 2011, for a review). These are based, for instance, on frames of Poisson wavelets (e.g., Maier and Mayer, 2003), slepian functions (e.g., Beggan et al., 2013), combination of band-limited spherical harmonics (e.g., Lesur, 2006), spherical tesselation (e.g., Stockman et al., 2009), or equivalent dipole source distributions (e.g., Langlais et al., 2004). We choose to rely on the Revised Spherical Cap Harmonic Analysis (R-SCHA, Thbault et al., 2006) technique as this approach was successfully used with the CHAMP magnetic data to build a global model of the Earths lithospheric eld (Thbault, 2006). We briey restate the basic properties of the modelling strategy we follow. A complete description of the Revised Spherical Cap Harmonic Analysis technique can otherwise be found in Thbault et al. (2004, 2006).

The solution is obtained within a conical domain bounded by two spherical caps of aperture 0 at altitudes 1 and 2; the surface being closed by the lateral boundary of the cone (see gure 1 in Thbault et al., 2006). The solution for the potential VSCHA can be expressed in a spherical coordinate system as a sum of two innite series

VSCHA(r, , ) = a

m=0

[parenleftBigg]

k=1 Gmkmk(r, , )

p=1 Gmp mp(r, , )

+

, (1)

where mk = Rk(r)eim Pmnk(cos ) are the complex spheri

cal cap harmonics of order m involving the Legendre function Pmnk(cos ) of real nk (see Haines, 1985) and Rk(r) their radial functions, while mp = Rp(r)eim K mp(cos ) are the

complex harmonics expressed in terms of the conical, or Mehler, functions and Rp(r) the real-valued radial functions (see Thbault et al., 2006, equations (18)(19)), a is the Earths mean radius and , , and r are respectively the colatitude, the longitude and the radius (in a reference frame

E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN 1259

whose axis is aligned with the axis of the cone), and Gm. the regional parameters for each order m written in their complex form. There are many ways of numerically dening the basis functions, one of which being to impose homogeneous boundary conditions on each of the three surfaces.In the present work, we apply the formalism described in Thbault et al. (2004), where the boundary value problem is solved with Neumann homogeneous boundary conditions on all three surfaces. In practice, the complete solution requires a second potential that is added to the one dened in Eq. (1)

VSH(r, , ) = a

n=1

1

2, (4)

which may be expressed in terms of the minimum length scale of the magnetic signal (Backus et al., 1996, section3.6.3)

2(a + )

nk + 1/2

a r

n+1 G0n0n, (2)

with 0n = P0n(cos ) the real-valued zonal spherical har

monic of integer degree n, the axis of which coincides with the cones axis, and G0n the zonal SH Gauss coefcients in the cones reference frame. This second potential is required to bypass a problem of incompleteness of the potential expressed in Eq. (1) when the geomagnetic eld is to be represented (see Thbault et al., 2004). The solution for the magnetic eld B is then obtained considering the gradient of the two potentials

B = (VSCHA + VSH) . (3)

Empirically, we observe that a solution dened by the Neumann boundary conditions converges fast for the vector magnetic eld components. This is a desired property when dealing with anomaly elds still possibly containing large spatial scales magnetic eld contaminations, which is usually the case when processing real satellite measurements or realistic synthetic Swarm datasets built from SH models. The spherical and spherical cap harmonics both allow an exact representation of any poloidal magnetic eld so that the two solutions converge in norm to each other. However, some representation errors invariably occur when the innite series (Eqs. (1) and (2)) are truncated to nite expansion indices. A regional solution converging fast can represent the magnetic eld B with fewer parameters than a slow converging series expansion. It guarantees a better comparison with the benchmark model during the synthetic validation of the algorithm as it mitigates the problem of truncation errors. There are some disadvantages, though, of combining the Eq. (1) with (2). In particular, the basis functions are not mutually orthogonal within the conical domain (Thbault et al., 2004). This leads to a non-diagonal design matrix where cross-talking between parameters implies some reorganization of their estimated values when changing the series truncation. In general, this is not ideal for carrying out spectral analysis of the magnetic eld signal or for isolating a single model parameter. For the purpose of this work, this is a minor issue because we do not aim at interpreting the values of each regional parameter independently and because we only deal with high spatial density satellite data. In such a case, the design matrix, although not diagonal, is never dense. The inverse problem is numerically stable.

Each type of functions (Mehler or Legendre) has specic properties (Thbault et al., 2006). In particular, the

maximum expansion index kmax of series in Eq. (1) can be asymptotically related to the degree nk of the Legendre functions (Haines, 1985; Thbault et al., 2006)

nk (2kmax m 1/2)

180

20

, (5)

where is the considered altitude. For kmax = 15 and

0 = 8, for example, we nd 130 km at 400 km alti

tude; these parameters are the ones selected for the present regional modeling.

3. The Dedicated Lithospheric Field Inversion (DLFI) Chain

The full chain is represented in Fig. 1. The Swarm Level 1b data are rst selected according to specic criteria and are corrected for various sources elds. As already discussed, this step is crucial for deriving the lithospheric model. It involves four processes (Process 1 to Process 4 in Fig. 2) that require auxiliary information and operators decisions. The selected and corrected measurements are then iteratively processed with the R-SCHA modeling. The data are distributed into geographical bins (Process. 5.1) corresponding to each spherical cap covering the Earths surface and transformed into the caps reference frame (Process 5.2). The inverse problem is solved, and the regional parameters stored, after building the design matrix (Process6.1) and solving the system by Singular Value Decomposition (SVD in Process 6.2). Once the regional inversions are completed for all caps, the regional parameters are used to predict the lithospheric vector eld values on the nodes of a Gauss-Legendre grid (Process 7) covering the Earths surface at a constant altitude. These grid values are nally converted into a unique set of spherical harmonics Gauss coefcients either by a fast spherical transform of the vertical magnetic eld component or by a linear inversion of the horizontal and vertical components (Process 8). The dedicated lithospheric eld model is the output product of the sequence of these numerical processes. Each process is described in more details in the following.3.1 Data selection

For this simulation, we use only the vector Level 1b data simulated for the lowest pair of satellites Swarm A and Swarm B (Olsen et al., 2013) since the third satellite, Swarm C, is expected to be too high to detect the ne lithospheric eld structures. Nevertheless, its data could prove useful for the chain in case of a failure of Swarm A or B and all data selection and correction procedures described below were also made operational for processing Swarm C magnetic eld measurements. Considering that vector magnetometers may suffer from malfunctions, the chain was also adapted to ingest scalar magnetic measurements if required.

More than four years of simulated vector Swarm A and B data are selected between a simulated time period extend-

1260 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

Fig. 1. General owchart of the six main processes run to generate a lithospheric eld model from the Level 1b magnetic eld measurements of the simulated Swarm satellite mission.

Fig. 2. This ow chart provides the sequence of the processing steps allowing to isolate the magnetic eld contributions from the Earths lithosphere.

E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN 1261

ing from July 1998 to December 2002. These data are sub-sampled every 10 seconds, corresponding to a spacing of about 70 km along the satellite track at 400 km altitude.This ensures that no aliasing occurs because of magnetic eld contributions at small wavelengths (recall that the minimum spatial wavelength detected by the regional model is 130 km with the current setting). Although the synthetic data do not contain toroidal external elds from eld aligned currents, we account for the different behavior of the external eld with latitudes by separating the vector data in two sets, one for mid-latitudes (52 to 52 magnetic

latitudes), and one for polar regions (>52 and <52). In

addition, we allow an overlap of 10 between mid and polar latitude regions (see also Maus et al., 2008, for instance) so that the provisional selected data for both sets range between 57 to 57 magnetic latitude for the mid-latitude set

and above |47| for the polar set. As discussed below, this

overlap is introduced to reduce correction errors caused by dedicated corrections applied along the satellite tracks. The mid-latitude data are further selected at 21:00 to 5:00 local time (LT) in order to minimize the contributions from the ionospheric Sq eld. All LT are considered for the polar data, provided they are outside the sunlit area (10 below

horizon). We use the Dst index to lower the effect of the magnetic eld generated by the ring current and select the data corresponding to Dst values lower than 20 nT and to

Dst variations smaller than 10 nT over the previous three hours. We also apply a selection criterion based on the Kp index to keep data corresponding to a global low level of geomagnetic activity. The data corresponding to a Kp index larger than 2, and to Kp variations larger than 2 over the previous three hours, are rejected. For this study, the selected Dst threshold is relatively large. We argue that rejecting many satellite data because of a high level of magnetic activity does not necessarily reduce the noise in the nal model, as the data error variance is proportional to the ratio v/M, where v is the variance of the external eld and M the number of available orbits over a given longitude (e.g., Lesur et al., 2013). It is important to dene an acceptable trade-off between the number of available satellite tracks and the level of external eld activity, more particularly when processing the data on regional scales because the local functions require a dense data coverage in space.

The selection procedure we applied is not exhaustive.

The process 1.2 (see Fig. 2) contains more data selection modules than the ones activated for this simulation.Among other selection criteria we included a processing block based on the measure of the Interplanetary Magnetic Field (IMF), whose value is related to the solar wind activity. We included also a selection block based on the ap index, which measures the general magnetic activity at the planetary scale, and a selection block based on quality ags such as the on/off status of the star camera. Other selection blocks are likely to be added as more and more experience is gained on the Swarm satellite data during the calibration/validation and operational phases.3.2 Data correction for the main and external elds

We developed two families of magnetic eld corrections for the selected vector data. The rst family relies on explicit models for the main eld, its secular variation, and the

magnetospheric eld, (processes 2 and 3 in Fig. 2), all expanded in spherical harmonics. These input models are provided by the other Swarm processing chains (see Olsen et al., 2013). Note that the data are not corrected for the ionospheric eld, which is assumed minimum at night-times (see Chulliat et al., 2013). The second family of correction is based on more empirical processes, performed either at the global scale, or along the satellite orbits (process 4). The settings of this suite of corrections depend on the analysis of the data residual after the completion of processes 2 and 3.

Within the overall process proposed by the consortium (see Olsen et al., 2013), the dedicated lithospheric eld chain is initiated as soon as the models for other source elds are available for use. The main eld correction can be performed using either the output of the comprehensive inversion chain (Sabaka et al., 2013) or that of the dedicated core eld chain (Rother et al., 2013). In principle, there is a better consistency between the data selection detailed in the previous section and the data selection applied in the dedicated chain for the core eld (Lesur et al., 2010; Rother et al., 2013), particularly at night times. This consistency is important if one keeps in mind that differences in data selection are acknowledged to be one of the major sources of discrepancies between main eld models yet parameterized in similar ways (e.g., Finlay et al., 2010). In the course of the Swarm preparation phase, we tested the inuence of each of the input main eld models and the results were similar. In this paper, the solution is obtained with the vector data corrected with the MCO SHAi2D (the output of the core eld dedicated inversion, see Fig. 2) for the main eld model and its secular variation up to SH degree 15. Since the synthetic Swarm vector measurements were generated using a secular variation model up to SH degree 20 (Olsen et al., 2013), the secular variation is only partly corrected for. The corrections for the external magnetospheric sources are made using the MMA SHAi2C model (Sabaka et al., 2013; see also Fig. 2). This model provides coefcients for the primary magnetospheric eld to SH degree 3 and a maximum order 1 with a time step of 1.5 hours. The induced counterpart is also dened every 1.5 hours but to maximum degree and order 5. This space and time resolution is lower compared to the benchmark magnetospheric eld model used to generate the synthetic data. This represents another source of signicant correction error.

Figure 3 displays the residual eld for the radial component after correcting explicitly for the main and magnetospheric elds (the known synthetic lithospheric eld are removed in this picture for the purpose of illustration). Thanks to the availability of the benchmark model used to compute the synthetic data, we know that these large-scale residuals are dominated by residual magnetospheric contributions. Another important source of errors comes from the ionospheric eld at night-times (particularly its induced counterpart), which was not corrected. The secular variation error of the main eld is of secondary importance. In total, the correction errors accumulated using SCARF scientic products reach almost 20 nT in two lobes located

near the South geographic pole.

We next identify the Swarm A and B data with the com-

1262 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

Fig. 3. Radial component of the residual eld after correction of the Swarm A and B simulated measurements for a main and a magnetospheric eld model. The known lithospheric eld is removed for this illustration. Projection is Mollweide.

mon reading times and built up a set of so-called gradient data. The terminology of gradient is misleading here as it refers to a point to point difference between the measured data across the two orbits of the lowest pair of satellites (Swarm A and B). This difference is expressed in the same units as the vector measurements (in nT), which avoids the unnecessarily complication of solving the inverse problem with different data types (the true gradients are in nT/km or nT/degree). Considering the nature of the residuals, and their variations in space and in time, taking the difference between the data from the lower pair of satellites in the East-West direction, spaced about 1.4 degree apart (at the geographic equator, see Friis-Christensen et al., 2006), appears to be an efcient way to mitigate the contributions of large-scale non-lithospheric elds. This gradient however is incomplete since the across difference between the two satellite measurements approximates only the three gradient components along the East-West direction. There is no redundancy with the vector measurements that are still needed to recover the lithospheric eld and to stabilize the inverse problem.

The gradient data are not corrected further and we remove the duplicate gradient data that were introduced by considering an overlap of 10 between the mid and the polar areas. In contrast, three additional dedicated corrections are applied to the vector data. Dedicated corrections applied only to portions of the sphere are necessary to lter out the large-scale residuals shown in Fig. 3. However, such along-track corrections require some care as they suffer from spectral leakage due to the impossible separation of the various source elds (see Sabaka and Olsen, 2006; Olsen et al., 2010c). This difculty is illustrated, for instance, by the expression of the zonal external eld coefcient

q0l that would be estimated from the magnetic eld vector data measured

only along a complete meridian and that would then be used to correct the data along this meridian. It can be shown that (Thbault et al., 2012)

q0l(, r) = q0l +

1

2l + 1

n=1n

m=0bm,ln(, r), (6)

where q0l is the true external coefcient, and bm,ln is a combination of all internal and external residual elds (expressed in SH) that contribute to the magnetic eld values along the meridian at longitude and radius r. The use of the coefcient

q0l estimated along a single orbit for correcting the magnetic eld measurements for the zonal external eld is inaccurate and risky depending on the magnitude of bm,ln(, r). A reduction of this bias is achievable by performing a maximum of prior analysis over the entire globe in order to minimize all internal and external eld contributions along all satellite orbits. The along-track correction on the residual data is then applied as a nal step.

The rst additional block process aims at correcting globally for the large spatial scale error of reduction, be it of internal or external origin. The internal and external Gauss coefcients are estimated to SH degree 8 in the inclined dipole reference frame with a one year knot-spacing over the full simulated mission duration. We also scale the measurements by the F10.7 solar radio ux values in order to account for the dependence of the ionospheric eld on the solar activity. The improvement brought by this scaling is rather marginal but reduces the offsets between close satellite encounters. Note that the internal and external Gauss coefcients estimated by the dedicated corrections have little scientic value as they contain numerous sources of magnetic elds. They are not intended to be interpreted but they remain useful for better understanding the origin of the

E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN 1263

large scale eld residuals. There is no major technical limit for this processing block to accept more sophisticated parameterization or improved time resolution. For instance, it could be ne-tuned to absorb most of the remaining large scale magnetic elds in case we lack the explicit models for the main and the magnetospheric elds described above.

A second dedicated correction in spherical harmonics is next applied along the satellite orbits separately for the mid and the polar latitude data. This along-track correction is an important step (e.g., Maus et al., 2006) because the external eld perturbations that contain relatively long wavelengths along a satellite track (during magnetically quiet times) change from one orbit to the next. It is worth stressing that some errors are due to uncertainties and changes in the Dst baseline (see Olsen (2002) and Thomson and Lesur (2007) for instance). These introduce lateral spatial oscillations that overlap with the high degree lithospheric eld and distort mostly the sectoral harmonics (degrees n = m).

Fortunately, one can further lessen such bias along satellite tracks (Eq. (6)) by applying the correction on large portions of orbits (Thbault et al., 2012). Therefore, all satellite tracks containing segments longer than about 2000 km without data are rejected. Then, an auxiliary lithospheric eld model (AUX LIT 2 in Fig. 2; taken here from the series of MF models, e.g., Maus et al., 2008) is subtracted from the remaining data. These precautions mitigate the problem of spectral leakage and avoid ltering out too much genuine lithospheric eld signals (although this is generally not guaranteed; see gure 5 in Maus et al., 2006). In a synthetic simulation, Sabaka and Olsen (2006; their gures 10 and 11) illustrated the difculty of recovering all magnetospheric eld coefcients up to degree 3 and order 1 at a high time frequency. We exploit this a priori information for tailoring the along-track correction module and we estimate the coefcients of an external eld to degree 3 and order 1 plus its induced internal part along each track using simultaneously Swarm A and B measurements. Note that the co-estimation for Swarm A and B represents a signicant improvement over estimating the correction along Swarm A and B separately (the bias bm,ln(, r) being longitude dependant in Eq. (6), it does not leak in a similar way along the Swarm A and B orbits). This estimated external eld is then applied as a second correction.

The nal dedicated correction is performed via a Fast Fourier Transform along the satellite tracks. The technique of high-pass ltering along the tracks dates back to the 1980s (Langel and Hinze, 1998, section 4.4). The inconvenience of this procedure is that there is no clear cutoff frequency in the signal to isolate the contributions of the lithospheric eld without ambiguities (see for instance Wang, 1987 for a discussion). Also, it can only be applied independently on the three magnetic eld components. This has much less physical justication than methods based on SH representation. However, the along-track Fourier ltering has interesting convergence properties. When the data are properly tapered, it signicantly lters the remaining intermediate wavelengths while producing little edge effect.We use this processing block to further reduce the magnetic eld contributions with length-scales larger than 10 in arc-length along the satellite orbit using a Butterworth lter

with a sharp cutoff (lter of order 20). The overlapping data at the junction between mid and polar areas are then used to conne possible edge effects outside the effective data geographic intervals. It is worth stressing that this processing block must be used with extreme care as it cannot discriminate between perturbation and lithospheric elds. In particular, the cut-off frequency must be adjusted depending on how well the AUX LIT 2 model values temporarily removed from the data agree with the measurements in the ltered waveband.

Finally, the redundant data that were kept to improve the continuity between the processing in mid and polar latitudes (see Subsection 3.1) are removed from the dataset. Then, the standard deviation and the mean of the residuals are computed for each vector component. All values outside the 3 range are considered as outliers and thus rejected and the auxiliary lithospheric eld model values subtracted before the series of dedicated corrections are added back.

Figure 4 shows the residual map for the vertical component after all corrections (where, again, the lithospheric eld contributions have been subtracted for the purpose of illustration). The processing applied to the vector data is heavy but successful (compare Figs. 3 and 4). However, the corrected measurements now contain a certain amount of artifacts caused by the ad hoc corrections we applied. This shows how important the gradient data are to balance the correction errors made when correcting the vector data for other contributions than those from the lithospheric eld. Still, the corrected data can be assumed to reect mainly the lithospheric eld contributions and are now ready to be used for modelling purposes.3.3 Regional iterative modelling of the pre-processed data

The selected and corrected vector data, with their East-West gradients, are processed regionally (Processes 5 and 6 in Fig. 1). The rst step consists in distributing the data within each of the 600 spherical cones that patch the Earth at the Swarm A and B altitudes. Then, the inverse problem is solved for each region. The cones have 8 aperture, are centered on the nodes of an equal area grid (Saff and Kuijlaars, 1997) about 8.5 spaced apart, and therefore overlap with each other. The lower and upper radii are automatically dened by the minimum and maximum altitudes of the available data within each cap. The estimated maximum likelihood regional parameters of each type Gm. (see

Eqs. (1) and (2)) are stored for each region. The maximum horizontal spatial resolution is about 130 km (see Eq. (5) with kmax = 15).

The R-SCHA technique, initially set up for processing the CHAMP satellite data (Thbault, 2006), was further developed to ingest the Swarm satellite vector and scalar data. This required the implementation of specic tools to deal also with the magnetic eld satellite track differences (referred to as gradient data), with the inversion procedure, its weighting, and its regularization.

Indeed, the gradient and vector data have (here) the same units but they do not have the same statistics. We implemented the possibility of incorporating data error covariance in the inversion. Formally, the noise distribution in the gradient data results from the addition (or subtraction) of

1264 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

Fig. 4. Vertical component of the residuals after the completion of process 4. The benchmark lithospheric eld values have been subtracted for the purpose of this illustration. Mollweide projection.

Fig. 5. Vertical component of the difference between the gradient (estimated by nite difference between Swarm A and Swarm B data, see the text for details) that would be estimated using lithosphere only data and the gradient obtained from the data imperfectly corrected for the main internal and external source elds. Projection is Mollweide.

noise in the Swarm A and B measurements. If this noise is Gaussian with a standard deviation of v, the gradient noise would also be Gaussian with a standard deviation of g = 2 v. For the synthetic simulation, the level of in

strumental noise added to the synthetic measurements was low enough to neglect this statistical difference. However, it illustrates why we should not expect noise in the gradient

and vector data to behave similarly. In addition, the longitudinal distance of 1.4 between the parallel orbits of Swarm

A and B was devised to maximize the gain of the gradient data for SH degrees 120130 (gure 4 in Friis-Christensen et al., 2006). The counterpart is that they have a low gain for low SH degrees. But since the lithospheric eld signal measured from space is dominated by the low SH degree contri-

E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN 1265

Fig. 6. Geographical distribution of the mist after completion of the block processing dedicated to the regional inversions. Mollweide projection.

Fig. 7. Spherical Cap Correlation Analysis. The correlation coefcient nk (Eq. (9)) between the regional parameters estimated with lithosphere magnetic eld data only and the parameters estimated with all source measurements for each spherical cap degree nk is shown for each band of geographic latitude.

butions (large spatial scales), the root mean square (RMS) of the gradient data is reduced compared to the RMS of the vector data. In terms of the inverse problem, this is an important statement since the functional to be minimized in a least-squares sense is

2 =

Nv

n=1B2 + Nv/2 n=1

B

2 , (7)

where B and B are the differences between the re

gional model and the vector and gradient data, respectively (in the east-west direction), and is a real positive parameter. The parameter can be tuned to increase the rela-

tive weight of gradient data in order to balance the fact that their RMS is lower than the RMS of the vector data and that their number is on average half the number of vector data. However, one must also keep in mind that in principle the vector measurements carry all the lithospheric eld information whereas the gradient, being incomplete (we only estimate the East-West gradient), do not. In a thorough prospective work investigating the benets of gradiometry from space, Kotsiaros and Olsen (2012) indeed illustrate that the East-West gradient data poorly constrain the low SH orders m at all SH degree n (gure 2 in Kotsiaros and Olsen, 2012) showing that vector data are essential to constrain the model. Too much down-weighting of the vector

1266 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

data can lead to an ill-posed inverse problem. From a practical viewpoint, it can be seen in Fig. 5 that the corrections for the explicit main and external eld models being imperfect, the gradient of the lithospheric eld is also affected by correction errors (note the similarity between the structures seen in Figs. 3 and 5). These correction errors are different from those seen in the vector data (Fig. 4). It thus appears that both corrected datasets have different and maybe contradictory information contents. The issue, difcult to address, is to assess which of the contamination affecting the gradient data, and of the reduction errors affecting the vector data, is the largest.

These formal and experimental differences will need to be considered carefully on a case by case basis when the Swarm mission will be in operation. For these reasons, the block processes 6.1 and 6.2 (Fig. 1) solve the least-squares inverse problem independently for each spherical cap

G =

Kmax

m=0

Gmk

Gmk + Hmk

Hmk

[parenrightbig]

=

[braceleftBigg][bracketleftBigg]K

[bracketrightBigg][bracketleftBigg]K

max

m=0

Gmk

[parenleftbig][tildewide]Gm

max

m=0

2

2+

Hmk

[parenleftbig] [tildewide]

k

2+

Hmk

2

[bracketrightBigg][bracerightBigg]1/2

(9)

Gmk, Hmk[parenrightbig]the estimated and true local parameters of the Legendre basis function. A similar correlation analysis could be performed with the two other sets of coefcients (The Mehler and the SH ones, see Eqs.(1) and (2)) and, strictly speaking, the quality information provided by this indicator is incomplete. However, only the coefcients of the Legendre basis functions can be related to an equivalent SH degree nk (Eq. (4)) and to a minimum horizontal spatial wavelength (Eq. (5)). The true local parameters are obtained by inverting perfect lithospheric magnetic eld measurements containing SH contributions to degree 180. Figure 7 shows the degree per degree correlation analysis sorted in geographic latitude. The correlation is often better than 0.8 up to equivalent SH degree ranging between 130 and 140 but it degrades more rapidly in the region where mid and high-latitude data were separated (between |40| and |60| latitude). Note also that the

correlations are in general lower for the R-SCHA degrees between 10 and 30 as a result of the large-scale residual external elds overlapping with the lithospheric eld (these large scale elds are visible in the gradient data, see Fig. 5). The message conveyed here is that the regional parameters could be in some places, mostly at mid-latitudes, reliably estimated to a high equivalent spatial degree nk.

We nally convert the regional models into a set of Gauss SH coefcients. The most efcient way to do this is to exploit the sampling theorem on the Earths sphere (see for instance Schaeffer, 2013 and references therein). We rst compute the forward problem for each regional model on the nodes of a Gauss-Legendre grid in the geocentric reference frame (see process 7 in Fig. 1) at 400 km altitude (this corresponds to the median altitude of the considered dataset). This grid satises the sampling theorem to SH degree 300. To ensure that the numerical solution is not affected too much by the lack of measurements in the polar gap, we rely on a nal module. Indeed, large errors, or lack of data, near the geographic pole have a clear adverse effect on the recovery of some constituents of the lithospheric eld. The zonal SH coefcients (of order m = 0)

are mostly constrained by the magnetic eld vector components (mostly the horizontal ones) at high latitudes. The Swarm orbital conguration tested in this study has a polar gap of about 5.2. This represents a diameter of about 615 km at 400 km altitude. Discontinuities caused by the polar gaps introduce ringing when local functions are used. This problem is addressed using a post-processing based on an iterative regional modelling of the vector values computed on the Gauss-Legendre grid above the two polar gaps. The vector data outside the Northern (respectively the Southern) polar gap are rst used to derive a low spatial resolution

[parenleftbig][tildewide]Gm

with

k ,

Hmk

[parenrightbig]

and

AT C1dA + R

[parenrightbig]

AT d (8)

where G is the column matrix of the regional parameters to be recovered (Eqs. (1) and (2)) A the matrix of basis functions for the vector and gradient data, C1d the covariance matrix of the vector and gradient data, d the vector of magnetic eld measurements and their gradients, and R a regularization matrix multiplied by a scalar Lagrange parameter . Finding the best trade-off between mitigating the data noise, rejecting the correction errors, and exploiting the information content of each dataset is not trivial. We did not try to tackle this problem in the present simulation. We set the covariance matrix to identity and neglect the fact the gradient and the vector data errors could be correlated. The regularization matrix R can take various forms, depending on the a priori information one would like to exploit regionally, but was left empty. A solution could also be reached in a L1-norm sense by applying an iteratively re-weighted least-squares algorithm with Huber weights in order to obtain a robust solution by iteratively weighting the effect of outlier data in the parameter estimation. We performed several tests and came to the conclusion that this option should be handled with care at the regional scale when processing lithospheric magnetic eld measurements. Applying one iteration appears to be sufcient in many cases to reach an optimum solution while not smoothing genuine lithospheric eld signal too much. Figure 6 shows the geographical distribution of the residual mean square after completion of the dedicated regional inversion block (Fig. 1). As can be seen most errors are to be found near the polar gaps and around the transitions between mid and polar latitudes that introduce discontinuities in the magnetic eld values.3.4 Conversion to a SH model and assessment

For comparison with other lithospheric eld models and a more stable continuation of the predictions to the Earths surface, it is convenient to represent the nal lithospheric eld model in the form of a set of SH Gauss coefcients. Before we convert the local models into a set of global Gauss coefcients, it is interesting to analyze statistically how close the estimated regional models are to their expected values. By analogy with the spherical harmonics, the degree correlation between two spherical cap harmonic analyses can be estimated by (e.g., Langel and Hinze, 1998;

section 4.5.2)

nk

E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN 1267

Fig. 8. Residual map of the vertical component between the true and the estimated lithospheric eld model at the Earths mean radius. The models are calculated to SH degree 150. Projection is Mollweide.

regional model (i.e., with a minimum spatial wavelength larger than the size of the polar gap). Synthetic data are then generated from this low resolution model on the Gauss-Legendre points located within the Northern (respectively the Southern) polar gap areas and a new set of model parameters is estimated at the next iteration to a higher spatial resolution. The process stops when the spatial resolution of the model over the polar gaps reaches the same resolution as in other regions. This iterative procedure could have been implemented directly in the block dedicated to the regional inversions (Fig. 1). However, implementing it as a post-processing step offers a better control because the ratio between the number of gridded data inside and outside the polar gaps is then the same from one simulation to the other.

Finally, we cast the spatial vector data computed on the Gauss-Legendre grid into the SH using a fast spherical transform of the vertical component (process 8 in Fig. 1). This is generally done up to SH degree 300 to provide a theoretical spatial resolution for the model in agreement with the one provided by the local basis functions. The transformation into SH Gauss coefcients could also be done by a slower module that computes the problem in a least-squares sense. This process would allow considering the three vector eld components, some global scale regularization to be introduced if needed, and to double check the model obtained by the fast spherical transform. The ofcial and nal results will then be provided as a set of Gauss coefcients up to SH degree 150 (MLI SHA 2D, see Fig. 1). However, recognizing that some regions are better represented than others (see Fig. 6), a similar set up to SH 200 (MLI SHA 2E) will also be provided. This latter model, however, should not be downward continued to the Earths surface unless otherwise stated (for this reason, we call it a map in Fig. 1).

The performance of the lithospheric eld model expressed in SH can be evaluated through four different criteria dened during the preparation of the Swarm mission (Swarm Level 2 Processing System Consortium, 2013). The rst one shown in Fig. 8, relies on a geographic map showing the residuals between the true and the estimated lithospheric eld models to SH degree 150 at the Earths mean radius. As can be seen most of the errors are found near the polar gaps and at the mid to high-latitude transitions. More quantitatively, the performance can be estimated using the three complementary criteria in the spectral domain (Fig. 9-left): the spatial power spectrum (Lowes, 1966) of the difference between the true and the estimated lithospheric eld model, the spherical harmonic correlation analysis between both models as dened by Eq. (9) for integer degrees n (Fig. 9-middle), and the azimuthal power spectra that provide information about the energy distribution along the orders m (Fig. 9-right). The power spectra and the correlation analysis are presented with the lithospheric eld model estimated prior and after the implementation of the post-processing in the polar gaps. As can be seen, this simple expedient signicantly improves the model around the order m = 0. Finally, the resolution ma

trix between the true and the model estimated with the post-processing for the spectral gaps can also be computed. It tells us that the error is not distributed along a preferred direction except along order m = 0 (Fig. 10). This is no sur

prise since neither the vector nor the gradient data are able to constrain these spherical harmonic contributions (see the discussion above). The synthesis of these quality indicators suggests that the model was reliably recovered at least up to SH degree 133, which was the target of this Swarm end-to-end simulation.

1268 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

Fig. 9. Power spectral analysis and degree coherence between the true and the estimated lithospheric eld models. Left: power spectrum of the true model (black), the estimated model without polar gap treatment (red), and the estimated model with the polar gap processing (blue). The spatial power spectrum of the difference between both types of models and the true one are also shown. Middle: Spherical Harmonic Correlation Analysis between the true model and the model estimated with (blue) and without (red) polar gap post-processing. The dashed vertical line indicate SH degree 133. Right: Azimuthal power spectrum of the true model (black) and the model estimated with (blue) and without (red) polar post-processing.

Fig. 10. Sensitivity matrix between the true model and the model estimated after a polar gap processing.

4. Conclusions

In this paper we presented the Swarm SCARF Dedicated Lithospheric Field Inversion chain, which will be used to process the magnetic data of the forthcoming Swarm satellite mission and to produce a model for the Earths lithospheric eld. This approach differs signicantly from the comprehensive approach (Sabaka et al., 2013) as it heavily relies on data selection, corrections, and visual inspections. It also uses a modeling scheme based on local functions that were designed to better detect the small scales of the lithospheric eld. During the operational phase, the criteria discussed in this paper will be subject to changes and will therefore be described in a companion report of the Swarm SCARF dedicated lithospheric eld model. With the current setting, the algorithm provides a model that meets the

required performance at least up to SH degree 133 when considering and the synthetic data provided to us for the full mission duration. This is a better spatial resolution than achieved so far by single satellite missions at comparable altitudes. This result conrms that we are able to develop strategies for recovering the lithospheric magnetic eld to a high spatial resolution in a favorable situation when high quality vector magnetic eld measurements are available for the two lower Swarm A and B satellites. This simulation was realistic in the sense that it was realized with a mathematical method different from the one used to synthesize the simulated data, with a separation into regions that inevitably introduces discontinuities, with imperfect main and external eld models for correcting the synthetic measurements for non-lithospheric magnetic elds, and with real

E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN 1269

proxies for selecting the data for quiet magnetic days along the entire simulated mission duration. However, in the operational phase, the performance of this chain will depend on our ability to estimate the real quality of the input magnetic eld measurements and the accuracy of the scientic products used to correct for the main and the external elds.For instance, complementary simulations (not shown) performed with the same parameters as the ones described in this paper suggest that a small bias of 0.3 nT affecting intermittently Swarm A or B vector measurements, would degrade the overall performance down to about SH degree 110120. It is worth stressing that even this degraded performance remains impressive compared to those achieved by single satellite missions, which based on the experience with CHAMP are of the order of SH degree 85 (this it the maximum SH degree for which models start being in signicant disagreement). An improvement from n = 85 to

n = 120 means that twice as many SH Gauss coefcients

would be recovered (the number of Gauss coefcients is n(n + 2)). This spatial resolution limit may not be a hard

bound but it shows that the recipe we followed may need to be revised and ne-tuned when real data will become available.

The Swarm calibration/validation phase will be a valuable learning experience to better control the intrinsically different information content and contamination errors of the vector and the gradient data. Some tools have been developed to account for the statistical differences between these two types of datasets. In this paper, we did not try to optimize the weights or to regularize the inverse problem but this possibility was left open on both the regional (Process 6.2) and the global scales (Process 8). More theoretical work will be needed to dene the most relevant physical norms to be used for the lithospheric eld and fully benet from this possibility. We may also select the spectral band corresponding to the highest information content for each type of data to constrain the model with vector data at the low SH degrees and with gradient data at high SH degrees (see Kotsiaros and Olsen, 2012 for a discussion). This is an option we did not consider at present because the gradient estimated from the Swarm lower pair of satellites is not complete. Its implementation is not necessarily complex but it adds another level of subjectivity concerning the respective spectral bands that are better constrained by each data type.

A major issue we did not address is that of the impact on the East-West gradient of undesired magnetic signals due to electric currents that would cross the two orbits of the lower pair of satellites. At present, the R-SCHA technique is efcient for processing poloidal magnetic eld measurements but cannot handle toroidal elds without some significant amount of development. This is the main shortcoming of the Dedicated Lithospheric Field Inversion chain in its present form. But we may still exploit the information provided by the East-West gradient to at least estimate the amount of current crossing the Swarm A and B mesocenter (see Kotsiaros and Olsen, 2012, equation (3.8); or Shen et al., 2012, for instance) and to use this as a selection criterion. We could then use the gradient at times (and location) of negligible estimated currents. We note that, as it stands,

the chain could also readily be used with scalar gradient data only that are in general less affected by toroidal elds.

One key advantage of the Dedicated Lithospheric Field Inversion chain is that it is based on a regional modeling approach. This allows the inverse problem to be solved independently for each region. Depending on the available computer resources, the process divided in as many independent jobs is intrinsically fast. For the present simulation, a few hours were necessary to select and process more than ve million out of more than 250.106 available vector synthetic measurements of the Swarm A and B satellites. A fast algorithm is advantageous considering the chosen philosophy. Isolating the lithospheric eld relies on a careful selection of data and on compromises made in correcting the data more or less empirically. Therefore, fast processes open an avenue for computing parades of statistically equivalent models in order to estimate the variability of each SH degree part depending on the applied selections and corrections. An immediate application of this principle would be to play with the boundary separating mid and polar regions as it introduces one of the largest source of errors in this simulation (see Fig. 8).

The algorithm is also highly modular. Complementary processing blocks could be included as more and more experience is gained with the real Swarm satellite measurements. Preprocessing blocks based on quality indicators, ags, and star camera switches could be useful to account for instrumental failures or malfunctions. The algorithm is therefore not frozen and will benet from the scientic advances prompted by the Swarm mission along its entire duration. Finally, we note that the median altitude of Swarm A and B measurements was about 400 km for the simulations we discussed. An extension of the mission to lower altitudes would amplify the small-scale signals, which are currently strongly attenuated at 400 km altitude, and would possibly allow us to derive lithospheric eld models to even higher resolution.

Acknowledgments. The authors gratefully acknowledge support from the Centre National des Etudes Spatiales (CNES) within the context of the Travaux prparatoires et exploitation de la mission Swarm project and from the European Space Agency (ESA) through the ESTEC contract No. 4000102140/10/NL/JA Development of the Swarm Level2 Algorithms and Associated Level2 Processing Facility. We also acknowledge the comments of two anonymous referees that helped to improve the manuscript. This is IPGP contribution number 3410.

References

Alken, P., S. Maus, P. Vigneron, O. Sirol, and G. Hulot, Swarm SCARF equatorial electric eld inversion chain, Earth Planets Space, 65, this issue, 13091317, 2013.

Backus, G., R. L. Parker, and C. Constable, Foundations of Geomagnetism,

Cambridge Univ. Press, Cambridge, U.K., 1996.

Beggan, C. D., J. Saarimaki, K. A. Whaler, and F. J. Simons, Spectral and spatial decomposition of lithospheric magnetic eld models using spherical Slepian functions, Geophys. J. Int., doi:10.1093/gji/ggs122, 2013.

Chulliat, A., P. Vigneron, E. Thbault, O. Sirol, and G. Hulot, Swarm

SCARF Dedicated Ionospheric Field Inversion chain, Earth Planets Space, 65, this issue, 12711283, 2013.

Finlay, C. C., S. Maus, C. Beggan, M. Hamoudi, V. Lesur, F. J. Lowes,N. Olsen, and E. Thbault, Evaluation of candidate geomagnetic eld models for IGRF-11, Earth Planets Space, 62(10), 787804, 2010. Fox Maule, C., M. Purucker, N. Olsen, and K. Mosegaard, Heat ux

1270 E. THBAULT et al.: SWARM SCARF DEDICATED LITHOSPHERIC FIELD INVERSION CHAIN

anomalies in Antarctica revealed from satellite magnetic data, Science, doi:10.1126, 2005.

Friis-Christensen, E., H. Lhr, and G. Hulot, Swarm: A constellation to study the Earths magnetic eld, Earth Planets Space, 58(4), 351358, 2006.

Haines, G. V., Spherical cap harmonic analysis, J. Geophys. Res., 90,25832591, 1985.

Hemant, K. and S. Maus, Geological modeling of the new CHAMP magnetic anomaly maps using a Geographical Information System (GIS) technique, J. Geophys. Res. B, 110, B12103, doi:10.1029/2005JB003837, 2005.

Hulot, G., N. Olsen, and T. J. Sabaka, The present eld, in Volume 5 Geomagnetism, Treatise on Geophysics, vol. 5 (Elsevier, Amsterdam, 2007), pp. 3372, 2007.

Kotsiaros, S. and N. Olsen, The geomagnetic eld gradient tensor Properties and parametrization in terms of spherical harmonics, Int. J.

Geomath., 3, 297314, 2012.

Langel, R. A. and W. J. Hinze, The Magnetic Field of the Earths Litho-sphere: The Satellite Perspective, 429 pp., Cambridge Univ. Press, New York, 1998.

Langlais, B., M. E. Purucker, and M. Mandea, Crustal magnetic eld of Mars, J. Geophys. Res., 109(E2), doi:10.1029/2003JE002048, issn: 0148-0227, 2004.

Lesur, V., Introducing localized constraints in global geomagnetic eld modelling, Earth Planets Space, 58, 477483, 2006.

Lesur, V. and S. Maus, A global lithospheric magnetic eld model with reduced noise level in the Polar Regions, Geophys. Res. Lett., 33, L13304, doi:10.1029/2006GL025826, 2006.

Lesur, V., I. Wardinski, M. Hamoudi, and M. Rother, The second generation of the GFZ Reference Internal Magnetic Model: GRIMM-2, Earth Planets Space, 62(10), 765773, 2010.

Lesur, V., M. Rother, F. Vervelidou, M. Hamoudi, and E. Thbault, Postprocessing scheme for modeling the lithospheric magnetic eld, Solid Earth, 4, 105118, doi:10.5194/se-4-105-2013, 2013.

Lowes, F. J., Mean square values on the sphere of spherical harmonic vector elds, J. Geophys. Res., 71, 2179, 1966.

Maier, T. and C. Mayer, Multiscale downward continuation of CHAMP FGM-data for crustal eld modelling, in First CHAMP Mission Results for Gravity, Magnetic and Atmospheric Studies, edited by C. Reigber,H. Lhr, and P. Schwintzer, pp. 288295, Springer Verlag, 2003.

Maus, S., An ellipsoidal harmonic representation of Earths lithospheric magnetic eld to degree and order 720, Geochem. Geophys. Geosyst., 11, Q06015, doi:10.1029/2010GC003026, 2010.

Maus, S., H. Lhr, and M. E. Purucker, Simulation of the high-degree lithospheric eld recovery for the Swarm constellation of satellites, Earth Planets Space, 58(4), 397407, 2006.

Maus, S., F. Yin, H. Lhr, C. Manoj, M. Rother, J. Rauberg, I. Michaelis,C. Stolle, and R. D. Mller, Resolution of direction of oceanic magnetic lineations by the sixth-generation lithospheric magnetic eld model from CHAMP satellite magnetic measurements, Geochem. Geophys. Geosyst., 9, Q07021, doi:10.1029/2008GC001949, 2008.

Olsen, N., A model of the geomagnetic eld and its secular variation for epoch 2000, Geophys. J. Int., 149(2), 454462, 2002.

Olsen, N., M. Mandea, T. J. Sabaka, and L. Toffner-Clausen, The CHAOS-3 Geomagnetic Field Model and Candidates for the 11th Generation of IGRF, Earth Planets Space, 62, 719727, 2010a.

Olsen, N., G. Hulot, and T. J. Sabaka, Sources of the geomagnetic eld and the modern data that enable their investigation, in Handbook of Geomathematics, edited by W. Freeden, M. Z. Nashed, and T. Sonar, chapter 5, pp. 106124, doi:10.1007/978-3-642-01546-5-5, 2010b. Olsen, N., K.-H. Glassmeier, and X. Jia, Separation of the Magnetic

Field into External and Internal Parts, Space Sci. Rev., 152, 135157, doi:10.1007/s11214-009-9563-0, 2010c.

Olsen, N., E. Friis-Christensen, R. Floberghagen, P. Alken, C. D Beggan,A. Chulliat, E. Doornbos, J. T. da Encarnaao, B. Hamilton, G. Hulot, J. van den IJssel, A. Kuvshinov, V. Lesur, H. Lhr, S. Macmillan, S. Maus,

M. Noja, P. E. H. Olsen, J. Park, G. Plank, C. Pthe, J. Rauberg, P. Ritter, M. Rother, T. J. Sabaka, R. Schachtschneider, O. Sirol, C. Stolle,E. Thbault, A. W. P. Thomson, L. Tffner-Clausen, J. Velmsk, P. Vigneron, and P. N. Visser, The Swarm Satellite Constellation Application and Research Facility (SCARF) and Swarm data products, Earth Planets Space, 65, this issue, 11891200, 2013.

Purucker, M., B. Langlais, N. Olsen, G. Hulot, and M. Mandea, The southern edge of cratonic North America: Evidence from new satellite magnetometer observations, Geophys. Res. Lett., 29(15), article number 013645, doi:10.1029/2001GL013645, 2002.

Reigber, C., H. Lhr, and T. Schwintzer, Champ Mission status, Adv. Space

Res., 30(2), 129134, doi:10.1016/S02731177(02)00276-4, 2002. Rother, M., V. Lesur, and R. Schachtschneider, An algorithm for deriving core magnetic eld models from the Swarm data set, Earth Planets Space, 65, this issue, 12231231, 2013.

Sabaka, T. J. and N. Olsen, Enhancing comprehensive inversions using the

Swarm constellation, Earth Planets Space, 58, 371395, 2006. Sabaka, T. J., L. Tffner-Clausen, and N. Olsen, Use of the Comprehensive Inversion method for Swarm satellite data analysis, Earth Planets Space, 65, this issue, 12011222, 2013.

Saff, E. B. and A. B. J. Kuijlaars, Distributing many points on a sphere,

The Mathematical Intelligencer, 19(1), 511, 1997.

Shen, C., Z. J. Rong, and M. Dunlop, Determining the full magnetic eld gradient from two spacecraft measurements under special constraints, J. Geophys. Res., 117, A10217, doi:10.1029/2012JA018063, 2012. Schaeffer, N., Efcient Spherical Harmonic Transforms aimed at pseudo-spectral numerical simulations, Geochem. Geophys. Geosyst., 14(3), 751758, 2013.

Schott, J. J and E. Thbault, Modelling the earths magnetic eld from global to regional scales, in Geomagnetic Observations and Models, IAGA Special Sopron Book Series, Volume 5, 229264, doi:10.1007/978-90-481-9858-0-9, 2011.

Stockman, R., C. C. Finlay, and A. Jackson, Imaging Earths crustal magnetic eld with satellite data: A regularized spherical triangle tessellation approach, Geophys. J. Int., 179, 929944, doi:10.1111/j.1365-246X.2009.04345.x, 2009.

Swarm Level 2 Processing System Consortium, Product specication for

L2 Products and Auxiliary Products, Doc.no: SW-DS-DTU-GS-0001, 2013.

Thbault, E., Global lithospheric magnetic eld modeling by successive regional analysis, Earth Planets Space, 58, 485495, 2006.

Thbault, E., J. J. Schott, M. Mandea, and J. P. Hoffbeck, A new proposal for Spherical Cap Harmonic Analysis, Geophys. J. Int., 159, 83105, doi:10.1111/j.1365-246X.2004.02361.x, 2004.

Thbault, E., J. J. Schott, and M. Mandea, Revised spherical cap harmonic analysis (R-SCHA): Validation and properties, J. Geophys. Res., 111, B01102, doi:10.1029/2005JB003836, 2006.

Thbault, E., K. Hemant, G. Hulot, and N. Olsen, On the geographical distribution of induced time-varying crustal magnetic elds, Geophys. Res. Lett., 36, L01307, doi:10.1029/2008GL036416, 2009.

Thbault, E., M. Purucker, K. Whaler, B. Langlais, and T. J. Sabaka, The

Magnetic eld of the Earths lithosphere, Space Sci. Rev., 155(14), 95 127, doi:10.1007/s11214-010-9667-6, 2010.

Thbault, E., F. Vervelidou, V. Lesur, and M. Hamoudi, The along-track satellite analysis in planetary magnetism, Geophys. J. Int., 188(3), 891 907, doi:10.1111/j.1365-246X.2011.05281.x, 2012.

Thomson, A. W. P. and V. Lesur, An improved geomagnetic data selection algorithm for global geomagnetic eld modelling, Geophys. J. Int., 169(3), 951963, doi:10.1111/j.1365-246X.2007.03354.x, 2007. Wang, Z., Understanding models of the geomagnetic eld by Fourier analysis, J. Geomag. Geoelectr., 39, 333347, 1987.

E. Thbault (e-mail: [email protected]), P. Vigneron, S. Maus, A. Chulliat, O. Sirol, and G. Hulot