RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 20, ES6013, doi:10.2205/2020ES000710, 2020
Andrei V. Khokhlov1,2, Valeriy P. Shcherbakov3, Florian Lhuillier4
1Institute of Earthquake Prediction Theory and Mathematical Geophysics RAS, Moscow, Russia
2Geophysical Center RAS, Moscow, Russia
3Geophysical Observatory "Borok" IPE RAS, Borok, Yaroslavl region, Russia
4Ludwig-Maximilian University of Munich, Munich, Germany
The description of the behavior of the geomagnetic field in the geological past is greatly hampered by the paucity of the data both in space and time. In paleomagnetic studies, this circumstance is partly overcome through the use of the Geomagnetic Axial Dipole (GAD) model, which however becomes unsatisfactory when multipolar terms with non-zero time averages are introduced. Under these conditions, the only way to describe the spatio-temporal evolution of the paleofield is to investigate the temporal evolution of its statistical characteristics. An applicable method for such a description is the use of the Giant Gaussian Process model which does not determine the individual parameters but the probabilistic compatibility of a given model with empirical data. However, the specificity of the data entails many technical difficulties for the implementation of this method when applied to paleomagnetic and paleointensity data. An example of such an analysis applied to the Brunhes epoch clearly revealed the necessity of introducing the significant non-dipole terms in the ancient geomagnetic field configuration for this epoch. In addition, the problem of the discrepancy between models constructed separately for paleointensity and paleoinclination data was discovered. Hypotheses are to explain this discrepancy.
First of all, the study of the Earth's main magnetic field requires the study of the available data. The corresponding theoretical descriptions may vary depending on the time scale we are interested in. Here we address the time scales of secular variation that exceed. Since, at smaller time scales, the magnetic field is also not constant in time, we assume that the value ${\mathbf B} ({\mathbf r}, t)$ (of the magnetic field at time $t$ at a given geographical location ${\mathbf r}$ is constant over smaller time scales. In practice we may set it equal to the time-average of the observable values over some reasonably small time interval. Such an approach guarantees that the annual variation of ${\mathbf B} ({\mathbf r}, t)$ is of the order of few nanoteslas – that is negligible when compared with the precision of paleomagnetic records. The complete description of the Earth's magnetic field should ideally provide (with an appropriate precision) the vectorial value of ${\mathbf B} ({\mathbf r}, t)$ at any epoch in the past and at any geographical location $\mathbf r$. Unfortunately, this desirable target description is hampered by the fact that paleomagnetic data rarely come as 3D vectors. Retrieved from sedimentary and igneous rocks, paleomagnetic observations most often comprise only directional records whereas intensity records are much more seldom. Moreover, such records are sparse in time and space. Understanding the limitations due to the intrinsic of paleomagnetic data is thus of utmost importance.
Mathematically speaking [Backus, 1968] the Earth's magnetic field can only be recovered in a unique way if its vertical component is known everywhere at the Earth's surface (This follows from the uniqueness of the harmonic magnetic potential under the external Neumann boundary conditions). Unfortunately, the vertical component itself is not represented in the existing databases and can hardly be accurately recovered from the remanent magnetization of paleomagnetic samples. Whether we can reconstruct the vector ${\mathbf B} ({\mathbf r}, t)$ from the non-vectorial data, even assuming perfect measurements, is an open question. Two cases are of the special interest: (i) the directional-only data and (ii) the intensity-only data. The following questions are of special interest.
The first two questions are known as "uniqueness problems"; the answers however depend on the additional constraints. For instance, for directional-only data we need to know the number $n$ of loci where the field vector is known to be either zero (no direction) or normal to the surface. The result [Hulot et al., 1997] states that the dimension of the space of solutions cannot exceed $n-1$ in the general case. For instance, the geomagnetic field with only two poles (say, South and North magnetic poles) can be recovered – uniquely up to a global numerical multiplier – from directional-only data gathered at the Earth's surface. Since, in case of perfect measurements, the resulting uncertainty is an unknown parameter, it is enough to additionally know the single value $|{\mathbf B} ({\mathbf r}, t_0)|$. In real case anyway we need not to use the rich intensity data $|{\mathbf B} ({\mathbf r}, t_0)|$ for the recovery the full vector field.
Consider now the case of intensity-only data. This is in some sense the complementary problem to the directional one. G. Backus [1970] showed the absence of uniqueness for the reconstruction, therefore we add to the intensity data the geographical locations where the magnetic field is tangent to the Earth's surface (i.e. dip-equator) as a supplementary condition. Then the main result of [Khokhlov et al., 1997, 1999] states the uniqueness of the reconstructed ${\mathbf B} ({\mathbf r}, t)$.
The conclusion is that the non-vectorial nature of most paleomagnetic records is not a dramatic obstacle if one assumes a GAD field geometry, but the problem of reconstruction still persists if one assumes a more complex geometry [Kaiser, 2012]. However, magnetic excursions and polarity reversals strongly suggests the existence of short-term departures from the GAD geometry. More surprisingly, recent palaeodirectional and paleointensity studies of the Devonian and the Ediacaran periods proved significantly non-dipolar geometry of the field on long-term time intervals. The detailed description of the field structure during these events/epochs is much more complicated than the one produced in the framework of the GAD hypothesis and needs to have more extensive full-vector paleomagnetic data. But even if this condition is fulfilled, this may be not enough as no complete theory of the uniqueness is known yet in the general case.
The assumption of the dipolar-like structure was long time taken for granted, in particular the documented human history of the last millennia proved no significant deviations from the modern structure of the main magnetic field. This dipolar assumption gives rise to many heuristic approximation methods to deal with geographically sparse data at a given moment. It is also crucial to all results in reconstruction of the ancient continents. In the following section, we will thus place a special emphasis on the dipolar. assumption.
From the seventeenth century the magnetic field at the surface of the Earth is approximated by a magnetic dipole placed at the Earth's center and tilted on the order of ten degrees with respect to the rotational axis of the Earth. The quality of such approximation can be measured in average over the surface and provide roughly 80 percents of the real Earth's magnetic field there. The dipole axis intersects the Earth's surface at two points, referred to as the geomagnetic poles – they are differ in general from loci where the real field vector ${\mathbf B} ({\mathbf r}, t)$ is normal to the surface. The quality of such a dipolar approximation and the direction of the dipolar axis are both non-constant in time – these are the known facts about the observable modern secular variations of the Earth's magnetic field. The rather strong assumption is that the similar approximation is valid for almost all time periods in the Earth's history, i.e. any measured paleomagnetic direction or intensity value must be considered as the remanence of the ancient magnetic field with approximately similar structure (except for the value of the tilting angle $\phi$). In particular we have the dispersion of geomagnetic poles at the Earth's surface Another very strong assumption is that the corresponding probability density function is bimodal with with maxima at the geographic poles. Moreover, paleomagnetists widely believe in the rotational symmetry of this probability density function.
In particular, if we consider the expectancy of unit vectors that correspond to the geomagnetic poles distributed in the northern (southern) hemisphere then the expectancy vector is aligned to the rotational axis of the Earth. Let this be a weak form of the famous Geocentric Axial Dipole hypothesis which actually claims much more: the time-averaged geomagnetic field is that of a geocentric axial dipole [Merrill et al., 1998]. Neither the weak form nor the strong form of the GAD statement say nothing about the scale of the time-averaging, it is widely believed (after [Hospers, 1954]) that several thousands of years would be enough to provide the good fit. However the realness of this depends upon the hypothesis related to the value of the tilting angle and the rate of secular variations at all scales. In the more rigorous fashion we have to assume something about the stationarity of the corresponding random process of paleosecular variation (PSV) and the time-scale of its autocorrelation function. Indeed the rate of convergence of the arithmetical averages to expectancy is subjected to the properties of the correlation function. It is not at all traditional to incorporate these statistical considerations to paleomagnetic studies. As a result a lot of heuristic arguments related to the structure of the ancient magnetic field coexist even though being even in a disagreement with each other.
The problem of GAD consistency obviously persists especially for the poor and small datasets related to the deep past (say a billion year-old records, see [Meert, 2009]), however even for the more recent times the corresponding arguments are highly suspicious. Indeed, the approximation of the geomagnetic poles scattering that resulted from the VGP transformation of the paleomagnetic directional data has no definitive precision, since the latter depends upon the unknown tilting angle $\phi$ used in the dipolar approximation. This does not mean that we must reject the GAD but a certain refinement of it is needed.
Before all we have to set the appropriate language of the descriptions and this would be of course the spectral representation of the Earth's main magnetic field (in magnetostatic approximation) at any given moment $t$ of time – that is the linear combination of the gradients of the spherical harmonic functions endowed with the the time-dependent coefficients (usually referenced as gaussian coefficients). In literature the traditional formula of this linear combination is using the magnetic potential in the spherical coordinate system $(r, \theta, \varphi)$ (distance from the Earth's center, colatitude and longitude):
\begin{eqnarray*} \mathbf{B} \:(r,\theta ,\varphi ,t) = - \nabla \Bigl [ a \sum_{l=1}^{\infty} \sum_{m=0}^{l} ( \frac{a}{r})^{l+1} (g_{l}^{m} (t) \cos m \varphi \:+ \end{eqnarray*} \begin{equation} \tag*{(1)} h_{l}^{m} (t) \sin m\varphi) P_{l}^{m} (\cos \theta )\Bigr ] \end{equation}where $g_l^m (t)$ and $h_l^m (t)$ are the time-varying coefficients (gaussian coefficients), $a$ denotes Earth's radius, $P_l^m $ is for Legendre polynomial function. Therefore the problem of reconstruction of the magnetic field is equal to the problem of reconstruction at a time $t$ (with a suitable precision) all gaussian coefficients from the paleomagnetic data that is sparse both spatially and in time. Unfortunately as a general rule the reconstruction of spectra from unevenly sampled data is not robust and in general has very poor precision (for instance this are well known facts for the Fourier spectrum of a real-valued function). Even though the problem of recovering gaussian coefficients from paleomagnetic data is highly nontrivial one can test the compatibility of GAD and paleomagnetic data. Indeed the dipolar approximation means that the magnetic field computed from the magnetic potential 1 with $g_k^l (t) = h_k^l (t) =0$, $k>1$ provides the good approximation of the real paleomagnetic data. Such a test was implemented several times on the timescale 0–5 Ma (where the continental drift is small): the world-wide paleomagnetic poles (recovered by VGP-transform of the normal and reversed paleomagnetic data) for this time span show the scatter around the geomagnetic poles see [Cox and Doell, 1960; Merrill et al., 1998] – this indicates support for GAD at least for its weakest form. Even for the weakest form some departures from the For the strong form of GAD we need the convergency rate of the time-averaged magnetic field – no exact global test is known yet because the accurate recovery of the time-dependency of gaussian coefficients is hardly achievable. For older than a few Ma past time the continental drift is not negligible but the position of ancient continents are usually determined by invoking the GAD: therefore the only argument in support of GAD is the absence of contradiction while positioning the extended regions of the continents during various geological epochs. Again here we use only the weak form of GAD.
The quantitative analysis of dipolar approximation therefore is not at all trivial and known attempts need some additional statistical settings and constraints; we consider the example of such a statistical approach below in more general fashion.
Even though the direct spectral transform of paleomagnetic data to the set of gaussian coefficients cannot be directly applied because of the sparse uneven distribution of the data one can apply mild conditions and consider certain finite approximations of the gaussian coefficients. The corresponding approach after the seminal papers [Bloxham, 1987; Bloxham et al., 1989] is implemented in sequence of geomagnetic models of Holocene, see the latest model in [Constable et al., 2016]. The idea is to find the truncated sequence of gaussian coefficients by means of the least square estimator of the differences between paleomagnetic data and the model field. The continuous in time functions $g_k^l (t)$, $h_k^l (t)$ are defined using finite set of nodes and the corresponding set of cubic B-splines. The robustness test of this construction in [Licht et al., 2013] show the wide enough error bars so this type of modelling provides valuable approximations of gaussian coefficients only for a few initial indices also with the very poor time resolution for times older than 1500 BC. No argument in support of GAD follow from such a model moreover these computations show the presence of the non-dipolar terms.
The conclusion therefore is that the geomagnetic field during the stable polarity epoch likely may have time-averages that are significantly distinct from the pure dipole aligned to the rotational axe of the Earth. This imply the study of the explicit description of the time-average geomagnetic field (i.e. averaged over the given time interval) together with the description of the deviations (paleosecular variations) from this average – these are so-called TAF and PSV studies, for the modern state of this research see [Johnson and McFadden, 2015]. We wouldn't repeat here the known facts from numerous TAF and PSV publications here, instead we recall the ordinary statistical approach and will describe below the corresponding implementation of it for GGP geomagnetic model.
We recall the well-known Mathematical Statistics construction that is key important for our paleomagnetic field studies. A standard statistical procedure tests the relationship between two statistical data sets, or a data set and synthetic data drawn from an idealized model. Such an idealized model typically describes the distribution laws of the random values involved, the non-parametric test then measures the probability to get the observed data in assumption that this theoretical laws are valid. For the statistical test we need the data records that we may assume statistically independent and obey the fixed theoretical distribution, for instance, the Kolmogorov–Smirnov test (KS-test) and the Anderson–Darling test (AD-test). For the most simple and robust case these two tests rely on the fact that, if a given data set $\{x_i\}$, $i=1, \ldots N$ is compatible with a uniform distribution over $[0,1]$, its empirical cumulative distribution function (cdf) $F_N(x)$ should fluctuate within known limits about the theoretical cdf value $F(x) = x$. The tested hypothesis should then be rejected if the empirical cdf $F_N(x)$ departs too far from the function $y=x$ and the measure of such a departure is the probability of the corresponding distance for the idealized dataset. The KS-test and AD-test differ in a way chosen to assess the distance $\| F_N (x) - x\|$ over $[0,1]$:
Because of the weight $[x(1-x)]^{-1}$, it is much more sensitive to the behaviour of $\{x_i\}$ at both extremes of the segment $[0,1]$.
Note, that the quantitative result of these test is always the probability $p$ to observe from the ideal dataset (that obey the uniform law) the distances that are equal the quantity computed for the real data. As a common rule for the natural sciences we declare the incompatibility of the real data and the theoretical law when this probability $p$ is less than some threshold (the given level of consistency), usually 0.05. In brief: we reject those theories that explain the real case as somewhat almost not real! The non-parametric test do not insist on the "proof" of the theory, we stay with conclusion "not enough arguments found in the dataset to reject this theory", however this conclusion is strongly meaningful when we consider the big datasets.
Coming back to the PSV and TAF studies we may try to implement this approach providing we have the explicit theory for the paleomagnetic observables such as directions or intensities. Such a theory was suggested in [Constable and Parker, 1988] and is known as the Giant Gaussian Process or just GGP-model. Moreover this GGP-model defines the special transformation of the paleomagnetic data that converts it into the (presumably uniformly distributed at $[0,1]$) sequences $\{x_i\}$, $i=1, \ldots N$ – therefore KS- and AD-tests are applicable without any significant changes.
The basic assumption of the GGP model claims that at times of stable polarity, the (finitely dimensional) vector with the components $g_{k}^{l} (t)$, $h_{k}^{l} (t)$, $l \le k \le M$
\begin{eqnarray*} \bigl \{ {\mathbf g}(t), {\mathbf h}(t)\bigr \} = \bigl \{ g_{1}^{0} (t), g_{1}^{1} (t), h_{1}^{1} (t), g_{2}^{0} (t), \ldots \end{eqnarray*} \begin{eqnarray*} g_{M}^{M} (t), h_{M}^{M} (t)\bigr\} \quad M \gg 1 \end{eqnarray*}is a single realization of a stationary random vectorial process. According to the general Kolmogorov consistency theorem (also independently discovered by British mathematician Percy John Daniell) [Tao, 2011] a suitable family of the probability distributions at any finite set $t_{1} ,t_{2} , \ldots t_{M} $ will define the random process. In particular, if all such distributions are multidimensional normal (i.e. have Gauss probability densities), the random process called Gaussian process (do not confuse this with the name of spectral coefficients $g_{k}^{l} (t)$, $h_{k}^{l} (t)$) and, therefore the explicit formula of the multidimensional probability densities needs only knowledge of expectations and covariation functions.
In other words, the GGP model assumes that the temporal evolution of the geomagnetic field during a given long period of time (say, for instance, between two reversals) can be described in terms of statistical fluctuations of the field about a mean field.
Because of representation (1) and the stationarity condition the temporal evolution can be described in terms of fluctuations of the random vector $\{ \mathbf{ g,h} \} $ (made of coefficients $\mathbf{g} = \{ g_{1}^{0} ,g_{1}^{1} ,\ldots,$ $g_{l}^{m} \ldots \} $, $\mathbf{h} =\{ h_{1}^{1} ,h_{2}^{1} , \ldots ,h_{l}^{m} \ldots \} $) about some average model $\{ E({\mathbf g}),E({\mathbf h})\}$. Assuming that these fluctuations can be described in terms of a short term memory stationary random Gaussian process (which is consistent with both historical geomagnetic and archeomagnetic variations – see e.g. [Hongre et al., 1998; Hulot and Le Mouël, 1994]), and that any two paleomagnetic observations are always separated by a period of time larger than the memory (short memory is described in terms of rapidly decreasing autocovariation functions) of the process (memory is of the order of a couple centuries), each paleomagnetic data can then be viewed as a local (both in time and space) independent realization of a random Gaussian drawing. Describing the paleomagnetic field in terms of generalized GGP simply consists in identifying the first and the second moments of the Gaussian statistics of ${\mathbf k}$ best predicting the observed statistics for the paleomagnetic data. The mean and fluctuating fields are then characterized by the expectations $E(g_{i}^{j})$, $E(h_{k}^{l} )$ and the covariance matrix.
This description of the paleomagnetic field generalizes the original GGP of [Constable and Parker, 1988] which was build upon two additional important simplifying assumptions: first, that the gaussian coefficients could be considered as being independent from one another, and second, that all Gauss coefficients sharing the same degree $n$, would share the same variance $\sigma_{n} $. In other words, the original GGP further assumed that
\begin{eqnarray*} \mathrm{cov} (h_{i}^{j} , h_{k}^{l} ) = \mathrm{cov}(g_{i}^{j} ,g_{k}^{l} ) = \delta_{ik} \delta_{jl} \sigma _{k}^{2} \end{eqnarray*} \begin{eqnarray*} \mathrm{cov} (g_{i}^{j} , h_{k}^{l} )\equiv 0 \end{eqnarray*}After the studies [Constable and Johnson, 1999; Johnson and Constable, 1997] of the past 5 Ma all recent TAF models for that period include some amount of anisotropy of the fluctuating field. This is indeed the refinement of the GAD conjecture for the given epoch of stable polarity; at a larger scale of several such epochs testing the GAD conjecture rely on the GGP-parameters at smaller scales.
At present all gaussian coefficients are usually assumed to be mutually independent i.e. the autocovariance matrix ${\mathbf C} (0)$ is diagonal. There exist several candidates (that differed in their statistical parameters) to be the best GGP-model and it is worth to compare them using the mentioned above non-parametric statistical approach and various types of paleomagnetic data: intensity-only data, directional data and available collections of 3D paleomagnetic measurements. Below we explain the necessary transformations of the corresponding types of data in order to apply KS- and AD-tests.
We have shown above, that the GGP model for the past geomagnetic field has at least clear statistical structure: for a given set of parameters (expectations and covariation matrices) the testing process of this set against the real data looks like as the routine non-parametric statistical test. This contrast with those approaches that use the GAD hypothesis: for instance computations with VADM and VGP-transform are in fact nonlinear on $g_k^l (t)$, $h_k^l (t)$ and therefore shadow in an uncontrolled way the non-dipole content in the geomagnetic data. Such a convenience of the statistical methods within GGP-approach is already sufficient argument but there is more fundamental consideration in support of GGP-approach in geomagnetic modelling: namely the study of the computational numerical solutions of so-called geodynamo differential equations. It appears reasonable (see e.g. [Bouligand et al., 2016]) to conduct the GGP simple tests using the modern calculations of the numerical geodynamo models: series of such calculations was carried out in the past two decades in a number of studies (the detailed review see [Christensen, 2011]).
In general a statistical analysis can only be correctly conducted when based on the ensemble of the solutions; meanwhile, the corresponding geodynamo data are unavailable for us as of now. Generally speaking, the appearance of non-Gaussian distributions in turbulent dynamical systems was noted in the analysis of the properties of fluid dynamic equations a fairly long time ago. The refined analysis of the hydrodynamic empirical data also showed that the corresponding distributions are similar to the Laplace distribution. The idea to interpret the statistical properties of these data in terms of mixtures of Gaussian laws emerged in the article [Barndorff-Nielsen, 1979]. Computational geodynamo models are also considered for turbulent regimes and, thus, the effects we observe in the numerical geodynamo solutions reasonably well agree with the properties of the data that were previously yielded by the physical experiments with turbulent systems; hence, these effects cannot be attributed to purely computational artifacts.
The cases when the geodynamo parameters correspond, in the context of the up-to-date notions, to the parameters of the magnetic field of the Earth are particularly interesting. In these cases, the key interest lies in the statistical characteristics of the time-dependent coefficients $g_k^l (t)$, $h_k^l (t)$, primarily in respect of elucidating the validity conditions of the GGP hypothesis. However, all hypotheses cannot be tested correctly and correct statistical estimates cannot be obtained with a single numerical solution; this can only be achieved with an ensemble of solutions (that are sufficiently extensive for statistics) with different initial conditions. As is well known, the substitution of ensemble statistics to time statistics is only possible when a random process is stationary and ergodic. In practice, the stationarity and ergodicity of data are frequently assumed without reasonable grounds for this; however, this assumption is not always true and its violation can result in the unexpected consequences.
For instance, in [Khokhlov et al., 2017] we studied the numerical solutions of [Lhuillier et al., 2013] with a resolution up to the spherical harmonics of order 44, which guaranteed the ratio of the minimal-to-maximal amplitudes of the spatial spectrum (averaged over the volume of the liquid core). The total volume of the obtained data contains 2,676,712 time points corresponding to the time interval of 33.45 Ma (It is worth to mention also the early attempt [Bouligand et al., 2005] to test the GGP approach related to the numerical solutions of the early Earth-like model of [Glatzmaier et al., 1999]).
We analyzed the statistics for coefficients $g_k^l (t)$, $h_k^l (t)$ on the intervals of constant polarity. The direct tests of the distribution parameters at reasonable time-scales show the zero-valued expectations for the coefficients of degree $> 1$, however, this may reveal just the nature of the particular computational geodynamo model but may be not relevant to real paleomagnetic data. Though the short memory statement seems with the high confidence to be true, there is also a feature that is for the first glance incompatible with the GGP-model statements.
|
| Figure 1 |
Namely, we constructed the distribution histograms of the coefficients (based on the 2,323,995 data points of the time series) up to the degree $n \leq 12$. These histograms (see for instance Figure 1 left) turned out to be noticeably different from the Gaussian and, moreover, the obtained shapes were much better described by the Laplace distribution formula $y(x) = \mathrm{const} \times \exp (- |b x|)$. On the other hand, from the GGP description, it follows that the corresponding $g_n^l (t)$, $h_n^l (t)$, distributions should be Gaussian and, hence, the appearance of non-Gaussian forms in the numerical solution of a particular geodynamo not only raises the question of GGP applicability limits but also challenges some estimates associated with the statistical processing of the magnetic data.
It should be born in mind that the shape of the histogram depends on the amount of data and the way they are binned. In the case of the applicability of statistical methods the shape of the histogram stabilizes as the amount of the data increases. When instead of an ensemble of independent realizations there is only a single realization of a certain random process, for the applicability of the statistics, it is required that the statistical properties of the data be the same irrespective of the time interval selected for the study. In formal language this property is termed the stationarity of the process and, of course, it is by no means always the case in practice. The statistical estimates based on short segments are usually inaccurate and may only describe the general properties of the distribution. With the increase of the segment, in the case of a stationary process, the corresponding pattern of the distribution becomes more refined; therefore, generally speaking, the characteristic of the stationary behavior for a single realization refers to a certain (particular) length of the time segment. Similarly, in a study of an ensemble of realizations on a fixed time segment, statistical estimates will require a large number of such realizations.
Therefore, we recovered from these numerical solutions that we must be careful when speaking about the stationary character of the geodynamo model and therefore geomagnetic field. For instance, Figure 1 on the right, shows a histogram for coefficient $h_2^2$ for a relatively short segment of the solution containing data samples at successive time points. This shape of the histogram on a semilog scale is reasonably accurately described by the convex curve close to the parabola. Hence, albeit, with some caution, we may speak of the Gaussian character of the process at reasonable time-scales (in geomagnetic units spanning $\approx 100,000$ years and therefore consistent with the time of the constant polarity) while the larger time intervals must be considered as the union of several stable processes however differ in their parameters. In such conditions the stabilization of the histogram about the non-Gaussian form with the increase of the analyzed time interval suggests that the randomization model (hence a non-gaussian model that is also referred to as a mixture of distributions, see [Barndorff-Nielsen, 1979]) is much more suitable here.
The study of the time evolution of coefficients $g_k^l (t)$, $h_k^l (t)$ in the geodynamo solutions [Lhuillier et al., 2013] shows that the secular variations on short time segments and those on intervals that are long in terms of the distribution of the probable values are dissimilar. Specifically, on short segments, the GGP model can be (with a certain degree of caution) assumed to be valid; i.e., secular variations can be considered as a stationary Gaussian process (although, strictly speaking, we should ascertain that a mixed multivariate distribution is also Gaussian). A sufficiently long time interval is however subdivided into several segments so that the parameters of the corresponding stationary Gaussian processes on these segments differ. Hence, the statistics of the global behavior of the process of secular variations will differ from the local statistics, which is reflected, e.g., in the form of the histograms for $g_k^l (t)$ and $h_k^l (t)$. A refined description of the global behavior of secular variations can be given in terms of the mixture of several Gaussian stationary processes, which can be informally described as switching between different modes of behavior. This property is by no means new for dynamical systems since it accompanies the fairly common phenomenon of intermittence in the observed solutions [Frisch, 1995]. The particular numerical geodynamo solutions [Lhuillier et al., 2013] were obtained in a model setting with particular parameters; therefore, the type of the observed intermittence (and even its probable absence) also depends on these parameters. The question on whether intermittence exists, and if so, which particular intermittence exists in the magnetism of the real Earth remains unclear. The answer to this question can be obtained by two approaches: by accumulating representative and accurate data for the intervals of stable polarity in the geological history of the Earth or by studying the question about the correspondence of these types of parameters to the real conditions in the liquid core of the Earth. As of now, the second approach appears to be more realistic, in terms of the available possibilities.
We may conclude that the GGP is indeed a good candidate to model the real geomagnetic field.
The local magnetic field probability distribution function can easily be derived from the formula for magnetic potential: obviously the local pdf is Gaussian. When the data is vectorial, vector errors can first be considered as independent gaussian vectorial increments added to the error-free vector value. As explained in, e.g., [Constable and Parker, 1988] the 1 formula implies that vector samples $ {\mathbf x} = (x_{1} ,x_{2} ,x_{3})$ of the modelled field at a given site (i.e. location $r,\theta ,\varphi $) at the Earth's surface will behave as if drawn from a 3D Gaussian distribution defined by a mean vector ${\mathbf m} = (m_{1}, m_{2}, m_{3})$ (in local cartesian coordinates) and a covariance matrix $\mathrm{Cov} ({\mathbf x}, {\mathbf x}) = \left[ \mathrm{cov} (x_{i} ,x_{j} )\right]$, the details of which depend on the site location, the mean field Gauss coefficients $\{ g_{l}^{m} , h_{l}^{m} \}$ and the covariance matrix $\mathrm{Cov} ({\mathbf k}, {\mathbf k})$ of the GGP model (see, e.g., [Khokhlov et al., 2001, 2006] for detailed formulae, which need not be made explicit here).
In principle, testing whether a GGP model is compatible with paleomagnetic data simply consists in testing such pdfs against data at each site where data have been collected. But these data are often sparse and only relatively few data can be tested against the corresponding distribution for a given site. In addition, these data are always measured and archived with some information about their errors, and this too needs to be taken into account.
When the data is vectorial, dealing with such issues is relatively straightforward. Vector errors can first be considered as independent gaussian vectorial increments added to the error-free vector value. The corresponding 3D-Gaussian error pdf can then be convolved with the 3D-Gaussian pdf to produce yet another 3D-Gaussian pdf to be tested against the data from a given site. At such a single site, and for such a classical comparison, numerous statistical tests are available [Press et al., 2002]. Simultaneously testing data from different sites (to test the regional or global compatibility of a GGP model against such data, assuming the data from different sites are independent) is then also possible. It just requires some preliminary data transformation to ensure that the local pdfs are reduced to a common standard isotropic 3D-Gaussian distribution. This transformation is a linear coordinate change in the local (site) cartesian frame. Regional or global tests can then easily be performed by comparing the transformed data against the common 3D-Gaussian pdf, again using standard tests. This possibility, however, is linked to the fact that all local data satisfy 3D-Gaussian pdfs.
Unfortunately, only a small part of palaeomagnetic records are 3D-vectorial but most of them are directional-only and some – are intensity-only.
In this section we follow our paper [Khokhlov and Shcherbakov, 2015]. The GGP-model prescribes all statistical characteristics of the model data $\mathbf{B} (r,\theta ,\varphi )$ at a given site $r,\theta ,\varphi $, note, that $|\mathbf{B}|$ has the distribution of length of the 3D-Gaussian vector and by no means is Gaussian – we denote the corresponding model probability density function as $f_{|\mathbf{B}|}(x)$. For each value $b$ the probability $u(b) = P\left\{ | \mathbf{B} | < b \right\}$ is the left $u(b)$ - quantile of random $|\mathbf{B}|$ value. On the other hand at the same site we may compare this theoretical quantile with the relative number of those $b_{i} $ from real paleointensity data $\left\{ b_{i} \right\}$ that are smaller than $b$, i.e. we compare the theoretical quantile with empirical quantile. When the model is true then the empirical quantile should be close to the theoretical one; that is the empirical cumulative distribution function (cdf) of $\left\{u(b_{i} )\right\}$ is close to the diagonal. In other words the set $\left\{u(b_{i} )\right\}$ would be close to uniformly distributed points at $[0,1]$. We call uniformization the transformation $b_{i} \mapsto u(b_{i} )=u_{i} $ because after this we expect the uniform population $\left\{u_{i} \right\}$ in case of compatibility between the GGP-model and initial data. For scalar random populations this approach is classical. In the more general case of several sites we start from GGP-model and produce local statistical behaviours at each site and deduce a local, regional or global measure of the adequacy of the model to the data. To test this final data set against a uniform distribution over $[0,1]$, we relied on two different tests already mentioned above: the KS-test and AD-test.
In practice, for each of the uniformized data set $\{ u_{i} \} $ we had to test, we computed the merit values $M_{N} $ and $I_{N} $, together with the probabilities $P(M_{N} )$ and $P(I_{N} )$ for the null hypothesis to have 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), 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.
For the uniformization transform the declared measurement error should be taken into account. In the error-free case we have for $u(b)$ and the model probability density $f_{|\mathbf{B}|} (s)$
\begin{eqnarray*} u(b) = P\left\{|\mathbf{B}| < b \right\} = \int_{-\infty }^{b} f_{|\mathbf{B}|} (s)ds \end{eqnarray*}and in presence of unbiased error $\alpha $ having pdf $g_{\alpha } $ we must substitute the probability density function $f_{|\mathbf{B}|} $ with the convolution $f_{|\mathbf{B}|} * g_{\alpha } $. The explicit analytic cumbersome expression for $f_{|\mathbf{B}|} $ in a given geographical location exists, but the use of it is not that easy (see Appendix in [Khokhlov and Shcherbakov, 2015]), in practice, therefore, we may use a Monte Carlo approach to synthesize a large number of model data with a declared error and compute with the necessary accuracy the uniformization $b_{i} \mapsto u_{i} $.
Similar reasoning can be applied also for the inclination-only paleomagnetic data providing the locations of geographical locations are known – this is rarely the case for the records older than 5 Ma because of the continental drift impact. In general we may apply the non-parametric test to 1D paleomagnetic data and moreover we may use simultaneously the different types of 1D data providing that these records are independent.
The directional-only data no longer consist of ${\mathbf x} = (x_{1} ,x_{2} ,x_{3} )$ local cartesian coordinate values but of unit vectors ${\mathbf u} = {\mathbf x}/|{\mathbf x}|$. Local non-parametric tests against a GGP-model require the explicit form of the local probability density function predicted by the GGP model in terms of the directional vector ${\mathbf u}$ on the unit sphere $S^{2} $. Note that this Angular Gaussian distribution ${\mathbf s}( {\mathbf u})$ thus results from integration over all lengths $\rho $ of the 3D-gaussian distribution. As shown by [Khokhlov et al., 2006] (see also [Bingham, 1983] for series expansions in the case of an Angular Gaussian distribution corresponding to an isotropic 3D-Gaussian distribution), this indeed gives the explicit formula for the Angular Gaussian distribution. Errors in paleomagnetic directional measurements are commonly treated as Fisherian (Tauxe, L., et al., 2018. Essentials of Paleomagnetism, 5th Web Edition. http://ltauxe.github.io/Essentials-of-Paleomagneti sm/WebBook3.html). To test a given GGP-model against a given directional data set with associated errors, one thus has to convolve a Fisher distribution with the local Angular Gaussian distribution, details see [Khokhlov et al., 2006]. He uniformizing transformation of the directional data however is not as simple as in the case of vectorial or 1D data, see details in [Khokhlov and Hulot, 2013; Khokhlov et al., 2006]. The final simultaneous tests of the uniformized directional data from one or several different sites by means of mentioned above non-parametric AD- and KS-tests are arranged exactly as before.
Taking advantage of the 2D nature of the unit vector ${\mathbf u} \in S^{2} $ we may apply even more sophisticated uniformization procedure that transforms the directional measurements $\left\{ {\mathbf u}_{1} , {\mathbf u}_{2} , \ldots \right\}$ to the vectors in $[0,1] \times [0,1] \subset\mathbb{R} ^{2} $ so that the statistical compatibility of directional measurements with the given GGP-model is equivalent to the statistical compatibility of the population $\left\{(t_{1} ,s_{1} ),(t_{2} ,s_{2} ), \ldots \right\}$ with a uniform distribution in the unit square $[0,1]\times [0,1]$. In that case, both KS- and AD-tests are applicable together with less known additional empirical and approximate tests [Khokhlov and Hulot, 2013].
Assume we are given the particular set $\mathrm k$ of GGP-model parameters and the paleomagnetic database. We may preselect $N$ paleomagnetic records from this database to have them uncorrelated and then using the described above procedure of uniformization (that is based on parameters $\mathrm k$) transform them into the population $\{ x_{i} \} $, $i=1, \ldots N$ from $[0,1]$. After that KS-test and AD-test produce the corresponding confidence values $0
 |
| Figure 2 |
We started this research of databases by using the databases of paleomagnetic directions [Cromwell et al., 2018; Johnson et al., 2008; Quidelleur et al., 1994]), results were partly published [Khokhlov and Hulot, 2013; Khokhlov et al., 2006] or presented at AGU [Hulot et al., 2012; Khokhlov et al., 2013] – in all those cases shown are the compatibilities only for GGP-models with statistically noticable non-zero expectations of non-dipole terms. Recently we published [Shcherbakov et al., 2019] the similar research that uses Brunhes chron intensity records: that is an attempt to recover the directional features of the ancient field from only field intensities. For our analysis, we use the world palaeointensity database http://wwwbrk.adm.yar.ru/palmag/database.html maintained by the Geophysical Observatory "Borok" that includes practically all known published paleointensity results. The result of [Shcherbakov et al., 2019] shown on left on Figure 2 crudely speaking does not contradict but not in a perfect agreement with the results [Hulot et al., 2012; Khokhlov et al., 2013]. Moreover, comparing results for the paleoinclinational and paleointensity data, we see that no bright areas in the left and right figures on Figure 2, corresponding high probabilities of their occurrence, fit each other all over the diagrams. It means that the results that use the intensity-only records differ radically from those that use paleoinclinations from the same database, see right plot on Figure 2 right.
Finding the reasons for this discrepancy is beyond the scope of this paper, because even a thorough statistical analysis provides only a formal assessment of the validity of the hypothesis, which is tested on data from a catalogue. In other words, it is a discrepancy is not explained in the GGP model. However, since the discussion of this problem is of great interest for paleomagnetism and rock magnetism, we present our current view on this issue below.
Differences between the results for paleointensity and paleoinclination are due to a significant underestimation of paleointensity in high latitudes, both in the northern and southern hemispheres [Shcherbakov et al., 2019]. A possible explanation is based on the influence of a significant quadrupole term. A similar observation was recently made by [Muxworthy, 2017].
Processing the database we see that the VADM values obtained for the low-latitudes sites are systematically higher than the values for the high-latitude sites so that the mean VADM values for these distributions are 8.02 and 6.7, respectively. It is intuitively clear that the difference between these values is too large to assume the data for latitudes $> 45\mbox{°}$ and $< 45\mbox{°}$ to be generated by a random process corresponding to the axial dipole hypothesis in which case the mean VADM values should not depend on the latitude of their determination. This consideration is also supported by the outcome of the KS-test which shows that the probability of this event is in this case close to zero.
We note by the way that a mirror explanation of this effect by the overestimation of the paleointensity at high latitudes is also probable; however, for definiteness, we hold the first concept. At the same time, the behavior of paleoinclinations is free of these peculiarities (remember that the statistical tests from above never give unique GGP-model, but only help to exclude those GGP-models that are incompatible with the data).
This naturally raises the question about a discrepant behavior of the paleointensity and paleoinclination data. There are two probable explanations.
The first variant calls for seeking a new model describing the geometry of the geomagnetic field in the past geological epochs and secular variations of the field, which would be able to explain the mentioned peculiarities in the distributions of the paleointensity and paleoinclinations. Clearly, this task requires much effort and goes far beyond the scope of this paper. However, we note one important point: the specifics of the spatiotemporal distribution of the field can be due to the nonstationarity of the geodynamo process. Indeed, according to [Bouligand et al., 2016; Khokhlov et al., 2017], secular variation can be considered as a stationary Gaussian process (as required by the GGP-model) only on a relatively short time interval because on sufficiently long intervals, the nonstationarity of the generation processes of the geomagnetic field becomes apparent. As was shown by the analysis of the statistical characteristics of the geomagnetic field generated in the numerical geodynamo models [Khokhlov et al., 2017], secular variations should rather be described in terms of the intermittence when the generation process is subdivided into the segments with a duration of about 100 ka because the parameters of the respective stationary Gaussian processes differ on these segments. However, the GGP scheme, generally speaking, ignores intermittence providing stationary distributions of the coefficients on all the time intervals despite the fact that the parameters can (and most likely do) vary.
At the same time, it must be recognized that in any case there will be serious doubt that we will manage to obtain the desired result in this way because both the intensity and the angular components of the field are derived from the same potential and are therefore correlated with each other. Hence, we can barely imagine a geometry of the field in which the distribution of the intensities across the globe would strongly differ from dipole while the distribution of the angular elements remain close to dipole (see notes on uniqueness from above).
In our opinion, a much more probable cause of the discussed controversy lies in the artifacts associated with erroneous determinations of the paleointensity [Khokhlov and Shcherbakov, 2015]. We are speaking here of paleointensity because the determinations of paleodirections, particularly from young rocks, are much less prone to the risk of obtaining a wrong result. Here we primarily need to mention the risk of misidentifying the nature of remanent magnetization (NRM). The point is that for the correct paleointensity determination by the Thellier method it is necessary that a rock carry thermoremanent magnetization (TRM). However, as hypothesized in [Draeger et al., 2006; Smirnov et al., 2005], TRM can have a similar temperature stability and, hence, close spectra of the blocking temperatures with chemical remanent magnetization (CRM). This hypothesis was subsequently supported in [Gribov et al., 2017; Shcherbakov et al., 2017], where it was shown by the experiments and calculations that, when applied to CRM, the Thellier method yields underestimated the intensity value. All these points allow us to associate the phenomenon of underestimating VADM with the fact that NRM was misidentified as TRM, whereas the tested rocks actually carried CRM. This interpretation was suggested in [Khokhlov and Shcherbakov, 2015] where, based on the analysis similar to the one described here, the presence of more than 40 underestimated VADM values for the Brunhes in the WDB were pointed out.
However, here we should make a reservation that TCRM could be yet a reliable source of palaeomagnetic information, providing unbiased palaeointensity determinations if TCRM was acquired during oxy-exsolution of TM grains below Curie temperature [Shcherbakov et al., 2019].
Thus, it must be recognized that as of now we cannot offer a self-consistent reasonable explanation for the discussed peculiarities of the Brunhes distribution of paleointensities, and this problem remains an enigma to be solved by future research studies. Undoubtedly, its solution is vitally important and highly relevant for the interpretation of the paleointensity data in all the geological epochs. Developing a technique for the joint analysis of the paleodirections and paleointensities instead of analyzing them separately could be a possible way of solving this problem.
Backus, G. E. (1968) , Application of a non-linear boundary-value problem for Laplace's equation to gravity and geomagnetic intensity surveys, The Quarterly Journal of Mechanics and Applied Mathematics, 21, no. 2, p. 195–221, https://doi.org/10.1093/qjmam/21.2.195.
Backus, G. E. (1970) , Non-uniqueness of the external geomagnetic field determined by surface intensity measurements, Journal of Geophysical Research, 75, no. 31, p. 6339–6341, https://doi.org/10.1029/ja075i031p06339.
Barndorff-Nielsen, O. (1979) , Models for non-Gaussian variation, with applications to turbulence, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 368, no. 1735, p. 501–520, https://doi.org/10.1098/rspa.1979.0144.
Bingham, C. (1983) , A series expansion for angular gaussian distribution, Statistics on Spheres, Watson G. S. (Ed.), p. 226–231, Wiley, US.
Bloxham, J. (1987) , Simultaneous stochastic inversion for geomagnetic main field and secular variation: 1. A large-scale inverse problem, Journal of Geophysical Research: Solid Earth, 92, no. B11, p. 11,597–11,608, https://doi.org/10.1029/jb092ib11p11597.
Bloxham, J., D. Gubbins, A. Jackson (1989) , Geomagnetic secular variation, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 329, no. 1606, p. 415–502, https://doi.org/10.1098/rsta.1989.0087.
Bouligand, C., N. Gillet, et al. (2016) , Frequency spectrum of the geomagnetic field harmonic coefficients from dynamo simulations, Geophysical Journal International, 207, no. 2, p. 1142–1157, https://doi.org/10.1093/gji/ggw326.
Bouligand, C., G. Hulot, et al. (2005) , Statistical palaeomagnetic field modelling and dynamo numerical simulation, Geophysical Journal International, 161, no. 3, p. 603–626, https://doi.org/10.1111/j.1365-246x.2005.02613.x.
Christensen, U. R. (2011) , Geodynamo models: Tools for understanding properties of Earth's magnetic field, Physics of the Earth and Planetary Interiors, 187, no. 3–4, p. 157–169, https://doi.org/10.1016/j.pepi.2011.03.012.
Constable, C. G., C. L. Johnson (1999) , Anisotropic paleosecular variation models: implications for geomagnetic field observables, Physics of the Earth and Planetary Interiors, 115, no. 1, p. 35–51, https://doi.org/10.1016/s0031-9201(99)00065-5.
Constable, C. G., R. L. Parker (1988) , Statistics of the geomagnetic secular variation for the past 5 m.y., Journal of Geophysical Research, 93, no. B10, p. 11,569, https://doi.org/10.1029/jb093ib10p11569.
Constable, C. G., M. Korte, S. Panovska (2016) , Persistent high paleosecular variation activity in southern hemisphere for at least 10,000 years, Earth and Planetary Science Letters, 453, p. 78–86, https://doi.org/10.1016/j.epsl.2016.08.015.
Cox, A., R. R. Doell (1960) , Review of paleomagnetism, Geological Society of America Bulletin, 71, no. 6, p. 645, https://doi.org/10.1130/0016-7606(1960)71[645:rop]2.0.co;2.
Cromwell, G., C. L. Johnson, et al. (2018) , PSV10: A Global Data Set for 0–10 Ma Time-Averaged Field and Paleosecular Variation Studies, Geochemistry, Geophysics, Geosystems, 19, no. 5, p. 1533–1558, https://doi.org/10.1002/2017gc007318.
Draeger, U., M. Prévot, et al. (2006) , Single-domain chemical, thermochemical and thermal remanences in a basaltic rock, Geophysical Journal International, 166, no. 1, p. 12–32, https://doi.org/10.1111/j.1365-246x.2006.02862.x.
Frisch, U. (1995) , Turbulence, Cambridge University Press, Cambridge, https://doi.org/10.1017/cbo9781139170666.
Glatzmaier, G. A., R. S. Coe, et al. (1999) , The role of the Earth's mantle in controlling the frequency of geomagnetic reversals, Nature, 401, no. 6756, p. 885–890, https://doi.org/10.1038/44776.
Gribov, S. K., A. V. Dolotov, V. P. Shcherbakov (2017) , Experimental modeling of the chemical remanent magnetization and Thellier procedure on titanomagnetite-bearing basalts, Izvestiya, Physics of the Solid Earth, 53, no. 2, p. 274–292, https://doi.org/10.1134/s1069351317010062.
Hongre, L., G. Hulot, A. Khokhlov (1998) , An analysis of the geomagnetic field over the past 2000 years, Physics of the Earth and Planetary Interiors, 106, no. 3–4, p. 311–335, https://doi.org/10.1016/s0031-9201(97)00115-5.
Hospers, J. (1954) , Rock Magnetism and Polar Wandering, Nature, 173, no. 4416, p. 1183–1184, https://doi.org/10.1038/1731183a0.
Hulot, G., J. L. Le Mouël (1994) , A statistical approach to the Earth's main magnetic field, Physics of the Earth and Planetary Interiors, 82, no. 3–4, p. 167–183, https://doi.org/10.1016/0031-9201(94)90070-1.
Hulot, G., A. Khokhlov, C. Johnson (2012) , How different is the time-averaged field from that of a geocentric axial dipole? Making the best of paleomagnetic directional data using the statistical giant gaussian process approach, Abstract GP24A-03 presented at 2012 Fall Meeting, AGU, San Francisco, Calif., 3–7 Dec., p. GP24A-03, AGU, US (https://ui.adsabs.harvard.edu/abs/2012AGUFMGP 24A..03H/abstract).
Hulot, G., A. Khokhlov, J. L. L. Le Mouël (1997) , Uniqueness of mainly dipolar magnetic fields recovered from directional data, Geophysical Journal International, 129, no. 2, p. 347–354, https://doi.org/10.1111/j.1365-246x.1997.tb01587.x.
Johnson, C. L., C. G. Constable (1997) , The time-averaged geomagnetic field: global and regional biases for 0–5 Ma, Geophysical Journal International, 131, no. 3, p. 643–666, https://doi.org/10.1111/j.1365-246x.1997.tb06604.x.
Johnson, C. L., P. McFadden (2015) , The Time-Averaged Field and Paleosecular Variation, Treatise on Geophysics, p. 385–417, Elsevier, Amsterdam, https://doi.org/10.1016/b978-0-444-53802-4.00105-6.
Johnson, C. L., C. G. Constable, et al. (2008) , Recent investigations of the 0-5 Ma geomagnetic field recorded by lava flows, Geochemistry, Geophysics, Geosystems, 9, no. 4, p. 1–31, https://doi.org/10.1029/2007gc001696.
Kaiser, R. (2012) , Uniqueness and non-uniqueness in the non-axisymmetric direction problem, The Quarterly Journal of Mechanics and Applied Mathematics, 65, no. 3, p. 347–360, https://doi.org/10.1093/qjmam/hbs005.
Khokhlov, A., G. Hulot (2013) , Probability uniformization and application to statistical palaeomagnetic field models and directional data, Geophysical Journal International, 193, no. 1, p. 110–121, https://doi.org/10.1093/gji/ggs118.
Khokhlov, A., V. Shcherbakov (2015) , Palaeointensity and Brunhes palaeomagnetic field models, Geophysical Journal International, 202, no. 2, p. 1419–1428, https://doi.org/10.1093/gji/ggv236.
Khokhlov, A., G. Hulot, C. Johnson (2013) , On the contribution of g20 and g30 in the time-averaged paleo- magnetic field: First results from a new giant gaussian process inverse modeling approach, Abstract GP53D-04 presented at 2013 Fall Meeting, AGU, San Francisco, Calif., 9–13 Dec., p. GP53D-04, AGU, US (http://abstract search.agu.org/meetings/2013/ FM/GP53D-04.html).
Khokhlov, A., G. Hulot, C. Bouligand (2006) , Testing statistical palaeomagnetic field models against directional data affected by measurement errors, Geophysical Journal International, 167, no. 2, p. 635–648, https://doi.org/10.1111/j.1365-246x.2006.03133.x.
Khokhlov, A., G. Hulot, J. Carlut (2001) , Towards a self-consistent approach to palaeomagnetic field modelling, Geophysical Journal International, 145, no. 1, p. 157–171, https://doi.org/10.1111/j.1365-246x.2001.01386.x.
Khokhlov, A., G. Hulot, J. L. Le Mouël (1999) , On the Backus Effect–II, Geophysical Journal International, 137, no. 3, p. 816–820, https://doi.org/10.1046/j.1365-246x.1999.00843.x.
Khokhlov, A., G. Hulot, J. L. L. Le Mouël (1997) , On the Backus Effect-I, Geophysical Journal International, 130, no. 3, p. 701–703, https://doi.org/10.1111/j.1365-246x.1997.tb01864.x.
Khokhlov, A. V., F. Lhuillier, V. P. Shcherbakov (2017) , Intermittence and peculiarities of a statistic characteristic of the geomagnetic field in geodynamo models, Izvestiya, Physics of the Solid Earth, 53, no. 5, p. 695–701, https://doi.org/10.1134/s106935131705007x.
Lhuillier, F., G. Hulot, Y. Gallet (2013) , Statistical properties of reversals and chrons in numerical dynamos and implications for the geodynamo, Physics of the Earth and Planetary Interiors, 220, p. 19–36, https://doi.org/10.1016/j.pepi.2013.04.005.
Licht, A., G. Hulot, et al. (2013) , Ensembles of low degree archeomagnetic field models for the past three millennia, Physics of the Earth and Planetary Interiors, 224, p. 38–67, https://doi.org/10.1016/j.pepi.2013.08.007.
Meert, J. G. (2009) , In GAD we trust, Nature Geoscience, 2, no. 10, p. 673–674, https://doi.org/10.1038/ngeo635.
Merrill, R., M. McElhinny, P. McFadden (1998) , The Magnetic Field of the Earth, Academic Press, San Diego.
Muxworthy, A. R. (2017) , Considerations for Latitudinal Time-Averaged-Field Palaeointensity Analysis of the Last Five Million Years, Frontiers in Earth Science, 5, https://doi.org/10.3389/feart.2017.00079.
Press, W., S. Teukolsky, et al. (2002) , Numerical Recipes in C++, Cambridge University Press, Cambridge.
Quidelleur, X., J.-P. Valet, et al. (1994) , Long-term geometry of the geomagnetic field for the last five million years: An updated secular variation database, Geophysical Research Letters, 21, no. 15, p. 1639–1642, https://doi.org/10.1029/94gl01105.
Shcherbakov, V. P., S. K. Gribov, et al. (2019) , On the Reliability of Absolute Palaeointensity Determinations on Basaltic Rocks Bearing a Thermochemical Remanence, Journal of Geophysical Research: Solid Earth, 124, no. 8, p. 7616–7632, https://doi.org/10.1029/2019jb017873.
Shcherbakov, V. P., A. V. Khokhlov, N. K. Sycheva (2019) , Analysis of the Hypothesis of a Giant Gaussian Process as a Means for Describing Secular Variations of the Geomagnetic Field Vector, Izvestiya, Physics of the Solid Earth, 55, no. 1, p. 182–194, https://doi.org/10.1134/s1069351319010099.
Shcherbakov, V. P., N. K. Sycheva, S. K. Gribov (2017) , Experimental and numerical simulation of the acquisition of chemical remanent magnetization and the Thellier procedure, Izvestiya, Physics of the Solid Earth, 53, no. 5, p. 645–657, https://doi.org/10.1134/s1069351317040085.
Smirnov, A. V. (2005) , Thermochemical remanent magnetization in Precambrian rocks: Are we sure the geomagnetic field was weak?, Journal of Geophysical Research, 110, no. B6, https://doi.org/10.1029/2004jb003445.
Tao, T. (2011) , An Introduction to Measure Theory, American Mathematical Society, Providence, Rhode Island.
Received 4 March 2020; accepted 6 April 2020; published 13 December 2020.
Citation: Khokhlov Andrei V., Valeriy P. Shcherbakov, Florian Lhuillier (2020), Using the Giant Gaussian Process model from paleodirectional and paleointensity data to investigate paleomagnetic secular variation, Russ. J. Earth Sci., 20, ES6013, doi:10.2205/2020ES000710.
Copyright 2020 by the Geophysical Center RAS.