Bouligand

A Khokhlov

G Hulot

‡ C Bouligand

Testing statistical paleomagnetic field models against directional data affected by measurement errors

Testing statistical palaeomagnetic field models against directional data affected by measurement errors

Introduction

Thanks to the large amount of magnetic data provided by satellite missions spherical harmonic models of the main magnetic field produced within the core can readily be computed, providing high resolution pictures of the way this field has been behaving over the past few decades [START_REF] Hulot | Smallscale structure of the geodynamo inferred from Oersted and Magsat satellite data[END_REF]. Additional information can also be recovered from observations carried out by many generations of observers, explorers and navigators worldwide to produce spherical harmonic models describing the main magnetic field over the past four centuries [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF]. Reconstructing similar, albeit much less accurate, spherical harmonic models of the Earth's main magnetic field further back in time is also possible [START_REF] Hongre | An analysis of the geomagnetic field over the past 2000 years[END_REF], Korte & Constable 2005]. This however requires that indirect measurements be used, only available through human artefacts, lava flows and sediments that have been magnetized in the ancient field. It also requires that a good age control of each sample is available to ensure a satisfying synchronization of the data used in computing spherical harmonic models for a given epoch.

Because changes in the non-dipole component of the main field occur on time scales all the shorter than the spatial scales considered are smaller [Hulot & Le Mouël 1994], this unfortunately limits our ability to produce such spherical harmonic models back in time (so far, only up to 7000 years in the past, [Korte & Constable 2005]).

To recover information about the Earth's main magnetic field further back in time, one must rely on paleomagnetic data and acknowledge the impossibility of accurately (within a few decades at most) synchronizing data acquired at different locations. This is especially true for the so-called paleosecular variation [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF], Johnson & Constable 1996, McElhinny & McFadden 1997] and paleointensity [START_REF] Tanaka | Some global features of paleointensity in geological time[END_REF][START_REF] Perrin | [END_REF] databases covering the past few million years.

Those databases encompass data recovered from volcanic samples and testifying for the instantaneous value of the direction and/or intensity of the field at well-known locations but poorly known times (typically within much more than a millennium). Fortunately, ages can nevertheless be measured with enough accuracy to identify the chron during which each of these samples acquired its magnetization.

Assuming that the geodynamo essentially remained in a stationary state at times of stable polarity, this then opens the possibility, recognized long ago (see e.g. [START_REF] Merrill | The magnetic field of the Earth[END_REF] for an extensive review of earlier work), that at least some statistical properties could be recovered, characterizing the ancient field at times of stable polarity over the past millions years.

Several approaches have been used to try and infer such statistical properties (again, see e.g. [START_REF] Merrill | The magnetic field of the Earth[END_REF]). One approach has drawn particular attention, the so-called Giant Gaussian Process (GGP) approach first introduced by [Constable & Parker 1988], and next generalized by [Hulot & Le Mouël 1994] and [Kono & Tanaka 1995]. In its most general form (see e.g., [Hulot & Bouligand 2005]), this approach consists in 1) using the classical spherical harmonic decomposition of the Main Field:

B( , Θ, Ψ, t) = -∇ a ∞ l=1 l m=0 a l+1 g m l (t)cos mΨ + h m l (t)sin mΨ P m l (cos Θ) (1)
where g m l (t) and h m l (t) are the time-varying Gauss coefficients, a denotes Earth's radius and { , Θ, Ψ} are the standard spherical coordinates (i.e., distance from the Earth's center, colatitude and longitude), and 2) assuming that the g m l (t) and h m l (t) are the coordinates of a vector k(t) in a multidimensional model space, which behaves as a single realization of a multidimensional stationary random Gaussian process characterized by a statistical mean (or mean model) E(k(t)) and a covariance matrix Cov(k(t), k(t )) (which may very well be non-diagonal).

In addition to being defined with the help of a fully consistent formalism, this GGP approach has the unique advantage of providing a common tool to analyse just as well the historical ( [Constable & Parker 1988, Hulot & Le Mouël 1994]), archeomagnetic [START_REF] Hongre | An analysis of the geomagnetic field over the past 2000 years[END_REF]) and paleomagnetic fields ( [Constable & Parker 1988] and many studies since, see e.g. [Kono et al 2000a, Tauxe & Kent 2004]). Even more, it can also be used to carry on similar analysis of the field produced by dynamo numerical simulations ( [START_REF] Kono | Dynamo simulations and paleosecular variation models[END_REF], McMillan et al 2001, Bouligand et al 2005]). In the latter case, considerable advantage can be taken of the fact that the Gauss coefficients (and very long time series) are readily available and can thus be analyzed in much detail without having to worry about any observational errors. Those studies have shown that under reasonable circumstances, the field produced by a numerical dynamo is indeed compatible with a GGP description, which can then be used to detect very interesting statistical properties, most notably, spontaneous and forced (because of the boundary conditions imposed at the core-mantle boundary) symmetry breaking properties ( [START_REF] Bouligand | Statistical paleomagnetic field modelling and dynamo numerical simulation[END_REF], Hulot & Bouligand 2005]). The GGP approach thus appears to potentially be a tool of choice to characterize the Earth's main magnetic field at various epochs (present and past), detect possible changes in the regime of the geodynamo, some of which could be due to changes in the boundary conditions, and decide whether a numerical dynamo simulation is "Earth-like" or not.

For this to be possible, we however need to have robust and sensitive numerical tools to actually decide which GGP model best fits a given data set. The issue is unfortunately not so trivial when the data to be analyzed are the paleomagnetic data from the paleosecular variation databases (the largest available paleomagnetic database for that type of study). In that case indeed, because the (directional) data are non-linearly related to the Gauss coefficients of the main field, no simple method is available to directly produce estimates of the mean and covariance matrix of the GGP process best describing the data. Furthermore, using approximate approaches can easily lead to biased estimates ( [START_REF] Khokhlov | Towards a self-consistent approach to paleomagnetic field modelling[END_REF], Hatakeyama & Kono 2001]). The purpose of the present paper, which expends on the earlier work of [START_REF] Khokhlov | Towards a self-consistent approach to paleomagnetic field modelling[END_REF] (hereafter paper I), is to show that in fact a rigorous method can be implemented which 1) is fully consistent with all the assumptions involved in the GGP approach, 2) takes into account the measurement errors provided with each data point (a major improvement with respect to paper I), 3) provides a rigorous quantified assessment of the compatibility of any GGP model with any directional data set, and 4) proves very discriminating.

Dealing with error-free paleodirectional data

In this section, we first recall the approach proposed in paper I to test any error-free paleodirectional data set. [Note however that we will make use of slightly different notations to avoid confusion among the various quantities we need to define. As a general rule, in the present paper, we will indeed always make use of gothic letters, namely g, s, k and p (possibly with indexes) for the various probability density functions we will need to introduce].

