# A proposal for regional modelling at the Earth's surface, R-SCHA2D

Erwan Thébault

► To cite this version:

Erwan Thébault. A proposal for regional modelling at the Earth's surface, R-SCHA2D. Geophysical Journal International, 2008, 174, pp.118-134. 10.1111/j.1365-246X.2008.03823.x. insu-03603717

HAL Id: insu-03603717

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

Submitted on 10 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
<figure><img src="image_2.png" /></figure>

# A proposal for regional modelling at the Earth's surface, R-SCHA2D

Erwan Thébault

Equipe de Géomagnétisme Institut de Physique du Globe de Paris et Université Paris Diderot CNRS-INSU, 4 place Jussieu, 75005 Paris, France.

E-mail: Ethebault@jpgp.jussieu.fr

Accepted 2008 April 14. Received 2008 April 11; in original form 2007 November 27

## SUMMARY

Modelling the magnetic field at a regional scale has always been challenging. Several techniques were proposed in the past but few solve Laplace equation and meet the requirements of potential theory. Among the regional techniques, the Spherical Cap Harmonic Analysis (SCHA) that solve Laplace equation has become widespread. However, recognizing the inherent difficulties of SCHA to represent the main field over a small spherical cap, a revision of this technique the R-SCHA was proposed to deal with multilevel magnetic measurements. For the reasons explained in this paper, R-SCHA is not suitable either when data are available at one surface only. The present paper proposes an alternative to SCHA and R-SCHA, named R-SCHA2D, that is suitable for modelling geomagnetic data available at a quasi-constant altitude. The efficiency of R-SCHA2D is successfully tested on a set of repeat station data. The magnetic field over France and its secular variation are reliably represented between 1965 and 2007.5. The regional parameters are obtained from the joint inversion in space and time of 29 metropolitan repeat stations and 1 observatory vector data using time power series as a temporal basis function. This approach does not require additional constraints like adding synthetic data, regularizing the inverse problem, or removing a trend prior to the inversion to be stable. It is shown that the geomagnetic jerks over France are detected in the regional parameters themselves and that it is possible to estimate the crustal bias. Considering the real data distribution in space and time, and the parameters used for this modelling, the R-SCHA2D model compared with the CM4 comprehensive model showed no indisputable evidence proving the existence of a secular variation anomaly over France. The main field and the secular variation are sketched for various epochs, in particular for epoch 2007.5. The present approach does neither address the problem of internal/external fields separation nor the upward continuation of the model to satellite altitudes.

**Key words:** Spatial analysis; Magnetic anomalies: modelling and interpretation; Satellite magnetics.

## 1 INTRODUCTION

There are many cases where magnetic field models with temporal variations are requested. The demand is usually related to data reduction to a common epoch, prediction of the magnetic field over a specific area, or updating of the magnetic charts for navigational purposes. To meet these aims, the International Association of Geomagnetism and Aeronomy (IAGA) publishes the International Geomagnetic Reference Field (IGRF) every 5 yr. At the same time, some individual countries measure the magnetic field at ground based repeat stations in order to publish their own national charts useful for airplane navigation. Moreover, since the beginning of the satellite era and the past decades of magnetic field measurements, ambitious efforts have been made to derive magnetic field models continuous in time from epochs 1960 to 2006 (Sabaka et al. 2004; Olsen et al. 2006; Maus et al. 2006, for instance).

Apart from complementing spherical harmonic models that are not well constrained in remote regions, there are lots of other scientific interests for representing the field at regional scales. By construction, a regional model captures the magnetic field features difficult to represent in spherical harmonics because of their small spatial extensions. This approach is also useful when a particularly dense data distribution is available in a region. Whatever the physical origin of the magnetic field sources, the regional modelling approach is able to represent the magnetic field in a spatial spectral band, typically between 4000 and 5 km wavelengths, that is difficult to reach in spherical harmonics. Such methods allowed, for instance, to describe the lithospheric field by exploiting all kinds of magnetic data available from ground to satellite altitudes over France (Thébault et al. 2006a). Rapid external field variations, together with the equivalent current density function, may
be monitored by taking advantage of a dense network of observatories, as was done in Europe (for instance, Haines & Torta 1994; Torta et al. 1997). Historical data, when processed at a regional scale, may also provide a high resolution model of the magnetic field for the past millennia (Pavón-Carrasco et al. 2007). This, in turns, could offer an archaeological dating tool. It could also help to better localize the time of geomagnetic jerks that do not seem to be global in occurrence (Chambodut & Mandea 2005). Other debated questions could potentially be addressed. Geomagnetic fields induced by ocean tides (Tyler et al. 2003) may have a regional contribution (Thébault & M. Mandea 2006). The existence of a continental scale secular variation induced by the lithospheric field could also be explored but so far failed to be evidenced (Korte & Holme 2003). For all these reasons, a well-established regional technique is required to scrutinize the small spatial structures of the magnetic field.

Among the regional approaches that guarantee the conditions imposed by the electromagnetic theory, which requires that the magnetic potential satisfies Laplace's equation in source-free regions, the Spherical Cap Harmonic Analysis (SCHA) has become widespread (Haines 1985a). This technique was first applied to study the secular variation over Canada (Haines 1985b) but later proved to suffer from various practical drawbacks that were better identified in many studies (De Santis & Falcone 1995; De Santis et al. 1999). It was observed that the modelling of the main field showed gross instabilities (De Santis et al. 1990; Torta et al. 1992) and sensitivity to data noise and gaps (Thébault & Gaya-Piqué 2008). These problems are usually reduced by adding synthetic data calculated with spherical harmonic models (Torta et al. 1992), or by regularizing the inverse problem (Korte & Holme 2003). In addition, common to almost all studies, SCHA models are obtained by artificially opening the cap's aperture and data require a 'detrending', a removal of a main field model (Bullard 1967). Regional models are thus not computed independently from global models. At last, modelling data measured at different altitudes is usually problematic as a result of an inappropriate radial variation of the basis functions.

In opposition to a common belief, but following the work of De Santis et al. (1999) and unpublished results of Lowes (in Langel & Hinze 1998, chapter 5.3), Thébault et al. (2004, 2006b) argued that SCHA is not a well founded technique because of its lack of generality. Although it may provide realistic results (see Torta et al. 2006, for a review), it must be used with caution if the data do not meet specific requirements, roughly if the magnetic field does not average out over the spherical cap (Thébault et al. 2006b).

In order to remedy some of these difficulties, the Revised-Spherical Cap Harmonic Analysis (R-SCHA) was proposed to allow the simultaneous inversion of data measured at different altitudes (Thébault et al. 2006b). However, for the reasons discussed in Appendix A of this work, R-SCHA is not convenient either for representing data measured at one altitude only. Continuing the quest towards the modelling at a regional scale, another approach inspired from SCHA and R-SCHA, named R-SCHA2D is presented in this paper. The steps leading to the construction of the regional basis functions are explained in the first section. This is followed by a discussion about the completeness of the proposed solution. The formalism is then applied to the French repeat station measurements between epochs 1965 and 2007.5. The regional model is obtained by a joint inversion in space and time using a polynomial temporal basis. It comprises the main field, the secular variation and, in a first version, the large-scale lithospheric field structures. It is thus demonstrated that a regional main field model, having the required potential field properties, can be obtained from ground-based measurements independently from global models.

## 2 METHOD

The SCHA technique is introduced by Haines (1985a) and some important properties are discussed by De Santis et al. (1999). An exhaustive reference list related to SCHA and its caveats is provided by Torta et al. (2006). For this reason, only the fundamental equations are written here.

The lateral boundaries of the spherical cap are sketched in Fig. 1. Let us assume that we wish to find the expression of $V_1$, a potential solving the following incomplete boundary value problem:

$$ \nabla^2 V_1 = 0 \qquad (1a) $$

$$ \left. \frac{\partial V_1}{\partial \theta} \right|_{\partial \theta_0 \Omega} = 0 \qquad (1b) $$

