# Wavelet frames: an alternative to spherical harmonic representation of potential fields

Aude Chambodut, I. Panet, M. Mandea, M. Diament, M. Holschneider, O.

Jamet

► To cite this version:

Aude Chambodut, I. Panet, M. Mandea, M. Diament, M. Holschneider, et al.. Wavelet frames: an alternative to spherical harmonic representation of potential fields. Geophysical Journal International, 2005, 163, pp.875-899. 10.1111/j.1365-246X.2005.02754.x. insu-03601100

HAL Id: insu-03601100

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

Submitted on 8 Mar 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
# Wavelet frames: an alternative to spherical harmonic representation of potential fields

A. Chambodut,<fnref n="1" />* I. Panet,<fnref n="2" />,<fnref n="3" /> M. Mandea,<fnref n="1" />† M. Diament,<fnref n="2" /> M. Holschneider<fnref n="4" /> and O. Jamet<fnref n="3" />

Accepted 2005 July 19. Received 2005 June 1; in original form 2004 September 14

## SUMMARY

Potential fields are classically represented on the sphere using spherical harmonics. However, this decomposition leads to numerical difficulties when data to be modelled are irregularly distributed or cover a regional zone. To overcome this drawback, we develop a new representation of the magnetic and the gravity fields based on wavelet frames.

In this paper, we first describe how to build wavelet frames on the sphere. The chosen frames are based on the Poisson multipole wavelets, which are of special interest for geophysical modelling, since their scaling parameter is linked to the multipole depth (Holschneider et al.). The implementation of wavelet frames results from a discretization of the continuous wavelet transform in space and scale. We also build different frames using two kinds of spherical meshes and various scale sequences. We then validate the mathematical method through simple fits of scalar functions on the sphere, named ‘scalar models’. Moreover, we propose magnetic and gravity models, referred to as ‘vectorial models’, taking into account geophysical constraints. We then discuss the representation of the Earth’s magnetic and gravity fields from data regularly or irregularly distributed. Comparisons of the obtained wavelet models with the initial spherical harmonic models point out the advantages of wavelet modelling when the used magnetic or gravity data are sparsely distributed or cover just a very local zone.

**Key words:** magnetic and gravity field, spherical harmonics, wavelets.

## 1 INTRODUCTION

Magnetic and gravity observations are of great importance for the understanding of geodynamic activity of our planet. Measurements of the Earth's magnetic and gravity fields undertaken by satellites (without forgetting those on land, sea and air) are of particular interest, as they provide a global and uniform survey of these fields and of their temporal evolution.

Models of the **magnetic field** have been derived by means of several Earth's satellite missions, which have been carrying magnetic sensors. Satellite-borne magnetometers provide information on strength and direction of the internal and external Earth's magnetic field and its time variations. The Earth is surrounded by a large and complicated field caused to a large extent by a dynamo operating

in the fluid core. Currents flowing in the ionosphere, magnetosphere and oceans and magnetized rocks also influenced the geomagnetic field.

Three magnetic missions (Ørsted—launched in 1999, CHAMP and SAC-C—launched in 2000) have collected measurements providing new insights into the composition and the processes in the interior, and surrounding of the planet. These observations are also used in a range of applications, including navigation systems, resource exploration drilling, spacecraft attitude control systems and assessments of the impact of space weather. The coming decade will see further missions planned for more in-depth, dedicated studies of magnetic field including DEMETER, 2004; ESPERIA, 2006; Swarm, 2008; etc.

**Gravity field** observations from space can advance our knowledge of the geoid and its time variations. The geoid is the surface of equal gravitational potential at mean sea level, and reflects the irregularities in the Earth's gravity field at its surface due to the inhomogeneous mass and density distribution in the Earth's interior. Such measurements are vital for quantitative determination, in

<footnote n="1">¹Laboratoire de Géomagnétisme et Paléomagnétisme, UMR 7577, Institut de Physique du Globe de Paris, France. E-mail: chambodu@ipgp.jussieu.fr</footnote>

<footnote n="2">²Laboratoire de Gravimétrie et Géodynamique, Département de Géophysique Spatiale et Planétaire, UMR 7096, Institut de Physique du Globe de Paris, France</footnote>

<footnote n="3">³Laboratoire de Recherche en Géodésie, Institut Géographique National, France</footnote>

<footnote n="4">⁴Department of Applied and Industrial Mathematics, University of Potsdam, Germany</footnote>

<footnote>*Now at: Department of Applied and Industrial Mathematics, University of Potsdam, Germany.</footnote>

<footnote>†Now at: GeoForschungsZentrum of Potsdam, Section 2.3: Geomagnetism, Germany.</footnote>
combination with satellite altimetry, of permanent ocean currents, for improvement of global height references, for the study of the Earth's internal structure, for estimates of the thickness of the polar ice sheets and its variations and for estimates of the mass/volume redistribution of freshwater in order to further understand the hydrological cycle.

Gravity field measurements often utilize combinations of different instrument types in order to derive the necessary information: single or multiple accelerometer, precise satellite orbit determination systems and satellite-to-satellite tracking systems. CHAMP (resp. GRACE) gravity data have been providing new information on the Earth's gravity field since 2000 (resp. 2002). The future gravity mission GOCE will provide data sets that are capable of deriving the geoid with 1 cm accuracy and gravity anomalies up to 1 mGal at 100 km resolution.

With the advent of space exploration and all of the subsequent technological advancements it is possible to systematically study the Earth as a whole entity. However, it is realized that to comprehend the myriad of interactions between Earth systems, we must utilize a multidisciplinary approach within which the mapping of the Earth's magnetic and gravity fields, to a high degree of accuracy and resolution, plays a crucial role. High-resolution models of the magnetic and gravity fields of the Earth help us to understand the structure and the driving forces behind plate tectonics, lithospheric motions, mantle convection and core fluid flows.

Until recently, magnetic and gravity models have been realized by applying the same technique: the spherical harmonic analysis (SHA). This method is well suited for global representations, but is very demanding for high-resolution models. This is the main reason why during the last years some other new methods have been investigated. Among them, a promising approach emerges: the wavelets technique. In the present paper the spherical wavelet models are introduced as an alternative to spherical harmonic models of the Earth's magnetic and gravity fields like IGRF, OIFM, CM4 and, respectively, EIGEN-1S, EIGEN-2, EIGEN-GRACE01S, GGM01S, UCPH2002.0.5, EGM96. Thereby the localizing properties of spherical wavelets and their approximating capacity are shown. A detailed description of the inverse problems is given. We validate the mathematical method through the simple fits of scalar functions on the sphere: this approach will be referred to as the 'scalar case' in the following, and the derived representations as 'scalar models'. Finally, we propose magnetic and gravity models: this approach will be referred to as the 'vectorial case', and the derived models as 'vectorial models'.

# 2 STATE OF THE PROBLEM: GLOBAL AND REGIONAL MODELLING

## 2.1 Magnetic field

Current geomagnetic field models include contributions from the core, crust, ionosphere and magnetosphere and are derived by a joint analysis of ground-based and satellite magnetic observations. In doing so, a separation of the magnetic field contributions at particular epoch is necessary, and their inadequate separation is one of the limiting factors for a more accurate determination of the core field, of its secular variation, and of the lithospheric field.

The standard method for modelling the three dimensional magnetic field is called SHA. In the following this method is summarized; for the technical details, see books by Jacobs (1987) and Merill et al. (1996). The internal part of the magnetic field is the negative

spatial gradient of a scalar potential $V_M(r, \theta, \phi, t)$, which satisfies Laplace's equation. Each internal field model comprises a set of spherical harmonics, each of which being a solution to Laplace's equation:

$$V_M(r, \theta, \phi, t) = R_E \sum_{\ell=1}^{N} \sum_{m=0}^{\ell} \left( \frac{R_E}{r} \right)^{\ell+1} (g_{\ell}^{m}(t) \cos m\phi + h_{\ell}^{m}(t) \sin m\phi) P_{\ell}^{m}(\cos \theta), \quad (1)$$

where $R_E = 6371.2$ km is the mean radius of the Earth, $r \ge R_E$ denotes the radial distance from the centre of the Earth, $\theta$ denotes the geocentric colatitude, $\phi$ denotes the east longitude, $t$ is the time, $P_{\ell}^{m}(\cos \theta)$ are the Schmidt semi-normalized associated Legendre functions of degree $\ell$ and order $m$ and $g_{\ell}^{m}(t)$ and $h_{\ell}^{m}(t)$ are the corresponding Gauss coefficients. The maximum spherical harmonic degree of the expansion is $N$, which leads to $N(N+2)$ real coefficients.

Since the magnetic field changes in space and time, measurements need to be continually performed and models are frequently generated to accurately represent the magnetic field as it is. The internal field at the Earth's surface contains clearly defined components from the core at least up to harmonics $N = 13$, beyond which they begin to become comparable to those from the crust (Mauersberger 1956; Lowes 1974). For higher $N$ they likely continue to decrease in power until they are submerged below the crustal component. The crustal component has a much flatter spectral fall-off, and dominates the observed field for $N \ge 16$, which corresponds to scale sizes smaller than 2500 km (Fig. 1). No method has yet been found to separate completely the two sources. The general practice has been to ignore the crustal contribution below $N \approx 13$, and core component above that number.

The described approach can also be applied for a regional magnetic field model known as the spherical cap harmonic analysis (SCHA) and firstly proposed by Haines (1985, 1990). SCHA permits the use of data from only a portion of the Earth while satisfying the constraints of potential field theory. Indeed, the region to be studied has to be defined by a spherical cap and satisfies the zero curl and divergence conditions. The method is claimed to be valid over any spherical cap at any altitude above the Earth's surface. On the basis of these assertions, SCHA has been widely used for getting regional magnetic maps (De Santis et al. 1997; Hwang & Chen 1997; Korte & Haak 2000). However, two kinds of difficulties are to be taken into account when SCHA is performed. The first one, crucial, stressed by De Santis & Falcone (1995) comes from the failure of the correctly modelling the radial dependence (i.e. that the models can not be continued and data acquired at different altitudes can not be simultaneously inverted). The second is that there is not a simple

<figure><img src="image_2.png" /><figcaption>Figure 1. The energy spectrum of the magnetic field at the Earth's surface (CO2 model, from Holme et al. (2003)).</figcaption></figure>
relation with the global spherical harmonics. A new approach for the spherical cap harmonic modelling has been proposed by Thébault et al. (2004) in order to solve these two difficulties.

## 2.2 Gravity field

The gravity field of the Earth reflects the internal structure of the solid Earth as well as the distribution of masses in the surrounding fluid envelops (oceans, atmosphere, ice caps, hydrology). Models of the static and time-varying gravity field lead to a better understanding of the internal geodynamical processes and of the superficial and external envelops. However, time-varying gravity effects, mostly due to the contribution of the fluid envelops but also to the solid Earth processes such as post-glacial rebound, are three to four order of magnitude smaller than static contributions.

For the past three decades, global models of the Earth's static gravity field were derived from the combination of high-altitude satellite data (including altimetry) with ground-based measurements. With the advent of low-altitude gravity missions, new global models are currently released, dramatically improving our knowledge of the static field and giving an insight into its temporal variations at large scale for the first time.