When dealing with such data, which consist in instantaneous values of the paleodirection recorded at times separated by much more than the memory of the GGP process, each paleomagnetic datum can be viewed as a local (both in time and space) independent realization of the GGP pro-cess. As a result, the covariance matrix of the process can be assumed to be of the simplified form Cov(k(t), k(t )) = Cov(k, k)δ(t -t ) and temporal correlations ignored altogether. Rewriting (1) in the simplified form B( , Θ, Ψ, t) = A( , Θ, Ψ)k(t), where A( , Θ, Ψ) is a matrix projecting the multidimensional vector k(t) into the 3D vector B( , Θ, Ψ, t), a given instantaneous value B of the field at a given location ( , Θ, Ψ) at the Earth's surface can then itself be viewed as a random drawing from a 3D Gaussian distribution of random vectors V characterized by E(V ) = A( , Θ, Ψ)E(k) and is then defined by:

Cov(V , V ) = A( , Θ, Ψ)Cov(k, k)A( , Θ, Ψ) T .
s(u) = ∞ 0 g(ρu)dρ. (2)
In a local Cartesian coordinate system, we may write

E(V ) = m = (m 1 , m 2 , m 3 ) and Cov(V , V ) = [cov(V i , V j )]. Then, let Λ = [Λ ij
] be the inverse (hence, also symmetric) matrix of Cov(V , V ). With respect to this local Cartesian coordinates, the probability density function of V is:

g(x 1 , x 2 , x 3 ) = det Λ (2π) 3 exp   - 1 2 3 i,j=1 Λ ij (x i -m i )(x j -m j )   , (3)
or

g(x) = det Λ (2π) 3 exp - 1 2 Λ(x -m), x -m . ( 4
)
Making use of the Λ-inner product (x, y) Λ = (Λx, y) = 3 i,j=1 Λ ij x i y j together with the corresponding Λ-norm |x| Λ = (x, x) Λ , and turning back to the spherical coordinates (u, ρ) in R 3 , we may then write:

g(ρu) = det Λ (2π) 3 e -1 2 |ρu-m| 2 Λ ρ 2 , (5)
and, after ρ-integration (2), finally get

s(u) = e -1 2 m 2 • √ det Λ 4π|u| 3 Λ z 2 π + e 1 2 z 2 (1 + z 2 ) 1 + Erf z √ 2 , (6)
where

z = (m, u) Λ |u| Λ , m = |m| Λ , (7)
correspond to respectively the Λ-projection of m on the direction u, and the Λ-norm of m. Since Λ is positive -m z m.

Equation ( 6) makes it possible to predict the pdf of the direction of the field at any location at the Earth's surface. Any given GGP model can then easily be tested against the data set corresponding to such a location.

Unfortunately, relatively few data are usually available at a given location. In paper I, we explained how this drawback could be overcome, and sparse data distributed over various sites combined into a uniformized data set to globally test the GGP model.

This consists in converting the directional data set {u i } under consideration (where the "i" index varies between 1 and N, the total amount of data) into a univariate data set {t i } with the help of:

t i = P {u|s i (u) s i (u i )} = {u|si(u) si(ui)} s i (u)dU ( 8
)
where dU is the elementary surface about u on Σ, and testing the distribution of the {t i } against a uniform probability density over the segment [0,1].

This uniformization procedure simply amounts, for each data u i , to identify on Σ the iso-probability line of the pdf s i (u) on which the data u i lie (note that this pdf depends on the location the data u i come from, which is why an"i" index is being added to s i (u)), and to assign to t i the value of the probability that u i could have been lying on a higher iso-probability line on Σ. Then, if the GGP model is compatible with the data set (hence, if the u i are compatible with the s i (u)), the {t i } should be compatible with a uniform distribution on [0, 1].

Whether this is indeed the case can finally be assessed with the help of various standard tests, (in paper I we used the Kolmogorov-Smirnov (KS) and the χ 2 tests). Note that this uniformization procedure is very general and does not rely on any properties of the Angular Gaussian distribution. It can be applied in more general situations.

Taking measurement error into account

A main issue that paper I did not address is that of data measurement errors. The tests we just described are appropriate only in the event these errors are negligible. Unfortunately this is far from being the case and we now need to take this into account. In the case of directional data, these measurement errors are classically described in terms of a Fisher distribution characterized by a socalled concentration (or precision) parameter K. Within this Fisherian frame-work, stating that the observed direction s reflects the actual direction u on Σ with an error characterized by a concentration parameter K, amounts to state that s is the result of a random drawing from a Fisher distribution centered on u with pdf:

k K (u, s) = K 2π e K -e -K e K•(u,s) . (9)
In spherical coordinates and when centered on u 0 = (θ = 0, ϕ = 0), this takes the more usual form:

k K (u 0 , s) = Ke K cos θ 2π e K -e -K . ( 10
)
This means that in order to test a given GGP model against a given data set {u i } with associated errors characterized by {K i }, we now need to test whether each u i can be considered the result of a realization s of the GGP process (distributed over Σ with the pdf s i ), shifted to u i as a result of an error with a Fisher distribution k Ki (s, u) centered on s and a concentration parameter K i .

The infinitesimal probability that the GGP process first produces a direction within the elementary surface dS about s, and that the error next independently shifts this direction to within the elementary surface dU about u i is:

s i (s)dS • k Ki (s, u)dU. (11)
Then, integrating over all possible intermediate s directions leads to the probability p i (u)dU of finding u i within dU about u, where:

p i (u) = Σ s i (s)k Ki (s, u)dS (12)
is the new pdf against which u i should be tested. Thus, taking data errors into account only amounts to use p i (u) in place of s i (u) in Eq.( 8). Tests can then again be used, either on a site-by-site basis or, more interestingly, on a regional and global scale after using the uniformization procedure, which now becomes:

t i = P {u|p i (u) p i (u i )} = {u|pi(u) pi(ui)} p i (u)dU (13)

Practical implementation

In practice, testing a given GGP model against a given directional data set thus involves four successive steps: 1) for each data u i , to compute the error-free pdf s i (u) predicted by the GGP model at the site where u i was collected (the analytical form of which is given by ( 6)); 2) to compute the error-included pdf p i (u) through the convolution (12); 3) to produce the uniformized data t i with the help of (13); 4) to test the uniformized data set {t i } against a uniform distribution over [0, 1].

Unfortunately, no exact analytical solutions of ( 12) and ( 13) are known to us. In principle, this is not too much of a problem, since both ( 12) and ( 13) can be computed numerically. In practice however, the numerical implementation of formula (12) requires quite some computational time. This drawback can be considered negligible if we simply test a single GGP model against a small database, or if all data u i coming from the same site share the same data error (in which case p i (u) only needs to be computed once for each of the typically thirty something sites). But if p i (u) needs to be computed for each data of the data set (i.e. typically a thousand times), this computational time can become very long (on the order of a few days on a current PC). For that kind of more realistic situations (and for the "forward" testing method we propose here to possibly become of any use for future much more computer-intensive "inverse" searches of best models), faster algorithms are clearly desirable. This prompted us to further look for an approximate, accurate enough, analytic solution to the convolution (12).