$$ \left. \frac{\partial V_1}{\partial r} \right|_{\partial a \Omega} = \tilde{Z}_a(a, \theta, \varphi), \qquad (1c) $$

where $\tilde{Z}_a$ is the vertical component of the magnetic field at the Earth's surface rotated in the cap's frame. The spherical cap surface domain is defined by $\partial_a \Omega = \{r = a, 0 \le \theta \le \theta_0, 0 \le \varphi \le 2\pi\}$ and $\partial_{\theta_0} \Omega = \{r > 0, \theta = \theta_0, 0 \le \varphi \le 2\pi\}$. In spherical coordinates, the Laplace operator is:

$$ \nabla^2 = \left( \frac{\partial^2}{\partial r^2} + \frac{2}{r} \frac{\partial}{\partial r} \right) + \frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial}{\partial \theta} \right) + \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2}{\partial \phi^2} \qquad (2) $$

with the radius $r$, the colatitude $\theta$ and the longitude $\varphi$ in the cap's reference frame (see Fig. 1). Laplace equation is solved by separating the variables:

$$ V_1 = R(r)\Theta(\theta)\Phi(\varphi). \qquad (3) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Repeat stations with their code (dots), Chambon-La-Forêt observatory (star) and the spherical cap (solid line) with its geometric characteristics. The three unit vectors used in this study are specified at the bottom-right-hand side of the spherical cap.</figcaption></figure>

This gives the well-known triplet of ordinary differential equations:

$$ \frac{d^2 \Phi (\varphi)}{d\varphi^2} + m^2 \Phi (\varphi) = 0 \qquad (4a) $$

$$ \frac{1}{\sin \theta} \frac{d}{d\theta} \left[ \sin \theta \frac{d\Theta(\theta)}{d\theta} \right] - \frac{m^2}{\sin^2 \theta} \Theta(\theta) = -\lambda \Theta(\theta) \qquad (4b) $$

$$ \frac{d}{dr} \left[ r^2 \frac{dR(r)}{dr} \right] = \lambda R(r). \qquad (4c) $$

The homogeneous eq. (1b) implies that the solution of (4b) are the associated Legendre functions of real degrees $n_k$ and integer order $m$, with $n_k \ge m$ (Haines 1985a). The solution may be written in terms of the hypergeometric function $\mathcal{F}$ (Hobson 1965) and we have the classical result (Abramowitz & Stegun 1965):

$$ \Theta(x) \equiv P_{n_k}^m(\cos \theta) = \frac{(-1)^m}{2^m m!} N_{n_k}^m \sin^m \theta \mathcal{F} \left( m - n_k, n_k + m + 1, 1 + m; \sin^2 \frac{\theta}{2} \right) \qquad (5) $$

with $N_{n_k}^m$ the normalization factor. The expression of the potential $V_1$ expands in terms of two infinite series:

$$ V_1(r, \theta, \varphi, t) = a \sum_{k \ge m} \sum_{m \ge 0} \left(\frac{a}{r}\right)^{n_k+1} \left[ G_{n_k}^{i,m}(t) \cos(m\varphi) + H_{n_k}^{i,m}(t) \sin(m\varphi) \right] P_{n_k}^m(\cos\theta) \\ + a \sum_{k \ge m} \sum_{m \ge 0} \left(\frac{r}{a}\right)^{n_k} \left[ G_{n_k}^{e,m}(t) \cos(m\varphi) + H_{n_k}^{e,m}(t) \sin(m\varphi) \right] P_{n_k}^m(\cos\theta) \qquad (6) $$

with $a$ the reference radius of the modelling, taken here as the mean Earth's radius, and $t$ the epoch. Note that the complete solution includes both sets of regional parameters $\{G_{n_k}^{m,i}, H_{n_k}^{m,i}\}$ and $\{G_{n_k}^{m,e}, H_{n_k}^{m,e}\}$. The omission of $\{G_{n_k}^{m,e}, H_{n_k}^{m,e}\}$ in the solution is illustrated in the next section.

The term $G_0^{i,0}$ is essential when dealing with a magnetic field that does not average out over the surface. De Santis et al. (1999) showed that $G_0^{i,0} = -\bar{Z}(a, \theta, \varphi)$ (considering only one basis function of SCHA as is the case here), where $\bar{Z}$ is the average of the vertical field over the spherical cap. This term is only defined for $n_k = 0$, a solution of the homogeneous condition (1b). Eq. (1b) was thus chosen instead of $V_1|_{\theta_0 \Omega} = 0$ because the solution comprises the term $G_0^{i,0}$.

In a very general situation, $V_1$ is not sufficient to fully represent the magnetic field. In the spherical cap reference frame defined by the unit vectors $(\mathbf{e}_\theta, \mathbf{e}_\varphi, \mathbf{e}_r)$ (Fig. 1), we have $\mathbf{B}_1 = B_{1,\theta} \mathbf{e}_\theta + B_{1,\varphi} \mathbf{e}_\varphi + B_{1,r} \mathbf{e}_r$ and eq. (1b) implicitly imposes that $\int_0^{2\pi} B_{1,\theta}(\theta_0, \varphi, r) \sin \theta_0 d\varphi = 0$. However, the component $\tilde{B}_\theta$, in the cap's reference frame, of an internal magnetic field $\mathbf{B} = \tilde{B}_\theta \mathbf{e}_\theta + \tilde{B}_\varphi \mathbf{e}_\varphi + \tilde{B}_r \mathbf{e}_r$ expanded in spherical harmonics writes along this closed contour:

$$ \int_0^{2\pi} \tilde{B}_\theta (\theta_0, \varphi, r) \sin \theta_0 d\varphi = 2\pi a \sin \theta_0 \sum_{n \ge 1} \tilde{g}_n^0 \left. \frac{dP_n^0}{d\theta} \right|_{\theta_0} \qquad (7) $$

with $\tilde{g}_n^0$ the spherical harmonic Gauss coefficients of order $m = 0$ rotated in the cap's reference frame. A complementary potential $V_2$ must thus be added to $V_1$ so that $\mathbf{B}_2 = -\nabla V_2$ satisfies (7). It is also important to choose a potential $V_2$ satisfying Laplace equation and we may
look for the solution to the following problem:

$$ \Delta V_2 = 0 \qquad (8a) $$

$$ \left. \frac{1}{a} \frac{\partial V_2}{\partial \theta} \right|_{\partial \theta_0 \Omega} = \tilde{X}(a, \varphi, \theta_0) \qquad (8b) $$

with $\tilde{X}$ the North component of the magnetic field in the cap's reference frame ($\tilde{X} = -B_0$). In principle, any potential $V_2$ satisfying (8a) and (8b) is suitable. However, eq. (1c) must hold when adding $V_2$ to $V_1$ and the further constrain is added:

$$ \left. \frac{\partial V_2}{\partial r} \right|_{\partial \theta_0 \Omega} = 0. \qquad (9) $$

The solution to this second boundary value problem is non-unique because boundary conditions were not set at infinity (see Appendix A for the condition at $r = \infty$ and the justification of our choice). Setting $r = ae^t$, eq. (4c) can be rewritten in the form:

$$ \frac{d^2 R_2(t)}{dt^2} + \frac{d R_2(t)}{dt} - \lambda_2 R_2(t) = 0. \qquad (10) $$

This gives two possibilities. For $\lambda_2 \neq -1/4$, with $\lambda_2 = n_2(n_2 + 1)$, the solution takes the form:

$$ R_2(r) = A \left(\frac{r}{a}\right)^{n_2} + B \left(\frac{a}{r}\right)^{n_2+1} \qquad (11) $$

and for $\lambda_2 = -1/4$, with $n = -1/2$:

$$ R_2(r) = \sqrt{\frac{a}{r}} \left[ A \ln \left( \frac{r}{a} \right) + B \right]. \qquad (12) $$