Those models are classically expressed as a series of fully normalized spherical harmonics :

$$V_G(r, \theta, \phi, t) = \frac{GM}{R_{E_{eq}}} \sum_{\ell=0}^{N} \sum_{m=0}^{\ell} \left( \frac{R_{E_{eq}}}{r} \right)^{\ell+1} (C_{\ell}^{m}(t) \cos m\phi + S_{\ell}^{m}(t) \sin m\phi) R_{\ell}^{m}(\cos \theta), \quad (2)$$

where $R_{E_{eq}} = 6378.1$ km is the mean equatorial radius of the Earth, $G = 6.67 \cdot 10^{-11}$ m³ kg⁻¹ s⁻² is the Newtonian gravitational constant, M the Earth's total mass including the atmosphere, $C_{\ell}^{m}(t)$ and $S_{\ell}^{m}(t)$ dimensionless coefficients and $R_{\ell}^{m}(\cos \theta)$ are the fully normalized associated Legendre functions of degree $\ell$ and order $m$. The spectrum of the gravity field decreases as shown on Fig. 2, reflecting the continuous distribution of masses inside the Earth and its fluid envelops. Major density discontinuities of the planet's interior can be recovered from a fine analysis of the spectrum (Hipkin 2001).

For regional gravity field modelling, the SCHA can be applied although references are mostly found in the geomagnetic literature. In the geodetic literature, Hwang & Chen (1997) introduced fully normalized SCHA to analyse sea-level data. Applications are also found for regional gravity field representation (De Santis & Torta 1997), for example, over China (Li et al. 1995). Another approach was developed to deal with the problem of off-polar orbits in satellite geodesy, leading to polar gaps in the data, and to study bounded

<figure><img src="image_3.png" /><figcaption>Figure 2. The energy spectrum of the gravity field, up to $N = 360$ (EGM96 model, from Lemoine et al. (1998)).</figcaption></figure>

domains such as the oceans. As known, spherical harmonics are global functions orthogonal over the sphere. Over bounded domain, they are no longer orthogonal. This is why a new set of functions is generated from the SHA, in order to be orthogonal over the limited domain. For example, Albertella et al. (1999) define a basis on a spherical belt to avoid the problem of the polar gaps in the data. Their approach can also be applied on a spherical cap. To represent an oceanic signal, spherical harmonics are ortho-normalized over the oceans by a Gram-Schmidt procedure (Hwang 1993).

## 2.3 Limitations of the spherical harmonic analysis

As it was shown before, at the global and regional scales, spherical harmonics are among the standard mathematical procedures for describing scalar and vector fields. However, with the advent of the new satellite missions, the geo-scientific community has a new challenge in defining better techniques to describe the very high-accuracy and high-resolution data sets provided by these satellite missions.

Spherical harmonics are well suited for regular distribution of data on the whole Earth. They form an orthonormal basis. This leads to the most compact representations at global scale. Furthermore, the spherical harmonics represent a complete set of eigenfunctions for a large set of observable functionals (Rummel & van Gelderen 1995). However, they have some drawbacks as soon as irregular or local distribution of data are considered. Thus, research on various aspects in new mathematical tools is ongoing on. In order to better define which methods have to be developed, some arguments on the advantages and disadvantages of spherical harmonic models are given below:

(i) global support is required for each harmonic term;

(ii) total number of terms may not be commensurate (either too few or too many) in some areas with that required to achieve requisite accuracy;

(iii) computation of error estimates implied by very high-degree expansions, if pursued via complete covariance matrix propagation, is extremely demanding computationally;

(iv) computationally cumbersome propagation of error statistics;

(v) models yield uniform global resolution: thus, very high-degree models generally do not reflect available data resolution everywhere, and the global and regional information are not found in the same set of coefficients.

For numerous applications in geomagnetism and geodesy/ gravimetry, one common strategy is to have a global spherical harmonic expansion to the highest resolution possible for the global data and then switch to a spatial representation for any further regional resolution. One crucial requirement is to ensure that no information is lost when refining a global spherical harmonic representation to a regional one, and the solution can be provided by the wavelet analysis.

# 3 WAVELETS FAMILIES

Wavelets were firstly introduced by Morlet, working on seismic data analysis (1985). Since this time, they are more and more widely used and have been spreading among many communities (signal processing in medicine, geophysics, finance; image or sound processing and compression; etc).

For geophysical purposes, the need has been expressed to design wavelets suitable for representing the potential fields. In particular
on the sphere, the used functions should satisfy the following properties:

(i) functions have to admit a physical interpretation;

(ii) harmonic prolongation must be easily computable;

(iii) function itself has to be numerically easy to compute;

(iv) functions must be localized on the sphere.

These requirements led to constructions based on the Poisson kernel of spherical functions. Let us notice that Poisson wavelets are also used in one and two dimensions (Sailhac et al. 2000; Martelet et al. 2001; Sailhac & Gibert 2003) to analyse potential fields. Spherical wavelet constructions of that kind are well known by now (Schroeder & Sweldens 1995; Freeden & Winterheuser 1996; Holschneider 1996; Freeden et al. 1998; Freeden & Schneider 1998; Dahlke et al. 2001).

## 3.1 Frame of wavelets

Special collection of functions called frames are of primary interest for representing potential fields on the sphere. The concept of frame is more general than the basis one. In fact, it is a complete set of functions but may include some redundancy, which makes frames much more flexible than bases. More precisely, a collection {$g_n$}$_{n=0,1,...}$ in a Hilbert space H is a frame if for all $s \in H$, with two constants called frame-bounds ($0 < A \le B < \infty$), the following inequalities exist:

$$A \|s\|^2 \le \sum_n |g_n \cdot s|^2 \le B \|s\|^2, \quad (3)$$

where the expression: $g_n \cdot s$ denotes the scalar product of $g_n$ with $s$. It is possible to build discrete frames based on wavelets as the constitutive functions, by properly sampling the continuous wavelet transform in space and frequency (Holschneider 1995; Freeden & Winterheuser 1996).

In this study, we use the Poisson wavelets. More precisely, a Poisson wavelet at point $\vec{x}$ on the sphere of radius $\|\vec{x}\| = R$ is expressed as (Holschneider et al. 2003):

$$\psi_a^n(\vec{x}) = N_a^n \sum_l (a\ell)^n e^{-a\ell} Q_l \left( \frac{\vec{e}}{\|\vec{e}\|} \cdot \frac{\vec{x}}{\|\vec{x}\|} \right), \quad (4)$$