As is described in the appendix, one such approximation can indeed be found. For any given data u i , this consists in approximating p i (u) by the Angular Gaussian distribution defined by ( 6) and ( 7), where m is then the mean field value predicted by the GGP model to be tested at the site associated with the data u i (i.e., predicted by g i (x) as defined by ( 4) for the site associated with the data u i ),

and

Λ is Cov(V , V ) + |m| 2 K I -1
, where Cov(V , V ) is the covariance matrix associated with g i (x) and K(= K i ) is the concentration parameter characterizing the Fisherian error associated with u i .

In what follows, whenever we use this approximation of p i (u) to implement the tests, we will refer to the "approximate" method, and whenever a direct numerical implementation of the exact convolution ( 12) is used, we will refer to the "exact" method. It turns out that, as we shall see, the approximate method is most often just as good as the exact method, and about 100 times faster to run.

Models, data set, and statistical tools used in this study

To illustrate the power of the method we propose, several published GGP models have been tested.

We will refer to these as:

• CP model, which is the preferred model of [Constable & Parker 1988];

• QC model, C1 (preferred) model of [START_REF] Quidelleur | [END_REF]];

• CJ model, which is the CJ98 model proposed by [Constable & Johnson 1999];

• JC model, which is the CJ98.nz model also proposed by [Constable & Johnson 1999];

• TK model, which is the TK03.GAD model recently proposed by [Tauxe & Kent 2004] and further discussed in [Tauxe 2005];

• HK model, which is the final model for the normal polarity of [Hatakeyama & Kono 2002].

Models CP, QC, CJ, JC and TK share many characteristics. They are defined by simple axisymmetric mean models for which only E(g 0 1 ) and E(g 0 2 ) can take non-zero values, and purely diagonal

covariance matrices with Cov(g m n , g m n ) = Cov(h m n , h m n ) = (σ m n ) 2 , except in the case of model JC, which assumes different values for Cov(g 1 2 , g 1 2 ) = (σ(g 1 2 )) 2 and Cov(h 1 2 , h 1 2 ) = (σ(h 1 2 )) 2
. In all but one case, for n 3, the σ m n are further assumed to be independent of m and defined by σ m n = σ n where 1 gives the values of the relevant parameters for each of those five models. Finally, model HK differs from the other models because of a more elaborate mean field, defined up to degree and order 4. But its covariance matrices are otherwise defined in much the same way. A full description of this model can be found in Table 2 of [Hatakeyama & Kono 2002].

