# Post-processing scheme for modelling 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 modelling the lithospheric magnetic field. Solid Earth, 2013, 4, pp.105-118. 10.5194/se-4-105-2013. insu-03581800

HAL Id: insu-03581800

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

Submitted on 20 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>

# Post-processing scheme for modelling the lithospheric magnetic field

V. Lesur¹, M. Rother¹, F. Vervelidou², M. Hamoudi³, and E. Thébault²

¹Helmholtz Centre Potsdam, GFZ German Research centre for Geosciences, Telegrafenberg 14473, Germany

²Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Univ. Paris Diderot, UMR7154 CNRS, 75005 Paris, France

³Université des sciences et de la technologie, 16111 Bab Ezzouar, El-Alia Alger, Algeria

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

Received: 9 October 2012 – Published in Solid Earth Discuss.: 30 October 2012

Revised: 29 January 2013 – Accepted: 31 January 2013 – Published: 7 March 2013

**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 is a random variable, normally distributed with zero mean. 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 degree 120. The model is in agreement with other existing models. The technique can, in principle, be extended to all sorts 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 properly handling 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 no 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 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 mostly 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 for example Langel et al. (1989); Holme and Bloxham (1996); Rygaard-Hjalsted et al. (1997); Holme (2000), 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 (see for a review, e.g. Hamoudi et al., 2010). The method has also been adapted to satellite magnetic data. In such cases a large-scale field of external origin is fitted to a data set made of only
a few tracks. This allows us to successfully derive magnetic field models of the lithosphere to a 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 Thébault 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 feature of this so-called “along-track filtering” approach is the impossibility to estimate how much the processing applied will distort 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 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. This therefore lead 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 obtained 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 that is not published but the MF6 is presented in Maus et al. (2008); CHAOS-4; Olsen et al., 2010b) as well as older models (see for a review Thébault 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 detail the two-step process towards the final lithospheric field model; the resulting model is then discussed. We conclude 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-time periods. For simplicity we will assume that a track follows a meridian – i.e. it corresponds to a single longitude value – which is a reasonable approximation for near-polar orbiting satellites. Several magnetic field sources are contributing to these data (Hulot et al., 2007), typically the core and lithospheric fields, the ionospheric and magnetospheric fields, and the fields generated by field aligned currents. Other contributions exist, as the field induced in the conductive layers of the Earth, but they are of much weaker amplitudes. Mathematical models are available for all these contributions and can be subtracted from the data, leaving residuals due mainly to the limited precision of these models. In particular, the description of the external field is not very accurate and the residuals obtained along that half-orbit track contain relatively long wavelengths. We assume that these residuals 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 can be written

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

where $\epsilon_n^k$ is the Gauss coefficient of degree $n$ and order $k$, $a = 6371.2$ km is the Earth’s reference radius, and $Y_n^k(\theta, \phi)$ are the Schmidt semi-normalized spherical harmonics (SHs). We use throughout 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 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 (Olsen et al., 2010a) – i.e. a single half-orbit – hence we can fit by least-squares the residuals defined in Eq. (1) with the lithospheric model given in Eq. (2) and find a non-zero solution. This least-squares solution is found by minimizing 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 *η*<sub>j</sub> 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

$$ \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 *Φ* for the Gauss coefficients $\tilde{g}_l^m$ leads to a system of equations:

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

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

$$ \{\mathbf{A}^t \mathbf{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). \quad (7) $$

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 (8) $$

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 (9) $$

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 \eta_i & \text{if } m, k \ge 0 \\ \frac{1}{M} \sum_{i=1}^{M} \cos m\phi_i \sin |k|\phi_i \eta_i & \text{if } m k < 0 \\ \frac{1}{M} \sum_{i=1}^{M} \sin |m|\phi_i \sin |k|\phi_i \eta_i & \text{if } m, k < 0. \end{cases} \quad (10) $$

As for $\Pi_{mm'}$, it is symmetric relative to its subscripts: $\chi_m^k = \chi_k^m$. We note at this point that it is important to have the $\eta_i$ constant along half-orbits; otherwise, the summations over latitudes and longitudes could not be separated.

For very large number *M* of orbits uniformly distributed along longitudes, the quantity $\Pi_{mm'}$ tends to a $\delta$-function – i.e. $\Pi_{mm'} \simeq (\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

<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 \neq 0$ and $m \neq k$</th></tr></thead><tbody><tr><td>$k=0$</td><td>$\frac{v^\eta}{M}$</td><td>-</td><td>$\frac{v^\eta}{2M}$</td></tr><tr><td>$k \neq 0$</td><td>$\frac{3v^\eta}{8M}$</td><td>$\frac{v^\eta}{8M}$</td><td>$\frac{v^\eta}{4M}$</td></tr></tbody></table>

$\langle P_l^{|m|}, P_{l'}^{|m'|} \rangle = \frac{4-2\delta_{m0}}{2l+1} \delta_{ll'}$. The value given to the weights $w_i$ is less important than insuring the orthogonality of the Legendre functions through the product $\langle P_l^{|m|}, P_{l'}^{|m'|} \rangle$. This is also what a modeler tries to achieve when building a lithospheric magnetic field model from real data. However, assuming that both $\Pi_{mm'}$ and $\langle P_l^{|m|}, P_{l'}^{|m'|} \rangle$ can be regarded as $\delta$-functions, it is easy to see from Eq. (6) that the product matrix $\mathbf{A}^t \mathbf{A}$ is diagonal. Now turning to Eqs. 9, 10, 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 obtained by combining Eqs. (5), (6) and (9):

$$ \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 (11) $$

They correspond to the noise in a lithospheric field model that would be generated by un-modelled 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 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-modelled induced fields generated in the conductive layers of the Earth by rapid variations of the external fields. It gives the following result:

$$ \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 t_n^k \chi_m^k, \quad (12) $$

where $t_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 $v^\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 $v^\chi$ that depends on $v^\eta$, the number of half-orbits *M*, and the orders *k* and *m*. Possible values of the variance $v^\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. They can also be estimated analytically as shown in Appendix B. Figure 1 presents the histograms for several values of *m* and *k*.

In the remaining parts of this section we consider only the noise model given by Eq. (12). The general behaviour of the noise characterized by Eqs. (12) and (11) 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 $t_n^k$ with $(2N+1)(2L+1) - 2N^2$ random variables $\chi_m^k$ through Eq. (12). For small values of $N$ – e.g. $N=10$ – there is a very significant reduction of number of parameters, but Eq. (12) is non-linear.

## 2.2 Examples

In order to understand the main characteristics of the noise model defined by Eqs. (2) and (12), we present in this section the results of forward modelling calculations for a given choice of Gauss coefficients $t_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 (Piessens et al., 1983; Kahaner et al., 1989) 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 $t_1^1 = 1$ nT and $t_n^k = 0$ nT 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. (12) 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 t_1^1 \chi_m^1, \quad (13) $$

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

$$ \tilde{B}_{ri}(\theta, \phi, r') = t_1^1 2 \left(\frac{a}{r'}\right)^3 \sum_{l,m} \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 (14) $$

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 wave-

lengths dominate the model due to the ratio $r/r'$ raised to the power $l-1$ in the right-hand side of Eq. (14).

In Fig. 2, the model defined by Eq. (14) is mapped for the model coefficient $t_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 $v^\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 on this sign in Eq. (14). The obtained symmetry of the model is due to the product $\langle P_1^{|1|}, P_l^{|m|} \rangle$, which vanishes if the Legendre function $P_l^{|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 $t_1^0 = 1$ nT would have been chosen in place of $t_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' = 6371.2$ km is also plotted in Fig. 2. It presents some variability due to the use of a single SH in Eq. (13). Nonetheless, the behaviour is generally along a $(r/r')^{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 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 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 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) t_n^k Y_n^k(\theta, \phi). \quad (15) $$

The model was built simply by defining a circular ridge in a geomagnetic system of coordinates, which was then rotated in the usual reference frame. 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, derived from Eq. (12), is mapped in Fig. 3, right side. This model is also fairly well
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Histograms of the random variable $\chi_m^k$ for several values of k and m. Also plotted is the dashed curve $\frac{M \cdot S}{\sqrt{2\pi v\chi}} \exp\{-e^2/(2v\chi)\}$ where S is the histogram step length and e is the error.</figcaption></figure>

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 as 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 from these

results that there is no need to describe precisely the longitudinal dependence of the perturbation field to obtain a realistic noise model. Therefore, in Eq. (12), 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.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Left: Mapping at $r' = 6371.2$ km of the model defined in Eq. (14) where $t_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 $v^\eta = M$. Right: Associated power spectrum. The dashed line is proportional to $(r/r')^2$.</figcaption></figure>

<figure><img src="image_5.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 $v_M^\eta$ set to 1.</figcaption></figure>

### 3 Application to magnetic models of the lithosphere

We are now using the results presented in the previous section to derive a model of the magnetic field generated in the lithosphere from real CHAMP satellite data. The process we applied to calculate such a model is in two stages. First, we estimate a rough lithospheric field model from satellite data using a usual approach and a straightforward least-squares process (e.g. Lesur et al., 2008). Second, in the post-processing stage, we co-estimate a new lithospheric field model and a model of the noise where the output model of the first stage is used as data. The final results depend on the

processes applied during the two stages and therefore both are described in independent subsections below.

In order to avoid confusion between the different models, we use the following notations for the fields and Gauss coefficients:

- the noisy lithospheric model, output of the first stage, is denoted using $\hat{a}$ – e.g. $\hat{B}_i$ for the magnetic field vector,

- the lithospheric field model output of the post-processing (second) stage does not have any distinctive sign – e.g. $B_i$.
<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>

– the noise model is denoted, in the same way as in the previous section, using a ~ – e.g. $\tilde{B}_i$.

## 3.1 First stage: Data set, data selection, model parameterization and model estimation

Three component vector magnetic readings acquired during the ten years of the German CHAMP satellite mission are used. The data are selected for night-time periods and magnetically quiet days, in the same way as data are selected for the GRIMM series of core field models (Lesur et al., 2008, 2010). However, here the three components of the vector data are used and data in single star camera mode are rejected, whereas in the GRIMM selection scheme only the X and Y SM components are selected at mid- and low-latitudes. 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., 2010b).

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 strongly affect the final lithospheric field model. At mid- and low-latitudes only few data are rejected, and those rejected 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 were adjusted such that the

lithospheric field model remained smooth while the misfit to the original satellite data stayed unchanged. 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. (16) below, was fitted through a simple least-squares process to the data. The data weights are set to mimic a homogeneous data repartition and therefore depend only on the inverse of the data density.

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

The power spectrum at the Earth surface of the resulting lithospheric field model is presented in Fig. 4, left side, 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 spectrum 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 side. 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 Second stage: 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 $l_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_l^m$ of the field model $B_i$ by the
<figure><img src="image_6.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 $\hat{B}_i$ 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_7.png" /></figure>

following 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 (17) $$

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_l^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 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 (18) $$

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 fields are 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. (12) 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. (12) 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 unknowns 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$ for the $\chi_m^k$ (i.e. 721 values for $L=120$). These numbers have to be compared with the number of unknowns 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_l^m$ by minimizing $\Phi$ (Eq. 17) imposing $\chi_m^k = 0$ for all possible $m$ and $k$ values.
<figure><img src="image_8.png" /></figure>

**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.

- Step 2: Keeping the $g_l^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_l^m$, $i_n^k$ and $\chi_m^k$ that minimize $\Phi$ in Eq. (17), starting from the output of step 2.

### 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 a power spectrum in its expected range. 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°.

The output of the step 1 described above 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.

The residuals to the fit to the data after the last step of the fitting process are mapped in Fig. 6, left side. 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 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 a 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.5 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 side, shows the power spectra of the field model $\mathbf{B}_i$ and of the noise model $\tilde{\mathbf{B}}_i$. 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 as 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, left side, maps the noise model $\tilde{\mathbf{B}}_i$, and, on the right side, maps 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 $\hat{\mathbf{B}}_i$; Fig. 5 is only a smoothed version of it. The perturbation model (Fig. 7, right side) is dominated by a dipole term consistent with unmodelled 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 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 resulting map 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 levelling (or micro-levelling) with mostly unknown consequences on the final map, and the only data used are the CHAMP satellite data. We have run numerous experiments, and it appears that the determinant step for the final quality of the map is the data selection used to build
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Left: Residuals map to the final model fit to (noisy) lithospheric field $\hat{B}_i$ at $r = 6671.2$ km. 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 and therefore does not appear in these residuals. 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_10.png" /><figcaption>**Fig. 7.** Left: Map of the vertical down components of (left) the noise model and (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>

the model $\hat{B}_i$. Out of all these trials, the maps presenting the lowest level of noise are systematically the outputs of step 2 of our iterative inversion process. 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 models 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, although applications to
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Map of the vertical down components of the final lithospheric field model $B_i$. 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>

vector gravimetry or gradiometry may be possible. 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_l^{|m|}, P_{l'}^{|m|} \rangle \propto \delta_{ll'}$ is verified. We also make the assumption that the lithospheric field model is derived through a simple un-regularized least-squares process. This latter approximation is well verified for our application, but the three former are rather rough – e.g. in our data set the altitudes varies between 480 km and 250 km. However, it appears that the final result does not suffer too much from these hypotheses. We insist here on the fact that the noise models do not represent the expected noise in the satellite data but the noise leaking into the derived lithospheric models.

It is interesting to note that the amplitude of the noise generated depends on the variance of the random variable $\chi_m^k$, which 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 questionable. 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 $v_M^\eta$ 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 regard 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. 11, 12), such a hypothesis would not make any difference.