$$\text{with:} \quad \left\{ \begin{array}{l} N_a^n = \left( \int_{S(r)} \psi_a^n(\vec{x}), \psi_a^n(\vec{x}) ds(\vec{x}) \right)^{-\frac{1}{2}} \\ Q_l \left( \frac{\vec{e}}{\|\vec{e}\|} \cdot \frac{\vec{x}}{\|\vec{x}\|} \right) = (2\ell + 1) P_l \left( \frac{\vec{e}}{\|\vec{e}\|} \cdot \frac{\vec{x}}{\|\vec{x}\|} \right). \end{array} \right. \quad (5)$$

The wavelet is defined by three parameters: its scale a, its position $\vec{e}$ and its order n. In other words, a is a measure of the wavelength of the wavelet without any dimension, $\vec{e}$ is the position of the wavelet and order n characterizes the global shape of the wavelets (number of oscillations). These wavelets on the sphere may also be considered as the restriction on the sphere of a function that is harmonic outside and that has singularities inside it. In the following, $\psi_a^n(\vec{x})$ will represent either the wavelet on the sphere, or its harmonic continuation.

In the above equations, $N_a^n$ is a $L^2$-normalization factor, $\ell$ an integer which corresponds to the degree of the Legendre polynomials $P_\ell$ and $Q_\ell$ are the related kernel functions.

In fact Poisson wavelets show many useful properties. They may be identified with multipoles of order n, located at points ($R \cdot e^{-a\vec{e}}$) inside the sphere: the depth of the multipole is thus linked to the scale parameter. Let us recall the main geometrical, mathematical and numerical properties of these functions.

**Geometrical.** The wavelets are zonal functions, both localized in space and frequency. As their scales decrease, they become more sensitive to the local features of the signal and their spectra are shifted towards the higher degrees of the spherical harmonics.

**Mathematical.** The wavelets are generated by the Legendre polynomials as the spherical harmonics are, and when considering functions in the whole space, they satisfy the Laplace equation except for a pole inside the sphere. The use of the same generating functions constitutes a 'bridge' between the two representations and allows easier comparisons between the two developments. The harmonic continuation of the wavelets in the space outside sources can be defined in a similar way as for the spherical harmonics.

**Numerical.** The wavelets can be calculated easily thanks to analytical solutions. Indeed, they can be computed using a finite superposition of multipoles at suitable positions related to their scales instead of summing an infinite series of spherical harmonics. Moreover, the scalar product between two wavelets can be expressed as another wavelet.

## 3.2 Discretization of the positions

In order to build the appropriate wavelet frames, we have to discretize the continuous wavelet transform in space and define scales interval in such a way that all frequencies are covered over the whole sphere. We choose to locate wavelets on the vertexes of a spherical mesh derived by subdividing the facets of a regular convex polyhedron centred with respect to the sphere and projecting the directions of the vertexes onto the sphere. The more the polyhedron's facets are subdivided, the finer the mesh: it is then possible to create a set of hierarchical meshes associated with wavelets at different scales.

Let us note that any of the polyhedrons (Cromwell 1997) might be used. Here it is chosen to discretize the positions of wavelets on the sphere within a cube and an icosahedron. An approach based on successive generations is used, the generation j corresponding to a given mesh. It means that a level of subdivision corresponds to a given wavelet scale.

### 3.2.1 The cube

The cube is initially chosen due to its easier implementation. A description of the method is given in (Holschneider et al. 2003).

Positions are defined by recursive subdivision of each facet of the cube into four squares. Denoting j the generation, the number of vertexes of the j-level mesh $V(j)$ is defined as:

$$V(j) = \begin{cases} 1 & \text{for: } j = 0, \\ 6 \cdot (4^{j-1}) & \text{for: } j > 0. \end{cases}$$

Fig. 3 (top) shows the implementation of vertexes from $j = 1$ to $j = 4$.

The coordinates of positions are determined in two steps. At first, the points on the cube are defined within a 3-D-space Cartesian system. The middle of each facet is taken as a wavelet position. Then, they are projected onto the sphere according to a straight line, which passes from the centre of the sphere, the point at the surface of the cube and up to the sphere. With such a method, the points obtained are on the same sphere and are not perfectly homogeneously distributed (see Subsection 3.2.3).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Facets of polyhedrons at each generation *j* for a cube (top) and an icosahedron (bottom). Black crosses represent positions of wavelets.</figcaption></figure>

### 3.2.2 The icosahedron

The second chosen polyhedron is the icosahedron. The description of its implementation is given in (Kenner 1976).

Positions are defined by recursive subdivision of each facet of the icosahedron into four triangles connecting the middle of the sides. As for the cube, the number of vertexes of the *j*-level mesh is denoted $V(j)$:

$$V(j) = \begin{cases} 1 & \text{for: } j = 0, \\ 10 \cdot (4^{j-1}) + 2 & \text{for: } j > 0. \end{cases}$$

Fig. 3 (bottom) shows the incrementation of vertexes from $j = 1$ to $j = 4$.

The coordinates of positions are determined in the same two steps as for the cube. Nevertheless, the method is slightly different in the sense that the vertexes defining positions are not the middle of facets as for the cube. More interesting, the points are more regularly distributed than for the cube.

### 3.2.3 Comparison of the two meshes

The number of points at each generation is depicted in Fig. 4. The icosahedric meshes comprise more vertexes than the cubical ones: the ratio tends to about 1.7 as *j* increases. Moreover, at level (*j* + 1), it includes all the vertexes of level *j*. On the contrary, vertexes from the cubical meshes never coincide between

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Meshes on the sphere at two generations *j*, calculated for a cube (top) and an icosahedron (bottom).</figcaption></figure>

*j* and (*j* + 1) generations. Both kinds of meshes (Fig. 5) show a good regularity even if a non-negligible dispersion of distances between points is observed. The dispersion of distances between points at different generation for cubical and icosahedrical meshes is roughly coming to, respectively, 30 and 10 per cent of the mean value. This phenomenon is not significative for the resolution obtained in this study. In this paper, both kind of meshes are used in order to provide two different examples of wavelet frame discretization.

The regularity of the meshes is a difficult problem, known as ‘*Le problème des dictateurs ennemis*’: how to distribute territories of sphere to several dictators so that they have all the same territory and that they are as distant as possible from each other. The interested reader may find more details in an abundant bibliography on the subject (Hicks & Wheeling 1959; Muller 1959).

## 3.3 Discretization of the scales

In the previous subsection, we described the discretization of wavelet positions ($\theta$, $\phi$) on the sphere. We have now to associate to each generation *j* of positions the corresponding scale $a_j$. This last parameter has to be carefully chosen in order to satisfy two main constraints: (i) the spectrum should be covered and (ii) the number of wavelets of each scale should be sufficient, but not too large, to generate the corresponding spherical harmonics.

The sequence of the scales is defined on a unit sphere $\Sigma_1$ and corresponds to a geometric progression:

$$a_j = a_0 \cdot \gamma^{j-1}, \quad (6)$$

where *j* is the generation, $a_0$ is a chosen initial scale and $\gamma$, a constant verifying: $0 < \gamma < 1$.

Then, the position ($r_j$, $\theta$, $\phi$)$_{\Sigma_1}$ of the corresponding multipole inside $\Sigma_1$ is:

$$\left\{ \begin{array}{l} r_j = e^{-a_j} \\ (\theta, \phi) \text{ given by the mesh} \end{array} \right.$$

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Number of points of the polyhedrons at each generation *j*.</figcaption></figure>
The bounds of scales and of positions for ‘sources’ within $\Sigma_1$ are:

$$ \begin{cases} 0 < a_j \le a_0 \\ 1 > r_j \ge e^{-a_0} \end{cases} $$

Considering the Earth’s surface noted $\Sigma_E$, we define new positions $R_j$ and their associated scales $A_j$:

$$ A_j = a_j - \ln\left(\frac{R_{\text{ref}}}{R_E}\right), \qquad (7) $$

$$ R_j = R_{\text{ref}} \cdot r_j. \qquad (8) $$

$R_E$ is the mean radius of the Earth. $R_{\text{ref}}$ allows to introduce a *priori* information in the incrementation of scale, as it corresponds to radii of known discontinuities of the Earth. Indeed, all multipoles located inside (resp. outside) a sphere of radius $R_{\text{ref}}$ have a scale parameter greater (resp. lower) than:

$$ A_{\text{ref}} = -\ln\left(\frac{R_{\text{ref}}}{R_E}\right). \qquad (9) $$

Then, we can choose the scales of wavelet frames in order to sample the desired layers of the Earth’s interiors, taking into account the *a priori* structures.

The mathematical relation between $R_j$ and $A_j$ is the same as for the unit sphere:

$$ R_j = R_E \cdot e^{-A_j}. \qquad (10) $$

In the following, we discuss two possible examples of frames to represent magnetic and gravity data.

### 3.3.1 Magnetic field

For the magnetic field, we choose to implement a frame based on the multipoles of order 2. This order allows a good coverage of low degrees of the spherical harmonics. Moreover, it yields a precise localization in both space and frequency. The scale $a_j$ associated to the j-level on the unit sphere verifies $a_0 = 2$ and $\gamma = (\frac{1}{2})$. Positions $(\theta, \phi)$ are discretized on the cubical mesh.

Different frames are then used for the scalar and vectorial cases. Indeed, the demanding constraints and the physical meanings of this last case are more important than the simple fit of a scalar function on the sphere.

For the scalar case, the wavelets regularly sample the wavelengths present in the modelled scalar function, which is the intensity of the magnetic field $|\vec{B}|$. We thus simply use: $R_{\text{ref}} = R_E$. Table 1 shows the multipole characteristics for different generations. Fig. 6 *top* shows that the spectrum is covered. In the Fig. 6 *bottom*, we compare the number of wavelets with the number of spherical harmonics. Each

<table><caption>**Table 1.** Sequence of scales chosen for the frame used in scalar case for magnetic field modelling.</caption><thead><tr><th>j</th><th>a<sub>j</sub> = A<sub>j</sub></th><th>r<sub>j</sub></th><th>R<sub>j</sub> = r<sub>j</sub>R<sub>E</sub> (km)</th><th></th></tr></thead><tbody><tr><td>1</td><td>2</td><td>0.135</td><td>860</td><td rowspan="5">Large scale</td></tr><tr><td>2</td><td>1</td><td>0.368</td><td>2345</td></tr><tr><td>3</td><td>0.5</td><td>0.607</td><td>3867</td></tr><tr><td>4</td><td>0.25</td><td>0.779</td><td>4963</td></tr><tr><td>5</td><td>0.125</td><td>0.883</td><td>5626</td></tr><tr><td>6</td><td>0.0625</td><td>0.939</td><td>5985</td><td rowspan="2">Small scale</td></tr><tr><td></td><td>0.0</td><td>1.0</td><td>6371.2</td></tr></tbody><tfoot><tr><td colspan="4"></td><td>Earth's surface</td></tr></tfoot></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Discretization of the scales with order 2 multipoles. The energy spectra for the generations 1 to 6 are computed (top panel). The number of wavelets (*areas defined by thin segments*) is compared to the number of spherical harmonics (*areas below thick line*) for each degree $\ell$ (bottom panel).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Discretization of the scales with order 2 multipoles. The energy spectra for the generations 1 to 5 on the core (*top panel, solid curves*) and 1 to 6 on the crust (*top panel, dashed curves*) are calculated (top). The squares define the limit spectrum of the core-wavelets. The number of wavelets from the core (*areas defined by thick segments*) and from the crust (*areas defined by thin segments*) is compared to the number of spherical harmonics (*areas below thick line*) for each degree $\ell$ (bottom panel). Note that (i) the spectra of core-wavelets cover spherical harmonic up to degree $\ell \simeq 13$; (ii) the core-wavelets are more redundant than the crustal ones below $\ell = 13$ and (iii) only crustal-wavelets represent the field when $\ell \ge 16$.</figcaption></figure>

area defined by thin segments represents the number of wavelets of a given scale. The limits of the areas depend on the spectral coverage of the wavelets. The wavelet spectrum has an infinite support. However, we only consider the part retaining most of the wavelet energy. The number of wavelets is large enough, comparing to the number of spherical harmonics.

For the vectorial case, the positions of multipoles correspond to the core of the Earth: $R_{\text{ref}} = R_{\text{CMB}} = 3485$ km, and to the crust:
<table><caption>**Table 2.** Sequence of scales chosen for the frame used in vectorial case magnetic field modelling.</caption><thead><tr><th>j</th><th>a<sub>j</sub></th><th>a<sub>j</sub></th><th>r<sub>j</sub></th><th>r<sub>j</sub> (km)</th><th></th></tr></thead><tbody><tr><td>1</td><td>2</td><td>2.6</td><td>0.135</td><td>472</td><td>Core: large-scale</td></tr><tr><td>2</td><td>1</td><td>1.6</td><td>0.368</td><td>1282</td><td></td></tr><tr><td>3</td><td>0.5</td><td>1.1</td><td>0.607</td><td>2114</td><td></td></tr><tr><td>4</td><td>0.25</td><td>0.85</td><td>0.779</td><td>2714</td><td>Core: small scale</td></tr><tr><td></td><td>0.0</td><td>0.60</td><td>1.0</td><td>3485</td><td>CMB</td></tr><tr><td>1</td><td>2</td><td>2</td><td>0.135</td><td>860</td><td>Crust: large scale</td></tr><tr><td>2</td><td>1</td><td>1</td><td>0.368</td><td>2345</td><td></td></tr><tr><td>3</td><td>0.5</td><td>0.5</td><td>0.607</td><td>3867</td><td></td></tr><tr><td>4</td><td>0.25</td><td>0.25</td><td>0.779</td><td>4963</td><td></td></tr><tr><td>5</td><td>0.125</td><td>0.125</td><td>0.883</td><td>5626</td><td></td></tr><tr><td>6</td><td>0.0625</td><td>0.0625</td><td>0.939</td><td>5985</td><td>Crust: small scale</td></tr><tr><td></td><td>0.0</td><td>0.0</td><td>1.0</td><td>6371.2</td><td>Earth's surface</td></tr></tbody></table>

$R_{\text{ref}} = R_E = 6371.2 \text{ km}$. Table 2 shows the multipole characteristics for different generations.

Fig. 7 shows that the spectrum is covered, and that the number of wavelets is large enough compared to the number of spherical harmonics. The core and crustal wavelets cover all the spherical harmonic degrees. Indeed, it is necessary even for low degrees to consider both fields as the synthetic data in spherical harmonics do not distinguish the two contributions. In an ideal case, with an infinite number of wavelets, it would be possible to consider a wavelet model with positions of multipoles that would purely correspond to physical sources (for example: $R_{\text{crust}} \ge 6341.2 \text{ km} \rightarrow A_{\text{crust}} \le 0.005$). The radial positions of the wavelets inside the Earth, noted $R_j$ in Table 2, are given as an equivalent representation by using the non-unicity of the solution given by wavelet frames. Thus, the $R_j$ constitute the positions of equivalent sources.

### 3.3.2 Gravity field

For the gravity field, we choose to implement a frame based on order 3 multipoles. The localization in space and frequency is still satisfactory: indeed, the wavelets only show one spatial undulation and their spectra narrow. Let us notice that the spectra of these wavelets shift towards the higher degrees and show less power in the lower degrees. Choosing a frame based on higher order multipoles would degrade the spatial localization of the wavelets.

The scale associated to the j-level verifies $a_0 = 3$ and $\gamma = \frac{1}{2}$. The sequence of multipole depths $R_j$ inside the Earth, regularly samples its successive concentric envelops, to reflect the distribution of masses: $R_{\text{ref}} = R_{E_{\text{eq}}} = 6378.1 \text{ km}$. Table 3 shows the multipole

<table><caption>**Table 3.** Sequence of scales chosen for the frame used in gravity modelling.</caption><thead><tr><th>j</th><th>a<sub>j</sub> = A<sub>j</sub></th><th>R<sub>j</sub> = r<sub>j</sub>R<sub>E</sub> (km)</th><th>Location</th></tr></thead><tbody><tr><td>1</td><td>3</td><td>318</td><td>Inner core</td></tr><tr><td>2</td><td>1.5</td><td>1423</td><td>Outer core</td></tr><tr><td>3</td><td>0.75</td><td>3013</td><td>Outer core</td></tr><tr><td>4</td><td>0.375</td><td>4384</td><td>Lower mantle</td></tr><tr><td>5</td><td>0.1875</td><td>5288</td><td>Lower mantle</td></tr><tr><td>6</td><td>0.09375</td><td>5807</td><td>Upper mantle</td></tr><tr><td>7</td><td>0.046875</td><td>6086</td><td>Upper mantle</td></tr><tr><td>8</td><td>0.023438</td><td>6230</td><td>Upper mantle</td></tr><tr><td>9</td><td>0.011719</td><td>6304</td><td>Upper mantle (lithosphere)</td></tr><tr><td>10</td><td>0.005859</td><td>6341</td><td>Crust</td></tr><tr><td></td><td>0.0</td><td>6378.1</td><td>Earth's surface</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Discretization of the scales with order 3 multipoles. The energy spectra for the generations 1 to 7 are calculated (top panel). The number of wavelets (areas defined by thin segments) is compared to the number of spherical harmonics (areas below thick line) for each degree $\ell$ (bottom panel).</figcaption></figure>

depths for different generations. Those multipoles are considered as equivalent sources when modelling the disturbing potential (vectorial case, see Section 4).

Positions $(\theta, \phi)$ are discretized on the icosahedric meshes since the gravity anomalies are modelled at a rather high resolution.

Fig. 8 shows that the spectrum is homogeneously covered, and that the number of wavelets is large enough as compared to the number of spherical harmonics. Let us notice that this frame is more redundant than the other one on the cube.

# 4 INVERSE PROBLEM

## 4.1 The least-squares method

The magnetic and gravity fields (resp. $\vec{B}$ and $\vec{g}$) can be expressed as a linear combination of wavelets. In the following, we are focusing on the two quantities used by magnetic and gravity communities: the magnetic field and the free air gravity anomaly.<fnref n="1" /> In the scalar case, the intensity of magnetic field and gravity anomaly are directly written as a sum of wavelets. In the vectorial case, only potentials are expressed as a sum of wavelets modelled as superposition of multipolar potentials. Thereafter, the magnetic field and the gravity anomaly are derived.

Denoting $E$ the function to be represented, $\alpha$ the vector of wavelet coefficients and $\psi$ the wavelet frame, the following general equality holds :

$$E = \sum_{i} \alpha_{i} \psi_{i} \qquad (11)$$

<footnote n="1">¹The free air gravity anomaly is defined as the difference between the intensity of the real gravity field at the point of measurement and the intensity of the normal gravity field at the same point: the resulting value is called the ‘gravity disturbance’ among the geodetic community, but here it will be called ‘gravity anomaly’. It can be linked to the disturbing potential $T$, which is defined as the difference between the real gravity potential of the Earth and the normal potential of the reference ellipsoid (Moritz 1989; Hackney & Featherstone 2003).</footnote>
In order to find the coefficients $\alpha_j$ of the wavelet development, we set a classical least-squares problem. This method consists in deriving the set of coefficients, which minimizes the residuals between the data and the model in a $L^2$ sense. However, many sets of coefficients lead to a good fit of the data since the problem is often underdetermined. Thus, we have to take into account additional constraints in order to eliminate the overfitted solutions that show large oscillations. Here, we introduced a smoothness constraint to regularize the problem. However, the smoother the solution is, the worse the measurement residuals are: we have to find a trade-off between a good fit of data and the global smoothness. We actually minimize the following quantity:

$$ (b - M\alpha)^t \cdot W \cdot (b - M\alpha) + \lambda\alpha^t \cdot L \cdot \alpha \quad (12) $$

leading to the normal system:

$$ (M^t \cdot W \cdot M + \lambda L) \cdot \alpha = M^t \cdot W \cdot b, \quad (13) $$

where $b$ is the vector of measurements. $M$ is a $(i \times j)$ matrix. The $j$th column of the matrix $M$ contains the $j$th wavelet sampled at the $i$th observation points. $\alpha$ is the vector of wavelet coefficients; $W$ is a matrix of data weighting; $L$ is the matrix of a quadratic form that controls the regularity, and $\lambda$ is a parameter balancing between fit and smoothness. The parameter $\lambda$ has to be chosen in such a way that, on the one hand, it avoids the model to fit the data with a precision better than their noise (case of overfitting, for which the solution is too oscillating). On the other hand, it avoids the case of underfitting, for which the solution is too smooth.

Eq. (13) is actually related to the general theory of inverse problems (see for instance Tarantola 1987). In this theory, on the one hand, we suppose that the observation errors have a Gaussian distribution with covariance matrix $W^{-1}$. On the other hand, the *a priori* probability of a model, of a set of coefficients $\alpha$, is again a Gaussian distribution with density $\exp(-\lambda (\alpha^t \cdot L \cdot \alpha))$ where $L$ is the matrix of a bilinear quadratic form. Typically, $L$ describes how the field energy decreases from the large scales to the small ones: coefficients at large scales show indeed larger variances than those at small scales.

Measurements bring additional information, so that an *a posteriori* probability on the coefficients can be computed, taking the measurements as well as the *a priori* knowledge into account. The vector of coefficients for which this *a posteriori* probability reaches a maximum is given by eq. (13).

In the following, we assume an uncorrelated noise on the data, so $W$ is diagonal with:

$$ W_{ij} = \frac{1}{\sigma_j^2} \cdot \delta_{ij}, \quad (14) $$

where $\sigma_j^2$ is the variance of the noise of $j$th measurement. When measurements are very precise, the diagonal terms of $W$ become very large, implying a stronger constraint on measurement residuals. We provide two examples of regularization, keeping in mind the notion of spectral decrease.

We now detail eq. (13) for both scalar and vectorial cases.

## 4.2 Scalar case

In the scalar case, the function $E$ (eq. 11) represents the intensity of the magnetic field $|\vec{B}|$ or the gravity anomaly $\Delta g$. Thus, $b$ is directly the vector of measurements of $|\vec{B}|$ or $\Delta g$.

Studying the scalar case is mainly a way to check the applicability of the method rather than to get a representation of potential fields. It allows us to apprehend the behaviour of the wavelets and the

influence of each parameter of the model (incrementations of scales, discretizations of positions, choices of the regularization parameter $\lambda$). The scalar case also allows to test the capability of the wavelets to represent a given function $E$ on the sphere.

An important task is how to define the $L$ matrix. In the following, we present two examples in how it can be chosen. The first example is mainly applied to the magnetic field and the second to the gravity field. However, it is possible to exchange the presented approaches between the two fields.

### 4.2.1 Magnetic field

In this example, we directly parametrize the regularization matrix $L$. It expresses the decreasing power spectrum of the magnetic potential (Lowes 1974). $L$ contains the scalar product of the wavelets on the sphere as:

$$ L_{i,i'} = \langle K \psi_i, \psi_{i'} \rangle_{R_E} \quad (15) $$

where $R_E$ is the mean radius of the Earth. The notation of the scalar product should be interpreted as an integral:

$$ \langle s, u \rangle_R = \int_{\Sigma_R} s \cdot u \, du, $$

where $\Sigma_R$ is the sphere of radius $R$. $K$ is an operator of weighting the spectrum, defined by its action on spherical harmonics:

$$ K : Y_{\ell,m} \rightarrow \ell^t Y_{\ell,m} \quad (16) $$

with $t$ a positive exponent. Thanks to the wavelets chosen this quadratic form is explicit and leads to an analytical solution, which corresponds to another wavelet of higher order and scale (see Holschneider et al. 2003).

Note that $K$ may be written as:

$$ (Ks)(\vec{r}) = \langle \Pi(r, \cdot), s \rangle, \quad (17) $$

$$ \text{where } \Pi(\vec{r}, \vec{r}_1) = \sum_{\ell} \ell^t (2\ell + 1) P_{\ell}(\vec{r}, \vec{r}_1). \quad (18) $$

### 4.2.2 Gravity field

In this example, we take matrix $L^{-1}$ as the *a priori* covariance matrix between coefficients and then invert it to get $L$.

The covariances of the coefficients are thus linked to the spatial correlations of the gravity anomalies. We denote $C(\vec{r}, \vec{r}_1)$ the covariance function of gravity anomalies at points $\vec{r}$ and $\vec{r}_1$ on the sphere. We make the assumption that $C(\vec{r}, \vec{r}_1)$ only depends on the spherical distance between $\vec{r}$ and $\vec{r}_1$. In this case, $C(\vec{r}, \vec{r}_1)$ can be written as a series of Legendre polynomials (Moritz 1989):

$$ C(\vec{r}, \vec{r}_1) = \sum_{\ell} c_{\ell} P_{\ell}(\vec{r} \cdot \vec{r}_1). \quad (19) $$

The coefficients $c_{\ell}$ are equal to the variance of gravity anomalies for degree $\ell$. We assume that the power spectrum of the gravity potential follows Kaula's rule of quadratic decrease (Kaula 1966). Thus, the power spectrum of the gravity anomaly at degree $\ell$ decreases as:

$$ c_{\ell} = (\ell + 1)^2 (2\ell + 1) \frac{\eta}{\ell^4}, \quad (20) $$

where $\eta$ is a real constant. This quadratic form may be compared to the one used for the magnetic field.
Let us now denote $K$ the operator associating to each square integrable function on the sphere $f$ its scalar product with $C(\vec{r}, \vec{r}_1)$:

$$ (Kf)(\vec{r}) = \langle C(\vec{r}, \cdot), f \rangle_{R_{E\text{eq}}} \quad (21) $$

where $R_{E\text{eq}}$ is the mean equatorial radius of the Earth. We derive the covariance between two coefficients $\alpha_i$ and $\alpha_{i'}$ as a scalar product between the corresponding wavelets. This comes from the formulae of covariance propagation by (Moritz 1989):

$$ L_{i,i'}^{-1} = \langle K \psi_i, \psi_{i'} \rangle_{R_{E\text{eq}}} \quad (22) $$

The last step consists in inverting the matrix given by eq. (22) to get $L$.

## 4.3 Vectorial case

In the vectorial case, the function $E$ represents the magnetic potential $V_M$ or the gravity disturbing potential $T$ (see footnote 1 in Section 4.1).

The measurement vector $b$ contains no values of $V_M$ or $T$, but the vectorial components of the magnetic field or the gravity anomaly in the radial spherical approximation (the gravity anomaly is a scalar but in the spherical approximation, it is oriented in the radial direction and then considered here as vectorial). Thus the matrix $M$ is different for each inverse problem. Sections 4.3.1 and 4.3.2 present the equation systems for the two fields.

### 4.3.1 Magnetic field

The magnetic field is represented via a superposition of the derivatives of the wavelets in the spherical coordinate system :

$$ \vec{B} = \begin{pmatrix} B_r \\ B_\theta \\ B_\phi \end{pmatrix} = \begin{pmatrix} -\frac{\partial V}{\partial r} \\ -\frac{1}{r} \frac{\partial V}{\partial \theta} \\ -\frac{1}{r \sin \theta} \frac{\partial V}{\partial \phi} \end{pmatrix}, \quad (23) $$

$$ \vec{B} = \sum_i \alpha_i \begin{pmatrix} -\frac{\partial \psi_i}{\partial r} \\ -\frac{1}{r} \frac{\partial \psi_i}{\partial \theta} \\ -\frac{1}{r \sin \theta} \frac{\partial \psi_i}{\partial \phi} \end{pmatrix} \quad (24) $$

In matricial notation, eq. (24) becomes:

$$ \vec{B} = M \alpha \quad (25) $$

The $L$ matrix is, as for the scalar case, chosen in order to allow a balance between a good fit and a global smoothness. We choose to implement it in order to assume a regularity on the potential through constraint on $B_r$ the radial component of the magnetic field. Thus:

$$ L_{i,i'} = \left\langle \frac{\partial}{\partial r} \psi_i, \frac{\partial}{\partial r} \psi_{i'} \right\rangle_{R_E \text{ or } R_{CMB}} \quad (26) $$

The value of the radial derivative increases when the scale decreases. Thus, small-scale wavelets are more expensive than large-scale ones. This matrix is applied on both groups of wavelets supposed to represent the core and the lithosphere.

### 4.3.2 Gravity field

In the spherical approximation, the free air gravity anomaly is related to the disturbing potential Moritz (1989):

$$ \Delta g = -\frac{\partial T}{\partial r}, \quad (27) $$

where $r$ is the spherical radius. The wavelet expansion of the gravity anomaly is then derived by linearity:

$$ \Delta g = - \sum_i \alpha_i \left( \frac{\partial \psi_i}{\partial r} \right). \quad (28) $$

In matricial notation, eq. (28) becomes:

$$ \Delta g = M \alpha \quad (29) $$

As for the scalar case, the matrix $L$ is defined as the inverse of the covariance matrix of the wavelet coefficients, but these coefficients correspond now to the wavelet approximation of the disturbing potential. This matrix is derived in the same way as described in Section 4.2.2. However, the power spectrum of the gravity anomalies has to be replaced by the power spectrum of the disturbing potential. The Kaula's rule for the potential reads :

$$ c_{\ell} = (2\ell + 1) \frac{\beta}{\ell^4}, \quad (30) $$

where $\beta$ is a real constant.

# 5 RESULTS

In the present paper, we only discuss tests obtained from synthetic data sets. Indeed, this choice allows: (1) to know exactly the spatial and spectral contents of the used information, (2) to test as many distributions of data as possible and (3) to directly compare with the initial model in spherical harmonics. The two fields, magnetic and gravity, differ in their spatial and spectral characteristics. Their behaviours help us to study the impact of different parameters of the used wavelet frames.

The results we obtained are presented for magnetic and gravity fields, separately. Moreover, we made a distinction between the global and regional representations, obtained with both possible regular and irregular distributions.

## 5.1 Magnetic models

### 5.1.1 Global representations

**Data.** We used a synthetic data set computed from the CO2 model (Holme et al. 2003). This model is obtained from the measurements

<figure><img src="image_10.png" /><figcaption>Figure 9. Map for the magnetic field intensity ($|\vec{B}|$) computed at the Earth's surface, from CO2 model (Holme et al. 2003) up to degree/order 13. Black dots represent the chosen irregular distribution of data, which directly correspond to the distribution of 670, past and present, magnetic observatory locations.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Same as Fig. 9 for Northern X (top), Eastern Y (middle) and Downward vertical Z (bottom) components.</figcaption></figure>

<table><caption>**Table 4.** Parameters of the global tests—magnetic scalar case.</caption><thead><tr><th>Parameters</th><th>Regular case<br>(629 data)</th><th>Irregular case<br>(670 data)</th></tr></thead><tbody><tr><td>Order of multipoles</td><td>2</td><td>2</td></tr><tr><td>Generations of the frame</td><td>1 to 4</td><td>1 to 4</td></tr><tr><td>Scales</td><td>See Table 1</td><td>See Table 1</td></tr><tr><td>Number of wavelets</td><td>510</td><td>510</td></tr><tr><td>W matrix</td><td>1||</td><td>1||</td></tr><tr><td>λ parameter</td><td>10<sup>-10</sup></td><td>10<sup>-9</sup></td></tr><tr><td>t exponent parameter</td><td>7</td><td>7</td></tr></tbody></table>

provided by three magnetic satellites: CHAMP, Ørsted and SAC-C, from 2000 August to 2001 December. This spherical harmonic model is developed up to degree/order 29 when describing the internal field. The magnetic field spectrum is clearly decreasing up to degree/order 13 (Fig. 1). We truncated the CO2 model at this degree,

<table><caption>**Table 5.** Parameters of the global tests—magnetic vectorial case.</caption><thead><tr><th>Parameters</th><th>Regular case<br>(629 data)</th><th>Irregular case<br>(670 data)</th></tr></thead><tbody><tr><td>Wavelets belong to:</td><td>Core/crust</td><td>Core/crust</td></tr><tr><td>Order of multipoles</td><td>2/2</td><td>2/2</td></tr><tr><td>Generations of the frame</td><td>1 to 4/1 to 3</td><td>1 to 3/1 to 3</td></tr><tr><td>Scales</td><td>See Table 2</td><td>See Table 2</td></tr><tr><td>Number of wavelets</td><td>510/126</td><td>126/126</td></tr><tr><td>W matrix</td><td>1||</td><td>1||</td></tr><tr><td>λ parameter</td><td>10<sup>-21</sup></td><td>10<sup>-14</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Map of residuals on the magnetic field intensity (|B̃|) at the Earth's surface (initial spherical harmonic model CO2–510 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 629 data.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Histograms of residuals on the magnetic field intensity |B̃| (initial synthetic data from spherical harmonics model CO2—reconstructed data from 510 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 629 data.</figcaption></figure>

for global representations of main internal contributions (Figs 9 and 10).

From this model, we constituted a regular and an irregular distribution of data. We did not apply a Gaussian filter on synthetic data as it was done for gravity modelling (see below). The regular distribution of synthetic data comprises 629 samples, with one data per bin of 10° × 10° on the whole Earth's surface. The irregular distribution contains 670 samples, located at all observatory
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Maps of residuals on the magnetic field Northern *X* (top), Eastern *Y* (middle) and Downward vertical *Z* (bottom) components at the Earth's surface (initial spherical harmonic model CO2 – 636 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 629 data.</figcaption></figure>

positions that have ever run on the Earth (past and present observational locations). We chose this realistic distribution because it allows to get dense (Europe) and sparse (Pacific) covered areas.

**Tests.** Tables 4 and 5 summarize the parameters used in global magnetic field modelling. We fixed a noise level at 1 nT for all tests presented below: matrix *W* is thus set to the identity matrix in all computations. This choice is motivated by the fact that data are supposed to be ‘perfect’. The order of multipoles, the number of generations of the frame and the scale sequences are fixed with respect to some geophysical constraints (see paragraph 3.3.1). Thus the residuals between the synthetic data and their wavelet representation are minimized by adjusting the following parameters: the exponent of regularization *t* in the scalar case, the quadratic form (matrix *L*) in the vectorial case and the regularization parameter λ in both cases.

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Histograms of residuals on the magnetic field vectorial components (initial synthetic data from spherical harmonic model CO2—reconstructed data from 636 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 629 data.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Map for the magnetic field intensity ($|\vec{B}|$) computed at the Earth's surface from 510 wavelet model. Black dots represent the chosen irregular distribution of 670 synthetic data used in the scalar case.</figcaption></figure>

**Results—regular case.** The obtained wavelet models are not mapped here. Indeed, they mimic the initial synthetic data provided by CO2 model (Figs 9 and 10). For the scalar case, the residuals between the CO2 and the wavelet models (Fig. 11) are rather important, up to ±1800 nT. Nevertheless, the residuals between the initial synthetic and reconstructed data (Fig. 12) are well centred on zero and reach limits of ±400 nT. The scalar case is presented just as a simple fit of the data: the regularization is rough and it is not based on geophysical considerations.

The results obtained for the vectorial case show the real advantages when applying the wavelets. First of all, the residuals between the CO2 and the wavelet models (Fig. 13) are smaller than the considered noise of 1 nT. Furthermore, the residuals between initial synthetic and reconstructed data (Fig. 14) are well centred on zero, being no larger than ±0.25 nT. Due to the number and the repetition of data, low regularizations (see λ parameter in Tables 5) allow an easy recovery of the initial spherical harmonic model up to degree/order 13.

**Results—irregular case.** Fig. 15 shows a map of the intensity of the magnetic field modelled by wavelets. When comparing with Fig. 9, one can see additional oscillations, without major changes on the global shape of the intensity. The residuals between the CO2 and the wavelet models (Fig. 16) reach some huge values, up to –10 000 nT in the South Atlantic area or 5000 nT in India Ocean.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Map of residuals on the magnetic field intensity ($|\vec{B}|$) at the Earth's surface (initial spherical harmonic model CO2–510 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 670 data.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Histograms of residuals on the magnetic field intensity $|\vec{B}|$ (initial synthetic data from spherical harmonics model CO2—reconstructed data from 510 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 670 data.</figcaption></figure>

In this case, the residuals on the data (Fig. 17) vary up to ±300 nT. This model points out the key role of the regularization. Indeed, the small-scale wavelets are not sufficiently constrained. If the regularization parameter $\lambda$ is increasing, the resolution decreases, leading to an increase of the residuals on the data. If $\lambda$ is decreasing, the residuals also decrease up to ±5 nT but the obtained model is less regular. Fig. 18 shows maps for X, Y and Z magnetic field components modelled by wavelets. When comparing these maps with those shown in Fig. 10, no remarkable differences appear between the initial spherical harmonic and the wavelet models. For all regions well covered by data the residuals for all three components are as small as a few tens of nT. As expected, the largest residuals (Fig. 19) are in areas without data (Pacific Ocean, Southern Atlantic region). For regions without data, these residuals are larger. These differences are due to the spherical harmonic model, as wavelet model does not introduce spurious artefacts when data are missing. Let us emphasize that a recalculated spherical harmonic model from the same data set represents these data in a less accurate way than the wavelet model we presented here. The residuals on data (Fig. 20) are larger

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Maps for the magnetic field Northern X (top), Eastern Y (mid-dle) and Downward vertical Z (bottom) components computed at the Earth's surface from 252 wavelet model. Black dots represent the chosen irregular distribution of 670 synthetic data used in the vectorial case.</figcaption></figure>

than for the regular case up to ±25 nT. These last results were obtained with an increased regularization parameter. We had to choose a trade-off between a regular field at Earth's surface and a reasonable fit of data. Again, decreasing the regularization parameter $\lambda$ up to $1 \cdot 10^{-21}$ leads to a better fit of initial data, with residuals no larger than ±10 nT.

### 5.1.2 Regional representation

**Data.** For the regional representation we chose an area centred over Europe, as it contains the largest number of magnetic observations. The intensity $|\vec{B}|$ and the three components X, Y and Z of the magnetic field computed up to degree/order 29 from the CO2 model centred over the considered area are shown in Figs (21) and (22).

From the same model, we also computed an regular and an irregular distribution of data centred over Europe (15° N/75° N lat.;
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Maps of residuals on the magnetic field Northern *X* (top), Eastern *Y* (middle) and Downward vertical *Z* (bottom) components at the Earth's surface (initial spherical harmonic model CO2 – 252 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 670 data.</figcaption></figure>

15° W/45°E lon.). The regular distribution of data contains 576 values with one data per bin of 3° × 3° and a sparse grid of data elsewhere with one data per bin of 20° × 20°. The regular sampling takes into account the Nyquist frequency constraint. The irregular distribution comprises synthetic data at all European and Northern African observatory locations, and a sparse distribution with one data per bin of 20° × 20° on the remaining surface. We used this sparse, global data set because the regional magnetic field contains large wavelengths that need to be constrained. Thus, to derive the regional model, we selected small-scale wavelets in the neighbourhood of the area under study, and large-scale wavelets all over the globe. The extension of the neighbourhood depends on the shape of the wavelet: for a given scale, we selected all wavelets whose radius of influence intersects the area under study. The radius of influence of a wavelet is defined as the distance between the wavelet location

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Histograms of residuals on the magnetic field vectorial components (initial synthetic data from spherical harmonic model CO2—reconstructed data from 252 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 670 data.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Map for the magnetic field intensity ($|\vec{B}|$) computed from CO2 model up to degree/order 29, above Europe. The black dots represent the chosen irregular distribution of data, which directly correspond to the distribution of past and present magnetic observatory locations.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Maps for the magnetic field Northern *X*, Eastern *Y* and Downward vertical *Z* components computed from CO2 model up to degree/order 29, above Europe. The black dots represent the chosen irregular distribution of data, which directly correspond to the distribution of past and present magnetic observatory locations.</figcaption></figure>

<table><caption>**Table 6.** Parameters of the regional tests—magnetic scalar case (Parameters: with stars as superscripts denote wavelets localized on the whole sphere; without superscript, wavelets are localized on the European area.).</caption><thead><tr><th>Parameters</th><th>Regular case<br>(576 data)</th><th>Irregular case<br>(200 data)</th></tr></thead><tbody><tr><td>Order of multipoles</td><td>2</td><td>2</td></tr><tr><td>Generations of the frame</td><td>(1 to 3)*, 4 and 5</td><td>(1 to 3)*, 4 and 5</td></tr><tr><td>Scales</td><td>See Table 1</td><td>See Table 1</td></tr><tr><td>Number of wavelets</td><td>127* and 218</td><td>127* and 218</td></tr><tr><td>W matrix</td><td>1Ⅱ</td><td>1Ⅱ</td></tr><tr><td>λ parameter</td><td>10<sup>-15</sup></td><td>10<sup>-13</sup></td></tr><tr><td>t exponent parameter</td><td>7</td><td>7</td></tr></tbody></table>

<table><caption>**Table 7.** Parameters of the regional tests—magnetic vectorial case (Parameters: with stars as superscripts denote wavelets localized on the whole sphere; without superscript, wavelets are localized on the European area.).</caption><thead><tr><th>Parameters</th><th>Regular case<br>(576 data)</th><th>Irregular case<br>(200 data)</th></tr></thead><tbody><tr><td>wavelets belong to:</td><td>Core/crust</td><td>Core/crust</td></tr><tr><td>Order of multipoles</td><td>2/2</td><td>2/2</td></tr><tr><td>Generations of the frame</td><td>(1 to 3)* and 4/3 to 5</td><td>(1 to 3)* and 4/1 to 4</td></tr><tr><td>Scales</td><td>See Table 2</td><td>See Table 2</td></tr><tr><td>Number of wavelets</td><td>127* and 44/216</td><td>126* and 44/63</td></tr><tr><td>W matrix</td><td>1Ⅱ</td><td>1Ⅱ</td></tr><tr><td>λ parameter</td><td>10<sup>-20</sup></td><td>10<sup>-17</sup></td></tr></tbody></table>

(where it gets its maximal value) and the position where the wavelet amplitude has decreased by 90 per cent.

**Tests.** Tables 6 and 7 summarize the parameters used in regional magnetic modelling, for the scalar and vectorial cases, respectively. Compared with the global tests, we used here a larger number of generations of the frame. Again, for the vectorial case we considered

that magnetic sources are located in the core and in the crust, and they are consequently modelled.

**Results—regular case.** As for the global tests, the maps obtained for the regular case are not plotted here. Indeed they are similar to those obtained from the initial spherical harmonic model (Figs 21 and 22).

The residuals between the CO2 and the wavelet models (Figs 23 and 25) do not exceed ±80 nT for the scalar case and ±12 nT for the vectorial one. Let us remark that in the regional case the range covered by residuals is one order of magnitude larger than in the global case (Fig. 13). This difference comes from the fact that the wavelet model is not truncated at degree/order 29 as the spherical harmonic model does. The reconstructed wavelet model contains small scales that do not exist in the initial data set. The residuals on the data (Figs 24 and 26) are of order ±70 nT for the scalar case and of order ±8 nT for the vectorial case.

**Results—irregular case.** The wavelet models (Figs 27 and 30) we obtained mimic the initial synthetic data provided by CO2 model (Figs 21 and 22) for both scalar and vectorial cases.

The residuals between the CO2 and the wavelet models (Figs 28 and 31) are larger in the Southern part of the area for both scalar and vectorial maps. This is the outcome of the very sparse data distribution in Northern Africa. In the upper part of the studied area, the two observatories situated in Northern Sea allow to better constrain the wavelet model. The residuals on the data (Figs 29 and 32), as expected, are larger for the scalar case (±50 nT) than for the vectorial (±10 nT).

## 5.2 Gravity models

### 5.2.1 Global representations

**Data.** To assess the ability of the wavelets to represent the gravity field at a global scale, we inverted synthetic gravity anomalies in
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Map of residuals on the magnetic field intensity (|B|) at the Earth's surface (initial spherical harmonic model CO2-345 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 576 data.</figcaption></figure>

both cases of regular and irregular distribution of data. We constructed the synthetic data sets using the first gravity model based on GRACE data: GGM01S Tapley et al. (2004). This model was established with 111 days of GRACE measurements, from 2002 April to November, and is developed up to degree and order 120. We truncated the gravity anomaly model at degree and order 12, and applied a Gaussian filter to the spherical harmonic coefficients in order to avoid artificial oscillations Sandwell & Renkin (1988). The applied filter is given by:

$$w(\ell) = e^{-\frac{(\ell-2)^2}{2(\sigma-2)^2}}, \quad (31)$$

where $\ell$ is the degree of the spherical harmonics. We chose $\sigma = 7$. This value corresponds to a characteristic attenuation of around 0.6.

From this model, we computed a regular and an irregular distribution of data. The regular distribution comprises 1106 samples, with one data per bin of 7.5° × 7.5° on the whole Earth's surface. The 7.5° interval was chosen in order to respect the Nyquist frequency constraint for the degree 12. The irregular distribution comprises 99 samples, with large gaps at low latitudes and areas of higher concentration. We chose an arbitrary distribution of data, without any preference for the continental areas. Indeed, gravity data are avail-

<figure><img src="image_25.png" /><figcaption>**Figure 24.** Histograms of residuals on the magnetic field intensity |B| (initial synthetic data from spherical harmonic model CO2—reconstructed data from 345 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 576 data.</figcaption></figure>

able all over the Earth: ground-based measurements are numerous over continents, and satellite altimetry provides a complete coverage of the oceans.

Fig. 33 shows the gravity anomaly model GGM01S filtered as previously explained, with the irregular distribution of data superimposed.

**Tests.** Table 8 summarizes the parameters used in gravity modelling at a global scale. Parameters are the same for both scalar and vectorial cases. We did not use the first generation of the frame: indeed, gravity anomalies have no component on the degrees 0 and 1 whereas the spectrum of the first generation is centred onto degree 1.

As the data are ‘perfect’ (we did not spoil them with any synthetic noise), matrix W of data weighting is arbitrary. The uniform weighting applied corresponds to an arbitrary data noise of $10^{-2}$ mGal for the irregular distribution of data, and $10^{-5}$ mGal for the regular one.

Lastly, we filtered the wavelet model in the same way as the gravity model, to insure their consistency. Indeed, real gravity anomalies, modelled as a sum of wavelets, have an infinite spectrum. In our examples, the synthetic data only constrain the low-frequency part: thus, they only give an access to the low-frequency part of the wavelet representation. The presented results are thus derived on a filtered wavelet frame.

**Results—regular case.** Wavelets succeed in representing the low harmonics of the gravity anomaly model. We do not show the wavelet models since there is visually no difference with the GGM01S model shown on Fig. 33.

The residuals between the initial GGM01S gravity anomaly model and the wavelet reconstruction from an regular distribution of data are of same magnitude as the data noise (Figs 34 and 36). Moreover, their aspect is rather isotropic. Measurement residuals (Figs 35 and 37) reach $5 \cdot 10^{-5}$ mGal: same order as the $10^{-5}$ mGal of assumed data noise.

**Results—irregular case.** Wavelet models obtained for the scalar and vectorial cases (Figs 38 and 41) show that the wavelets well handle the gaps in the data set: the wavelet model reproduces the
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Maps of residuals on the magnetic field Northern *X* (left), Eastern *Y* (middle) and Downward vertical *Z* (right) components at the Earth's surface (initial spherical harmonic model CO2 – 387 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 576 data.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 26.** Histograms of residuals on the magnetic field vectorial components (initial synthetic data from spherical harmonic model CO2—reconstructed data from 387 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 576 data.</figcaption></figure>

main structures of the gravity anomaly model even in poorly constrained areas. Residuals between wavelet and GGM01S models increase in the equatorial area, due to the large gaps in the spatial coverage of the irregular data set (Figs 39 and 42): the wavelets nicely predict the features in the gaps, but their amplitudes and localization are not perfect because of lack of constraints. Moreover, spherical harmonics would lead to strong oscillations with such a distribution of data, a phenomenon that the wavelets avoid. Let us notice that the number of measurements is much smaller than the number of wavelets: that is the reason why residuals can be high. However, in areas of higher density of data, residuals do not exceed 1 mGal. Residuals between measurements and reconstructed data from wavelet model are mostly of comparable magnitude as assumed data noise (Figs 40 and 43).

Lastly, scalar and vectorial cases yield very similar results in both regular and irregular cases. Indeed, the wavelets can represent any spherical function. The difference between both cases is that the wavelet coefficients have a physical meaning in the vectorial case only, the scalar case only being a mathematical fit of a function on a sphere.

### 5.2.2 Regional representations

**Data.** Here, we show the ability of the wavelet frames to represent the gravity field at a higher resolution for regional studies. To get a maximum resolution, we derived the synthetic data from the EGM96 gravity model (Lemoine et al. 1998), up to degree and order 360 (this corresponds to a resolution of 55 km). As previously, we applied a Gaussian filtering to the coefficients up to degree 360 with $\sigma = 250$.

From this model, we computed a regular and an irregular distribution of data. The studied area is located in South America, in the northern part of the Andes (5°/14° lat. S; 73°/82° lon. W). The regular distribution contains 1369 samples, with one data per bin of 0.25° at the Earth's surface. Again, the bin takes into account the Nyquist frequency constraint. The irregular distribution comprises 576 samples, in a larger area than regular case (0°/9° lat. S; 67°/85° lon. W). There is a higher concentration of data in the northern half (394 data) than in the southern one (182 data). Fig. 44 shows the gravity anomaly model EGM96 filtered as previously explained, in the whole area, with the irregular distribution of data superimposed.

**Tests.** Table 9 gives the used parameters. They are the same for both scalar and vectorial cases. We selected the wavelets located inside the area under study and in an appropriate neighbourhood. We defined this neighbourhood in the following way: first, we calculated the radius of influence of the wavelets at the different scales. This radius is equal to the square root of the spherical variance of the function. Then, for a given scale, we selected all the wavelets whose
<figure><img src="image_28.png" /><figcaption>**Figure 27.** Map for the magnetic field intensity ($|\vec{B}|$) computed at the Earth's surface from 345 wavelet model. Black dots represent the chosen irregular distribution of 200 synthetic data used in the scalar case.</figcaption></figure>

radius of influence intersects the area under study. The wavelet model is filtered in the same way as the data are. We applied the arbitrary noise of $10^{-2}$ mGal for the irregular distribution of data and $10^{-3}$ mGal for the regular one. This last value is larger than for the global case: indeed, generation 11 of the frame could also contribute to gravity anomaly modelling since its power spectrum still contains energy below degree 360. Nevertheless, we only took 10 generations in the frame so that the number of wavelets should not be too large. Moreover, the spatial selection of wavelets could also degrade the approximation. For all these reasons, we increased the arbitrary data noise.

**Results—regular case.** The subset of the wavelet frame succeeds in representing the regional gravity anomaly model with a good precision, but we could not avoid edge effects. We do not show the wavelet models since there is again no visual difference with the EGM96 model.

Residuals between the initial gravity model and the wavelet reconstruction are rather homogeneous in the central area (around $10^{-3}$ mGal), and they increase to $10^{-2}$ mGal when approaching the boundaries, a value still very acceptable considering the usual precision of regional gravity maps (data precision is around 0.1 mGal). The residuals reach the extrema of -0.3 mGal and 0.25 mGal at

<figure><img src="image_29.png" /><figcaption>**Figure 28.** Map of residuals on the magnetic field intensity ($|\vec{B}|$) at the Earth's surface (initial spherical harmonic model CO2 – 345 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 200 data.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 29.** Histogram of residuals on the magnetic field intensity ($|\vec{B}|$) at the Earth's surface (initial spherical harmonic model CO2 – 345 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 200 data.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Figure 30.** Maps for the magnetic field Northern *X* (top), Eastern *Y* (middle) and Downward vertical *Z* (bottom) components computed at the Earth's surface from 233 wavelet model. Black dots represent the chosen irregular distribution of 200 synthetic data used in the vectorial case.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 31.** Maps of residuals on the magnetic field Northern *X* (top), Eastern *Y* (middle) and Downward vertical *Z* (bottom) components at the Earth's surface (initial spherical harmonic model CO2 – 233 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 200 data.</figcaption></figure>

the edges (Figs 45 and 47). These effects are mainly localized on the northwest and the southeast edges. They are due to the spatial selection of the wavelets, and to the distribution of vertices of the icosahedral mesh. Indeed, the vertices do not follow lines of equal

latitude and longitude, but their two main directions show a slight obliquity with respect to the directions of parallels and meridians. As a consequence, when selecting wavelets in a neighbourhood of the area, the coverage is not perfectly isotropic and two corners over
<figure><img src="image_33.png" /><figcaption>**Figure 32.** Histograms of residuals on the magnetic field vectorial components (initial synthetic data from spherical harmonic model CO2—reconstructed data from 233 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 200 data.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 33.** Map for the gravity anomaly $\Delta g$ computed at the Earth's surface from model GGM01S up to degree/order 12, with Gaussian filtering. Black dots represent the irregular distribution of 99 synthetic data.</figcaption></figure>

the four are privileged. Let us notice that the two other corners are modelled with an accuracy of $10^{-3}$ mGal. A full treatment of the edge effects is beyond the scope of this paper and will necessitate further studies. Residuals between synthetic data and reconstructed data from wavelet model are of comparable magnitude as assumed data noise (Figs 46 and 48).

**Results—irregular case.** Wavelet models (Figs 49 and 52) for both scalar and vectorial cases again underline the ability of wavelets to model dispersed data without generating awkward large oscillations. Similar to the global case, the wavelets reproduce the main features of the field. Residuals between the wavelet model and the EGM 96 model (Figs 50 and 53) are larger in the southern part of the area than in the northern one: this is due to the unequal

<table><caption>**Table 8.** Parameters of the global tests—gravity case.</caption><thead><tr><th>Parameters</th><th>Regular case<br>(1106 data)</th><th>Irregular case<br>(99 data)</th></tr></thead><tbody><tr><td>Order of multipoles</td><td>3</td><td>3</td></tr><tr><td>Generations of the frame</td><td>2 to 4</td><td>2 to 4</td></tr><tr><td>Scales</td><td>See Table 3</td><td>See Table 3</td></tr><tr><td>Number of wavelets</td><td>846</td><td>846</td></tr><tr><td>W matrix</td><td>10<sup>10</sup>Π</td><td>10<sup>4</sup>Π</td></tr><tr><td>λ parameter</td><td>0.005</td><td>5</td></tr></tbody></table>

<figure><img src="image_35.png" /><figcaption>**Figure 34.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model GGM01S – 846 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 1106 data.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>**Figure 35.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model GGM01S—approximated data from 846 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 1106 data.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Figure 36.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model GGM01S – 846 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 1106 data.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Figure 37.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model GGM01S—approximated data from 846 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 1106 data.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>**Figure 38.** Map of the gravity anomaly $\Delta g$ computed at the Earth's surface from 846 wavelet model, with Gaussian filtering. Black dots represent the irregular distribution of 99 synthetic data used in scalar case.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>**Figure 39.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model GGM01S–846 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 99 data.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Figure 40.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model GGM01S—approximated data from 846 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 99 data.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>**Figure 41.** Map of the gravity anomaly $\Delta g$ computed at the Earth's surface from 846 wavelet model, with Gaussian filtering. Black dots represent the irregular distribution of 99 synthetic data used in vectorial case.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Figure 42.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model GGM01S–846 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 99 data.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>**Figure 43.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model GGM01S—approximated data from 846 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 99 data.</figcaption></figure>

<table><caption>**Table 9.** Parameters of the regional tests—gravity case.</caption><thead><tr><th>Parameters</th><th>Regular case<br>(1369 data)</th><th>Irregular case<br>(576 data)</th></tr></thead><tbody><tr><td>Order of multipoles</td><td>3</td><td>3</td></tr><tr><td>Generations of the frame</td><td>2 to 10</td><td>2 to 8</td></tr><tr><td>Scales</td><td>See Table 3</td><td>See Table 3</td></tr><tr><td>Number of wavelets</td><td>7189</td><td>1962</td></tr><tr><td>W matrix</td><td>10<sup>6</sup>||</td><td>10<sup>4</sup>||</td></tr><tr><td>λ parameter</td><td>1000</td><td>50</td></tr></tbody></table>

<figure><img src="image_45.png" /><figcaption>**Figure 44.** Map of the gravity anomaly $\Delta g$ on Northern Andean Cordillera computed at the Earth's surface from model EGM 96 up to degree/order 360, with Gaussian filtering. Black dots represent the irregular distribution of 576 synthetic data.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 45.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model EGM 96 – 7189 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 1369 data.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>**Figure 46.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model EGM96—approximated data from 7189 wavelet model). The wavelet model is computed applying the scalar case on the regular distribution of 1369 data.</figcaption></figure>

repartition of data, with a larger density in the northern half. Residuals are larger in areas with large variations of the gravity anomaly, where a larger amount of data should be used to constrain the wavelet model. Contrary to the regular data distribution case, we do not notice any edge effect: indeed, those effects are of too small amplitude to be detected when the limitation of precision comes from the number of available data. Residuals between synthetic data and reconstructed data from wavelet model (Figs 51 and 54) are mostly of same order of magnitude as the assumed data noise. Some residuals can reach slightly larger values, but still very acceptable considering the usual noise on ground measurements.
<figure><img src="image_48.png" /><figcaption>**Figure 47.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model EGM 96 — 7189 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 1369 data.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>**Figure 48.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model EGM96—approximated data from 7189 wavelet model). The wavelet model is computed applying the vectorial case on the regular distribution of 1369 data.</figcaption></figure>

Finally, scalar and vectorial cases give again very similar results as for the global tests, in both cases of regular and irregular distributions of data.

## 6 CONCLUSIONS

In this paper, we investigated an appealing alternative to the spherical harmonic representation of potential fields: the Poisson multipole wavelet frames. We took into account the specificities of the magnetic and the gravity fields. Indeed, their spectral and spatial characteristics are quite different. This reflects the different distributions

<figure><img src="image_50.png" /><figcaption>**Figure 49.** Map of the gravity anomaly $\Delta g$ computed at the Earth's surface from 1962 wavelet model, with Gaussian filtering. Black dots represent the irregular distribution of 576 synthetic data used in scalar case.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**Figure 50.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model EGM 96 – 1962 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 576 data.</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>**Figure 51.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model EGM96—approximated data from 1962 wavelet model). The wavelet model is computed applying the scalar case on the irregular distribution of 576 data.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>**Figure 52.** Map of the gravity anomaly $\Delta g$ computed at the Earth's surface from 1962 wavelet model, with Gaussian filtering. Black dots represent the irregular distribution of 576 synthetic data used in vectorial case.</figcaption></figure>

of magnetic and gravity sources inside the Earth. In particular, the magnetic field is dominated by an important dipolar component. Moreover, geomagneticians and gravimetricians are used to work with different functionals of the field: intensity and vectorial field for the former, anomaly for the latter. That is the reason why we proposed two different implementations of a wavelet representation, well suited for each kind of data. The magnetic field was thus modelled with order 2 multipoles on a cubical mesh, whereas the

<figure><img src="image_54.png" /><figcaption>**Figure 53.** Map of residuals on the gravity anomaly $\Delta g$ at the Earth's surface (initial spherical harmonic model EGM 96 – 1962 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 576 data.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>**Figure 54.** Histogram of residuals on the gravity anomaly $\Delta g$ (initial synthetic data from spherical harmonic model EGM96—approximated data from 1962 wavelet model). The wavelet model is computed applying the vectorial case on the irregular distribution of 576 data.</figcaption></figure>

gravity anomaly was modelled with order 3 multipoles on more regular meshes, based on the icosahedron, so that the finest scales could be correctly represented. Thus, the frame used for gravity modelling is more redundant than the one used for magnetic modelling. As a consequence of this redundancy and of the importance of the small scales, the inverse problem needed more regularization for the case of the gravity field than for the case of the magnetic field.
We proposed two different kinds of modelling, the one which we named scalar and the other vectorial. Whereas the scalar models constitute a simple fit of a function defined on a sphere, the vectorial ones take into account a fundamental characteristic of the used wavelets, namely their multipolar nature. The wavelets can thus be considered as equivalent sources of the magnetic—respectively, gravity and field.

Our calculations proved the interest of the Poisson multipole wavelets for global and regional modelling of the magnetic and the gravity fields. We wrote the normal equations in two different ways, one taking into account in the wavelet model the spectral truncation in the data sets, and the other not. We studied the impact of the different parameters involved in the computations of various magnetic and gravity models: order of the multipoles, meshes, scales, regularization, data distribution and spectral characteristic, etc. The results presented at a global scale and derived from regularly distributed data sets showed that the wavelets well represent an harmonic field. On the one hand, the number of wavelets involved is larger than the number of spherical harmonics. Thus, the spherical harmonics are more adapted to deal with global, regularly distributed data sets. On the other hand, wavelets are considerably more interesting than spherical harmonics as soon as the data sets are dispersed or regional. Indeed, the wavelet models based on data sets with gaps do not show the strong oscillations that are generated by spherical harmonics. On the contrary, they represent the data in a faithfully and realistic way, be they globally or regionally distributed. The main structures of the fields are correctly reproduced, with the right amplitudes and locations. The interpolation errors are due to the lack of constraints in the large data gaps. They are higher for the gravity anomaly, due to the importance of the small-scale components. The results obtained on regional, regularly distributed data sets stressed the ability of a subset of the frames to represent an harmonic field with a high precision, despite some edge effects easy to bypass. The numeric gain is quite impressive, since only a few thousands of wavelets are necessary to model a field up to degree and order 360, instead of about 130 000 spherical harmonics. Moreover, the number of wavelets may be decreased by adjusting the frame parameters.

Therefore, if spherical harmonics, as already mentioned, are well suited for global and regular distributions of data, wavelets are more appropriate for local and irregular distributions of data. Thus, spherical harmonics and wavelets are complementary tools for geomagneticians, gravimetricians and geodesists. By the way, hybrid models combining for low degrees the spherical harmonics and for higher degrees the wavelets can also be considered (Freeden & Winterheuser 1997). As far as wavelets are concerned, the next step will be to apply our approach to ground-based measurements localized in a specific region, in order to derive both magnetic and gravity models that will be useful for geodynamical purposes.

## ACKNOWLEDGMENTS

The authors wish to thank the editor, Prof Fernando Sansò and two anonymous referees for their very useful comments. The authors would like to thank Lars Amsel (University of Potsdam) for his computational help, and David Renault (Université Bordeaux 1 and ENSEIRB) for fruitful discussions about spherical meshes.

All maps were plotted using the GMT software Wessel & Smith (1991). This is IPGP contribution num 2083.

## REFERENCES

Albertella, A., Sanso, F. & Sneeuw, N., 1999. Band-limited functions on a bounded spherical domain: the Slepian problem on the sphere. *Journal of Geodesy*, **73**, 436–447.

Cromwell, P.R., 1997. *Polyhedra*, p. 460, Cambridge University Press, New York.

Dahlke, St., Steidl, G. & Teschke, G., 2001. Coorbit Spaces and Banach Frames on Homogeneous Spaces with Applications to Analyzing Functions on Spheres, Zentrum für Technomathematik Technical Report, 01-13, Bremen University.

De Santis, A. & Falcone, C., 1995. Special cap models of Laplacian potentials and general fields, in, *Geodetic Theory Today*, ed. F. Sansò, p. 141–150, Springer, Berlin.

De Santis, A., De Franceschi, G. & Kerridge, D.J., 1994. Regional spherical modelling of 2D functions: the case of the critical frequency of the F2 ionospheric layer, *Computers & Geosciences*, **20**, 849–871.

De Santis, A., Falcone, C. & Torta, J.M., 1997. SHA vs. SCHA for modelling secular variation in a small region such as Italy, *J. Geomag. Geoelectr.*, **49**, 359–371.

De Santis, A. & Torta, J.M., 1997. Spherical cap harmonic analysis: a comment on its proper use for local gravity field representation, *Journal of Geodesy*, **71**, 526–532.

Freeden, W. & Winterheuser, U., 1996. Spherical wavelet transform and its discretization, *Advances in Computational Mathematics*, **5**, 51–94.

Freeden, W. & Winterheuser, U., 1997. Combined spherical harmonic and wavelet expansion—a future concept in Earth gravitational determination, *Applied and Computational Harmonic Analysis*, **4**, 1–37.

Freeden, W. & Schneider, F., 1998. An integrated wavelet concept of physical geodesy, *Journal of Geodesy*, **72**, 259–281.

Freeden, W., Gervens, T. & Schreiner, M., 1998. *Constructive Approximation on the Sphere (With Applications to Geomathematics)*, Oxford Science Publication, Clarendon Press, Oxford.

Grossman, A. & Morlet, J., 1985. *Mathematics and Physics*, Vol. 2, Streit, L. ed., World Scientific Publishing, Singapore.

Hackney, R.I. & Featherstone, W.E., 2003. Geodetic versus geophysical perspectives of the gravity anomaly, *Geophys. J. Int.*, **154**, 35–43.

Haines, G.V., 1985. Spherical cap harmonic analysis, *J. geophys. Res.*, **141**, 2583–2591.

Haines, G.V., 1990. Regional magnetic field modelling: a review, *J. Geomag. Geoelectr.*, **141**, 1001–1018.

Hicks, J.S. & Wheeling, R.F., 1959. An efficient method for generating uniformly distributed points on the surface of an n-dimensional sphere, *Communications of the ACM*, **2**, 4, 17–19.

Hipkin, R.G., 2001. The statistics of pink noise on a sphere: applications to mantle density anomalies, *Geophys. J. Int.*, **144**, 259–270.

Holme, R., Olsen, N., Rother, M. & Lühr, H., 2003. CO2—A CHAMP Magnetic Field Model, in *First CHAMP Mission Results for Gravity, Magnetic and Atmospheric Studies*, pp. 220–225, eds Reigber, C., Lühr, L. & Schwintzer, P., Springer Verlag, New York.

Holschneider, M., 1995. *Wavelets: an Analysis Tool*. Oxford Sciences Publications, Oxford.

Holschneider, M., 1996. Wavelet analysis on the sphere. *J. Math. Phys.*, **37**(8), 4156–4165.

Holschneider, M., Chambodut, A. & Mandea, M., 2003. From global to regional analysis of the magnetic field on the sphere using wavelet frames, *Phys. Earth Planet. Inter.*, **135**, 107–124.

Hwang, C., 1993. Spectral analysis using orthonormal functions with a case study on the sea surface topography, *Geophys. J. Int.*, **115**, 1148–1160.

Hwang, C. & Chen, S.K., 1997. Fully normalized spherical cap harmonics: application to the analysis of sea-level data from TOPEX/POSEIDON and ERS-1, *Geophys. J. Int.*, **129**, 450–460.

Jacobs, J.A. (Ed), 1987. *Geomagnetism*, Vol. 2, p. 627, Academic Press, London, Orlando.

Kaula, W.M., 1966. *Theory of satellite geodesy*, Waltham, Blaisdell.

Kenner, H., 1976. *Geodesic math and how to use it*, University of California Press, Berkeley.
Korte, M. & Haak, V., 2000. Modelling European Magnetic Repeat Station Data by SCHA in search of time-varying anomalies, *Phys. Earth Planet. Inter.*, **122**, 205–220.

Lemoine, F.G. et al., 1998. The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96. NASA/TP – 1998 – 206861, Greenbelt, Maryland.

Li, J., Chao, D. & Ning, J., 1995. Spherical cap harmonic expansion for local gravity field representation. *Manuscripta Geodaetica*, **20**, 265–277.

Lowes, F.J., 1974. Spatial power spectrum of the main geomagnetic field and extrapolation to the core, *Geophys. J. R. astr. Soc.*, **36**, 717–730.

Martelet, G., Sailhac, P., Moreau, F. & Diament, M., 2001. Characterization of geological boundaries using 1D-wavelet transform on gravity data: Theory and application to the Himalaya, *Geophysics*, **66**(4), 1116–1129.

Mauersberger, P., 1956. Das Mittel der Energiedichte des geomagnetischen Hauptfeldes an der Erdoberfläche und seine säkulare Aenderung. *Gerlands Beitr. Geophys.*, **65**, 207–215.

Merrill, R.T., McElhinny, M.W. & McFadden, P.L., 1996. *The Magnetic Field of the Earth*, Academic Press, San Diego.

Moritz, H., 1989. *Advanced Physical Geodesy*, 2nd edn, Wichmann, Karlsruhe.

Muller, M.E., 1959. A note on a method for generating points uniformly on N-dimensional spheres, *Comm. Assoc. Comput. Mach.*, **2**, 19–20.

Rummel, R. & van Gelderen, M., 1995. Meissl scheme—spectral character-

istics of physical geodesy. *Manuscripta Geodaetica*, **20**, 379–385.

Sailhac, P., Galdeano, A., Gibert, D., Moreau, F. & Delor, C., 2000. Identification of sources of potential fields with the continuous wavelet transform: complex wavelets and application to aeromagnetic profiles in French Guiana. *J. geophys. Res.*, **105**, B8, 19455–19475.

Sailhac, P. & Gibert, D., 2003. Identification of sources of potential fields with the continuous wavelet transform: Two-dimensional wavelets and multipolar approximations. *J. geophys. Res.*, **108**, B5, 2262, doi:10.1029/2002JB002021.

Sandwell, D.T. & Renkin, M.L., 1988. Compensation of swells and plateaus in the North Pacific: no direct evidence for mantle convection. *J. geophys. Res.*, **93**, 2775–2783.

Schroeder, P. & Sweldens, W., 1995. Spherical wavelets: efficiently representing functions on the sphere. *Computer Graphics Proceedings (SIGGRAPH 95)*, 161–172.

Tapley, B.D., Bettadpur, S., Watkins, M. & Reigber, C., 2004. The gravity recovery and climate experiment: mission overview and early results. *Geophys. Res. Lett.*, **31**, L09607:10.1029/2004GL019920.

Tarantola, A., 1987. *Inverse problem theory*, Elsevier, New York.

Thébault, E., Schott, J.J., Mandea, M. & Hoffbeck, J.P., 2004. A new proposal for spherical cap harmonic analysis. *Geophys. J. Int.*, **159**, 83–103.

Wessel, P. & Smith, W.H.F., 1991. Free software helps map and display data. *Eos Trans. Am. geophys. Un.*, **729**, 441.