σ n = α(c/a) n /((n + 1)(2n + 1)) 1/2 . The only exception is model TK which distinguishes σ m n = σ n for (n -m) even from σ m n = βσ n for (n -m) odd. Table
All those models have been constructed in the hope that they would properly describe the statistical properties of the paleomagnetic field at times of normal polarity over the past five million years. In particular, they have been constructed with the help of databases including a large number (if not a majority) of directional data corresponding to the current Bruhnes chron (e.g. [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF], Johnson & Constable 1996]). They are therefore good candidates for an example test against a well-controlled data set covering the Bruhnes chron and corresponding to volcanic directional data acquired at various sites distributed worldwide. Ideally, such a test data set would have to be built by extracting data from the most recent and permanently updated IAGA paleomagnetic reference database for paleofield direction available at the National Geophysical Data Center (http://www.ngdc.noaa.gov/seg/geomag/paleo.shtml), and by relying on a set of fine-tuned criteria agreeable to the community. This however is not a trivial matter since, as noted by one of the reviewers, there currently is no general agreement among investigators on what comprises a satisfactory database for this kind of study. Such a substantial endeavour, which we plan to carry on at a later stage, is therefore clearly beyond the scope of the present paper, which only intends to introduce, test, and illustrate a new methodology. For such a purpose, we felt that a simpler, more readily available data set (and one that at least already went through some type of selection processes relevant to the present study), would be appropriate enough.

We therefore decided to use a test data set extracted from the [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] database,

originally used by [START_REF] Quidelleur | [END_REF] to build their QC model last updated in January 1998, and currently available at http://www.ipgp.jussieu.fr/rech/paleomag/var-secu/. We note however, even before proceeding further, that because of this choice, a close fit of the QC model to our test data set can be anticipated, so that test results reported here will necessarily not test the QC model as stringently as the other models. This test data set consists of a total of 990 independent estimates of the local direction of the paleomagnetic field at 36 sites (Figure 1, Table 2). Each site is located with the help of its latitude and longitude. At each site, each estimate is based on the direction of the resultant vector R of n ( 3) volcanic samples (unit vectors) and is given in the form of a declination D and an inclination I together with both n and the norm R = |R| of the resultant vector. Although R is usually not published as such in the original papers, it can be accurately recomputed from the published material. This then makes it possible to present all the data in an homogeneous and complete form in the database. Indeed, and as is well known (e.g. [McElhinny & McFadden 2000]), an estimate of the local direction of the paleomagnetic field is completely characterized (with is Fisherian error) as soon as the four parameters D, I, n and R are given.

For each estimate, the corresponding vector u i could then be located on the unit sphere Σ, and an estimate of the concentration parameter K i defining the Fisherian error computed [with the help of

K = R (n-1) (n-R)
; note that this is not the concentration parameter k = (n-1) (n-R) classically published with the data, which characterizes the dispersion of the n individual samples about their average direction, but the concentration parameter characterizing the error on u i ]. Starting from u i and K i , and for each GGP model to be tested, a uniformized value t i could finally be computed with the help of Eq. ( 13)

to produce the final uniformized data set {t i }.

To test this final data set against a uniform distribution over [0, 1], we relied on two different well-known tests: the Kolmogorov-Smirnov test (KS-test) and the Anderson-Darling test (AD-test).

These two tests rely on the fact that, if a given uniformized data set {t i } (with i = 1 to N ) is indeed compatible with a uniform distribution over [0, 1], its empirical cumulative distribution function (cdf)

F N (x) should fluctuate within predictable limits about the theoretical cdf value F (x) = x. The null hypothesis should then be rejected if the empirical cdf F N (x) either departs too much from x, or remains to close to it (the latter revealing the presence of suspicious regularities within the {t i }). The KS-test and AD-test differ in the measure chosen to assess how distant F N (x) is from x over [0, 1]:

• the KS-test (as defined by e.g. [START_REF] Press | Numerical Recipes in C. Second edition[END_REF]) uses the maximum value

M N of |F N (x) -x| over [0, 1],
and is therefore most sensitive to departures of the {t i } from a uniform distribution towards the middle of the segment [0, 1];

• the AD-test (as defined by e.g. [Marsaglia & Marsaglia 2004]) uses the integral quantity

I N = N 1 0 (FN (x)-x) 2 x(1-x) dx. Because of the weight [x(1 -x)]
-1 , it is much more sensitive to the behaviour of {t i } at both extremes of the segment [0, 1].

The KS and AD tests are both excellent tests for small samples and can also be used for large samples. Neither require any additional parameter and both can therefore be considered as completely objective. Also, they complement each other very well. In practice, for each of the uniformized data set {t i } we had to test, we therefore computed the values of M N and I N , together with (from known softwares, e.g. [START_REF] Press | Numerical Recipes in C. Second edition[END_REF] and [Marsaglia & Marsaglia 2004]) the probabilities P (M N )

and P (I N ) for the null hypothesis to have possibly produced such large, or even larger, values for respectively M N and I N . Whenever P (M N ) and P (I N ) were found to take values very close to 0 (typically 0.05 or less) or to 1 (typically 0.95 or above), the null hypothesis had to be rejected and the GGP model under consideration had to be considered incompatible with the data-set at this level of confidence.

Validation and discriminating power of the method

Before actually implementing the testing method we propose, we ran a number of instructive validating tests. For that purpose, various artificial data sets have been produced in the following way. Starting from a given GGP model, we first produced random error-free directional data values. Exactly the same number of data values were produced for each of the 36 sites, as available in the (real) test data set (thus amounting to a total of 990 independent artificial values). This could easily be done by just randomly drawing 990 independent sets of Gauss coefficients (g m l , h m l ) from a multidimensional gaussian generator with mean values and covariances matrices chosen as specified by the GGP model under consideration, and computing the predicted values at each site. A random directional error was next added to each of those 990 error-free directional data values. Because, as shown in the appendix, a

Fisher distribution is extremely close to an Angular Gaussian distribution (for realistic Fisherian errors such as those we deal with here), this was actually done by introducing a random Angular Gaussian error (an easier procedure to control than that involved in producing random Fisherian errors). Using this procedure and starting from various specified GGP models, a number of artificial directional data sets affected by specified errors could be produced, which we used to validate our method.

We first validated the "exact" method and checked that for artificial data produced from any GGP models and perturbed by some known error, this method would always conclude to the compatibility of the data set with the starting GGP model, provided the data are assumed to be affected by the correct amount of error. This first series of tests not only allowed us to validate the method. It also allowed us to check our software and the level of numerical accuracy needed to produce accurate results. It turns out that from a numerical point of view, the most sensitive step of the entire method is the numerical implementation of the uniformization procedure (i.e. the computation of Eq. ( 13)). This indeed requires the numerical computation of p i (u) on a grid over the unit sphere Σ. The most demanding computations, in terms of grid size, happen to be those corresponding to p i (u) predicting the largest dispersions of the u i on the unit sphere Σ. In that case, when the grid is not tight enough, numerical errors happen to mainly affect the computation of the t i lying close to 0 or 1. This in turn mainly affects the results of the AD-test which produces too large values of I N , and therefore too small values of P (I N ). As a result, the AD-test tends to more easily produce negative results than it should.

For the purpose of the present paper, we used a grid corresponding to a decent compromise, tight enough for the KS-test to always lead to an accurate result, but not too tight to make it possible to run all tests in a reasonable amount of time. This had the drawback that, for the reasons just outlined, the AD-test produced slightly erroneous results when testing the CJ and JC models. For those two models, we therefore concluded that only the KS-test could be considered reliable. For all other models (CP, QC , TK and HK), both the KS-test and the AD-test were otherwise found to be reliable. Those conclusions have been reached when testing artificial data affected by errors characterised by α 95 of up to 12.5 • . [Note also that to produce those artificial data sets and run all the tests reported in this study, all GGP models have been used only up to degree and order 7, since, as we shall later see, higher degrees contribute negligibly to the final results for all the GGP models we tested].

We next proceeded to compare the results of the "approximate" method to those of the "exact" method, using exactly the same artificial data sets in both cases. Again and as anticipated (see appendix), significant disagreements between the two methods were found to arise only when p i (u) predicts large dispersions of the u i on the unit sphere Σ. In practice, the approximate method was found to lead to accurate enough results for all the models we tested (testing artificial data with all levels of errors up to α 95 = 12.5 • ), except when testing models CJ and JC against data affected by errors with α 95 greater than 7.5 • . Table 3 summarises the domain of validity of each type of test as established from all those preliminary tests.

We then went on to test the ability of the method to properly discriminate the source GGP model and the level of error affecting the data. Because, as we shall soon see, the real data turn out to be compatible only with model QC (from the point of view of the tests described in the present paper), and since the average α 95 for the 990 data in the real database is of 4.7 • , this test was carried out on an artificial data set produced from model QC and affected by errors of α 95 = 4.7 • . This artificial data set was used to test all GGP models, assuming various possible levels of errors (0 • , 2.5 • , 5 • , 7.5 • , 10 • , 12.5 • ). An additional test was also run to check that the method would properly conclude that this data set is compatible with the starting model QC when assuming an error of α 95 = 4.7 • . Since both the exact and approximate methods provide reliable results for errors of that magnitude, both methods have been used with the KS-test and AD-test (except for the the CJ and JC models for which only the KS test was used, for the reasons outlined above). Table 4 reports the results of those tests.

[Note that Table 4 also provides a good illustration of the accuracy of the results of the approximate method, when compared to the results of the exact method].

This table shows that, as expected, the artificial data set is compatible with the starting model (QC model) provided the correct level of error affecting the data is being assumed (α 95 = 4.7 • ). It also shows that by contrast, making a wrong a priori assumption with respect to the errors affecting the data can lead to a negative result. A minimum error of α 95 = 2.5 • and a maximum error of α 95 = 7.5 • is indeed required for the test with model QC to produce positive results (with no rejection at a 95% level of confidence by either of the KS or AD tests). Finally, Table 4 shows that the method also appears to be able to discriminate the various GGP models. The data set produced by QC with α 95 = 4.7 • cannot be attributed to a CP model, unless (but marginally so) one wrongly assumes α 95 = 10 • . It can no more be attributed to a CJ, JC or HK model, unless one again wrongly assumes α 95 = 12.5 • . And it can hardly be attributed to a TK model, whatever the value assumed for α 95 .

Altogether those results thus show that the tests we introduced are sensitive to both the choice of the model to be tested and the level of error assumed to affect the data. In particular they show that if we have a reasonable knowledge of the error affecting the data, the method applied to the type of data set we deal with here (i.e. 990 data distributed over the 36 sites of Figure 1) makes it possible to fully discriminate the models. Indeed Table 4 shows that if we know a priori that the error affecting the data is somewhere between 2.5 • and 7.5 • , the method can detect that all of the CP, CJ, JC, HK and TK models are incompatible with the data set produced from QC model with α 95 = 4.7 • , at more than a 95% level of confidence. Those encouraging results show that the method proposed in the present paper has enough discriminating power to decide which, if any, of the six GGP models can best explain the test data set we have assembled.

Testing the real data

Having validated the method, and assessed its limits and discriminating power, we then went on to test the various GGP models against the real test data set.

Global tests

In a first series of tests, we did not rely on the actual errors provided with the data. Rather, for each model tested, we successively assumed values of α 95 = 0 • , 5 • , 7.5 • , 10 • , and 12.5 • . This made it possible to again run both the "exact" and "approximate" methods in parallel for further validation of the "approximate" method. It also provided us with a quick way of assessing which GGP model could possibly be considered compatible with the test data set. Table 5 shows the results of those tests.

Interestingly, we first note some common features with the results of the previous test (recall Table 4). Model QC is rejected at more than a 95% confidence level for α 95 = 0 • , 10 • , and 12.5 • , but not for α 95 = 5 • and 7.5 • . Again we see that model CP shows a tendency to be compatible with the data if one assumes α 95 = 10 • , though this possibility is now clearly less likely. Finally we see that only by assuming α 95 = 12.5 • can models CJ, JC, HK possibly be consistent with the data, while model TK is always rejected at a 99% confidence level. Based on these tests, we therefore reach the conclusion that if we were to ignore the actual level of error affecting the data in the test data set, all models would be rejected at more than a 95% confidence level except model QC when assuming α 95 = 5 • or 7.5 • , or model CJ and JC when assuming α 95 = 12.5 • .

The key point, now, is that we do have some knowledge of the errors affecting the data in the test data set. We already noted that on average, this error is of α 95 = 4.7 • . This value is clearly incompatible with the α 95 = 12.5 • value needed for models CJ and JC to possibly be compatible with the test data set. But it seems compatible with the results for model QC.

To finally assess the issue, we eventually ran the test for all six models, no longer assuming the same a priori error for all the data of the test data set, but the true error associated (i.e. published)

with each of the 990 data. Since, as already noted, this requires considerably more computational time, only the approximate method was used for this comparative assessment [Note that all previous tests have anyway shown that for errors of the order of 5 • , the approximate method always leads to accurate enough results]. Table 6 gives the results of this final test, and Figure 2 shows the cdfs of the corresponding uniformized data. Those results clearly show that only model QC is not rejected by the test data set. All other models can be rejected at more than a 99% confidence level, TK being the second, CP the third, HK the fourth, CJ the fifth and JC the sixth best models in terms of both the M N and I N measures (note that when both quantities are available in Table 6, M N and I N lead to the same relative ranking).

As an ultimate check and only for model QC, we finally duplicated this last test using the "exact" method. The result of this test (also shown in Table 6) confirms the previous results. Moreover, the difference found between those two results for model QC clearly appears to be much smaller than the differences among the results of the approximate method for all models. We therefore feel quite confident that these results can be trusted.

It is interesting to take advantage of model QC, the only model not rejected by the test data set so far, to also quickly address the issue of the number of Gauss coefficients that should be taken into account to carry on the tests. As stated earlier, we used all Gauss coefficients up to degree and order 7. But all models, including QC, assume non zero variances for all degrees up to infinity. However, the higher the degree n, the weaker the contribution of the coefficients to the local pdfs p i (u), and thus to the final result. Table 7 illustrates this point and shows that the calculation is indeed already converged when taking into account all degrees up to 7. Note that this also shows that the test (applied to the test data set used in the present study), is therefore hardly sensitive to Gauss coefficients of the GGP models with degrees above seven.

Local tests

It is also important to keep in mind that whereas the results presented so far show that models CP, CJ, JC TK and HK should be rejected, they do not imply that model QC is fully compatible with the data.

Those results only correspond to a global test, for which all the 990 data u i have been uniformized into a single global distribution of {t i } over [0, 1], and it may be that model QC, not rejected by this global test, could be rejected by a more local test. Table 8 shows results of such tests. The same testing procedure has been used as before (with both the exact and approximate method, using both the KS and AD tests and assuming the true error), except for the fact that for each test, only the data from one of the 36 sites of the test data set have been used at a time, to produce the (then much less numerous) uniformized {t i }.

Those results show that the KS-test reveals four sites (sites 17, 18, 30 and 35) rejecting model QC at a 99% confidence level because of too low values of P (M N ). This amounts to about 10% of the sites, an uncomfortably high number. By contrast, the AD-test only reveals two such sites (sites 3 and 18). But it confirms the strong level of disagreement of site 18 with model QC. This site clearly needs to be double-checked. Indeed, a closer look at sites 17 and 18, reveals some interesting features.

Those sites are geographically very close to one another, at latitude 50.0 • and longitude 7.0 • . Site 17 (data selected from [START_REF] Böhnel | Paleomagnetic investigation of Quaternary West Eifel volcanics (Germany): indication for increased volcanic activity during geomagnetic excursion/event?[END_REF]), is in the quaternary East-Eifel volcanic field, whereas site 18 (data selected from [START_REF] Böhnel | Palaeomagnetism of Quaternary volcanics of the East-Eifel, Germany[END_REF]), is in the quaternary West-Eifel volcanic field. For those two sites, we show in Figure 3 how the selected data u i plot on the unit sphere Σ. On the same plots, we also show the isovalues of the p i (u), which have been computed assuming α 95 = 2.76 • for site 17 and α 95 = 4.68 • for site 18. Those values correspond to the average errors for all data from respectively sites 17 and 18. (Using such averages is done here mainly for illustrative purpose, the plotted isovalues being then simultaneously meaningful for all the data shown on the same plot). Also shown in Figure 3 are the cdfs for the {t i } corresponding to this calculation, and the analogous cdfs corresponding to the exact calculation (i.e. when the true, rather than the average, error is being used for each individual data). The reasons sites 17 and 18 reject model QC are clear: for site 17, the data tend to display too little dispersion compared to the expected distribution (producing too many low values for the {t i }), whereas for site 18, the data show a strong tendency to lie East of the expected distribution (producing too many high values for the {t i }).

The interesting observation here is that beyond the disagreement of each site with model QC, we see that both sites also disagree with each other. Since they share virtually the same location, this would suggest that two different stable polarity data sets acquired at the same location but at different times during the Bruhnes chron, could disagree. Quite an embarrassing finding in view of the stationarity assumption underlying the GGP modelling approach. Fortunately, a closer look at the [START_REF] Böhnel | Palaeomagnetism of Quaternary volcanics of the East-Eifel, Germany[END_REF], Böhnel et al 1987] papers provides us with an explanation for this disagreement.