In the application to real data, the noise models were estimated 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 model 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 Eq. (A12). Such a covariance matrix could 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 Fig. 8 and Fig. 5).

One can question if parts of the lithospheric field can be removed by our post-processing steps and contribute to the noise model. The lithospheric noise model derived is only a combination of spherical harmonics with some strong correlations between the Gauss coefficients. Therefore, there is no doubt that part of the true lithospheric magnetic field model can contribute to the noise model. It is, however, not possible to estimate a priori what this part is because it clearly depends on both the noisy lithospheric field model on which the post-processing is applied and on the true lithospheric field we want to estimate. In order to test our scheme, we have first applied the processing on a synthetic data set built on a Gauss-Legendre grid where both the lithospheric field model and the noise are known. We used only the radial component of the field and verified that the noisy lithospheric field model derived from these data was contaminated by a noise corresponding exactly to Eq. (12). However, the full inversion process revealed that part of the lithospheric field was seen as noise. We also applied step 2 of our processing using a noise-free synthetic lithospheric field model and the noise model defined by Eq. (A12). Here again, despite a noise-free data set, the lithospheric field is partly interpreted as noise.
The amplitude of the obtained noise model is relatively large where there are strong anomalies, more or less aligned along orbits. This is the case for the Kursk anomaly (51° N, 37° E), whereas the Bangui anomaly (4° N, 16° E) is apparently not affected by the processing. Outside a few localized areas, the noise model remained relatively small. This impossibility to properly separate the noise from the lithospheric field is a common limitation of all existing processing methods. In our specific scheme, the only way we can reduce this effect is by constraining the perturbation model. We therefore recommend that the post-processing is used only when the noise is clearly identifiable at the smallest wavelengths, constraining this way the perturbation model. Overall, the proposed post-processing probably performs better than other smoothing techniques but that should be tested on a case-by-case basis.