The condition (9) imposes that $n_2A - (n_2 + 1)B = 0$ in eq. (11). Within this infinite set of solutions, we may choose $n_2 = 0$ and $B = 0$ (i.e. $\lambda_2 = 0$) or $n_2 = -1$ and $A = 0$ [i.e. $\lambda_2 < -1/4$]. Note that the case where $n_2$ is complex is discussed in Appendix A. In both situations, the radial function $R_2$ is constant and since $P_0^m = P_{-1}^m$ (Robin 1957), both solutions are identical. The potential $V_2$ writes in that case (Thébault et al. 2006b):

$$ V_2(r, \theta, \varphi, t) = a \sum_{m \ge 0} R_2 \left[ G_0^m \cos(m\varphi) + H_0^m \sin(m\varphi) \right] P_0^m(\cos\theta) $$

with $P_0^m = N_0^m \tan^m \frac{\theta}{2}$. However, for $m = 0$, $P_0^0 = N_0^0$ implies $dP_0^0/d\theta = 0$. Therefore, eq. (8b) and condition (7) are not satisfied and we thus discard this solution.

To the contrary, the condition (9) with eq. (12) imposes that $B = 2A$ and the radial function writes:

$$ R_2(r) = A \sqrt{\frac{a}{r}} \left[ \ln \left( \frac{r}{a} \right) + 2 \right]. \qquad (13) $$

Replacing $n_2 = -1/2$ in (5) gives:

$$ P_{-1/2}^{m}(\cos \theta) = \frac{(-1)^m}{2^m m!} N_{-1/2}^{m} \sin^m \theta \mathcal{F} \left( m + 1/2, m + 1/2, 1 + m; \sin^2 \frac{\theta}{2} \right). $$

The functions $P_{-1/2}^m$ are a special case of the Mehler functions obtained in Thébault et al. (2006b). They have the required property that $dP_{-1/2}^m/d\theta \neq 0$ for any $m \in \mathbb{N}$.

Finally the additional potential $V_2$ may be chosen as:

$$ V_2(r, \theta, \varphi, t) = a \sum_{m \ge 0} R_2(r) \left[ G_{-1/2}^{m}(t) \cos(m\varphi) + H_{-1/2}^{m}(t) \sin(m\varphi) \right] P_{-1/2}^{m}(\cos\theta). \qquad (14) $$

The complete solution $V$ resulting from the superposition of potentials $V_1$ and $V_2$ solves the boundary value problem:

$$ \nabla^2 V = 0 \qquad (15a) $$

$$ \left. \frac{1}{a} \frac{\partial V}{\partial \theta} \right|_{\partial \theta_0 \Omega} = \tilde{X}(r, \theta_0, \varphi) \qquad (15b) $$

$$ \left. \frac{\partial V}{\partial r} \right|_{\partial \theta_0 \Omega} = \tilde{Z}(a, \theta, \varphi). \qquad (15c) $$