As far as site 17 is concerned, it simply appears that the 44 data selected in the [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] database are not fully independent. In seven instances, several data points appear to have basically recorded the same field (each time reflecting a single volcanic event). Those events correspond to low values of t i and bias the distribution of the {t i } towards low values. Site 18 is a slightly different story.

In that case, [START_REF] Böhnel | Paleomagnetic investigation of Quaternary West Eifel volcanics (Germany): indication for increased volcanic activity during geomagnetic excursion/event?[END_REF]] make a good case that a subset of the data they published must have recorded a relatively short-lived excursion or event during the Bruhnes chron, biasing the data set towards the East. Although [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] rejected most of those transitional data (based on the selection criteria that the data must correspond to a period of stable polarity) some of the data biased towards that transitional direction have clearly been included within the database. Thus the issue with sites 17 and 18 appears to be more one of ill-sampling than one of overall disagreement with model QC. Interestingly, we further note that if we actually bin the two sites into a single one, as is also shown in Figure 3, the local test no longer rejects model QC. This then leaves only two sites (30 and 35) out of 36 truly rejecting model QC at a 99% confidence level because of too low values of P (M N ), and only one site (site 3) rejecting QC because of a too low value of P (I N ). These are no longer unexpected proportions and model QC thus indeed appears to be compatible with the [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] Bruhnes stable polarity data set.