The work presented here opens numerous possibilities for processing data acquired along linear paths, such 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 can be written

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

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{\mathbf{B}}_i(\theta_i, \phi_j, r) - \eta_j \cdot \mathbf{B}_p(\theta_i, \phi_j, r)|^2, \quad (A2) $$

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

$$ \tilde{\mathbf{B}}_i(\theta, \phi, r) = -\nabla \left[ a \sum_{l,m} \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 can be written:

$$ \begin{aligned} \{\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_j, r) Y_{l'}^{m'}(\theta_i, \phi_j, r)\} + \left(\frac{a}{r}\right)^{l+l'+4} \\ & \sum_{i,j} \{w_i \nabla_h Y_l^m(\theta_i, \phi_j, r) \cdot \nabla_h Y_{l'}^{m'}(\theta_i, \phi_j, r)\}. \end{aligned} \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 term we use the following 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

$$ \begin{aligned} \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. \end{aligned} \quad (A6) $$

Equation (A4) then becomes

$$ \begin{aligned} \{\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_j, r) Y_{l'}^{m'}(\theta_i, \phi_j, 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_j, r). \end{aligned} \quad (A7) $$

Defining $\Pi_{mm'}$ as in Eq. (8) gives in the limit of a large number $M$ of orbits

$$ \Pi_{mm'} \approx \left(\frac{1}{2} + \frac{1}{2}\delta_{m0}\right)\delta_{mm'}. $$

Further, the weights $w_i$ are chosen such that

$$ \sum_i w_i P_l^m(\cos\theta_i) P_{l'}^{m'}(\cos\theta_i) = \frac{4 - 2\delta_{m0}}{2l + 1} \delta_{ll'}, $$

which reduces for $l' = 0$ to $\sum_i w_i P_l^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}{r}\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 can be written

$$ \begin{aligned} \{\mathbf{A}^t \mathbf{b}\}_{l,m} &= M \sum_{n,k}^{N} \left(\frac{a}{r}\right)^{l+n+4} \\ & \{ (l+1)(n+1) \sum_{i,j} \eta_j w_i Y_l^m(\theta_i, \phi_j, r) Y_n^k(\theta_i, \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) \}. \end{aligned} \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 $\chi_m^k$ is given in Eq. (10). Equation (A9) becomes

$$ \begin{aligned} \{\mathbf{A}^t \mathbf{b}\}_{l,m} = & M \sum_{n,k}^{N} \iota_n^k \left(\frac{a}{r}\right)^{l+n+4} \\ & \{ (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 \} \quad (\text{A11}) \\ & + \left( \frac{\langle |k|P_n^{[k]} |, \langle |m|P_l^{[m]} | \rangle \dot{\chi}_m^k}{\sin\theta}, \frac{\langle |m|P_l^{[m]} |, \sin\theta \rangle \dot{\chi}_m^k}{\sin\theta} \right) \end{aligned} $$

leading, when combined with Eq. (A7), to

$$ \begin{aligned} \tilde{g}_l^m = & \sum_{n,k}^N \iota_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 \right. \\ & \left. + \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 (\text{A12}) $$

The $L(L+2)$ Gauss coefficients $\tilde{g}_l^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.

# Appendix B

## Estimating the variance of $\chi_m^k$

The random variable $\chi_m^k$ has been defined above by

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

where $\eta_i$ is a random variable with zero mean and normally distributed with variance $v^\eta$. Let us consider first the simplest case where $m, k \ge 0$. If we introduce the two vectors $\mathbf{D} = [\cos m\phi_i \cos k\phi_i]_{i=1,M}$ and $\mathbf{N} = [\eta_i]_{i=1,M}$, then Eq. (B1) can be written

$$ \chi_m^k = \frac{1}{M} \mathbf{D}^t \cdot \mathbf{N}. \quad (\text{B2}) $$

With this definition it is clear that

$$ v^\chi = \frac{1}{M^2} \mathbf{D}^t \cdot \mathbf{C}_N \cdot \mathbf{D}, \quad (\text{B3}) $$

where $v^\chi$ is the variance of $\chi_m^k$ and $\mathbf{C}_N = v^\eta \mathbf{I}_d$, with $\mathbf{I}_d$ being the identity matrix. We therefore obtain

$$ v^\chi = \frac{v^\eta}{M^2} \left\{ \sum_{i=1}^{M} (\cos m\phi_i \cos k\phi_i)^2 \right\}, \quad (\text{B4}) $$

and using the identity $\cos a \cos b = \frac{1}{2}(\cos(a-b)+\cos(a+b))$ it follows that

$$ v^\chi = \frac{v^\eta}{4M^2} \sum_{i=1}^{M} \left\{ \begin{aligned} & (\cos(m+k)\phi_i)^2 + (\cos(m-k)\phi_i)^2 \\ & + 2(\cos(m+k)\phi_i \cos(m-k)\phi_i) \end{aligned} \right\} \quad (\text{B5}) $$

Over the CHAMP mission there is a large number of orbits and the $\phi_i$ are uniformly distributed between $[0; \pi]$. It follows that the last term on the right-hand side vanishes unless $k=0$ or $m=k=0$, and it can be verified that

$$ \sum_{i=1}^{M} (\cos n\phi_i)^2 = \sum_{i=1}^{M} (\sin n\phi_i)^2 = \frac{M}{2} \quad (\text{B6}) $$

as long as $n$ is not too large. Therefore comes the following results:

$$ - \text{ if } k=m=0, v^\chi = \frac{v^\eta}{M} $$

$$ - \text{ if } k=0 \text{ and } m \neq 0, v^\chi = \frac{v^\eta}{2M} $$

$$ - \text{ if } m=k \neq 0, v^\chi = \frac{3v^\eta}{8M} $$

$$ - \text{ if } m, k \neq 0, m \neq k, v^\chi = \frac{v^\eta}{4M} $$

Extending these results to the two other cases (i.e. $km < 0$ and $m, k < 0$) is straightforward.

*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.

Edited by: F. Speranza

# 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., IAGA book series, 5, 57–105, Springer, doi:10.1007/978-90-481-9858-0, 2010.

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

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

Hulot, G., Olsen, N., and Sabaka, T. J.: The present field, Treatise in Geophysics, Geomagnetism, Elsevier Ltd., Amsterdam, 2007.

Kahaner, D., Moler, C., and Nash, S.: Numerical Methods and Software, Prentice-Hall, ISBN 978-0-13-627258-8, 1989.

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.

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

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, doi:10.1111/j.1365-246X.2008.03724.x, 2008.
Lesur, V., Wardinski, I., Hamoudi, M., and Rother, M.: The second generation of the GFZ Reference Internal Magnetic field Model: GRIMM-2, Earth Planet. Space, 62, 765–773, doi:10.5047/eps.2010.07.007, 2010.

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.

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.

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. Geophys. Geosyst., 9, 949, doi:10.1029/2008GC001, 2008.

Olsen, N., Glassmeier, K.-H., and Jia, X.: Separation of the Magnetic Field into External and Internal Parts., Space Science Reviews, 152, 135–157, doi:10.1007/s11214-009-9563-0, 2010a.

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, 2010b.

Piessens, R., Doncker-Kapenga, E., Überhuber, C., and Kahaner, D.: Quadpack, A Subroutine Package for Automatic Integration, Springer Series in Computational Mathematics, 1, Springer-Verlag, ISBN 978-3-540-12553-2, 1983.

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

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.

Thébault, E., Purucker, M., Whaler, K., Langlais, B., and Sabaka, T.: The Magnetic field of the Earth’s lithosphere, Space Sci. Rev., 155, 95–127, doi:10.1007/s11214-010-9667-6, 2010.

Thébault, E., Vervelidou, F., Lesur, V., and Hamoudi, M.: The shortcomings of along-track satellite analysis in planetary magnetism, Geophys. J. Int., doi:10.1111/j.1365-246X.2011.05281.x, 2012.