Considering (6), (14) and the boundary condition (1b) that defines the real degrees $n_k$, there is no longer any implicit restriction on $X$ and $Z$ when considering data at one surface only. Any magnetic field can be represented at a regional scale over a spherical cap, in particular at altitude $r = a$ with these basis functions. In theory, the solution developed here is not adapted to problems involving magnetic data available at different altitudes. Even though eq. (14) is part of the solution within a semi-infinite cone, by neglecting the boundary condition at infinity for potential $V_2$ we are missing the part of the solution that drives the correct upward/downward continuation of the model (Appendix A). The behaviour of the R-SCHA2D technique with altitude was not systematically tested, but there are cases for which the missing functions are indeed negligible. These cases were encountered when dealing with the lithospheric field, a topic that extends the purpose of this paper. For the main field, the solution given by (6) and (14) provides a reasonable approximation when altitude differences are weak (say several tens of kilometres depending on the cap's size).
The revised spherical cap method proposed here is named R-SCHA2D in the following. Contrary to R-SCHA, which was purposefully designed to work well for multilevel data (Thébault et al. 2006b), R-SCHA2D is more adapted to a single surface distribution of data. Compared to SCHA the expression of V involves only one set of Legendre basis function, the k-m even (Haines 1985a) with all parameters, including {$G_{nk}^{m,e}$, $H_{nk}^{m,e}$}. The other set of SCHA basis functions, the k-m odd, is replaced by (14) that only requires 2m + 1 parameters. A synthetic comparison between SCHA and R-SCHA2D is given in Thébault & Gaya-Piqué (2008).

## 2.1 Orthogonality and completeness of the solution

The basis functions of the potential V and its gradient **B** = -∇V are not orthogonal since ∇V₁ and ∇V₂ are not mutually orthogonal. In terms of the variable θ, the functions $P_{nk}^m$ are orthogonal with respect to the inner product $\int_0^{\theta_0} f(\theta)g(\theta) \sin\theta d\theta$ at each m and form a complete basis in the Hilbert space L² (0, θ₀). The resolution matrix of the inverse problem, calculated from a set of discrete values of the magnetic field **B** over the spherical cap, is block-diagonal but the non-diagonal terms are weak compared to diagonal elements (Appendix B).

It is not straightforward to establish a full set of relationships between regional parameters and spherical harmonic Gauss coefficients. This is due to the fact that the boundary conditions are not set on a closed surface but on an open surface instead. As a result, the coefficients {$G_{nk}^{m,i}$, $H_{nk}^{m,i}$} and {$G_{nk}^{m,e}$, $H_{nk}^{m,e}$} require the knowledge of both $\tilde{X}$ and $\tilde{Z}$ over the spherical cap to be uniquely determined. This situation suggests that dealing with intensity data will raise some difficulties and that probably {$G_{nk}^{m,e}$, $H_{nk}^{m,e}$} will be difficult to estimate in such situation if no vector data are available. In order to obtain the analytical expressions of the regional parameters directly from the boundary conditions, we must either close the conical volume and set boundary conditions on each surface (Thébault et al. 2006b), or set boundary conditions at infinity (see Appendix A). In both cases, this changes the mathematical expression of potential V₂ (potential V₁ remains unchanged unless we choose another homogeneous boundary conditions in eq. 1b).

To the contrary, the boundary condition for the second potential V₂ is set on a closed contour and a unique relationship between regional parameters and spherical harmonic Gauss coefficients can be established. From (8b) we have, integrating between $\varphi = 0$ and $2\pi$ and setting $R_2(a) = 1$ in (13), the simple expression at each order m:

$$ \left\{ \begin{matrix} G_{-1/2}^m \\ H_{-1/2}^m \end{matrix} \right\} = \sum_{l \ge m} \left\{ \begin{matrix} \tilde{g}_l^m \\ \tilde{h}_l^m \end{matrix} \right\} \frac{dP_l^m(\theta_0)}{d\theta} \frac{dP_{-1/2}^m(\theta_0)}{d\theta} \quad (16) $$

with {$\tilde{g}_l^m$, $\tilde{h}_l^m$} the spherical harmonic Gauss coefficients of a global model rotated in the cap's reference frame. Note that the Gauss coefficients can include internal and external contributions in this expression (writing for instance, $\tilde{g}_l^m = \tilde{g}_l^{m,i} + \tilde{g}_l^{m,e}$, the sum of internal and external spherical harmonic Gauss coefficients). This relationship may be very useful for constraining ill-posed inverse problem by imposing a continuity between spherical harmonics and R-SCHA2D main field or secular variation models at the cap's edges.

The completeness of the solution was tested numerically. In order to verify that the basis functions are suitable for representing the magnetic field, a numerical analysis is performed on a set of synthetic vector data. The synthetic vector data are distributed homogeneously at the Earth's surface within a spherical cap of 5° aperture. At each location, the core and crustal fields are calculated for spherical harmonic degrees 1-90 using Pomme 3.1 model (Maus et al. 2006). In practice, infinite series in (6) and (14) are truncated to maximum indices $K_{max}$ and $M_{max}$ (with $M_{max} = K_{max}$ here). The residual mean square (RMS) between the Pomme 3.1 model and the R-SCHA2D model obtained by least-squares is calculated for different series expansion. Table 1 shows that the RMS indeed converges to zero and demonstrates the completeness of R-SCHA2D basis functions for single surface data. For $K_{max} = 20$ and a cap half angle of 5°, the model represents the total field up to 0.4 nT. For a spatial index as low as $K_{max} = 4$, the main field is resolved with an approximation of 15 nT, with a slightly higher RMS for Z than for the horizontal component. When the set of parameters {$G_{nk}^{m,e}$, $H_{nk}^{m,e}$} are excluded from the modelling, the model converges with an error of about 280 nT. These coefficients are thus necessary even in the case of internal magnetic field only, which also means that they cannot be employed for separating the internal and the external fields. The main result here is that the regional modelling is able to capture the main magnetic field with a few parameters only (44 parameters for $K_{max} = 4$).

# 3 MAGNETIC FIELD MODEL OVER FRANCE BETWEEN 1965 AND 2007.5

## 3.1 Data

From 1965 onward, the French repeat stations network comprises 29 metropolitan stations located 150-200 km apart. The spatial data distribution remained essentially the same during the five decades of measurements and the temporal distribution is nearly homogeneous. Repeat station data were visited every 5 yr with two extra visits of subnetworks in 1985.5 and 2004.5 (Table 2). The observatory annual means of the Chambon-La-Forêt observatory are used every year between 1965.5 and 2007.5. The latest repeat station survey, started in 2007 July, is provisional and incomplete. Two repeat station measurements are missing in 2007.5 at Garchy (GAR) and Delettes (DEL) (see Fig. 1). One measurement in 2002.5 for the repeat station SAI proved to be a possible outlier and was removed. The altitude difference between the lowest and the highest repeat station is less than 1km in the geocentric reference frame and R-SCHA2D is thus valid.
<table><caption>**Table 1.** Convergence of R-SCHA2D model towards Pomme3.1 (spherical harmonic degrees 1–90 inclusive) according to the maximum spatial index $K_{max}$. RMS 1 gives the error in nT between both models considering all sets of regional parameters. RMS 2 gives the error when the so-called external field parameters are withdrawn.</caption><thead><tr><th>KMAX</th><th>RMS 1</th><th>RMS 2</th></tr></thead><tbody><tr><td>1</td><td>340.3</td><td>811.1</td></tr><tr><td>2</td><td>53.9</td><td>467.9</td></tr><tr><td>3</td><td>25.3</td><td>416.4</td></tr><tr><td>4</td><td>15.7</td><td>388.5</td></tr><tr><td>5</td><td>10.9</td><td>370.6</td></tr><tr><td>6</td><td>7.7</td><td>357.3</td></tr><tr><td>7</td><td>5.7</td><td>347.3</td></tr><tr><td>8</td><td>4.4</td><td>339.2</td></tr><tr><td>9</td><td>3.4</td><td>332.3</td></tr><tr><td>10</td><td>2.7</td><td>326.3</td></tr><tr><td>11</td><td>2.2</td><td>321</td></tr><tr><td>12</td><td>1.7</td><td>316.2</td></tr><tr><td>13</td><td>1.4</td><td>311.9</td></tr><tr><td>14</td><td>1.2</td><td>307.9</td></tr><tr><td>15</td><td>1</td><td>304.2</td></tr><tr><td>16</td><td>0.8</td><td>300.8</td></tr><tr><td>17</td><td>0.7</td><td>297.5</td></tr><tr><td>18</td><td>0.6</td><td>294.3</td></tr><tr><td>19</td><td>0.5</td><td>291.3</td></tr><tr><td>20</td><td>0.4</td><td>288.3</td></tr></tbody></table>

<table><caption>**Table 2.** Temporal distribution of the French repeat station measurements posterior to 1965.</caption><thead><tr><th></th><th>1965</th><th>1967.5</th><th>1972.5</th><th>1977.5</th><th>1982.5</th><th>1985.5</th><th>1987.5</th><th>1992.5</th><th>1997.5</th><th>2002.5</th><th>2004.5</th><th>2007.5</th></tr></thead><tbody><tr><td>PLO</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>SEG</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>NAN</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>LAM</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>LEO</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>MIG</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>SAI</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X*</td><td>-</td><td>X</td></tr><tr><td>AUV</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>CLA</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>LEF</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>NUR</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>NEU</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>SAL</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>DEL</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>-</td></tr><tr><td>RIE</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>TUC</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>LIZ</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>GAR</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>-</td></tr><tr><td>LEM</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>MAR</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>CLC</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>MAI</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>VIL</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>FRA</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>DRA</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>AIL</td><td>-</td><td>-</td><td>-</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>BIO</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>HAG</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td></tr><tr><td>MAH</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>-</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr></tbody></table>

Note: The * symbol is used when the measurement is considered as an outlier.

Prior to the real data inversion, R-SCHA2D was shown to be stable with regard to a normally distributed error of 10 nT added to a set of synthetic data calculated at repeat station locations (Fig. 1). For the purpose of discussion, no correction for external or crustal field contamination is applied to the repeat station data in a first iteration of the modelling. A second and definitive spatiotemporal model will be obtained after crust and external field corrections.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** R-SCHA2D model in 2000 calculated with uncorrected data. Oscillations in the centre of the figure are due to the strongest lithospheric field anomaly in France: the Paris Basin anomaly. Units are in degree for the Declination D and in nT for Z.</figcaption></figure>

The basis functions developed above are used to model the magnetic field and its secular variation over France between 1965 and 2007.5. The time variation of the regional parameters $g(t) = \{G_{n_k}^{m,i}, H_{n_k}^{m,i}, G_{n_k}^{m,e}, H_{n_k}^{m,e}, G_{-1/2}^{m}, H_{-1/2}^{m}\}$ are represented by the power series for each parameter $g_n^m(t)$:

$$g_n^m(t) = \sum_{q=0}^{Q} \alpha_{q,n}^m t^q. \qquad (17)$$

The coefficients $\alpha_{q,n}^m$ are coestimated in a least-squares sense in space and time. The maximum temporal series expansion is $Q = 6$. In the case of a more sparse distribution in space and time, it may be necessary to coestimate the temporal parameters with cubic B-splines (see for instance, Korte & Holme 2003) in order to mitigate the problem of wild oscillations that may occur when using the temporal polynomials. The general least-squares solution is:

$$g(t) = (\mathbf{A}^T \mathbf{C}_d^{-1} \mathbf{A} + \gamma \mathbf{E}) \mathbf{A}^T \mathbf{C}_d^{-1} \mathbf{b}(t)$$

with $\mathbf{A}$ the operator relating the regional parameters $g(t)$ and the data $\mathbf{b}(t)$ using the gradient of (6) and (14). $\mathbf{C}_d^{-1}$ is the data error covariance matrix set to identity in the absence of objective criteria for weighting the data. $\mathbf{E}$ is the damping matrix that can be defined, for instance, as the norm of the regional field with $\gamma$ the Lagrange multiplier. In this application, $\gamma = 0$. The maximum spatial index $K_{max}$ in eq. (6) is chosen according to the spatial data distribution. A mean distance of 200 km between two repeat stations imposes a minimum resolvable wavelength of 400 km. For the spherical cap used here, centred on the geographic longitude $\Phi_0 = 2^\circ$ and latitude $\Theta = 46.5^\circ$, with an open half angle of $5^\circ$, a model up to 400 km wavelength corresponds to $K_{max} = 4$. The data distribution does not bear higher series expansion if the inverse problem is not regularized.

## 3.2 Initial model and corrections

The total inversion standard deviation is 28 nT and a stable residual is obtained for each component for each epoch. The average RMS over the whole time-series is 23 nT for X, 20 nT for Y and 24 nT for Z. The comprehensive CM4 model (Sabaka et al. 2004) including core and lithospheric field contributions gives an RMS of 43 nT for X, 41 nT for Y and 53 nT for Z between 1965 and 2002.5. To that respect, the regional model is superior to CM4 because it takes the regional crustal sources and the residual external field into account. According to Table 1, an error of 15 nT is expected when a perfect set of main field data is modelled. The obtained RMS in the real situation is twice as large. Removing a known main field model (i.e. detrending the data) before carrying out the inverse problem does not improve the misfit. The misfit does not come either from a convergence problem because of the large wavelengths or from a systematic high noise level in the repeat station measurements. As it is discussed below, repeat stations lie on local crustal anomalies. This relatively high RMS reveals that there is power in the data at wavelengths shorter than those included in the model. A sketch of the main field is provided in Fig. 2 for epoch 2000 that shows oscillations as a result of the lithospheric field contamination.

### 3.2.1 Crustal bias

Theoretically, it should be possible to identify the parameters that do not vary with time. They should account for the lithospheric static field. However, it is difficult to isolate the parameters that exhibit nearly no time variation because of various sources of error like the modelling errors (due to the low spatial index in the series expansion) and external and crustal fields contaminations. When plotting the difference between the regional main field models and CM4 every year between 1965 and 2002.5, quasi-static features, shown in Fig. 3 for the Z component,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Average difference between R-SCHA2D and CM4 models over all epochs. The Z component shown here is quasi-static (filled contours) and correlates well with the Z component of the lithospheric field at 400 km wavelength (black countour lines). Solid lines represent positive values, dashed lines represent negative values, contour spacing is 5 nT.</figcaption></figure>

become visible. A magnetic field positive difference is visible in the northern part of France and two negative lobes in the southern part. There is also a positive difference in the Eastern part over the Marlenheim (MAH) repeat station (Fig. 1). These features are reminiscent of the lithospheric field model over France (Thébault *et al.* 2006a). The R-SCHA lithospheric field model calculated for wavelengths larger than 400 km compares well with the quasi-static features obtained in this study. This comparison, shown in Fig. 3 shows that the regional model also captures the essence of the lithospheric field over France for each epoch. The main magnetic anomaly over France at 5 km altitude, the Paris Basin Anomaly, has an amplitude comprised between -150 and 300 nT that is consistent with the crustal bias observed at the repeat stations. This crustal bias is responsible for the spatial oscillations in Fig. 2. These oscillations are genuine over central France but spread over large distances in regions where no data are available. They hamper the robust estimation of the secular variation from the regional parameter derivatives.

### 3.2.2 External field contamination

Another effect overlaps with the secular variation over France. The original data were reduced to the Chambon-La-Forêt observatory annual means and there is a good evidence that the parameter oscillations are partly correlated with the 11 yr solar cycle activity. This appears, for instance, when comparing measurements with the CM4 model. A magnetic jerk is conceived as a sudden change in the slope of the secular variation. Unfortunately, the solar cycle overlap with the apparent periodicity of jerk occurrence over France, roughly every 10 yr since 1969. Repeat station data are smeared with this small effect that contaminate the regional parameter time derivatives. The problem concerning the internal or external origin of the magnetic field sudden changes remains to be clarified at regional scales and the overlap between 'fast' internal and 'slow' external variations will be problematic until a good method for separating internal and external fields is found. Till then, the regional model will strongly depend on data pre-processing and on our ability to identify and filter out the unwanted contributions.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The external field template is calculated by comparing CM4 and CLF annual means. The procedure is described by Verbanac *et al.* (2007). The residual is anticorrelated with the Ap solar index showing that CLF annual mean observatory are contaminated with long-term external fields.</figcaption></figure>

## 3.3 A ‘clean’ regional field model

A second iteration of the modelling is performed using data corrected for the crustal and external fields. The crustal bias at the repeat station is estimated from the R-SCHA lithospheric field model over France (Thébault *et al.* 2006a). The principle of the *ad hoc* external correction applied is the same as in Korte & Thébault (2007) and relies on the CM4 model calculated between 1965 and 2002.5 (Sabaka *et al.* 2004). The CM4 main field annual means calculated at Chambon-La-Forêt are subtracted from the Chambon-La-Forêt annual means, then a constant value for each magnetic field components is removed. This leaves us with an estimate of the solar cycle effect at the observatory between 1965 and 2002.5 that is assumed to influence the repeat station data reduction to the same accuracy. This template is then subtracted from repeat station data. The external field correction template, shown in Fig. 4, is comparable to Fig. 5(a) of Verbanac *et al.* (2007). The external field contamination is particularly prominent for the North component of the magnetic field and the correction is more efficient on this component (Fig. 5). Note that the external field correction is not made for data posterior to 2002.5. The total inversion standard deviation is now 13.8 nT with RMS of 11, 10 and 12 nT, respectively for X, Y and Z. These values are consistent with the error predicted in Table 1.

A few magnetic field models between epochs 1970 and 2007.5 are sketched in Figs 6 and 7. For simplicity the declination *D* and the vertical component *Z* are represented only. *D* is directly computed from *X* and *Y*. The R-SCHA2D main field model is stable and is equivalent to the comprehensive model after crustal and external field correction. The deviations between CM4 and R-SCHA2D are notable in oceanic regions where no repeat station data were available. Note that spurious oscillations are noted in 2007.5 as a result of the two missing data for Garchy (GAR) and Delettes (DEL). Moreover, the provisional data are potentially noisier as no external field correction was applied. The model is thus less constrained by the data than for other epochs. In any case, this ‘clean’ main field model no more exhibits oscillation (compare Figs 2 and 7 in 2000.0) thanks to the crustal and external field corrections.

## 3.4 Secular variation modelling

Compared to a model obtained with homogeneously distributed data in space, the situation is stable for the main field but suffers from edge effects. These edge effects are enhanced when computing the temporal derivative of the magnetic field. In order to remedy this difficulty the secular variation models are calculated with the regional parameters to the spatial index $K_{max} = 2$ only. The purpose is to show the ability of R-SCHA2D to directly process the data. However, it is important to keep in mind that the secular variation model is more stable when modelling the time differences of repeat station measurements (see Haines, 1993 for a discussion). According to the modelling method, the $G_0^0$ Gauss coefficient reveals the time variation of the component *Z* averaged over the spherical cap. In the same vein, {$G_{-1/2}^m, H_{-1/2}^m$} are specific to the magnetic field horizontal component along the closed contour $\partial\theta_0\Omega$ at $r = a$. It is thus interesting to compare the time evolution of both types of coefficients. Fig. 8 shows the time variation of $G_0^0$ and $G_{-1/2}^0$ and their time derivative. The time derivative of $G_{-1/2}^0$, for instance, has points of inflexion occurring at the times when geomagnetic jerks are recorded at the Chambon-La-Forêt observatory in 1969, 1981, 1991 and 2000. The occurrence for these events are not exact in the R-SCHA2D secular variation model as a result of the poor
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The template of Fig. 4 is subtracted at each repeat station for epochs between 1965 and 2002.5. The scatter at each repeat station before (in black) and after (in grey) the correction is shown for X, Y and Z. The correction is particularly efficient on the X component.</figcaption></figure>

temporal resolution of repeat station measurements. The analysis of $G_0^0$ shows that the vertical components varied from 40 900 nT in 1965 to about 41 700 nT in 2007 in average over France. The variation of Z per year ($20\,\text{nT}\,\text{yr}^{-1}$ in average) is slightly lower than for $G_{0\,1/2}^0$. The temporal derivative of $G_0^0$ also seems to 'miss' the 1991 geomagnetic jerk confirming that Z is less sensitive to geomagnetic jerks. The other parameter variations are also relatively smooth with time but not linear. Changes in slope appear occasionally that show that the geomagnetic jerks are detected in the regional parameters themselves. This kind of result is difficult to reach when modelling the magnetic field at a global scale because geomagnetic jerks are usually not worldwide in occurrence (Chambodut & Mandea 2005). The magnetic power of a jerk may be strong over a specific region but weak at the global scale and the overall error on the coefficients is usually too high to isolate the real statistical variations of a spherical harmonic coefficient (Olsen & Mandea 2007). Chambodut & Mandea (2005) showed that the dates of jerk occurrence are relatively homogeneous over the entire European continent. This is why clear spikes are observed in the regional parameters derivatives over a smaller region like France.

The secular variation obtained directly from the parameter derivatives also suffers from errors near the edges (in time and space). Some of the problem are identified as coming from the modelling error as a result of the low spatial and series expansion. The absolute error estimated in the case of perfect synthetic data is $2\,\text{nT}\,\text{yr}^{-1}$ for the secular variation. The expected relative error thus does not exceed 10 per cent for the three components for a model interpolated on a regular grid. Figs 9 and 10 show the regional secular variation model at various epochs. When the difference between R-SCHA2D and CM4 is less than $2\,\text{nT}\,\text{yr}^{-1}$, the contour lines of CM4 are not labelled. In general, the R-SCHA2D model is consistent with CM4 but the difference is slightly increasing from 1990. R-SCH2D and CM4 secular variation have sometimes different magnitudes as the regional model is smoother in time than the CM4 model. In 2005.5, the comparison is done with the POMME 4 model, an update of POMME 3 (Maus *et al.* 2006) that spans epochs from 2000 to 2007.5. From 2005.5, the model seems to be less stable than for previous epochs because of the two mission repeat station measurements (GAR and DEL) and, again, because other sources of noise were not corrected for at this epoch. The power series expansion may also cause this behaviour near the edges of the time range 1965 2007.5.

For the whole time span, the secular variation model is calculated at the Chambon-La-Forêt observatory (Fig. 11). The secular variation model is directly compared with CLF annual means (not corrected for the 11 yr solar cycle) and the CM4 secular variation. The R-SCHA2D secular variation model is in good agreement even though it is evident from Fig. 11 that R-SCHA2D suffers from a lack of temporal resolution as a result of the average 5 yr interval between two repeat station measurements. This poor temporal resolution explains the differences seen in Figs 9 and 10 between R-SCHA2D and CM4. It is worth noting that R-SCHA2D secular variation model has a possible instability near the epoch 2007 for Z. The change of slope is too steep after 2005 to be realistic and, in the future, cubic splines will be used for the time basis functions. Note also that fast secular acceleration occurring within 5 yr intervals could induce a shift in the R-SCHA2D modelling, or a wrong localization of jerks, due to the poor 5 yr temporal resolution or the measurements. For this reason, the 2000 geomagnetic jerk is perhaps not really seen in the R-SCHA2D model (the acceleration does not have the correct sign in Fig. 11 for Y: while it is a local maximum in R-SCHA2D, it is a local minimum in the CLF annual means) but the change of slope probably results from the average of the 1996 and the 2003 geomagnetic impulses clearly visible in the CLF observatory data. This shift is also visible in the parameters (Fig. 8). To the contrary,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** R-SCHA2D main field model for epochs 1970 (top panel), 1980 (middle panel) and 1990 (bottom panel) in solid line. The CM4 model is sketched in dashed line. The *X*, *Y* and *Z* R-SCHA2D components are extrapolated on the regular grid (dots on the top left-hand side) and the declination is calculated from *X* and *Y*.</figcaption></figure>

the 2003 jerk, which is not worldwide in occurrence but mostly located in a rather limited area near 90°E (Olsen & Mandea 2007) and so weak over France, is correctly located in the curve showing $\partial Z / \partial t$ because there is no fast variation neither before nor after this event on this component (and obviously because the nearest repeat station measurements were done in 2002.5).

# 4 CONCLUSION

R-SCHA2D is an approach designed to work with data recorded at nearly constant altitudes although its behaviour with altitude was not thoroughly tested. So far, it seems that the reliability of the upward/downward continuation of R-SCHA2D models depends on their spectral content; upward/downward continuation being more reliable when dealing with smaller wavelengths than the cap's size. The R-SCHA2D is the solution of an incomplete boundary value problem. This comes from the fact that potential theory requires a closed surface to uniquely determine the solutions of a boundary value-problem. However, as explained in the first section of this paper, and in Appendix A, the chosen basis functions are arguably the most convenient for modelling the main geomagnetic field and its secular variation. This approach suffers from
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Fig. 6 for epochs 2000, 2005 and 2007.5.</figcaption></figure>

some of the drawbacks existing in SCHA such as the theoretical inability to separate internal and external fields and poor downward/upward continuation outside the altitude of the dominion of data. The first problem seems inherent to any regional modelling (Thébault *et al.* 2006b) and the second is imposed by the assumption of data availability at one surface only. Knowing these limitations, the main advantages are that the total set of basis functions is complete and that we were able to represent the main geomagnetic field and its secular variation at a regional scale. Taking all care to generalize the results of this paper (because regional models depend on the spherical cap dimensions), it thus means that any magnetic field recorded at the Earth's surface in a region can be modelled with R-SCHA2D without apparent restriction. In addition, the resolution matrix of the inverse problem is block-diagonal but the non-orthogonality terms are rather weak. This makes R-SCHA2D less sensitive to data gaps and noise and it is possible to represent large magnetic field structures with few parameters without regularization. These assertions were verified numerically using a set of noisy synthetic vector data equally and sparsely distributed calculated with a spherical harmonic model from degree 1 to 90. All the more, some regional parameters have a physical meaning, either related to the average value of the vertical magnetic field, the coefficient $G_0^0$ for instance, or more related to the horizontal geomagnetic field, like $G_{-1/2}^0$ (that is the average of the $\tilde{X}$ component along the closed spherical cap contour). This offers new perspectives towards the elaboration of a suitable spectral analysis from the regional parameters, an important tool difficult to establish at a regional scale as was shown by Haines (1991).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Two regional parameters and their temporal derivatives are sketched between epochs 1965 and 2008. The derivatives of $G_0^0$ show the 1970, 1981 and 2000 geomagnetic jerks (right panel). The temporal derivative of $G_{-1/2}^0$ also detects the 1991 geomagnetic jerk.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** R-SCHA2D secular variation model for epochs 1967.5, 1975 and 1985. Solid lines are R-SCHA2D SV model and dashed lines are CM4 SV model. Contour line spacing is 2 nT yr⁻¹.</figcaption></figure>

R-SCHA2D was applied to repeat station data over France between epochs 1965 and 2007.5 with a minimum of data pre-processing. The regional parameters were found to smoothly vary between epochs. In particular, it was possible to derive the secular variation directly from the model coefficients and to identify the geomagnetic jerks that were recorded over France from 1965 to 2007.5 in the parameters themselves. The importance of data pre-processing for the derivation of secular variation model came out from this first analysis because systematic noise related to the 11 yr solar cycle activity overlapped with the geomagnetic jerk occurrence. This was due to repeat data reduction to annual mean observatory. The regional model also showed that it was possible to capture the large structures of the lithospheric field from repeat station measurements. This result is two-sided, it is satisfactory because it proves the ability of R-SCHA2D to highlight the smallest possible structures but, at the same time, it demonstrates how difficult it is to disentangle the magnetic signals at regional scales. These small structures also have the ability to spoil the model time derivative and create oscillations in regions where no data were available. With these
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Fig. 9 but for epochs 1995, 2000 and 2005.5. When dashed lines are labelled, it means that the difference between CM4 and R-SCHA2D is greater than 2 nT yr⁻¹. Contour line spacing is 4 nT yr⁻¹ for X and 2 nT yr⁻¹ for Y and Z.</figcaption></figure>

very simple corrections, the definitive regional model and its secular variation are very consistent with the comprehensive CM4 model even though the temporal resolution is coarser. Another approach, more efficient for representing the secular variation, will consist in performing a fit by taking the first central difference between two successive measurements. This so-called 'time-derivative' fit removes the crustal field contamination and thus stabilizes the secular variation model. It was not presented here because the purpose of this work was to test the efficiency of R-SCHA2D in a general case. More sophisticated data processing and inverse modelling method may be envisaged.

An aside of this work was to investigate for possible secular variation anomaly over France. The results are not conclusive as it is not possible to find any convincing temporal regional variations. However, there is an increasing difference between CM4 and the R-SCHA2D regional model between 1990 and 2007.5 that is tenuous but larger than the estimated modelling error. The source of this difference remains to be investigated by a 'time-derivative' fit and requires a specific strategy that exceeds the scope of this work.

It was possible to derive a spatiotemporal magnetic field model over France with a method obeying Laplace equation independently from global models but the regional modelling technique did not bring improvements over a spherical harmonic approach over France for the study of the main magnetic field. In fact, the comprehensive model CM4 is doing a good job in that areas between 1960 and 2002.5. It is important to stress that the CM4 is well constrained over Europe but it is reckoned to be less reliable, particularly in the Southern hemisphere or isolated area, like in Antarctica (De Santis et al. 2002) or South Africa (Kotzé 2002). Moreover, global models have non-negligible crustal bias that can be corrected for at a regional scale considering jointly aeromagnetic and satellite data over a given region (Thébault et al. 2006a). The same lithospheric field model can then be added back to R-SCHA2D in order to derive national magnetic field models updated for navigational purposes. The method is in its infancy but it is reasonable to foresee various applications involving any kind of spatiotemporal magnetic field variations.

## ACKNOWLEDGMENTS

I would like to record my gratitude to the team of the Observatoire de Chambon-La-Forêt, especially to Arnaud Chulliat, Kader Telali and François Truong, who provided me with the provisional repeat station data. The result of Appendix A, concerning the limit at infinity of
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The R-SCHA2D secular variation model is calculated at the Chambon-la-Forêt observatory between 1965 and 2008 (pink line). The blue line is the CM4 secular variation model and the dashed line is the secular variation derived from CLF annual means.</figcaption></figure>

R-SCHA, benefited in 2004 from discussions with J.P Hoffbeck, a graduate student from ENS, and with Jean-Jacques Schott from EOST. J.M Torta and an anonymous reviewer are acknowledged for their constructive comments. This is the IPGP contribution number 2344.

## REFERENCES

* Abramowitz, M. & Stegun, A., 1965. *Handbook of Mathematical Functions*, Dover, Mineola, NY.

* Bullard, E.C., 1967. The removal of trend from magnetic surveys, *Earth planet. Sci. Lett.*, **2**, 293–300.

* Chambodut, A. & Mandea, M., 2005. Evidence for geomagnetic jerks in comprehensive models, *Earth Planets Space*, **57**(2), 139–149.

* De Santis, A. & Falcone, C., 1995. Spherical cap models of Laplacian potentials and general fields, in *Geodetic Theory Today*, pp. 141–150, Springer, New York.
De Santis, A., Battelli, O. & Kerridge, D.J., 1990. Spherical cap harmonic analysis applied to regional field modelling for Italy, *J. Geomag. Geoelectr.*, **42**, 1019–1036.

De Santis, A., Torta, J.M. & Lowes, F.J., 1999. Spherical cap harmonics revisited and their relationship to ordinary spherical harmonics, *Phys. Chem. Earth*, **24**(11–12), 935–941.

De Santis, A., Torta, J.M. & Gaya-Piqué, L.R., 2002. The first Antarctic geomagnetic Reference Model (ARM), *Geophys. Res. Lett.*, **29**(8), 1192, doi:10.1029/2002GL014675.

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

Haines, G.V., 1985b. Magsat vertical field anomalies above 40° spherical cap harmonic analysis, *J. geophys. Res.*, **90**, 2593–2598.

Haines, G.V., 1991. Power spectra of sub-periodic functions, *Phys. Earth planet. Inter.*, **65**(3), 231–247.

Haines, G.V., 1993. Modelling geomagnetic secular variation by main-field differences, *Geophys. J. Int.*, **114**, 490–500.

Haines, G.V. & Torta, J.M., 1994. Determination of equivalent current sources from spherical cap harmonic models of geomagnetic field variations, *Geophys. J. Int.*, **118**, 499–514.

Hobson, E.W., 1965. *The Theory of Spherical and Ellipsoidal Harmonics*, 2nd edn, Chelsea, New York.

Kim, H.R., von Frese, R.R.B., Golynsky, A.V., Taylor, P.T. & Kim, J.W., 2004. Application of satellite magnetic observations for estimating near-surface magnetic anomalies, *Earth, Planets Space*, **56**, 955–966.

Korte, M. & Holme, R., 2003. Regularisation of spherical cap harmonics, *Geophys. J. Int.*, **153**(1), 253–262.

Korte, M. & Thébault, E., 2007. Geomagnetic repeat station crustal biases and vectorial anomaly maps for Germany, *Geophys. J. Int.*, doi:10.1111/j.1365-246X.2007.03387.x

Kotzé, P.B., 2002. Modelling and analysis of Oersted total field data over South Africa, *Geophys. Res. Lett.*, **29**(15), doi:10.1029/2001GL013868.

Langel, R.A. & Hinze, W.J., 1998. *The Magnetic Field of the Earth's Lithosphere*, Cambridge Univ. Press, New York.

Maus, S., Rother, M., Stolle, C., Mai, W., Choi, S., Lühr, H., Cooke, D. & Roth, C., 2006, Third generation of the Potsdam Magnetic Model of the Earth (POMME), *Geochem. Geophys. Geosyst.*, **7**, Q07008, doi:10.1029/2006GC001269.

Olsen, N. & Mandea, M., 2007. Investigation of a secular variation impulse using satellite data: the 2003 geomagnetic jerk, *Earth planet. Sci. Lett.*, **255**, 94–105.

Olsen, N., Lühr, H., Sabaka, T.J., Mandea, M., Rother, M., Tøffner-Clausen,

L. & Choi, S., 2006. CHAOS—a model of the Earth’s magnetic field derived from CHAMP, Ørsted, and SAC-C magnetic satellite data, *Geophys. J. Int.*, **166**(1), 67–75, doi:10.1111/j.1365-246X.2006.02959.x

Pavón-Carrasco, F.J., Osete, M.L., Torta, J.M., Gaya-Piqué, L.R. & Lanos, Ph., 2007. Initial SCHA.DI.00 regional Archaeomagnetic model for Europe for the last 2000 years, *Phys. Chem. Earth*, in press.

Robin, L., 1957–59, *Fonctions sphériques de Legendre et fonctions sphéroïdales*, 3 Vols., Gauthier-Villars, Paris.

Sabaka, T., Olsen, N. & Purucker, M., 2004. Extending comprehensive models of the Earth’s magnetic field with Ørsted and CHAMP data, *Geophys. J. Int.*, **159**, 521–547.

Thébault, E. & Gaya-Piqué, L., 2008. Applied comparisons between SCHA and R-SCHA regional modelling techniques, *Geochem. Geophys. Geosyst.*, in press.

Thébault, E. & Mandea, M., 2006. Contribution of regional modelling techniques to the Swarm mission, ESA special publication 615, 4 pages.

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

Thébault, E., Mandea, M. & Schott, J.J., 2006a. modelling the lithospheric magnetic field over France by means of revised spherical cap harmonic analysis (R-SCHA), *J. geophys. Res.*, **111**, B05102, doi:10.1029/2005JB004110.

Thébault, E., Schott, J.J. & Mandea, M., 2006b. Revised spherical cap harmonic analysis (R-SCHA): validation and properties, *J. geophys. Res.*, **111**, B01102, doi:10.1029/2005JB003836.

Torta, J.M., Garcia, A., Curto, J.J. & De, Santis A., 1992. New representation of geomagnetic secular variation over restricted regions by means of spherical cap harmonic analysis: application to the case of Spain, *Phys. Earth. planet. Inter.*, **74**, 209–217.

Torta, J.M., Curto, J.J. & Bencze, P., 1997. Behavior of the quiet day ionospheric current system in the European region, *J. geophys. Res.*, **102**(A2), 2483–2494.

Torta, J.M., Gaya-Piqué, L.R. & De Santis, A., 2006. Spherical cap harmonic analysis of the geomagnetic field with application for aeronautical mapping, in *Geomagnetics for Aeronautical Safety: A Case Study in and Around the Balkans, NATO Security Through Science Series—C*, pp. 291–307, eds Rasson, J.L. & Delipetrov, T., Springer, Dordrecht.

Tyler, R., Maus, S. & Lühr, H., 2003. Satellite observations of magnetic fields due to ocean tidal flow, *Science*, **299**, 239–241.

Verbanac, G., Lühr, H., Rother, M., Korte, M. & Mandea, M., 2007. Contributions of the external field to the observatory annual means and a proposal for their corrections, *Earth Planets Space*, **59**, 251–257.

# APPENDIX A

This appendix summarizes the difficulty of using R-SCHA (Thébault *et al.* 2006b) or applying boundary conditions at infinity when data are available at one surface only.

## A.1 The finite but very thin cone

When considering a finite cone, an extra boundary condition is applied at an upper surface $\partial_b \Omega = \{b > a, 0 \le \theta \le \theta_0, 0 \le \varphi \le 2\pi\}$ such that $(\frac{\partial V_2}{\partial r})_{r=b} = 0$. This implies that the solution of equation (4b) is (Thébault *et al.* 2006a):

$$ \Theta \equiv K_p^m = \frac{(-1)^m}{2^m m!} N_p^m \sin^m \theta F \left( m + \frac{1}{2} - i \frac{p\pi}{\log \frac{b}{a}}, m + \frac{1}{2} + i \frac{p\pi}{\log \frac{b}{a}}, m + 1; \sin^2 \frac{\theta}{2} \right). \quad (A1) $$

It is clear from (A1) that $\lim_{b \to a} K_p^m (\cos \theta) = \infty$ and this solution is no more defined. It is therefore, not convenient to use R-SCHA to model data available at one surface or within a thin lens. This justifies the changes presented in this paper.

## A.2 The infinite cone

Let us now assume that we want to describe the magnetic potential within an semi infinite cone in order to allow our models to be upward or downward continued correctly. This boundary value problem writes:

$$ \nabla^2 (V_1 + V_2) = 0 \quad (A2) $$
$$ \left. \frac{1}{r} \frac{\partial V_2}{\partial \theta} \right|_{\partial \theta_0 \Omega} = \tilde{X}(r, \theta_0, \varphi) \qquad (A3) $$

$$ \left. \frac{\partial V_1}{\partial r} \right|_{\partial a \Omega} = \tilde{Z}(a, \theta, \varphi) \qquad (A4) $$

$$ \lim_{r \to \infty} \frac{\partial V_2}{\partial r} = 0 \qquad (A5) $$

$$ \lim_{r \to \infty} V_2 = 0. \qquad (A6) $$

The domain of interest, the semi-infinite cone, is regular at infinity and the boundary value problem is complete. The extra set of conditions (A5) and (A6) impose that {$G_{n_k}^{m,e}, H_{n_k}^{m,e}$} = 0 in (6). However, this also implies that $R_2(r) = A \left(\frac{r}{a}\right)^{n_2} + B \left(\frac{a}{r}\right)^{n_2+1}$ with $\lambda_2 < -1/4$. As a result, $n_2 = -1/2 + iy$ with $y \in \mathbb{R}^+$. The solution within a semi-infinite cone is the limit of the finite cone when the upper surface goes to infinity. It directly follows that $V_2$ writes:

$$ V_2(r, \theta, \varphi, t) = a \sum_{m \ge 0} \int_{y \in \mathbb{R}^+} R_y(r) [G_y^m(t) \cos(m\varphi) + H_y^m(t) \sin(m\varphi)] P_y^m(\theta). $$

In practice, there is thus an indefinite number of parameters {$G_y^m, H_y^m$} to be estimated by an inverse problem. This is why no conditions were set at infinity in eqs (15). Note that eq. (14) is an incomplete solution of the semi-infinite cone boundary value problem. Note also that only in the case of the infinite cone can the coefficients {$G_{n_k}^{m,e}, H_{n_k}^{m,e}$} in (6) be set to zero. The solution of Laplace equation with boundary conditions at one surface, within a closed and finite domain or a semi-infinite domain are thus different. For this reason, R-SCHA2D was developed to deal with data available at one surface only.

# APPENDIX B

Let us write the norm of the regional magnetic field over the spherical cap domain:

$$ \| \nabla (V_1 + V_2) \|_{\partial \Omega_a}^2 = \int_{\partial a \Omega} \nabla V_1 \cdot \nabla V_1 d \partial \Omega_a + \int_{\partial a \Omega} \nabla V_2 \cdot \nabla V_2 d \partial \Omega_a + 2 \int_{\partial a \Omega} \nabla V_1 \cdot \nabla V_2 d \partial \Omega_a \qquad (B1) $$

with $\partial_a \Omega = \{r=a, 0 \le \theta \le \theta_0, 0 \le \varphi \le 2\pi\}$. Let us focus on the cross-product between $\nabla V_1$ and $\nabla V_2$. The gradient operator may be split in its horizontal and radial part, $\nabla = \nabla_H + \nabla_r$ such that:

$$ \int_{\partial \Omega_a} \nabla V_1 \cdot \nabla V_2 d \partial_a \Omega = \int_{\partial \Omega_a} \nabla_H V_1 \cdot \nabla_H V_2 d \partial_a \Omega + \int_{\partial \Omega_a} \nabla_r V_1 \cdot \nabla_r V_2 d \partial_a \Omega \qquad (B2) $$

According to (9) the second integral on the right hand-side is zero. The first Green's identity extended to the surface case writes:

$$ \int_{\partial \Omega_a} \nabla_H V_1 \cdot \nabla_H V_2 d \partial_a \Omega = \int_0^{2\pi} V_1 \theta_0 \frac{\partial V_2}{\partial \theta} \bigg|_{\theta_0} \sin \theta_0 d\varphi - \int_0^{2\pi} \int_0^{\theta_0} V_2 \nabla_H^2 V_1 d \partial_a \Omega. \qquad (B3) $$

In addition, $\nabla_H^2$ also writes at the Earth's surface $\nabla_H^2 = -L^2/a^2$, with $L^2$ the angular momentum operator with the property that $L^2 V_2 = \lambda_2 V_2$. The cross-product between $\nabla V_1$ and $\nabla V_2$ is:

$$ \int_{\partial \Omega_a} \nabla V_1 \cdot \nabla V_2 d \partial_a \Omega = \int_0^{2\pi} V_1 \theta_0 \frac{\partial V_2}{\partial \theta} \bigg|_{\theta_0} \sin \theta_0 d\varphi + \lambda_2 \int_0^{2\pi} \int_0^{\theta_0} V_2 V_1 \sin \theta d\theta d\varphi \qquad (B4) $$

which gives:

$$ \int_{\partial \Omega_a} \nabla V_1 \cdot \nabla V_2 d \partial_a \Omega = \sum_{k \ge m \ge 0} \sum_{m \ge 0} C_{k, -1/2}^m \left\{ \sin \theta_0 P_{n_k}^m(\theta_0) \frac{d P_{-1/2}^m(\theta_0)}{d\theta} + \lambda_2 \int_0^{\theta_0} P_{n_k}^m P_{-1/2}^m \sin \theta d\theta \right\} \qquad (B5) $$

with $C_{k, -1/2}^m = (1 + \delta_{m,0})\pi a^2 \{(G_{n_k}^{m,i} + G_{n_k}^{m,e})G_{-1/2}^m + (H_{n_k}^{m,i} + H_{n_k}^{m,e})H_{-1/2}^m\}$, the product between regional parameters. Furthermore, it can be shown that:

$$ \int_0^{\theta_0} P_{n_k}^m P_{-1/2}^m \sin \theta d\theta = \frac{\sin \theta_0}{\lambda_1 - \lambda_2} \left[ P_{n_k}^m \frac{d P_{-1/2}^m}{d\theta} - P_{-1/2}^m \frac{d P_{n_k}^m}{d\theta} \right]_0^{\theta_0} $$

with $\lambda_1 = n_k(n_k + 1)$ and $\lambda_2 = -1/4$ as defined in Section 2. In order to have orthogonality, the term into brackets in (B5) must vanish. This happens only if $\lambda_2 = 0$ and $P_{n_k}^m(\theta_0) = 0$. This case was discarded in Section 2. In any other case, the cross-product is non-zero. It depends on the ratio between $\lambda_1$ and $\lambda_2$ but also on the value of the associated Legendre and Mehler functions near the edges.