Conclusion

In the present paper we introduced the first quantitative method capable of assessing the compatibility of a GGP model with a given paleodirectional database in a way consistent with both the statistical assumptions underlying the GGP approach, and the nature of the measurement error affecting paleodirectional data.

The method was successfully implemented and tested, in its exact form, but also in an approximate form, relying on some useful approximations which we introduced to significantly reduce the computational burden of the tests. The limits of the approximate method have been assessed. Those limits are reached when both the GGP model and the measurement errors introduce strong dispersions of the directional data on the unit sphere. For realistic measurement errors, this approximate method however turns out to be just as good as the exact method. This brought many encouraging results.

First, it was shown that the method is usefully discriminating. A GGP model is compatible with an artificial data set affected by a given (known) level of measurement error only if the model tested is the same as the model used to produce the artificial data and if the correct level of error affecting the data is assumed when carrying the test. But if one assumes either too small or too large an error affecting the data, or tests a GGP model different from the starting model, this conclusion may no longer be reached. The method we propose thus appears to be relevant for deciding which GGP model, if any, best describes the statistical behaviour of the past geomagnetic field.

We applied the method to test six published GGP models against paleodirectional Bruhnes stable polarity data from the [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] database. All but one model (model QC) had to be rejected. This shows that in practice also, our method can successfully discriminate GGP models attempting to describe a given paleomagnetic data set.

The method we propose further offers the possibility of assessing a GGP model at any scale, local, regional or global. This makes it possible to better scrutinise a GGP model which first passed tests at a global scale. Any disagreement between the prediction of the GGP model and the data set at a local level can then be used to double-check not only the GGP model, but also the data set itself.

Applying this checking procedure to model QC allowed us to identify such problems with two sites of the [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] database and to confirm the compatibility of this model with the rest of the test data set.

This success of at least one GGP model is motivating. It shows, for the first time, that the statistical behaviour of the geomagnetic field at times of stable polarity can indeed be described in a consistent way in terms of a GGP model. However the specific success of model QC should not be overemphasised. This model was inferred in an empirical way by [START_REF] Quidelleur | [END_REF] from the [START_REF] Quidelleur | Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database[END_REF] database from which our test data set was also extracted. All other models were inferred (also in empirical ways) from different databases, and even though all those databases share many common data, they may very well not be fully compatible with each other and with our own data set. Part of the failure of those GGP models at being compatible with our test data set might originate from this. With this respect, it is quite clear that finally confirming the success of the QC model (and the failure of other models) would still require a more thorough study, involving a more recent, extensive and independently assessed data set, extracted from databases such as the IAGA paleomagnetic reference database we mentioned early on.

As a matter of fact, it is interesting to note that such a data set could also include paleointensity data extracted from databases such as the IAGA paleointensity reference database (PINT, [START_REF] Perrin | [END_REF]). Indeed, although the present paper mainly described the way to deal with paleodirectional data, paleointensity data could easily be taken into account [if such data provide the full 3D paleofield at a given location, the way those data plot in 3D can be compared with the local 3D

Gaussian statistics predicted by the GGP model (i.e., g) by using just the same kind of uniformization procedure as the one we described, see paper I; and taking measurement error into account would just consist in making use of a modified, but still Gaussian, 3D statistics of the type g a as defined in the appendix, where g e is then directly the Gaussian error assumed to affect the data]. In fact such tests might eventually show that model QC also is incompatible with such a more complete data set. But if such is the case, the method we propose here could be used to seek yet another model. Several strategies could then be pursued. One could involve using the same kind of empirical procedure as those used by [START_REF] Quidelleur | [END_REF] to first guess which parameter of the GGP model should be changed, and using our method to assess the improvement brought. Another strategy, more advanced and much more far fetched at this point, could involve designing a more systematic "inverse" approach based on the present "forward" testing method. In either case however, several issues would still have to be faced. In particular, we mentioned the fact that the method is not sensitive to parameters corresponding to degrees above 7, which reflects the weak contribution of those high degrees to the geomagnetic field observed at the Earth's surface. But we did not investigate in detail the sensitivity of our test to the various parameters defining the mean field (E(k)) and the covariance matrix (Cov(k, k))

of a GGP model. Assessing this would help us better understand how complex a GGP model, and in particular its mean field, really needs to be to explain the data. For the time being, and as far as the present test study can suggest, the Bruhnes stable polarity data does not seem to call for more complex a model than model QC: a mean field with only a g 0 1 and a weak g 0 2 contribution, and a simple diagonal covariance matrix only involving an enhanced σ 1 2 contribution along the lines first suggested by [Kono & Tanaka 1995] and [Hulot & Gallet 1996].

We plan to make our software available upon request (please be in touch with the corresponding author: gh@ipgp.jussieu.fr).

Appendix: Introducing useful analytic approximations

Fisher and Angular Gaussian distributions

Let us first note that a Fisher distribution k K (u, s) centered on u ∈ Σ with concentration parameter K is very close (when K is large enough) to the Angular Gaussian distribution derived from a 3D

Gaussian distribution characterized by m and Λ ij = σ -2 δ i j , provided u = m/|m| and σ = |m|K -1 2 .

This fairly intuitive result (see also [START_REF] Love | Gaussian statistics for paleomagnetic vectors[END_REF]), is in the same line as the better known result that a Fisher distribution is very close to a 2D Gaussian distribution within the plane tangent to Σ at u (see e.g., [START_REF] Merrill | The magnetic field of the Earth[END_REF]).

In spherical coordinates, and when assuming the Fisher distribution centered on u 0 = (θ = 0, ϕ = 0), from ( 7) we infer z = m cos θ and m = K 1 2 . Then, the (axisymmetric) Angular Gaussian distribution (with θ as natural argument) takes the form:

s(θ) = K 1 2 e -K 2 cos θ 2 √ 2π 3/2 + + 1 4π e -1 2 K(sin θ) 2 1 + K cos 2 θ 1 + Erf K 1 2 cos θ √ 2 (14)
To leading order when K becomes large, this becomes:

s(θ) ≈ K (cos θ) 2 2π e -1 2 K(sin θ) 2 (15)
which, to leading order when θ is small enough, further becomes:

s(θ) ≈ K 2π e -1 2 Kθ 2 (16)
This indeed is the approximate form k K (u 0 , s) (as defined by ( 10)) takes when K is large enough, and θ small enough.

