# Post-processing scheme for modeling the lithospheric magnetic field

V. Lesur, M. Rother, F. Vervelidou, M. Hamoudi, E. Thébault

► To cite this version:

V. Lesur, M. Rother, F. Vervelidou, M. Hamoudi, E. Thébault. Post-processing scheme for modeling the lithospheric magnetic field. Solid Earth Discussions, 2012, 4, pp.1345-1378. 10.5194/sed-4-1345-2012. insu-03583014

HAL Id: insu-03583014

https://insu.hal.science/insu-03583014

Submitted on 22 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

This discussion paper is/has been under review for the journal Solid Earth (SE).
Please refer to the corresponding final paper in SE if available.

# Post-processing scheme for modeling the lithospheric magnetic field

V. Lesur<fnref n="1" />, M. Rother<fnref n="1" />, F. Vervelidou<fnref n="2" />, M. Hamoudi<fnref n="3" />, and E. Thébault<fnref n="2" />

Received: 9 October 2012 – Accepted: 12 October 2012 – Published: 30 October 2012

Correspondence to: V. Lesur (lesur@gfz-potsdam.de)

Published by Copernicus Publications on behalf of the European Geosciences Union.

Title Page

Abstract

Introduction

Conclusions

References

Tables

Figures

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

Back

Close

Full Screen / Esc

Printer-friendly Version

Interactive Discussion

<figure><img src="image_7.png" /></figure>

<footnote n="1">¹Helmholtz Centre Potsdam, GFZ German Research Centre for Geosciences, Telegrafenberg, 14473, Germany</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Sorbonne, Paris Cité, Univ Paris Diderot, UMR7154, CNRS, F-75005 Paris, France</footnote>

<footnote n="3">³Université des Sciences et de la Technologie, USTHB, 16123 Bab Ezzouar, Alger, Algeria</footnote>
## Abstract

We investigated how the noise in satellite magnetic data affects magnetic lithospheric field models derived from these data in the special case where this noise is correlated along satellite orbit tracks. For this we describe the satellite data noise as a perturbation magnetic field scaled independently for each orbit, where the scaling factor has a specified variance. Under this assumption, we have been able to derive a model for errors in lithospheric models generated by the correlated satellite data noise. Unless the perturbation field is known, estimating the noise in the lithospheric field model is a non-linear inverse problem. We therefore proposed an iterative post-processing technique to estimate both the lithospheric field model and its associated noise model. The technique has been successfully applied to derive a lithospheric field model from CHAMP satellite data up to spherical harmonic 120. The model is in agreement with other existing models. The technique can be in principal extended for all kind of potential field data with "along track" correlated errors.

## 1 Introduction

All geophysical data are contaminated by signals that cannot be easily described by models. These poorly parameterized contributions are often treated as errors and they most of the time exceed the pure instrumental noise. These kind of errors are particularly difficult to deal with because they are often correlated in space and/or time. Further they may not follow a gaussian distribution. Yet, handling properly the data errors is at the heart of the data interpretation process and it usually requires their full statistical description – i.e. for a set of discrete measurements, the knowledge of the full covariance matrix of the data errors.

Geopotential data – i.e. gravity and magnetic measurements – are not an exception. For these types of data, the inverse problem that consists in finding the sources of the signals is particularly ill posed, and the proper statistical description of the data

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

5

10

15

20

25
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

errors is necessary. Failing to do so may lead to false conclusions about the signal sources. From a practical point of view, scientists have been relatively successful in estimating a priori the noise in gravity or magnetic data sets, however correlations between errors have been most of the time ignored. This is partly because, when known, the full covariance matrix for the data errors is generally so large that it cannot be handled easily, even on modern computers (but see Langel et al., 1989; Holme and Bloxham, 1996; Rygaard-Hjalsted et al., 1997; Holme, 2000 as examples where correlated errors are accounted for in geomagnetism).

The effects of these correlation errors are obvious in airborne, marine and satellite data. Typically, in all these type of surveys, the data are collected along linear paths and, after processing, the correlation errors become apparent as offsets between adjacent tracks. They then appear in maps and models as spurious anomalies, elongated in the direction of the tracks. An example of such an effect is shown in this manuscript for magnetic models derived from satellite data. The traditional way of dealing with this noise has been to perform a “leveling” of the data. In airborne geophysics, the approach mainly consists in deriving for each track a polynomial expression that is subtracted from the data such as to minimize data differences at the cross-over points (Hamoudi et al., 2010). The method has been also adapted to satellite magnetic data. In that case a large-scale field of external origin is fitted to a data set made of only few tracks. This allows to successfully derived magnetic field models of the lithosphere to relatively high degree. A well known example is the MF series of models – e.g. Maus et al. (2008). However the method, as applied to satellite data, has its drawbacks. The effects of its application have been carefully studied in Thebault et al. (2012) and it appears that, depending on the way the method is applied, it can lead to significant distortions of the final model. However, the weakest point of this so called “along track filtering” approach is the impossibility to estimate how much the processing applied distorts the model. For this aspect, post-processing techniques are preferable.

So far post-processing techniques have been developed and applied only to models derived from satellite gravity data – e.g. Kusche (2007). To the authors knowledge such
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

techniques have never been applied to magnetic models although we should note the attempt to estimate the model covariance matrix in Lowes and Olsen (2004). In this manuscript we present and apply such a post-processing scheme for a model of the magnetic lithospheric field derived from ten years of CHAMP satellite data (Reigber et al., 2005). Although we are presenting this work from its application side, it has deeper roots: We investigated how typical noise correlated patterns leak, through a least squares fitting process, inside a magnetic model of the lithospheric field. Therefore this leads to a model of the noise inside the lithospheric model. Once such a noise model is available, numerous post-processing schemes are possible; we just applied one specific approach to show that the noise model we obtain is relevant. The final resulting model of the lithospheric field is nonetheless of high quality and compares well with other recently released models e.g. MF7 (not published, but see MF6 in Maus et al., 2008), CHAOS-4 (Olsen et al., 2010).

The manuscript is organized as follows. In the next section we set the hypothesis and approximations, derive the general expression for the noise model and give examples of possible noise depending on the characteristic of the perturbation magnetic field in the data. In the third section we describe in details the two steps process towards the final lithospheric field model. The resulting model is then discussed in the last section.

## 2 The lithospheric noise model

In this section we present a noise model for a lithospheric model estimated from a set of radial magnetic data. We choose to present this case only in the main part of this manuscript as the equations are relatively simple to derive. The description for the usual case where the lithospheric model is obtained from the three components of a magnetic data set is given in Appendix A.
## 2.1 Theory

We consider a magnetic data set made of radial component readings along a single CHAMP satellite half-orbit during night times. For simplicity we will assume that a track follows a meridian – i.e. it corresponds to a single longitude value. Models of the core field, of the lithospheric field and of the large-scale external field are then subtracted from these data. Because the description of the external field is not very accurate, the residuals obtained along that track contain relatively long wavelengths. We assume that these residuals along this single orbit are well approximated by the radial component of an external magnetic field model that does not present time dependencies. It is hereafter named as the perturbation field and writes:

$$B_{rp}(\theta, \phi, r) = - \sum_{n,k}^{N} \left(\frac{r}{a}\right)^{n-1} n e_n^k Y_n^k(\theta, \phi), \quad (1)$$

where $e_n^k$ are the Gauss coefficients of degree $n$ and order $k$, $a = 6371.2$ km is the Earth's reference radius, $Y_n^k(\theta, \phi)$ are the Schmidt semi-normalized Spherical Harmonics (SHs). We use along this manuscript the convention that negative orders, $k < 0$, are associated with $\sin(|k|\phi)$ terms whereas null or positive orders, $k \ge 0$, are associated with $\cos(k\phi)$ terms. We will also often use the index $l$ for the degree and $m$ for the order.

We consider also a model of the radial component of a magnetic field of internal origin with no temporal dependencies. This model becomes below the lithospheric noise model we want to derive:

$$\tilde{B}_{ri}(\theta, \phi, r) = \sum_{l,m}^{L} \left(\frac{a}{r}\right)^{l+2} (l+1) \tilde{g}_{l}^{m} Y_{l}^{m}(\theta, \phi). \quad (2)$$

It is not possible to separate external field contributions from internal field contributions for data collected along a single meridian – i.e. a single half-orbit, hence we can fit

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

by least-squares the residuals defined in Eq. (1) with the lithospheric model given in Eq. (2). We simply need to minimize the functional:

$$ \Phi_j = \sum_i w_i |\tilde{B}_{ri}(\theta_i, \phi_j, r) - B_{rp}(\theta_i, \phi_j, r)|^2, \quad (3) $$

where $\theta_i$ are sampling points along the half-orbit, $\phi_j$ is the longitude of the meridian that we labeled with the subscript $j$ and $w_i$ are weights that are defined below.

Over 10 yr, the CHAMP satellite has collected data along a large number $M$ of half-orbits. We assume now that for each orbit the perturbation field model defined by Eq. (1) is scaled by a number $\eta_j$ and that all orbits are at the same radius $r$. This latter point is clearly a strong approximation but there is no obvious way to avoid it. Again, these external field contributions can be interpreted as a field of internal origin. To estimate this field, the functional we have to minimize is then:

$$ \Phi = \sum_{i,j} w_i |\tilde{B}_{ri}(\theta_i, \phi_j, r) - \eta_j \cdot B_{rp}(\theta_i, \phi_j, r)|^2. \quad (4) $$

Minimizing $\Phi$ for the Gauss coefficients $\tilde{g}_i^m$ leads to a system of equations:

$$ A^t A \tilde{g} = A^t b \quad (5) $$

where $\tilde{g} = [\tilde{g}_i^m]_{\{l,m\}}$. The matrix product $A^t A$ is derived from Eqs. (4) and (2) and the elements of this product associated with the degrees and orders $l, l', m, m'$ writes:

$$ \{A^t A\}_{l,m,l',m'} = M \left(\frac{a}{r}\right)^{l+l'+4} (l+1)(l'+1) \langle P_l^{|m|}, P_{l'}^{|m'|} \rangle \Pi_{mm'} \quad (6) $$

where the product $\langle P_l^m, P_{l'}^{m'} \rangle$ is defined by $\langle P_l^m, P_{l'}^{m'} \rangle = \sum_i w_i P_l^m(\cos \theta_i) P_{l'}^{m'}(\cos \theta_i)$. The variable $\Pi_{mm'}$ has been introduced to cover three cases:

$$ \Pi_{mm'} = \begin{cases} \frac{1}{M} \sum_{i=1}^{M} \cos m\phi_i \sin|m'|\phi_i & \text{if } mm' < 0 \\ \frac{1}{M} \sum_{i=1}^{M} \cos m\phi_i \cos m'\phi_i & \text{if } m \ge 0, m' \ge 0 \\ \frac{1}{M} \sum_{i=1}^{M} \sin|m|\phi_i \sin|m'|\phi_i & \text{if } m < 0, m' < 0, \end{cases} \quad (7) $$
and is symmetric relative to its subscripts – i.e. $\Pi_{mm'} = \Pi_{m'm}$. The elements of the right hand side vector of Eq. (5) are:

$$ \{\mathbf{A}^t \mathbf{b}\}_{l,m} = -M \sum_{n,k}^{N} \left(\frac{a}{r}\right)^{l-n+3} n(l+1) \langle P_n^{|k|}, P_l^{|m|} \rangle \epsilon_n^k \chi_m^k. \quad (8) $$

Depending on the sign of the orders $m$ and $k$, $\chi_m^k$ takes the following values:

$$ \chi_m^k = \begin{cases} \frac{1}{M} \sum_{i=1}^{M} \cos m\phi_i \cos k\phi_i \quad \eta_i & \text{if } m, k \ge 0 \\ \frac{1}{M} \sum_{i=1}^{M} \cos m\phi_i \sin k\phi_i \quad \eta_i & \text{if } mk < 0 \\ \frac{1}{M} \sum_{i=1}^{M} \sin m\phi_i \sin k\phi_i \quad \eta_i & \text{if } m, k < 0, \end{cases} \quad (9) $$

and, as for $\Pi_{mm'}$, it is symmetric relative to its subscripts: $\chi_m^k = \chi_k^m$.

For very large number $M$ of orbits uniformly distributed along longitudes, the quantity $\Pi_{mm'}$ tends to a $\delta$-function – i.e. $\Pi_{mm'} \approx (\frac{1}{2} + \frac{1}{2}\delta_{m0})\delta_{mm'}$. Further, by setting the weights $w_i$ to $w_i = \sin\theta_i$ and assuming that the sampling points are evenly spaced over the full meridian, we have $\langle P_l^{|m|}, P_{l'}^{|m|} \rangle = \frac{4-2\delta_{m0}}{2l+1}\delta_{ll'}$ and therefore, from Eq. (6), the product matrix $\mathbf{A}^t\mathbf{A}$ is diagonal. Regarding Eqs. (8) and (9), if the $\eta_i$ form a set of uncorrelated random variables, the $\chi_m^k$ are also random variables with zero mean.

The Gauss coefficients for the lithospheric noise model in Eq. (2) are then:

$$ \tilde{g}_l^m = - \sum_{n,k}^{N} \left(\frac{r}{a}\right)^{l+n+1} n \frac{2l+1}{2l+2} \langle P_n^{|k|}, P_{l'}^{|m|} \rangle \epsilon_n^k \chi_m^k, \quad (10) $$

and they correspond to the noise in a lithospheric field model that would be generated by un-modeled external fields in the radial component of magnetic data. Similarly, it is straightforward to find the noise in a lithospheric field model (i.e. static internal field

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>
model) generated by a perturbation field of internal origin. This case is relevant for signals generated in the lower E-region ionosphere (e.g. at 110 km altitude) when data are acquired at satellite altitudes. Other possible sources for this type of noise are the un-modeled induced fields generated in the conductive layers of the Earth by rapid variations of the external fields. It gives:

$$ \tilde{g}_l^m = \sum_{n,k}^{N} \left(\frac{r}{a}\right)^{l-n} (n+1) \frac{2l+1}{2l+2} \langle P_n^{|k|}, P_l^{|m|} \rangle i_n^k \chi_m^k \quad (11) $$

where $i_n^k$ are the Gauss coefficients for the ionospheric and/or induced field models.

In order to understand the behaviour of the lithospheric noise model, it is important to have an estimate of the probability density function of the random variable $\chi_m^k$. Assuming the random variable $\eta$ is normally distributed with variance $\nu^\eta$ then $\chi_m^k$ appears to be also normally distributed. The set of $\chi_m^k$ are uncorrelated with the exception that $\chi_m^k = \chi_k^m$. Further the $\chi_m^k$ have a variance $\nu^\chi$ that depends on $\nu^\eta$, the number of half-orbits $M$, the orders $k$ and $m$. Possible values of the variance $\nu^\chi$ are given in Table 1. These variances have been derived from numerical experiments involving 20 000 independent realizations of the random variables $\chi_m^k$ calculated from the same number of uniformly distributed orbits. Figure 1 presents the histograms for few values of $m$ and $k$.

In the following we consider only the noise model given by Eq. (11). The general behaviour of the noise characterized by Eqs. (11) and (10) is basically the same. In particular they have the same dependence relative to the degree $l$. These two noise models are only relevant for the cases where the radial components of vector data are used. The way the noise propagates in a lithospheric model is different if the three components of the vector data are fitted. The corresponding equations for that case are relatively complex and given in Appendix A.

The noise model defined in Eq. (2) has $L(L+2)$ parameters – i.e. $L(L+2)$ Gauss coefficients. This number reduces to $N(N+2)$ Gauss coefficients $i_n^k$ with

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>
<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

$(2N + 1)(2L + 1) - 2N^2$ random variables $\chi_m^k$ through Eq. (11). For small values of $N$ – e.g. $N = 10$, there is a very significant reduction of number of parameters, but the Eq. (11) is non-linear.

## 2.2 Examples

In order to understand the main characteristics of the noise model defined by Eqs. (2) and (11), we present in this section the results of forward modelling calculations for a given choice of Gauss coefficients $i_n^k$ and one realization of the set of random variables $\chi_m^k$. The products $\langle P_n^{|k|}, P_l^{|m|} \rangle$ are calculated numerically. These products are relatively difficult to estimate accurately as the $P_l^m(x)$ functions are oscillatory. However, an adaptive Gaussian quadrature was ultimately chosen as it gave the best results.

### 2.2.1 Dipole perturbation field

For this first example we use a simple model for the perturbation field of internal origin made of a single spherical harmonic $n = 1, k = 1$. Specifically, we set $i_1^1 = 1$ and $i_n^k = 0$ for $\{n, k\} \neq \{1, 1\}$. This type of noise in satellite data could result from a poor modelling of the field induced by a large-scale external field in the conductive layers of the Earth. In that case Eq. (11) reduces to:

$$ \tilde{g}_l^m = \left(\frac{r}{a}\right)^{l-1} 2 \frac{2l + 1}{2l + 2} \langle P_1^{|1|}, P_l^{|m|} \rangle i_1^1 \chi_m^1, \quad (12) $$

and the noise in the radial component of the field of internal origin is:

$$ \tilde{B}_{ri}(\theta, \phi, r') = i_1^1 2 \left(\frac{a}{r'}\right)^3 \sum_{l,m}^{L} \left(\frac{r}{r'}\right)^{l-1} \frac{(l+1)(2l+1)}{(2l+2)} \langle P_1^{|1|}, P_l^{|m|} \rangle \chi_m^1 Y_l^m(\theta, \phi), \quad (13) $$

where $r'$ is the modelling radius that is set to $r' = a = 6371.2$ km in this example. As the observation radius $r$ is expected to be larger than the modelling radius, the short
wavelengths dominate the model due to the ratio $r/r_0$ raised to the power $l-1$ in the right hand side of Eq. (13).

In Fig. 2, the model defined by Eq. (13) is mapped for the model coefficient $l_1^1 = 1$ nT, an observation radius at 300 km altitude ($r = 6671.2$ km) and the random variables $\chi_m^1$ with variances defined in Table 1 using $\nu^\eta = M$. The maximum SH degree involved is $L=120$. We observe that the noise model is symmetric relative to the equator, vanishes at the poles, and is made of East-West oscillating anomalies typical of the noise in lithospheric field model derived from satellite data. We note that these characteristics are independent from the sign of the SH order $k$ as only the random variable $\chi_m^k$ depends of this sign in Eq. (13). The obtained symmetry of the model is due to the product $\langle P_1^{[1]}, P_1^{[m]} \rangle$ that vanishes if the Legendre function $P_1^{[m]}$ is anti-symmetric – i.e. $l-|m|$ is odd. An anti-symmetric model, vanishing at the equator but not at the poles, would have been obtained if $l_1^0 = 1$ nT would have been chosen in place of $l_1^1 = 1$ nT. These symmetry/anti-symmetry characteristics are specific to models derived from the radial component alone. It can be seen in Appendix A that these characteristics are lost when a noise model is obtained from the three vector components.

The power spectrum of the model calculated at $r' = 6671.2$ km is also plotted in Fig. 2. It presents some variability due to the use of a single SH in Eq. (12). Nonetheless, the behaviour is generally along a $(r/r_0)^{2l}$ trend as it would be expected for a white noise at satellite altitude. Although the small wavelengths overshadow the larger wavelengths, the latter are also present in the noise model. It is clear that any magnetic field model derived from satellite data is contaminated by such a noise at all wavelengths unless pertinent processing steps are applied.

### 2.2.2 Auroral electrojet and field aligned currents

Another expected source of noise in satellite data is associated with the auroral electrojet and/or associated field-aligned currents. We do not aim at a precise description of the disturbance field but just consider the radial component of a perturbation field of

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>
internal origin, mapped in Fig. 3 (left), and defined by:

$$B_{rp}(\theta, \phi, r) = \sum_{n,k}^{N} \left(\frac{a}{r}\right)^{n+2} (n+1) l_n^k Y_n^k(\theta, \phi). \quad (14)$$

We recall that in our approach this field is scaled by a random variable with zero mean for each orbit. Therefore it is more the geometry of the field that is important here than its true value. We see that the perturbation field model is centred on the geomagnetic North pole and takes relatively large values up to 60° colatitudes. The lithospheric noise model we obtain is mapped in Fig. 3 (right). This model is also fairly well localized in latitudes as it basically vanishes in the southern hemisphere. However, it seems that the noise is propagating over all longitudes. The power spectrum of the model has essentially the same characteristic than in the previous example.

The results of this example have to be analysed with some caution since real satellite orbits deviate from the exact polar direction at high latitudes. Nonetheless, we take out from these results that there is no need to describe precisely the longitudinal dependence of the rapidly varying field to obtain realistic noise model. Therefore, in Eq. (11), the range of SH order k can be restricted to small values – e.g. $k_{max} = 2$, even if the maximum SH degree in the model remains large – e.g. N = 30. This will reduce even further the number of parameters needed to describe the noise model.

### 3 Application to magnetic models of the lithosphere

The process we applied to generate an accurate field model is in two steps. First we estimate a rough model from satellite data using a straightforward least-squares approach. Second, in the post-processing step, we co-estimate a new lithospheric field model and a model of the noise where the output model of the first step is used as data. The final results depend on the processes applied during the two steps and therefore both are described in independent subsections below.

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>
<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

In order to avoid confusion between the different models we use the following notations:

1. the noise model is denoted, in the same way as in the previous section, using a “~” – e.g. $\tilde{B}_i$ for the magnetic field vector,

5

2. the noisy lithospheric model, output of the first step, is denoted using a “^” – e.g. $\hat{B}_i$,

3. the lithospheric field model output of the post processing step does not have any distinctive sign – e.g. $B_i$.

## 3.1 Data set, data selection, model parameterization and model estimation

10

Three component vector magnetic readings acquired during the ten years of the German CHAMP satellite mission are used. The data are selected for night-times and magnetically quiet days, in the same way as data are selected for the GRIMM series of core field model (Lesur et al., 2008, 2010). However, in the present case the three components of the vector data are used and data in single star camera mode are rejected. In the usual GRIMM selection scheme, at mid and low latitudes, only the X and Y SM components are selected. A core field model and a model of the large-scale external field with its internally induced counterpart are subtracted from these data, leaving mainly the contributions from the lithosphere and the noise. The core field model and external field models used are resulting from the derivation of GRIMM-3 (Lesur et al., 2011), but this is not seen as an important point in the processing: another core field model would have been possible – e.g. CHAOS-4 (Olsen et al., 2010).

15

20

Next a first lithospheric field model up to SH degree 60 is derived, but our aim here is to reject outliers. The data corresponding to residuals larger than 3 times the standard deviation are rejected. The value of the threshold, for each data type, is given in table 2. This selection process is known to potentially affect strongly the final lithospheric field model. At mid and low latitudes only few data are rejected, and those rejected

25
<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>

data do not present clusters: No major difficulties are therefore expected there. At high latitudes however, a large amount of data are rejected and it is not possible to assess at this point if magnetic anomalies are erased or minimized there. We checked however, that outside the polar gaps due to the satellite orbits, the final data density at satellite altitude is everywhere large enough to allow for a lithospheric field model to be estimated up to SH degree 120.

In order to avoid spurious oscillations of the lithospheric model, further vertical down component data values were added over the polar gaps at an altitude of 6371.2 km. The data values were arbitrarily set to zero and the associated weights for the inversion process adjusted such that the lithospheric field model looks acceptable. We have tested other possible approaches, but the one we used gave the best results. One could alternatively use vertical down component values derived from aeromagnetic maps.

The data set resulting from this selection process still consists in some 5014325 data values. A model of the lithosphere magnetic field, defined by Eq. (15) below, was fitted through a simple least-squares process to the data. The data weights used depend only on the data density.

$$ \hat{B}_i(\theta, \phi, r) = -\nabla \left[ a \sum_{l,m} \left(\frac{a}{r}\right)^{l+1} \hat{g}_l^m Y_l^m(\theta, \phi) \right]. \quad (15) $$

The power spectrum at the Earth surface of the resulting lithospheric field model is presented in Fig. 4, left, together with the power spectrum of the CHAOS-4 model. Both models present very similar spectra up to SH degree 60 or 65. Our model presents slightly less power around degree 70, possibly due to the selection technique used. Above SH degree 85, CHAOS-4 is strongly minimized, whereas our model presents a spectrum rising to high values, evidence of the predominance of noise in the model at these SH degrees. The final misfits to the data are given in Table 2.

The vertical down component of the model – i.e. the Z component, is mapped at 300 km above the Earth surface in Fig. 4, right. At this altitude the long wavelength
lithospheric signal dominates but the noise is clearly visible, mainly over oceans, as elongated anomalies in the north-south direction – e.g. to the south of Australia. We point out that there are strong correlations between the estimated Gauss coefficients of the model and therefore the model cannot be truncated at an arbitrary degree without introducing artefacts.

## 3.2 Model post-processing

The post-processing part consists in fitting a model of the magnetic field generated in the lithosphere together with the model of noise, to a 300 km altitude map of the vertical down component of the field model $\hat{B}_i(\theta, \phi, r)$ (see Fig. 4). The noise model $\tilde{B}_i$ we used is derived in Appendix A and is parameterized by the variable $\chi_m^k$ and the Gauss coefficients of the perturbation model $i_n^k$. This inverse problem that consists in fitting the noise model and the lithospheric field model to $\hat{B}_{Zi}$ values presents some difficulties that are described first, results are given in a second subsection.

### 3.2.1 Inverse problem

We map the vertical down component of the magnetic field model $\hat{B}_i(\theta, \phi, r)$ at 29161 positions on a Gauss-Legendre grid at $r = 300$ km altitude. These data values are related to the Gauss coefficients $g_i^m$ of the field model $B_i$ by the relation:

$$ \hat{B}_{Zi}(\theta_i, \phi_i, r) = - \sum_{l,m}^{L=120} (l+1) \left(\frac{a}{r}\right)^{l+2} g_l^m Y_l^m(\theta_i, \phi_i) + \tilde{B}_{Zi}(\theta_i, \phi_i, r) + \epsilon_i, \quad (16) $$

where $\tilde{B}_{Zi}(\theta_i, \phi_i, r)$ is the vertical down component of the noise model derived in Appendix A, and $\epsilon_i$ is an unknown noise. As the maximum SH degrees in $\hat{B}_i$ and $B_i$ are the same, it is clear that the $g_i^m$ can be estimated such as $B_i$ fits exactly the values of $\hat{B}_{Zi}(\theta_i, \phi_i, r)$, with the noise model and the $\epsilon_i$ not contributing to the problem. These latter contributions become necessary only when a priori smoothness requirements

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>
are introduced on $B_i$. Hence the inverse problem consists in minimizing the functional $\Phi$ defined by:

$$ \Phi = \sum_i \{ \hat{B}_{Zi}(\theta_i, \phi_i, r) - B_{Zi}(\theta_i, \phi_i, r) - \tilde{B}_{Zi}(\theta_i, \phi_i, r) \}^2 + \lambda \sum_{l,m} \frac{l(l+1)^3}{2l+1} (g_l^m)^2 \quad (17) $$

The first term insures the fit to the data $\hat{B}_{Zi}(\theta_i, \phi_i, r)$ whereas the second minimizes the integral of the squared horizontal gradient of the radial component of $B_i$ over a sphere of radius $a=6371.2$ km. The parameter $\lambda$ controls the smoothness constraint applied on $B_i$.

As stated above, the noise model $\tilde{B}_i$ (Eqs. A3 and A12) is parameterized by the variable $\chi_m^k$ and the Gauss coefficients of the perturbation model $i_n^k$. A possibility is to set the perturbation model coefficients $i_n^k$, such that the model corresponds to a dipole field, and to try to estimate the $\chi_m^k$. The inverse problem is then linear. However, for such a choice the derived lithospheric field model $B_i$ appears to be still contaminated by noise, probably because the perturbation model has to be more complex than a simple dipole. Therefore, there is no other option than to co-estimate the $\chi_m^k$ and $i_n^k$ values. As these quantities enter as products in Eq. (A12), the inverse problem is non-linear and must be solved iteratively. We want to point out that finding the $\chi_m^k$ and $i_n^k$ values in Eq. (A12) or in Eq. (11) are two different problems with their own specific null-space and difficulties. In particular, $i_n^k$ and $i_n^{-k}$ values cannot be estimated independently if Eq. (11) is used. With Eq. (A12) this estimation becomes possible solely because of the way the Y component data affect the noise model. However in both cases the maximum value for $n$ can be relatively large, whereas the maximum value of $k$ has to be small. We used in this work a maximum value of $n: N = 20$ and a maximum value for $k: K = 1$. As noted in Sect. 2.2.2 most of the complexity in longitude of the noise model is carried by the $\chi_m^k$; there is no need for a large longitudinal complexity of the perturbation model. With such settings, the number of unknown describing the noise model in Eq. (A12) is reduced to $N(2K+1)+K-K^2$ for the $i_n^k$ (i.e. 60 values) and $(2K+1)(2L+1)-2K^2$.

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>
<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /></figure>

for the $\chi_m^k$ (i.e. 721 values for $L = 120$). These numbers have to be compared with the number of unknown in the lithospheric field model $L(L+2) = 14640$.

The iterative inversion process we followed to reach the solution presented in the next subsection is described in three steps:

step-1 Find the $g_i^m$ in Eq. (16) minimizing $\Phi$ (Eq. (17)) imposing $\chi_m^k = 0$ for all possible $m$ and $k$ values.

step-2 Keeping the $g_i^m$ unchanged, and starting with $i_n^k = 1$ for all possible $n$ and $k$ values, find iteratively the $i_n^k$ and $\chi_m^k$ that minimize $\Phi$ in Eq. (17).

step-3 Iteratively find the $g_i^m$, $i_n^k$ and $\chi_m^k$ that minimize $\Phi$ in Eq. (17), starting from the output of the step 2.

10

## 3.2.2 Results

The results were obtained by iteratively minimizing the functional defined in Eq. (17) following the process described above, with the parameter $\lambda$ set to $\lambda = 4.0 \times 10^{-5}$ such that the resulting field model has an acceptable power spectrum. The level of noise is larger at high latitudes in the $\hat{B}_{Zi}(\theta_i, \phi_i, r)$, we therefore weight the data by $1/6$ for magnetic latitudes higher than 50°.

15

The output of the step-1 of the iterative process is a smoothed model obtained without co-estimation of the noise model. The map of this model vertical down component at radius 6371.2 km is shown in Fig. 5. The perturbation due to the along track noise in the satellite data are strong, particularly over Antarctica, and in the Indian, Atlantic and eastern Pacific oceans. This map is given here as reference for comparison with our final model obtained by co-estimation with the noise model.

20

The residuals to the fit to the data after the last step of the fitting process are mapped in Fig. 6, left. The largest anomalies, as the Bangui anomaly in central Africa or the Kursk anomaly in western Russia, are clearly identifiable on this residual map, although

25
<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /></figure>

<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>

<figure><img src="image_77.png" /></figure>

they are not associated with too large residuals. There are very large clusters of residuals at high latitudes, and some of these residuals obviously correspond to lithospheric magnetic anomalies – e.g. North America, Southern tip of Greenland, Northern Europe. This is an incentive to work with localized system of representation and to define local constraints. Here we want to keep the processing as simple as possible and did not follow such approaches. It should be noted, however, that the amplitude of the residuals are clearly smaller than 1 nT and that there is only few traces of the “along track” noise in these residuals. The effect of the smoothing on the model remains acceptable.

Figure 6 right, shows the power spectra of the field model $B_j$ and of the noise model $\tilde{B}_j$. Also plotted is the spectrum from MF7. The damping parameter $\lambda$ in Eq. (17) has been adjusted to $\lambda = 4.0 \times 10^{-5}$ such that the power spectrum does not present excessively high values at high degrees. Overall, the derived map has the same level of energy than MF7 up to degree 100. Above that degree the spectra is clearly decreasing. Our opinion is that we are reaching at these SH degrees the maximum “global” resolution of the CHAMP data selected and processed following the technique described above. Improvements are probably still possible locally, particularly above the largest anomalies seen as Bangui and Kursk anomalies.

Figure 7 maps, on the left, the noise model $\tilde{B}_j$, and, on the right the perturbation model defined in Eq. (A1). The noise model presents the expected East-West high frequency oscillations. The map cannot be directly compared with Fig. 5 because the patterns of the oscillations in Fig. 7 correspond to the noise present in $\tilde{B}_j$: the Fig. 5 is only a smoothed version of it. The perturbation model (Fig. 7, right) is dominated by a dipole term consistent with un-modelled contributions generated in 1-D conductive layers of the Earth by a large-scale, rapidly varying external field. Although this large-scale field is dominant, higher spherical harmonic contributions exist in the perturbation model and are determinant for the success of the post-processing.

Our final result is a map of the vertical down component of the lithospheric field calculated at the Earth’s surface (see Fig. 8). The map includes all SH degrees of the lithospheric field model. The model is displayed with two different central meridians for
<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /></figure>

a better view of the anomaly patterns. The anomaly patterns are not as clearly defined as in MF7, but in numerous areas – e.g. the northern pacific, the map resulting from our processing is remarkably detailed. However, in the present case, the only difference with regards to a straightforward least-squares approach is the co-estimation of the noise model. In particular, there are no pre-processing steps such as data leveling (or micro-levelling) with mostly unknown concequences on the final map, and the only data used are the CHAMP satellite data. We have made numerous experiments, and it appears that the determinant step for the final quality of the map is the data selection used to build the model $\hat{B}_i$. Out of all these tries, the maps presenting the lowest level of noise are systematically the outputs of the step-2 of our processing. We decided not to show these results here because they are not consistent with the noise model presented in Appendix A that assumes a model derived through a non-regularized scheme. It is however an approach worth studying: There are no major difficulties in estimating what the noise model should be for a lithospheric model built using a regularized least-squares process.

## 4 Conclusions

We have calculated the gauss coefficients describing the noise leaking in lithospheric magnetic field model when derived from satellite data. The noise models were derived to cover two cases: first when exclusively the radial components of the satellite data are used and second when all three components are used. The first case would be primarily applicable to gravity data, whereas the second, as we used it here, is better suited for magnetic data. We made several strong hypotheses to obtain these results. Particularly, we consider that the orbits are exactly polar, that they are at constant radius and that the sampling rate along an orbit is “ideal” – i.e. the relation $\langle P_i^{[m]} \rangle, P_i^{[m]} \rangle \propto \delta_{jj'}$ is verified. We also make the assumption that the lithospheric field model is derived through a simple un-regularized least-squares process. We insist here on the fact that
<figure><img src="image_83.png" /></figure>

<figure><img src="image_84.png" /></figure>

<figure><img src="image_85.png" /></figure>

<figure><img src="image_86.png" /></figure>

<figure><img src="image_87.png" /></figure>

the noise models do not represent the expected noise in the satellite data but the noise
in the derived lithospheric models.

It is interesting to notice that the amplitude of the noise generated depends on the
variance of the random variable $\chi_m^k$, that itself depends on the variance of the external
field scaling factor $\eta$ and the number of orbits $M$ (see Table 1). Therefore the usual
choice of rejecting a significant part of the data because of its level of noise is ques-
tionable. For example when dealing with magnetic data, rejecting a full year of satellite
data because of the high level of magnetic activity is unlikely to reduce the noise level
in the model since the ratio $\frac{\nu^\eta}{M}$ generally does not get smaller. We cannot comment
however on a data rejection criteria based on the satellite altitude.

Another remarkable property of the noise models is their weak dependence with re-
gard to the source of the noise. We used here perturbation models either from internal
or external origin, but both lead to similar noise models. The same developments could
be done for a noise described by spherical harmonics without reference to any specific
source. For the case where only radial component data are used (Eqs. 10 and 11),
such an hypothesis would not make any difference.

In the application to real data, the noise models were used in a post-processing
scheme. The reason for this choice is that we did not know what kind of perturbation
model $B_p(\theta, \phi, r)$ should be used. We have seen that for deriving a lithospheric field
model, a dipole perturbation field is not leading to the best results. In an ideal case
where the perturbation model is known, the best approach to the problem would be
to build a covariance matrix $C_n$ for the noise from the variances given in Table 1 and
the Eq. (A12). Such a covariance matrix can then be used as a regularization matrix
in the least-squares fit of the lithospheric field model to the satellite data. However,
even if the information provided by the estimated variances has not been used in our
post-processing scheme, the resulting lithospheric field model is nonetheless much
improved compared to what can be obtained through a simple smoothing (see the
differences between Figs. 8 and 5).
One can question if parts of the lithospheric field can be removed by our post-processing steps and contribute to the noise model. To test such a hypothesis we have applied the step-2 of our processing using a synthetic lithospheric field model defined in the framework of the end-to-end Swarm Level-2 simulator. The amplitude of the obtained noise model was less than 0.8 nT anywhere on the Earth surface in all three magnetic components. This does not stand as a proof, because results may be different for other lithospheric models. Nonetheless, it seems that only a very small part of the lithospheric field model can be interpreted as noise.

The work presented here opens numerous possibilities for processing data acquired along linear paths, as satellite data. The major difficulty when dealing with large data set, is to handle the correlated errors. Facing this problem we have here simply calculated how this correlated noise affects the derived model through a least-squares process. Extending this to regularized least-squares approaches is certainly possible. The same technique can be applied for calculating small-scale secular variations from satellite data, or to process yearly estimates of the core field. The technique is also applicable for airborne data using any local system of representation rather than spherical harmonics. Interesting developments are possible through the design of local filters. The link with oriented wavelets on the sphere is also promising.

# Appendix A

## Noise model for three component vector data

We follow here the same developments as in Sect. 2 but consider the case where the perturbation field is of internal origin and the three magnetic vector components are used. The perturbation field writes:

$$B_p(\theta, \phi, r) = -\nabla \left[ a \sum_{l,m}^{L} \left(\frac{a}{r}\right)^{l+1} I_l^m Y_l^m(\theta, \phi) \right]. \quad (A1)$$

<figure><img src="image_88.png" /></figure>

<figure><img src="image_89.png" /></figure>

<figure><img src="image_90.png" /></figure>

<figure><img src="image_91.png" /></figure>
<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /></figure>

<figure><img src="image_94.png" /></figure>

<figure><img src="image_95.png" /></figure>

It is scaled at each orbit by a factor $\eta$ and is fitted by least-squares with a field of internal origin constant in time. Therefore we minimize the functional:

$$ \Phi = \sum_{i,j} w_i |\tilde{B}_i(\theta_i, \phi_i, r) - \eta_j \cdot B_p(\theta_i, \phi_i, r)|^2, \quad (A2) $$

where the noise model $\tilde{B}_i$ is defined by:

$$ \tilde{B}_i(\theta, \phi, r) = -\nabla \left[ a \sum_{l,m}^{L=120} \left(\frac{a}{r}\right)^{l+1} \tilde{g}_l^m Y_l^m(\theta, \phi) \right]. \quad (A3) $$

This leads to a linear system equivalent to Eq. (5), where the left hand side writes:

$$ \{\mathbf{A}^t \mathbf{A}\}_{l,m,l',m'} = \left(\frac{a}{r}\right)^{l+l'+4} (l+1)(l'+1) \sum_{i,j} w_i Y_l^m(\theta_i, \phi_i, r) Y_{l'}^{m'}(\theta_i, \phi_i, r) \\ + \left(\frac{a}{r}\right)^{l+l'+4} \sum_{i,j} w_i \nabla_h Y_l^m(\theta_i, \phi_i, r) \cdot \nabla_h Y_{l'}^{m'}(\theta_i, \phi_i, r). \quad (A4) $$

The operator $\nabla_h$ is the horizontal gradient on a sphere of unit radius. The first term in the right hand side does not present difficulties. For the second we use the identity:

$$ Y_l^m Y_{l'}^{m'} = \sum_{L=|l-l'|}^{|l+l'|} \sum_M C_{l,l',m,m'}^{L,M} Y_L^M. \quad (A5) $$

Applying twice the gradient operator gives:

$$ \nabla_h Y_l^m \cdot \nabla_h Y_{l'}^{m'} = \frac{l(l+1) + l'(l'+1)}{2} Y_l^m Y_{l'}^{m'} - \frac{1}{2} \sum_{L=|l-l'|}^{|l+l'|} \sum_M C_{l,l',m,m'}^{L,M} L(L+1) Y_L^M. \quad (A6) $$

The Eq. (A4) becomes:

$$ \{\mathbf{A}^t \mathbf{A}\}_{l,m,l',m'} = \left(\frac{a}{r}\right)^{l+l'+4} \frac{(l+l'+1)(l+l'+2)}{2} \sum_{i,j} w_i Y_l^m(\theta_i, \phi_i, r) Y_{l'}^{m'}(\theta_i, \phi_i, r) \\ - \left(\frac{a}{r}\right)^{l+l'+4} \sum_{L=|l-l'|}^{|l+l'|} \sum_M C_{l,l',m,m'}^{L,M} \frac{L(L+1)}{2} \sum_{i,j} w_i Y_L^M(\theta_i, \phi_i, r). \quad (A7) $$
<figure><img src="image_96.png" /></figure>

<figure><img src="image_97.png" /></figure>

<figure><img src="image_98.png" /></figure>

<figure><img src="image_99.png" /></figure>

Defining $\Pi_{mm'}$ as in Eq. (7) gives in the limit of a large number $M$ of orbits $\Pi_{mm'} \approx (\frac{1}{2} + \frac{1}{2}\delta_{m0})\delta_{mm'}$.

Further, the weights $w_i$ are chosen such that $\sum_i w_i P_i^m(\cos\theta_i)P_{i'}^{m'}(\cos\theta_{i'}) = \frac{4-2\delta_{m0}}{2l+1}\delta_{ll'}$, which reduces for $l' = m' = 0$ to $\sum_i w_i P_i^m(\cos\theta_i) = 2\delta_{m0}\delta_{l0}$.

As a consequence, for the second term in the right hand side of Eq. (A7) only the term $L=0$ remains. It therefore vanishes because of the factor $L(L+1)$ and we obtain:

$$ \{\mathbf{A}^t \mathbf{A}\}_{l,m,l',m'} = 2M \left(\frac{a}{l}\right)^{2l+4} (l+1) \delta_{ll'} \delta_{mm'} \quad (A8) $$

The matrix $\mathbf{A}^t \mathbf{A}$ is therefore diagonal: The discrete summations in Eq. (A4) are equivalent to continuous integrations.

The product $\{\mathbf{A}^t \mathbf{b}\}_{l,m}$ in the right hand side of Eq. (5) now writes:

$$ \{\mathbf{A}^t \mathbf{b}\}_{l,m} = M \sum_{n,k} \left(\frac{a}{r}\right)^{l+n+4} \left\{ (l+1)(n+1) \sum_{i,j} \eta_j w_i Y_l^m(\theta_i, \phi_j, r) Y_n^k(\theta_j, \phi_j, r) + \sum_{i,j} \eta_j w_i \nabla_h Y_l^m(\theta_i, \phi_j, r) \cdot \nabla_h Y_{l'}^{m'}(\theta_i, \phi_j, r) \right\}. \quad (A9) $$

We further introduce the variable $\dot{\chi}_m^k$ defined by:

$$ \dot{\chi}_m^k = \begin{cases} \chi_{-m}^{-k} & \text{if } mk > 0 \\ -\chi_{-m}^{-k} & \text{if } mk < 0 \\ 0 & \text{if } mk = 0, \end{cases} \quad (A10) $$

where the expression of $\dot{\chi}_m^k$ is given in Eq. (9). The Eq. (A9) becomes:

$$ \{\mathbf{A}^t \mathbf{b}\}_{l,m} = M \sum_{n,k} \dot{\chi}_n^k \left(\frac{a}{r}\right)^{l+n+4} \left\{ (l+1)(n+1) \langle P_n^{[k]}, P_l^{[m]} \rangle \chi_m^k + \langle \partial_\theta P_n^{[k]}, \partial_\theta P_l^{[m]} \rangle \chi_m^k + \langle \frac{|k|P_n^{[k]}}{\sin\theta}, \frac{|m|P_l^{[m]}}{\sin\theta} \rangle \dot{\chi}_m^k \right\}, \quad (A11) $$

leading when combined with Eq. (A7) to:

$$ \begin{aligned} \tilde{g}_l^m = & \sum_{n,k} \dot{\chi}_n^k \left(\frac{a}{r}\right)^{n-l} \left\{ \frac{n+1}{2} \langle P_n^{[k]}, P_l^{[m]} \rangle \chi_m^k + \frac{1}{2l+2} \langle \partial_\theta P_n^{[k]}, \partial_\theta P_l^{[m]} \rangle \chi_m^k + \frac{1}{2l+2} \langle \frac{|k|P_n^{[k]}}{\sin\theta}, \frac{|m|P_l^{[m]}}{\sin\theta} \rangle \dot{\chi}_m^k \right\}. \end{aligned} \quad (A12) $$
<figure><img src="image_100.png" /></figure>

<figure><img src="image_101.png" /></figure>

<figure><img src="image_102.png" /></figure>

<figure><img src="image_103.png" /></figure>

<figure><img src="image_104.png" /></figure>

The $L(L+2)$ gauss coefficients $\tilde{g}_i^m$ of the noise model can be represented by only $N(N+2)$ coefficients $\iota_n^k$ of the perturbation model and $(2N+1)(2L+1)-2N^2$ independent random variables $\chi_m^k$ where all symmetry properties have been accounted for.

***Acknowledgements.*** We would like to acknowledge the work of CHAMP satellite data processing team. M. R. was supported through the DFG grant: LE 2477/3-1 in the framework of the SPP-1488 “Planetary magnetism”. IPGP contribution number: 3340.

## References

Hamoudi, M., Quesnel, Y., Dyment, J., and Lesur, V.: Aeromagnetic and Marine Measurements, in: Geomagnetic Observations and Models, edited by: Mandea, M. and Korte, M., vol. 5, IAGA book series, chap. 4, 57–105, Springer, doi:10.1007/978-90-481-9858-0, 2010. 1347

10

Holme, R.: Modelling of attitude error in vector magnetic data: application to Ørsted data, Earth Planets Space, 52, 1187–1197, 2000. 1347

Holme, R. and Bloxham, J.: The treatment of attitude errors in satellite geomagnetic data, Phys. Earth Planet. In., 98, 221–233, 1996. 1347

15

Kusche, J.: Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models, J. Geodesy, 81, 733–749, doi:10.1007/s00190-007-0143-3, 2007. 1347

Langel, R. A., Estes, R. H., and Sabaka, T. J.: Uncertainty estimates in geomagnetic field modelling, J. Geophys. Res., 94, 12281–12299, 1989. 1347

20

Lesur, V., Wardinski, I., Rother, M., and Mandea, M.: GRIMM – The GFZ Reference Internal Magnetic Model based on vector satellite and observatory data, Geophys. J. Int., 173, 382–394, doi:10.1111/j.1365-246X.2008.03724.x, 2008. 1356

Lesur, V., Wardinski, I., Hamoudi, M., and Rother, M.: The second generation of the GFZ Reference Internal Magnetic field Model: GRIMM-2, Earth Planets Space, 62, 765–773, doi:10.5047/eps.2010.07.007, 2010. 1356

25

Lesur, V., Wardinski, I., Hamoudi, M., Rother, M., and Kunagu, P.: Third version of the GFZ Reference Internal Magnetic Model: GRIMM-3, in: IUGG2011, MR207, Melbourne, Australia, 2011. 1356
<figure><img src="image_105.png" /></figure>

<figure><img src="image_106.png" /></figure>

<figure><img src="image_107.png" /></figure>

<figure><img src="image_108.png" /></figure>

<figure><img src="image_109.png" /></figure>

Lowes, F. J. and Olsen, N.: A more realistic estimate of the variances and systematic errors in spherical harmonic geomagnetic field models, Geophys. J. Int., 157, 1027–1044, 2004. 1348

Maus, S., Yin, F., Lühr, H., Manoj, C., Rother, M., Rauberg, J., Michaelis, I., Stolle, C., and Müller, R. D.: Resolution of direction of oceanic magnetic lineations by the sixth-generation lithospheric magnetic field model from CHAMP satellite magnetic measurements, Geochem. Geophy. Geosy., 9, Q07021, doi:10.1029/2008GC001949, 2008. 1347, 1348

Olsen, N., Lühr, H., T. Sabaka, T., Michaelis, I., Rauberg, J., and Tøffner-Clausen, L.: CHAOS-4 – A high-resolution geomagnetic field model derived from low-altitude CHAMP data, in: AGU Fall Meeting, Poster GP21A-0992, San Francisco, 2010. 1348, 1356

Reigber, C., Lühr, H., Schwintzer, P., and Wickert, J., eds.: Earth Observation with CHAMP Results from three years in orbit, Springer-Verlag, 2005. 1348

Rygaard-Hjalsted, C., Constable, C. G., and Parker, R. L.: The influence of correlated crustal signals in modelling the main geomagnetic field, Geophys. J. Int., 130, 717–726, 1997. 1347

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

5

10
<figure><img src="image_110.png" /></figure>

<figure><img src="image_111.png" /></figure>

<figure><img src="image_112.png" /></figure>

<figure><img src="image_113.png" /></figure>

<table><caption>**Table 1.** Estimated variance of $\chi_m^k$</caption><thead><tr><th></th><th>m = k</th><th>m = -k</th><th>m ≠ 0 and m ≠ k</th></tr></thead><tbody><tr><th>k = 0</th><td>$\frac{v^k}{M^k}$</td><td>-</td><td>$\frac{v^k}{2M^k}$</td></tr><tr><th>k ≠ 0</th><td>$\frac{3v^k}{8M}$</td><td>$\frac{v^k}{8M}$</td><td>$\frac{v^k}{4M}$</td></tr></tbody></table>
<table><caption>**Table 2.** Thresholds and misfit values obtained when estimating $\hat{B}_i$ from selected satellite data. Mid and low latitudes are defined by magnetic latitudes in between ±55 deg. Values are given in nT.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Mid and low latitudes</th><th colspan="3">High latitudes</th></tr><tr><th>X<sub>SM</sub></th><th>Y<sub>SM</sub></th><th>Z<sub>SM</sub></th><th>X<sub>SP</sub></th><th>Y<sub>SP</sub></th><th>Z<sub>SP</sub></th></tr></thead><tbody><tr><th>Threshold</th><td>9.0</td><td>8.5</td><td>10.5</td><td>36</td><td>27</td><td>36</td></tr><tr><th>Misfit</th><td>2.47</td><td>2.30</td><td>2.53</td><td>13.52</td><td>10.79</td><td>10.49</td></tr></tbody></table>

<figure><img src="image_114.png" /></figure>

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

<figure><img src="image_117.png" /></figure>
<figure><img src="image_118.png" /><figcaption>**Fig. 1.** Histograms of the random variable $\chi_m^k$ for several values of $k$ and $m$. Is also plotted the dashed curve $\frac{M-S}{\sqrt{2\pi\nu^2}}\exp\{-e^2/(2\nu^2)\}$ where $S$ is the histogram step length and $e$ the error.</figcaption></figure>

<figure><img src="image_119.png" /></figure>

<figure><img src="image_120.png" /></figure>

<figure><img src="image_121.png" /></figure>

<figure><img src="image_122.png" /></figure>
<figure><img src="image_123.png" /></figure>

<figure><img src="image_124.png" /></figure>

<figure><img src="image_125.png" /></figure>

<figure><img src="image_126.png" /></figure>

<figure><img src="image_127.png" /></figure>

<figure><img src="image_128.png" /><figcaption>**Fig. 2.** Left: mapping at $r' = 6371.2$ km of the model defined in Eq. (13) where $l_1^1 = 1$, $r = 6671.2$ km – i.e. 300 km altitude, and the random variables $\chi_m^1$ have a variance defined in table 1 using $\nu^\eta = M$. Right: associated power spectrum. The dashed line is proportional to $(r/r^2)^2$.</figcaption></figure>
<figure><img src="image_129.png" /></figure>

<figure><img src="image_130.png" /></figure>

<figure><img src="image_131.png" /></figure>

<figure><img src="image_132.png" /></figure>

<figure><img src="image_133.png" /></figure>

<figure><img src="image_134.png" /><figcaption>Fig. 3. Mapping at $r' = 6371.2$ km of the field of internal origin (left) and the resulting noise model (right). The data acquisition radius has been set to $r = 6671.2$ km – i.e. 300 km altitude, and the random variables $\chi_m^k$ have the variances defined in Table 1 with the ratio $\sqrt{M}$ set to 1.</figcaption></figure>
<figure><img src="image_135.png" /></figure>

<figure><img src="image_136.png" /></figure>

<figure><img src="image_137.png" /></figure>

<figure><img src="image_138.png" /></figure>

<figure><img src="image_139.png" /></figure>

<figure><img src="image_140.png" /><figcaption>**Fig. 4.** Left: power spectra of the lithospheric field model (solid line) and of CHAOS-4b model (dashed line) calculated at the Earth's surface (i.e. *r* = 6371.2 km). Right: mapping of the vertical down component of the (noisy) lithospheric field *B̂*ᵢ at *r* = 6671.2 km. The largest magnetic anomalies dominate, but the “along track” noise is nonetheless visible over oceanic areas.</figcaption></figure>
<figure><img src="image_141.png" /></figure>

<figure><img src="image_142.png" /></figure>

<figure><img src="image_143.png" /></figure>

<figure><img src="image_144.png" /></figure>

<figure><img src="image_145.png" /></figure>

<figure><img src="image_146.png" /><figcaption>**Fig. 5.** Map of the vertical down component of the lithosphere magnetic field model at *r* = 6371.2 km radius derived after the step-1 of the processing chain. This corresponds to a smoothed model without co-estimation of the noise model. It is given here as a reference to be compared with Fig. 8. Along track noise is particularly visible around Antarctica, and in the Indian, Atlantic and eastern Pacific oceans.</figcaption></figure>
<figure><img src="image_147.png" /></figure>

<figure><img src="image_148.png" /></figure>

<figure><img src="image_149.png" /></figure>

<figure><img src="image_150.png" /></figure>

<figure><img src="image_151.png" /><figcaption>**Fig. 6.** Left: residuals map to the final model fit to (noisy) lithospheric field $\hat{B}_i$ at $r = 6671.2$ km. The residuals have been scaled by a factor 10. At mid-latitudes the largest residuals are associated with the strong magnetic anomalies. The along track noise has been fitted by the noise model. Right: power spectra of the lithospheric field model $B_i$ (solid line), of MF7 (dashed line), and the noise model (dotted line) calculated at the Earth's surface (i.e. $r = 6371.2$ km radius).</figcaption></figure>
<figure><img src="image_152.png" /></figure>

<figure><img src="image_153.png" /></figure>

<figure><img src="image_154.png" /></figure>

<figure><img src="image_155.png" /></figure>

<figure><img src="image_156.png" /></figure>

<figure><img src="image_157.png" /><figcaption>Fig. 7. Map of the vertical down components of, left, the noise model, right, the perturbation model. Both maps have been calculated at *r* = 6371.2 km radius. By definition the perturbation model is very smooth in longitude, but that does not preclude a large complexity for the noise model.</figcaption></figure>
<figure><img src="image_158.png" /></figure>

<figure><img src="image_159.png" /></figure>

<figure><img src="image_160.png" /></figure>

<figure><img src="image_161.png" /></figure>

<figure><img src="image_162.png" /></figure>

<figure><img src="image_163.png" /><figcaption>**Fig. 8.** Map of the vertical down components of the final lithospheric field model $B_j$. The map has been calculated at the Earth's surface (6371.2 km). Although some noise is still visible in the northern Atlantic and over the southern polar cap, the noise level over mid latitudes has been greatly reduced. Anomalies are particularly well defined over continents, and Indian and Pacific oceans.</figcaption></figure>