To numerically show that for useful values of the concentration parameter K, s(θ) is indeed very close to the k K (u 0 , s) Fisher distribution, we computed the P -probability angles α(P ) (defined as the maximum angle the direction may depart from the mean direction u 0 with probability P ) predicted by k K (u 0 , s):

α {P } = arccos 1 K log e K + P e -K -P e K (17)
and checked that for those angles, the Angular Gaussian distribution s(θ) given by ( 14), leads to a correct prediction of P . For various relevant values of K and P , Table (9) shows values of the departure ε(P ) of the probability predicted by the Angular Gaussian distribution from the expected probability P :

ε(P ) = 2π {θ α {P }} s(θ) sin θdθ -P (18)
As can be seen, ε(P ) never exceeds 0.01 and gets all the smaller that K increases.

This then suggests that for our purposes, data error could be taken into account with the help of an equivalent Angular Gaussian distribution in place of the more conventional Fisher distribution.

Expanding the directional error into a 3D error

The near-equivalence we just established between the Fisher distribution and an Angular Gaussian distribution shows that a good approximation of the convolution (12) defining p i (u), could be obtained by computing the convolution between two Angular Gaussian distributions. Unfortunately there is again no simple analytic form for the result of such a convolution. But this result suggests that yet another approximation could usefully be made (in this appendix, we will drop the "i" indexes for simplicity, since only a single data u i is to be considered at a time).

To see this, first recall that s(u) is in fact the projection on Σ of the 3D Gaussian distribution g(x)

predicted by the GGP process at the site associated with the data to be tested. This then suggests that we:

• seek a 3D Gaussian distribution g e (x) expending in 3D the directional error,

• compute the convolution of g(x) with g e (x) in R 3 , which would then analytically lead to yet another 3D Gaussian distribution g a (x) and,

• analytically compute the Angular Gaussian distribution s a (u) associated with this new 3D Gaussian distribution (thanks to (6)).

Provided g e (x) is properly chosen, we may then hope that the final Angular Gaussian distribution s a (u) provides an accurate enough approximation of p(u) as defined by the convolution (12).

In the limit case the GGP process predicts a g(x) distribution characterized by E(V ) = m and Cov(V , V ) = 0, while the directional error is characterized by the concentration parameter K, we already know that choosing E(V e ) = 0 and Cov(V e , V e ) = |m| 2 K I for g e (x) would do the trick. In that case indeed, s e (u) reduces to the single direction m/|m| so that p(u) given by ( 12) reduces to a Fisher distribution centered on m/|m| with a concentration parameter K. In the same time, g a (x) would be characterized by E(V a ) = m and Cov(V a , V a ) = |m| 2 K I which leads to s a (u). This, we just saw, is very close to a Fisher distribution centered on m/|m| with a concentration parameter K.

In the more general case when, for g(x), Cov(V , V ) is non-zero but remains small enough, we can next anticipate that sticking to the same 3D expansion g e (x) would also do the trick. Then g a (x) would be characterized by E(V a ) = m and Cov(V a , V a ) = Cov(V , V ) + |m| 2 K I, and hopefully, s a (u)

would again be close enough to p(u).

Note that this again can only be an approximation, and that in fact, s a (u) is not even be the exact result of the convolution of s(u) with s e (u) (the near-Fisherian Angular Gaussian distribution associated with g e (x)). This is qualitatively illustrated in Figure ( 4), which shows that the same 3D Gaussian error g e (x) applied to two different realizations V 1 and V 2 of g(x) with |V 1 | larger than

|V 2 |, would lead to an Angular Gaussian error on Σ smaller for u

1 = V 1 /|V 1 | than for u 2 = V 2 /|V 2 |.
Figure (4) however also shows that provided |m| is large enough compared to the dispersion implied by g(x) and g e (x), the distortion introduced should remain small. For the GGP models we have tested, this fortunately happens to be the case most of the time.

CP QC CJ JC TK E(g 0 1 ) -30.0 -30.0 -30.0 -30.0 -18.0 E(g CP QC CJ JC TK HK KS-Test with exact method 12.5 12.5 12.5 12.5 12.5 12.5 KS-Test with app. method 12.5 12.5 7.5 7.5 12.5 12.5 AD-Test with exact method 12.5 12.5 N.A. N.A. 12.5 12.5 AD-Test with app. method 12.5 12.5 N.A. N.A. 12.5 12.5 Figure 2: Cumulative distribution functions (cdf) of the N=990 uniformized data {t i } for each of the six models, when testing them against the real test data set, assuming the true errors and using the "approximate" method. Those cdfs lead to the M N , P (M N ), I N , P (I N ) shown in Table 6. In the case of model QC, the cdf obtained when using the "exact" method leads to a cdf almost indistinguishable from the one shown here. Color code: CP(red), QC(green), CJ(dark blue), JC(pink), TK(light blue), HK(black). shows the way the directional data plots on the unit sphere Σ which is Lambert-projected, the centre point (pointing downwards through the page) corresponding to the direction produced by a pure axial dipole at the site; North is towards the top of the page, East towards the right. Angular distances with respect to the centre of the plot are shown every 5 • (green circles). Isovalues of the p i (u) (computed with α 95 = 2.76 • for site 17 (a), α 95 = 4.68 • for site 18 (b), and α 95 = 3.72 • for both sites combined (c)) are also shown as red ellipses. Those correspond to the bounds within which respectively 10%, 30%, 50%, 70%, 90% of the data points should plot. Central column shows the cdf of the corresponding uniformized data {t i }. Right column shows the analogous cdf when the true errors are used.

O Σ V 1 V 2
Figure 4: Applying the same 3D Gaussian error g e (x) (as symbolised by the two small-dash circles) to two different realisations V 1 and V 2 of a general 3D Gaussian distribution g(x) (symbolised by the dashed ellipses), leads to different Angular Gaussian errors on the unit sphere Σ if |V 1 | is larger than |V 2 | (as symbolised by the thick arcs on Σ).

  Denote g, the corresponding probability density function (pdf) in R 3 .Next, introduce the spherical coordinates (u, ρ) of the vector V (where u = V /|V | ∈ Σ, the unit sphere in R 3 , and ρ = |V | denotes the distance from the origin O of the sphere Σ). Denote s the probability density function associated with the direction u on Σ. This so-called Angular Gaussian distribution ([Bingham 1983], in paper I, we also referred to it as the Gaussian Directional distribution)

  error of 2.76°: Test with true errors: MN=2.30E-01 P(MN)=0.02 IN=2.84E+00 P(IN)=0.03 MN=2.34E-01 P(MN)=0.01 IN=2.81E+00 P(IN)=0.03 Test with true errors: MN=2.87E-01 P(MN)=0.00 IN=8.67E+00 P(IN)=0.00 Test with averaged error of 4.68°: MN=2.82E-01 P(MN)=0.00 IN=8.56E+00 P(IN)=0.00 Test with averaged error of 3.72°: MN=1.15E-01 P(MN)=0.18 IN=3.36E+00 P(IN)=0.02 Test with true errors: MN=1.16E-01 P(MN)=0.17 IN=3.20E+00 P(IN)=0.

Figure 3 :

 3 Figure3: Local tests for model QC at sites 17 (a), 18 (b), and both sites combined (c). Left column shows the way the directional data plots on the unit sphere Σ which is Lambert-projected, the centre point (pointing downwards through the page) corresponding to the direction produced by a pure axial dipole at the site; North is towards the top of the page, East towards the right. Angular distances with respect to the centre of the plot are shown every 5 • (green circles). Isovalues of the p i (u) (computed with α 95 = 2.76 • for site 17 (a), α 95 = 4.68 • for site 18 (b), and α 95 = 3.72 • for both sites combined (c)) are also shown as red ellipses. Those correspond to the bounds within which respectively 10%, 30%, 50%, 70%, 90% of the data points should plot. Central column shows the cdf of the corresponding uniformized data {t i }. Right column shows the analogous cdf when the true errors are used.

Table 1 :

 1 Parameters defining models CP, QC, CJ, JC and TK. Units are in µT , except for β only defined for model TK and which is dimensionless (see main text for further details). Parameters defining model HK are to be found in Table2of[Hatakeyama & Kono 2002].

	0 2 ) σ g 1 σ 0 1 σ 1 1 σ 0 2 2 σ h 1 2 σ 2 2	-1.8 3.0 3.0 2.14 2.14 2.14 2.14	-1.2 3.0 11.72 11.72 -1.5 -1.5 3.0 1.67 1.67 1.3 1.16 1.16 4.3 4.06 1.16 4.3 4.06 8.12 1.3 1.16 1.16	0.0 6.36 1.67 0.58 2.20 2.20 0.58
	α	27.7	27.7	15.0	15.0	7.5
	β					3.8

Table 2 :

 2 Site number, latitude, longitude and number of directional estimates available for each of the 36 sites of the test data set used in the present study.

Table 3 :

 3 Domain of validity of the tests used in the present study. The test can be considered reliable if the error affecting the data is less than the α 95 value shown in the table (units in degrees).

	degree		KS		AD
		M N	P (M N )	I N	P (I N )
	1	0.362	0.00	645.	0.00
	2	0.0808	0.00	21.3	0.00
	3	0.0467	0.03	4.36	0.01
	4	0.0356	0.16	2.13	0.08
	5	0.0330	0.23	1.64	0.15
	6	0.0323	0.25	1.52	0.17
	7	0.0319	0.26	1.48	0.18
	8	0.0319	0.26	1.47	0.18

Table 7 :

 7 Convergence test. Model QC is tested here against the real test data set, using the true errors and the "approximate" method. The test is repeated eight times, taking into account more and more of the Gauss coefficients, up to successively degree 1, 2, 3, 4, 5, 6, 7, and finally 8.

			KS			AD
	site	approx.	exact	approx.	exact
		M N	P (M N )	M N	P (M N )	I N	P (I N )	I N	P (I N )
	1	0.106	0.15	0.105	0.16	1.81	0.12	1.73	0.13
	2	0.149	0.84	0.163	0.75	0.626	0.62	0.598	0.65
	3	0.284	0.05	0.286	0.05	3.63	0.01	3.63	0.01
	4	0.106	0.86	0.105	0.86	0.485	0.76	0.478	0.77
	5	0.122	0.90	0.115	0.93	0.576	0.67	0.542	0.70
	6	0.162	0.82	0.162	0.82	1.01	0.35	0.991	0.36
	7	0.0745	0.98	0.0762	0.98	0.401	0.85	0.413	0.83
	8	0.288	0.16	0.287	0.16	2.96	0.03	2.92	0.03
	9	0.265	0.15	0.264	0.16	1.28	0.24	1.27	0.24
	10 0.208	0.05	0.211	0.05	3.46	0.02	3.33	0.02
	11 0.167	0.27	0.165	0.29	1.46	0.19	1.46	0.19
	12 0.199	0.50	0.198	0.51	0.469	0.78	0.458	0.79
	13 0.156	0.32	0.158	0.30	2.08	0.08	2.23	0.07
	14 0.0867	0.80	0.0878	0.79	0.686	0.57	0.625	0.62
	15 0.304	0.38	0.325	0.30	2.71	0.04	2.88	0.03
	16 0.116	0.66	0.108	0.74	0.876	0.43	0.789	0.49
	17 0.234	0.01	0.233	0.01	2.81	0.03	2.83	0.03
	18 0.287	0.00	0.288	0.00	8.67	0.00	8.76	0.00
	19 0.185	0.12	0.195	0.09	2.31	0.06	2.49	0.05
	20 0.136	0.55	0.136	0.54	0.618	0.63	0.614	0.63
	21 0.212	0.07	0.201	0.10	1.88	0.11	1.67	0.14
	22 0.200	0.28	0.198	0.30	1.77	0.12	1.85	0.11
	23 0.199	0.51	0.199	0.50	1.08	0.32	1.08	0.32
	24 0.147	0.62	0.158	0.52	1.79	0.12	1.91	0.10
	25 0.198	0.21	0.199	0.21	2.41	0.06	2.48	0.05
	26 0.0941	1.00	0.0916	1.00	0.301	0.94	0.315	0.93
	27 0.222	0.65	0.239	0.56	0.796	0.48	0.846	0.45
	28 0.368	0.08	0.373	0.07	2.64	0.04	2.79	0.04
	29 0.140	0.50	0.141	0.50	0.908	0.41	0.994	0.36
	30 0.364	0.01	0.357	0.01	2.52	0.05	2.43	0.05
	31 0.327	0.30	0.321	0.32	1.64	0.15	1.64	0.15
	32 0.179	0.76	0.179	0.76	1.45	0.19	1.40	0.20
	33 0.209	0.78	0.210	0.77	0.925	0.40	0.927	0.40
	34 0.273	0.38	0.274	0.38	1.47	0.19	1.45	0.19
	35 0.248	0.01	0.240	0.01	3.42	0.02	3.23	0.02
	36 0.119	0.91	0.112	0.94	0.480	0.77	0.455	0.79

Table 8 :

 8 Testing model QC against the real test data set on a site by site basis. The true errors and both the "approximate" and "exact" methods are being used. For each site, results are provided for both the KS and AD tests. Figure 1: Locations of the 36 sites of the test data set used in the present study.

	K	α {0.95}	ε(0.1)	ε(0.2)	ε(0.3)	ε(0.4)	ε(0.5)	ε(0.6)	ε(0.7)	ε(0.8)	ε(0.9)
	25	28.3 •	0.004	0.006	0.009	0.009	0.009	0.009	0.006	0.003	0.000

Acknowledgements

We thank Richard Holme and two anonymous reviewers for their constructive comments. IPGP contribution 2147.