Skip to main content

Geomagnetic observatory monthly series, 1930 to 2010: empirical analysis and unmodeled signal estimation

Abstract

Ground-based magnetic observatory series are the main source of information for constructing time-dependent spherical harmonic geomagnetic field models from sub-annual to pluri-decadal time scales. Assessing the reliability of such models requires accurate estimation of the data errors. We propose an analysis of observatory monthly means over the period 1930 to 2010, where we sequentially isolate (i) a stochastic regression for the main field at every site, performed in the framework of Gaussian processes, (ii) a local fit to annual and semiannual signals, (iii) a month by month estimate of global, large length-scale external and induced fields. We then estimate the unmodeled signal level (UMSL, which refers to the instrumental noise plus extra signals not captured by the above data treatment) from the standard deviation of the residuals to the sequential analysis. This may be used to estimate data error covariances in future field modeling studies. Mainly a function of the geomagnetic latitude, the UMSL is larger towards auroral regions and carries the temporal signature of solar activity. While the UMSL shows rather similar magnitudes in all three components in recent epochs (typically a few nT), a significant decrease is found in the downward component of the field around 1960, which correlates with the introduction of proton magnetometers. We detail the geographic distribution of the periodic signals and confirm the variation of their amplitude at pluri-decadal time scales. From the spherical harmonic description of horizontal and vertical fields, we isolate the main patterns of the inducing field in Z. These are dominated by a zonal structure of degree 1 (and to a lesser extent, of degree 3) in dipole coordinates. We nevertheless isolate secondary, non-zonal sources that are most active during the 1960s and around 1990, periods of particularly large solar activity, denoting an unusual morphology of the inducing field.

Background

Dissociating external (ionospheric, magnetospheric) and internal (induced, crustal, core dynamo) magnetic sources is a major limitation to the spherical harmonic global modeling of geomagnetic records (e.g. Olsen et al. 2010a). This has repercussions on studies concerning both core physics and external fields. Indeed, the recovery of core motions at periods shorter than a few years is hindered by the domination of external sources towards high frequencies, and the reconstruction of changes in the magnetospheric ring current over decadal periods (e.g. Gillet et al. 2013; McLeod 1996; Sabaka et al. 1997) is made ambiguous because of non-global coverage of observatory sites.

Observatory series are the only records that permit the analysis of the geomagnetic field from sub-annual to pluri-decadal periods. Several global modeling avenues may be employed using such data. In a comprehensive approach, Sabaka et al. (2002, 2004) simultaneously inverted all sources from hourly values spanning 1962 onward. Alternatively, the complex signature from external fields is considered as contributing to the data error budget when recovering main field changes from the annual means series starting in 1840 (Jackson et al. 2000). Intermediate strategies considering a simple external field parameterization (Gillet et al. 2013), or a sequential elimination of the sources, may also be used.

However, assessing the robustness of the output models is made difficult by the crude estimates of the observation errors. Accurate magnitudes of the data error covariances are required in order to provide realistic a posteriori confidence levels for the field models. For instance, by allocating constant error bars on observatory annual means series, Jackson et al. (2000) and Gillet et al. (2013) have relatively underweighted the information contained in the most accurate data when building, respectively, the gufm1 and COV-OBS field models. With such an approach, a posteriori uncertainties tend to be overestimated towards recent epochs (when the data are of better quality), and time changes in the variance of some unmodeled signals are ignored (when they should be accounted for through the data error covariance matrix).

One could, in principle, consider the best documented period (the satellite era) to characterize unmodeled sources at ancient epochs. However, this approach is limited by the existence of external field variations for periods longer than a decade. Furthermore, there is no consensus on the variance one should consider to characterize satellite data errors (Finlay et al. 2012; Lesur et al. 2010; Olsen et al. 2010b). In a situation where there exists a frequency overlap within the different sources (given the recent estimates of the mantle conductivity, there is possibly a signature in magnetic series of core processes at periods shorter than a year, see Velímski and Finlaý (2011)), it is tempting to analyze the observatory monthly mean (OMM) series to try to better isolate the several contributions at periods ranging from a few months to a few decades.

When allocating error variances to observatory series, residuals from smooth fits to individual series, or from an initial time-dependent global field model prediction, have sometimes been considered. As an example, the generalized cross-validation method was applied to individual observatory annual mean series (see Bloxham and Jackson 1992) for the construction of the gufm1 and COV-OBS models, and Wardinski and Holme (2006) and Wardinski and Lesur (2012) calculated error cross-covariances between residuals from three-component OMM series in constructing the C 3 FM models from 1957 onward. This pragmatic avenue nevertheless implicitly relies on the modeling assumptions, through a regularization process involving free damping parameters and specific choices of temporal basis functions. Those govern the temporal smoothness through their differentiability properties and density and thus implicitly state the cutoff between the noise and the signal.

Using principal component analysis, Wardinski and Holme (2011) found large correlations between the residuals from a main field model and the Dst magnetic index. This led them to remove from the data series a statistical proxy for unmodeled signals of external origin. Such an approach is promising in the quest for detailed secular variation patterns at observatories. However, indices such as Dst do not necessarily represent well the external activity worldwide (especially at high latitudes). Furthermore, it may be biased at periods longer than a year because of the time changes of the Dst baseline (see for instance Lühr and Maus 2010; Olsen et al. 2005b). The alternative RC-index (Olsen et al. 2014) has recently been introduced to avoid this specific issue.

In an ideal situation, in the context of global spherical harmonic modeling, one should follow a comprehensive approach where all signal sources are simultaneously inverted for, in order to see the unmodeled signals vanish as much as possible. This is unfortunately not possible when looking back to the early observatory era before 1960, as common measurement protocols had not yet been adopted via the INTERMAGNET global network of magnetic observatories (e.g., Love and Chulliat 2013), and the measurement quality was poorer (e.g., Matzka et al. 2010). In this configuration, modeling errors have to be considered along with measurement uncertainties.

The purpose of the present work is to present an empirical analysis of OMM series, by which we provide an estimate of the unmodeled signal level (UMSL). It relies on a sequential cleaning of the contributions from several magnetic sources in the monthly series of the northward (X), eastward (Y), and downward (Z) components of the geomagnetic field. These compose the magnetic field vector B in geographic coordinates. We isolate three main contributions as follows:

  • First, a stochastic fit for the main field signal is obtained separately for each observatory using a Gaussian process regression that includes our prior knowledge of the core field temporal cross-covariances.

  • Next, signals with annual and semiannual periods are fit to the residuals between the OMM series and the main field regression, performed separately for each observatory.

  • Finally, we perform a month-by-month large length-scale global inversion of the remaining residuals. This approximates the main external source together with its associated induced fields.

Using a monthly sampling rate, we have to consider the large signals of periods of 1 year and half a year. If these have been isolated in observatory series for more than a century (e.g., Chapman and Bartels 1940), their mechanisms are still under debate. For this reason and because of their complex geographic distribution, we must fit them separately from the global, large length-scale external fields.

After presenting the OMM series used throughout this study (Section ‘Data’), we detail in Section ‘Methods’ the three steps of our sequential analysis. The several contributions to the OMM series are presented in Section ‘Results and discussion’, together with the variations of the UMSL in space and time. Finally, we summarize our analysis in Section ‘Conclusions’.

Data

We present in this section the OMM database employed throughout this study. We refer for instance to Matzka et al. (2010) for a more complete description of observatory data and the several signals they contain. We benefit from the strategy of Chulliat and Telali (2007), who built the worldwide OMM archive from hourly mean values provided by the world data center (WDC), checking for discontinuities (baseline jumps) in magnetic series. We build our own OMM series from hourly means provided by the Edinburgh WDC Catalog, for two main reasons:

  • Some series for which hourly means are available are not included in the Chulliat and Telali (2007) database (either because of lack of continuity, too low quality, or the presence of incompatibilities with the annual means provided by the WDC). Instead, we consider that those data may be helpful to infer some information on the core signal at ancient epochs and find it interesting to estimate their associated UMSL.

  • We recalculate the OMM in a slightly different way. In anticipation of a future purpose of main field modeling, we use only the midnight 3-h data in order to avoid the diurnal variation from the ionospheric Sq current system, which is large (weak) in the local daytime (midnight). We are aware that induced fields, with amplitudes of a few nT, are associated with Sq variations (Olsen et al. 2005a), but these essentially consist of large length-scale patterns that should be accounted for through the global large length-scale fit. To retrieve the UMSL of the OMM, we thus do not consider the impact from Sq variations.

We paid much attention to the correction for discontinuities in the monthly series. Since a large number of these ‘jumps’ are neither documented nor integrated in the WDC archives, we follow the strategy of Chulliat and Telali (2007) to diagnose and fix the baseline change issues in the series of daily values of midnight to 3-h means. It is unfortunately not possible to recognize all these offsets, especially the small ones, which will be hidden by the signature of natural rapid variations from external or induced origin. We thus consider that this type of undiscovered discontinuities contribute to the unmodeled signals.

We finally end up with a total of 197 OMM series. The total number of sites changes with time (see Figure 1). Additionally, we introduce some monthly data from Chinese observatories for the period 1995 to 2008 that are not published in the WDC and are not part of INTERMAGNET. These data are also calculated from daily values and are collected and administrated by the Institute of Geology and Geophysics, Chinese Academy of Sciences. All observatories used throughout this study are listed in Table 1.

Figure 1
figure 1

Number of observatories as a function of time for horizontal and vertical components.

Table 1 List of the observatories used throughout this study

Methods

Local stochastic regression of the main field

We wish to avoid the use of a previously built time-dependent geomagnetic field model when removing the internal contribution to the observed series, because most spherical harmonic internal models rely on temporal basis functions presenting specific properties that are not justified in a physical or statistical manner. Along the same lines, we wish to avoid the use of global (e.g., polynomial) or local (e.g., splines) functions to directly fit the series. We thus employ the Gaussian processes modeling framework (e.g., Rasmussen and William 2006), where the a posteriori expected value p of the process, calculated at epochs stored in a vector t p (of size N p ), is entirely determined by the knowledge of (i) an observation vector d recorded at epochs stored in a vector t d (of size N d ), (ii) the covariance matrix C ee that describes the statistics of the observation errors (of size N d ×N d ), and (iii) the a priori covariance function (τ)=E(φ(t)φ(t+τ)) for the series φ(t) of time t (where E stands for ‘expectation’). It comes down to calculating the best linear unbiased estimate (BLUE):

$$\begin{array}{@{}rcl@{}} \hat{\mathbf{p}} = {\sf C_{pd}}\left({\sf C_{dd}} +{\sf C_{ee}}\right)^{-1}\textbf{d}\,. \end{array} $$
((1))

Here, the element in the ith row and jth column of the matrix C pd =E(p d T), of size N p ×N d , is C pd ij =(t p i t d j ); and similarly for C dd =E(d d T), of size N d ×N d , we have C dd ij =(t d i t d j ).

Our choice for (τ) is motivated by the analysis of observatory series. The auto-covariance function of a series defines its power spectrum (f), where f is the frequency. As illustrated in Figure 2, that of observatory records evolves approximately as (f)f −4 at periods longer than 5 years (Currie 1968; De Santis et al. 2003). As highlighted by Gillet et al. (2013), it is convenient to model such series in the framework of stochastic differential equations (SDEs), for which a −4 slope spectrum at high frequencies corresponds to an SDE of order 2. If there is no evidence for the continuation of the −4 spectrum at periods as short as 1 month, Gillet et al. (2013) nevertheless show that there is no need for a higher order process (and thus a steeper spectrum slope) at such short periods. It is indeed difficult to define the a priori information about the core field for interannual to monthly periods, due to the ambiguity between internal and external magnetic sources.

Figure 2
figure 2

Power spectral density (log10, in nT 2 /yr −1 ) obtained from OMM series. For the three components of the magnetic field at the Kakioka observatory, superimposed with −4 and −1 slopes (black). The series have been corrected for referenced jumps (Chulliat and Telali 2007). A multitaper analysis has been performed. For each sub-series, we removed the end-to-end linear trend before applying a Hanning window (see De Santis et al. 2003). Generating many realizations of 100-year synthetic AR-2 series (for which we know the theoretical −4 slope of the spectrum), we retrieve with the above method a slope of −4.2±0.3 (i.e., a typical error of about 5% to 10%) in the period band from 2 to 25 years.

We choose the particular a priori auto-covariance function

$$\begin{array}{@{}rcl@{}} {\cal C}(\tau)=\sigma^{2}\left(1+\frac{\tau}{\tau_{c}}\right)\exp\left(-\frac{\tau}{\tau_{c}}\right)\, \end{array} $$
((2))

within the Matérn family of functions. It is characteristic of the second-order statistics for the stationary, auto-regressive (AR) process of order 2 (Yaglom 1962)

$$\begin{array}{@{}rcl@{}} \frac{d^{2}\varphi}{dt^{2}}+\frac{2}{\tau_{c}}\frac{d\varphi}{dt}+\frac{\varphi}{{\tau_{c}^{2}}} = \epsilon(t)\,, \end{array} $$
((3))

where ε(t) is a white noise processa. We use the correlation time τ c =1,730 years and a typical standard deviation of σ=40,000 nT. This corresponds to the a priori function used for the dipole only series by Gillet et al. (2013) - see their Equations 14 to 16. Of course, observatory records contain more than the dipole signal. However, we wish here to keep the definition of (τ) as simple as possible, still allowing the predicted series to potentially display the same continuity properties at monthly periods as they do at interannual periods. We performed a test using a more complex a priori function (such as a sum of AR-2 processes, to account for higher order spherical harmonic coefficients), and found that it did not significantly change the local estimate of residuals to the AR-2 fit (this is because we consider all observatories separately, and not all together through a global model). We estimate, at all sites, the BLUE (1) separately for all three components of B. We note B A R2 the stochastic regression of the observed series B. We discuss in Section ‘On the sensitivity of the residuals to the Gaussian process regression’ the sensitivity of the variances of the residuals (BB A R2) to the a priori error variances entering C ee .

Local analysis of annual and semiannual signals

Annual and semiannual signals clearly appear in the OMM series and indices (e.g., Lyatsky and Tan 2003). This is illustrated by the peaks in the power spectral density shown in Figure 2. They essentially result from ionospheric and magnetospheric fields (plus their induced counterparts) and exist worldwide in observatory records (see for instance Wardinski and Mandea 2006). Their short latitudinal extent at high latitude nevertheless prevents them from projecting entirely onto the large length-scale model estimated in Section ‘Spherical harmonic modeling of remaining large length-scale fields’, which is why we treat them separately at this stage of our empirical analysis.

In the construction of comprehensive field models, periodic signals originating from ionospheric fields are accounted for through a fit of global functions in quasi-dipole coordinates, calibrated a priori in time on the solar radiation flux index F10.7 (see Sabaka et al. 2002). However, this index does not account for the whole modulation of periodic signals. Indeed, the underlying mechanisms responsible for annual variations are subject to discussion and may vary depending on the geographical site (Malin and Winch 1996). Malin and Isikara (1976) suggested that they are generated by seasonal changes in the morphology of the auroral electrojet and the ring current system and imply variations of conductivity in the ionosphere (see for instance Liu et al. 2009).

Variations with a period of half a year are related to a larger occurrence and intensity of geomagnetic storms and substorms at equinoxes than at solstices (e.g., Cliver et al. 2000). For a long time, semiannual variations have been linked to the largest geomagnetic activity, occurring when the southward component of interplanetary magnetic activity is at a maximum in geocentric solar magnetospheric coordinates (Russell and McPherron 1973). A more complex picture has emerged since then. The latter mechanism is now believed to explain only part of the observed signal, and a more important role is given to the equinoxial hypothesis (see for instance Cliver et al. 2000; Lyatsky et al. 2001), by which the angle between the Earth dipole axis and the solar wind direction governs the efficiency in the response of the magnetosphere to solar wind flow.

Given the complex origin of periodic signals, we favor an empirical approach where the amplitude and phase of the periodic signal are kept completely free, along the lines of the studies by Le Mouël et al. (2004a,b). These authors calculated the Fourier coefficients of annual and semiannual periodic signals from band-pass filtered observatory hourly means. Their analysis indicates that both the amplitude and phase of these signals are changing with time. A similar result was found from wavelet analysis of the semiannual signal in the aa magnetic index (Elias et al. 2011).

Here, we fit simple sine/cosine functions with periods of 6 and 12 months to the remaining residuals δ B=BB A R2. Since neither the amplitude nor the phase are constant, these periodic signals are fitted using a sliding window with a restricted length δ t=5 years. Minimizing

$${} \begin{aligned} &\sum_{|t_{j}-t_{m}|\le \frac{\delta t}{2}} \left[ \delta X(t_{j})-\hat{X}_{1}(t_{m})\cos \left(\omega_{1} t_{j}+\phi_{1}(t_{m})\right)\right.\\[-6pt] &\left.\qquad\qquad\;\;-\,\hat{X}_{2}(t_{m})\cos(\omega_{2}t_{j}+\phi_{2}(t_{m}))\right]^{2} \,, \end{aligned} $$
((4))

we obtain the slowly varying functions ϕ 1(t), ϕ 2(t), \(\hat {X}_{1}(t)\), and \(\hat {X}_{2}(t)\), with ω 1=ω 2/2=2π rad/year (we use similar notations for the Y and Z components). This defines the annual and semiannual signals B 1(t) and B 2(t), whose characteristics are presented in Section ‘Worldwide description of annual and semiannual signals’.

Spherical harmonic modeling of remaining large length-scale fields

Some of the remaining residuals are due to resolvable external and induced magnetic fields at large length-scales. We approximate them by means of global modeling, constructing a discrete set of monthly models. These are considered independent from one epoch to another. Indeed, the power spectrum observed for external contributions, excluding annual and semiannual peaks, is only slightly red (with a slope about −1, see Figure 2).

In principle, external and internal sources can be separated from vector observations on a sphere. However, this separation is ambiguous because of the non-global coverage of observatory sites. Nevertheless, horizontal and vertical components of the field do probe complementary combinations of internal and external sources. Indeed, let us assume that the remaining large length-scale magnetic field B LLS derives from a potential,

$${} \begin{aligned} V_{LLS}(r,\theta,\phi,t)&=a\sum_{n=1}^{N}\sum_{m=0}^{n} \left(\! \left(\frac{a}{r}\right)^{n+1} {g}_{n}^{m}(t) +\left(\frac{r}{a}\right)^{n} {q_{n}^{m}}(t)\!\right)\\ &\quad\times{P_{n}^{m}}(\theta) e^{im\phi} + c.\,c.\,, \end{aligned} $$
((5))

where \({g}_{n}^{m}\) and \({q}_{n}^{m}\) are the complex internal and external Gauss coefficients, N is the truncation level, \({P_{n}^{m}}\) is the Schmidt semi-normalized Legendre polynomials of degree n and order m, and ‘ c. c.’ stands for complex conjugate. Then at the Earth’s surface, of radius r=a, one has (see Section 2.1 in Olsen et al. 2010a)

$$ \left\{ \begin{array}{rl} X_{LLS} = &\sum\limits_{n,m} {h}_{n}^{m}(t) \frac{d{P_{n}^{m}}}{d\theta} e^{im\phi} + c.\,c. \\ Y_{LLS} = &\sum\limits_{n,m} {h}_{n}^{m}(t) \frac{i m}{\sin\theta} {P_{n}^{m}} e^{im\phi} + c.\,c.\\ Z_{LLS} = &\sum\limits_{n,m} {v}_{n}^{m}(t) {P_{n}^{m}} e^{im\phi} + c.\,c. \end{array} \right.\,, $$
((6))

with the complex coefficients

$$ \left\{ \begin{array}{rl} {h_{n}^{m}}(t) = & {g}_{n}^{m}(t) + {q_{n}^{m}}(t)\\ {v_{n}^{m}}(t) = & (n+1) {g}_{n}^{m}(t) - n {q_{n}^{m}}(t) \end{array} \right.\,. $$
((7))

Below, \(h_{n}^{mc}\) and \(h_{n}^{ms}\) stand for the real and imaginary part of the complex coefficient \({h_{n}^{m}}\) (with similar notation for \(v_{n}^{mc}\) and \(v_{n}^{ms}\)). We will thus model X and Y data on the one hand, and Z data on the other, resorting to a least-squares fit. We emphasize that in this context, what will be isolated is a combination of global-scale fields that may arise from sources both external (e.g., magnetospheric) and internal (e.g., induced or remaining core signals) to the observation sites.

At each epoch t spanning 1930 to 2010 every month, we store residuals B−(B A R2+B 1+B 2) at all available observatories in two data vectors y h (t)={X(t),Y(t)} and y v (t)={Z(t)}. The number of available observatories (N x ,N y ,N z ) for the three components are a function of time (cf. Figure 1). Then, from Equations 5 and 6, we write at every epoch t:

$$ \left\{ \begin{array}{rl} \textbf{y}_{h} =& {\sf A}_{h} \textbf{h}+\textbf{e}_{h}\\ \textbf{y}_{v} =& {\sf A}_{v} \textbf{v}+\textbf{e}_{v} \end{array} \right.\,. $$
((8))

The unknown vectors h(t) and v(t), of size P = N(N + 2), contain, respectively, the coefficients \({h_{n}^{m}}(t)\) and \({v_{n}^{m}}(t)\), whose statistics are defined by the a priori model covariance matrices R h =E(h h T) and R v =E(v v T). Matrices A h and A v , of size, respectively, (N x +N y P and N z ×P, translate into matrix form the linear relation (Equation 6) between the data vectors (y h and y v ) and the vectors of model coefficients (h and v). Vectors e h and e v stand for the errors in the horizontal and vertical components, whose statistics are given by the covariance matrices \({\sf C}_{h} = E\left (\textbf {e}_{h} \textbf {e}_{h}^{T}\right)\) and \({\sf C}_{v} = E\left (\textbf {e}_{v} \textbf {e}_{v}^{T}\right)\). The Bayesian average solution is then

$$\begin{array}{@{}rcl@{}} \hat{\textbf h} = {\sf H}_{h} {\sf A}_{h}^{T}{\sf C}_{h}^{-1} \textbf{y}_{h}\,, \end{array} $$
((9))

with \({\sf H}_{h} = \left [{\sf A}_{h}^{T}{\sf C}_{h}^{-1} {\sf A}_{h} + {\sf R}_{h}^{-1}\right ]^{-1}\) the covariance matrix defining the a posteriori error statistics of the model h (with similar definitions for v). From \(\hat {\textbf h}\) and \(\hat {\textbf v}\), we build the large length-scale field vector B LLS .

The quantities entering the a priori data error and model error covariance matrices are a priori unknown. For the sake of simplicity, we consider uncorrelated errors with stationary statistics, i.e., \({\sf C}_{h}={e^{2}_{h}}{\sf I}_{h}\) at every epoch (I h stands for the identity matrix of rank N x +N y ), with similar notation for C v . Of course, errors are probably increasing towards the past, when records were less accurate. However, we tested several values for \({e^{2}_{h}}\) and \({e^{2}_{v}}\) and found that within a reasonable range, the final estimate of B LLS , and then of the UMSL, was not significantly affected. The results presented in Sections ‘Main patterns of remaining external and induced fields’ to ‘The UMSL throughout 1930 to 2010’ are obtained using e h =e v =6 nT. In practice, we fix the truncation level to N=3, which results in 15 unknowns per epoch. This is much smaller than the number of observatories after 1960 but not necessarily in the first half of the twentieth century.

If we were retrieving B LLS from Equation 8 in geographic coordinates, the model cross-covariances would vary with time depending on the dipole orientation. Thus, in order to solve Equation 8 using stationary variance properties, we rotate at each epoch the forward equations in dipole coordinates (for which we use the COV-OBS internal dipole field). Thus, we build matrices R h and R v from the variances of the coefficient series \({h_{n}^{m}}(t)\) and \({v_{n}^{m}}(t)\) in dipole coordinates, first obtained over the period 1960 to 2010 for which the model is weakly sensitive to the a priori matrices (we would underestimate the variances of the model parameters by considering the era before 1960, where fewer, less accurate data are available). Coefficients are then rotated back into geographic coordinates.

To finally obtain the UMSL presented in Section ‘The UMSL throughout 1930 to 2010’, we consider (i) the residuals \(\delta \textbf {y}_{h} = \textbf {y}_{h}-{\sf A}_{h} \hat {\textbf h}\) and \(\delta \textbf {y}_{v} = \textbf {y}_{h}-{\sf A}_{v} \hat {\textbf v}\), i.e., the final residual series B−(B A R2+B LLS +B 1+B 2), and (ii) the modeling errors arising from the uncertainties of the large length-scale field parameters. We thus estimate, for each epoch t, the UMSL at all observatories as the square root of the diagonal elements of the covariance matrices

$$\begin{array}{@{}rcl@{}} {\sf K}_{h}(t) = {\sf A}_{h}(t){{\sf H}_{h}(t)}{\sf A}_{h}(t)^{T} + {\delta\textbf{y}_{h}}(t){\delta\textbf{y}_{h}}(t)^{T}\, \end{array} $$
((10))

and K v (t) (with similar notation).

Results and discussion

On the sensitivity of the residuals to the Gaussian process regression

In this section, we test the sensitivity of the variance of the residuals δ B=BB A R2, obtained using the BLUE (Section ‘Local stochastic regression of the main field’), to the a priori noise amplitude that enters the covariance matrix C ee in Equation 1. We first study the case of synthetic data, before we analyze the geophysical series.

We generate monthly values of a noisy AR-2 process over a time window Δ T=50 years. We consider first the case of stationary white noise, varying the a priori noise variance that enters C ee from 1/100 to 100 times the true noise variance. We tested two values for the dimensional r.m.s. of the true noise, 10 and 100 nT, corresponding to a noise over signal ratio of 0.33% and 0.033%, respectively. The ratio between the magnitudes of the residuals to the BLUE and of the true noise, summarized in Table 2, remains close to unity in all cases. Given the range of a priori values, we conclude that we are able to retrieve well the true noise amplitude in this synthetic case. We checked that this result is not sensitive to the window length Δ T. We also consider a situation that mimics a quality jump (such as the advent of proton magnetometers), where the noise magnitude is ten times larger over the first half of the window (case C in Table 2). Again, we correctly retrieve the magnitude of the true noise, over both halves of the time window.

Table 2 Ratio between the true and retrieved noise magnitude when varying the prior noise variance

We now apply the same type of fit to the geophysical series. Examples of residual magnitudes obtained for the X, Y, and Z components for several observatories are displayed in Table 3. We used values for the a priori error variance η (the diagonal elements of the matrix C ee ) of 10, 100, and 1,000 nT 2. Given the standard error bars generally associated with observatory series, by using the smaller (larger) values of η, we supposedly under- (over-) estimate the signal that cannot enter the AR-2 fit. We observe that the a posteriori standard deviation of the residuals is only weakly sensitive to the a priori choice for η in the range of 10 to 1,000 nT 2. In comparison, when using a polynomial fit instead of the BLUE, we obtain (not shown) residual variances that are much more sensitive to the cutoff frequency (the number of polynomials) for both synthetic and geophysical series. We conclude that the method presented in Section ‘Local stochastic regression of the main field’ is relatively robust for estimating the variance of the residuals of the main field regression.

Table 3 Dimensional r.m.s. of the residuals to an AR-2 fit at several observatory sites

As expected, the amplitude of the residuals is weaker for the Y component, which is less affected by the magnetospheric ring current. It is also larger for X in the equatorial area (see for instance Alibag or Guam in Table 3). We observe larger values towards higher geomagnetic colatitudes (for instance at Sitka), where the magnetic fields of ionospheric origin associated with the auroral ring currents are known to produce signals of particularly large amplitude. From now on, we use an a priori noise variance η=36 nT 2. Examples of series of the residuals between the observed field and the AR2 fit are presented in Figures 3, 4, and 5 for all three components at equatorial, mid-, and high latitude sites.

Figure 3
figure 3

Residuals to OMM series and the UMSL: example of a low-latitude site. Top to bottom: OMM series superimposed with AR-2 fit at the Alibag observatory; residuals between OMM series and the AR-2 proxy, superimposed with the LLS field model with and without periodic signals; final residuals after removing the LLS field; our estimate of the UMSL. Left to right: X, Y, and Z components.

Figure 4
figure 4

Residuals to OMM series and the UMSL: example of a mid-latitude site. Same as Figure 3, at Chambon-la-forêt observatory.

Figure 5
figure 5

Residuals to OMM series and the UMSL: example of a high-latitude site. Same as Figure 3, at Qeqertarsuaq (Godhavn) observatory.

Worldwide description of annual and semiannual signals

We concentrate first on the time evolution of the amplitude and phase of the periodic signals. These are displayed in Figure 6, where we filter out the effect of the solar cycle using an 11-year running mean. We find trends similar to those computed from a Fourier analysis of filtered hourly means (see Figures four to five and eight to nine in Le Mouël et al. 2004a). In particular, we retrieve the two large variations, at pluri-decadal time scales, in the amplitude of the semiannual signal. The epochs at which these extrema occur are similar to those found from the analysis of the aa index by Elias et al. (2011). These variations are in phase at all observatories. The phase of the semiannual signal also shows decadal changes, in agreement with those found from aa (see Figure seven in Elias et al. 2011), with small delays from one site to another.

Figure 6
figure 6

Time changes of the amplitude and phase of periodic signals at various sites. Amplitude (top) and phase (bottom) obtained for the annual (right) and semiannual (left) periodic signals for the X component at several observatories, for a window length of δ t=5 years, after filtering with an 11-year running mean. The phase of the periodic signal (originally in radians) is translated into the date at which occurs the first maximum of the annual (or semiannual) fit to the data (for each epoch, a phase of zero would correspond to January 1).

We also see decadal changes in the phase and amplitude of the annual signal. These variations vary depending on the location, the phase showing, for instance, shifts as large as a couple of months. Several origins have been proposed for annual changes, including local effects of induction in the ocean, temperature dependency of the soil magnetization around sensors (Mishima et al. 2013), or spurious temperature effects in the Z component before the introduction of fluxgate sensors (Le Mouël et al. 2004b). Some effects related to measurement issues cannot be discarded from early epochs. Indeed, Figure 6 indicates that annual signals display less dispersion among observatories towards more recent periods.

We also observe differences when comparing our results to those of Le Mouël et al. (2004a), which can be attributed to several factors. We present in Figure 7 the impact of data selection (night only versus all-hours values) and of the window length δ t, for the Sitka and Tucson observatories (some differences may also be related to the filtering strategy used by Le Mouël et al. 2004a). At high geomagnetic latitudes, the amplitude of semiannual oscillations (and to a lesser extent, of annual ones) is much increased by the use of a shorter window length δ t (with weak impact on the phase), whereas this parameter has a minor impact at low to mid-latitude sites. Selecting nighttime data generates a downward (upward), almost constant bias of about 2 nT for the 6 (12) months signal amplitude at mid-latitudes, and a positive phase delay of about 15 days for the annual signal, when the phase of the semiannual signal is not much affected. Compared to all-hours data, using nighttime data also enhances the amplitude of the peaks in the 6-month variation at high latitudes and generates a time-dependent phase shift between zero and about ten days. It also shifts the amplitude of the annual signal by a few nT, due to aliasing, while its phase is not significantly affected.

Figure 7
figure 7

Sensitivity to the data treatment of the amplitude and phase of periodic signals. Amplitude and phase obtained for annual (right) and semiannual (left) periodic signals for the X component at Sitka (57°N, 225°E, four top panels) and Tucson (32°N, 249°E, four bottom panels), for different data selection (OMM constructed from night only data or from all-hours values) and window lengths δ t. The phase of the periodic signal (originally in radians) is translated into the date at which occurs the first maximum of the annual (or semiannual) fit to the data (for each epoch, a phase of zero would correspond to January 1).

The geographical dependence of the amplitude of annual and semiannual periodic signals can essentially be summarized as a function of geomagnetic colatitude, as shown in Figure 8. The amplitude of semiannual changes shows a large maximum at geomagnetic latitudes between 60° and 70° in all three components. It is particularly strong in X, which also displays a secondary maximum towards the geomagnetic equator. In the polar regions, the semiannual signal is larger in Z and weaker in X (note that the amplitude in Y can be significantly larger than that found in X at some sites). Such features could be associated with an increase of the ring currents and the westward auroral electrojet at equinoxes (see Lyatsky and Tan 2003). The latitudinal distributions observed in the X and Z components recall those found by Fujii and Schultz (2002) for the residuals between observatory data and the induction response to a \({P_{1}^{0}}\) source at periods ranging from a few days to a few months (see their Figure nine). These authors associate such features mainly with the complementary effect of auroral ring currents (cf. their Figure fourteen).

Figure 8
figure 8

Latitudinal variation of the annual and semiannual periodic signals amplitude. Amplitude of the annual (bottom) and semiannual (top) periodic signals for all three components of the magnetic field as a function of the geomagnetic latitude, averaged over the period 1960 to 2010, obtained for a window length of δ t=5 years.

The annual signal is generally larger in X, and its amplitude strongly increases at latitudes above 60°. A secondary maximum clearly appears in X at mid-latitudes. Also at mid-latitudes, non-zonal contributions are observed in the X and Y components, which may be associated with regional induction effects, as suggested by Wardinski and Mandea (2006). We note some discrepancies between the latitudinal dependence of annual and semiannual signals (see for instance in X). This may indicate that the source fields responsible for annual changes are different from those proposed by Fujii and Schultz (2002).

Main patterns of remaining external and induced fields

Examples of the contributions from periodic and large length-scale fields to the three components of observatory series are shown, from low to high latitudes, in Figures 3,4 and 5. In the horizontal components, modulations with the solar cycle are observed for the large length-scale field at mid- to low latitudes (more obvious in the X direction), and for periodic signals, in horizontal components at high latitudes. The modulation of periodic signals at longer time scales, highlighted in the previous section, also clearly appears.

We use below the notations \(\tilde {h}_{n}^{m}\) and \(\tilde {v}_{n}^{m}\) for the coefficients in dipole coordinates. We present in Figure 9 their associated (normalized) series, together with the standard deviation obtained for each of the coefficients. In association with the magnetospheric ring current and its induced counterpart, the largest contribution to the X and Y components comes from \(\tilde {h}_{1}^{0}\), while to Z, it is shared between \(\tilde {v}_{1}^{0}\) and \(\tilde {v}_{3}^{0}\). The important contribution of \(\tilde {v}_{3}^{0}\) to the vertical component is associated with the well-known larger impact of externally induced fields on Z (see for instance the discussion of Figure ten in Hulot et al. 2007). Zonal coefficients are generally better resolved than the non-zonal ones, and sectorial coefficients are associated with relatively larger error bars (see the \(\tilde {h}_{2}^{2}\) and \(\tilde {v}_{3}^{3}\) series). We also notice weaker uncertainties from 1960 onward (for instance \(\tilde {h}_{2}^{1c}\), \(\tilde {v}_{1}^{1s}\), \(\tilde {v}_{2}^{1s}\), and \(\tilde {v}_{3}^{2s}\)). The a posteriori errors we obtain exclude the occurrence of large variations for some coefficients at particular epochs (see the \(\tilde {h}_{n}^{0}\), \(\tilde {h}_{1}^{1}\), or \(\tilde {h}_{3}^{1}\) coefficients in the early 1930s, a period of relatively weak solar activity).

Figure 9
figure 9

Series of the LLS field coefficients. Normalized times series of the coefficients \(\tilde {h}_{n}^{m}\) (top) and \(\tilde {v}_{n}^{m}(t)\) (bottom). Error bars (grey-shaded area) are the square root of the diagonal elements of the posterior covariance matrices \({\sf H}_{\sf h}^{-1}\) and \({\sf H}_{\sf v}^{-1}\). Series are normalized by the standard deviation of the mean model coefficients (in black) calculated over 1960 to 2010, whose value (in nT) is given to the right of each series.

The signature of the solar activity is clear for most coefficients (see the K p index in Figure 10). In the 1960s and around 1990, periods of particularly intense solar cycles, we see unusually large fluctuations of coefficients, showing at other epochs a relatively modest dynamics (see the \(\tilde {v}_{1}^{1}\) and \(\tilde {v}_{2}^{1}\) coefficients after 1960). This suggests a particular response of the external (and then induced) fields to the solar flux at those epochs. We note that many observatories started with the International Geophysical Year 1957 to 1958 (cf. Figure 1), some of them having displayed unstable magnetic series. We corrected as much as possible for anomalous data during this period as well (see Section ‘Data’). However, if the remaining outliers partly enter the UMSL, they may also induce biases in the global model coefficients of the least-squares fit (a drawback which may be reduced by using a Huber norm for the measure of residuals, see e.g., Olsen 2002).

Figure 10
figure 10

Projection of the coefficients z( t ) on the first eigen modes, and comparison with the K p index. Top: monthly averages of the K p index (see e.g., Menvielle and Berthelier 1991). Bottom: normalized projection of the coefficient vectors \(\tilde {z}_{n}^{m}\) on the first four eigen modes \({\mathcal {Z}}_{1}\) to \({\mathcal {Z}}_{4}\), which represent 70.6% (top left), 9.4% (top right), 5.2% (bottom left), and 4.5% (bottom right) of the \(\tilde {z}_{n}^{m}\) variances.

The combination of \(\tilde {h}_{n}^{m}\) and \(\tilde {v}_{n}^{m}\) series is used to separate internal and external contributions; we obtain from Equation 7 the complex coefficients

$$ \left\{ \begin{array}{rl} \tilde{g}_{n}^{m} = & \frac{ n\tilde{h}_{n}^{m}+\tilde{v}_{n}^{m}}{ 2n+1}\\ \tilde{q}_{n}^{m} = & \frac{ (n+1)\tilde{h}_{n}^{m}-\tilde{v}_{n}^{m}}{ 2n+1} \end{array} \right. $$
((11))

(here in dipole coordinates). Their analysis as a function of frequency may be used to constrain C responses and the mantle conductivity profile (see e.g., Olsen 1999), which is out of the scope of the present study. We concentrate below on the description of the inducing field. From Equations 6 and 7, we obtain the coefficients describing the vertical component Z from external origin, \(\tilde {z}_{n}^{m}=-n\tilde {q}_{n}^{m}\), which compose the time-dependent vector z(t). From S=[z(1930.0)…z(2010.0)], we build the matrix M=S S T and perform on it a principal component analysis. We show in Figure 11 the vertical component obtained at the Earth’s surface for the first four eigen modes 1 to 4. These account for, respectively, 70.6%, 9.4%, 5.2%, and 4.5% of the variance of z. Most of the variability is carried by 1, which is associated with the axial dipole field (in dipole coordinates), although it is also composed of some non-axisymmetric features. 2 essentially summarizes the role played by the zonal coefficient of degree 3. We present in Figure 10 the projection of the \(\tilde {z}_{n}^{m}\) coefficients onto the first four eigen modes. 3 (with large contributions from \(\tilde {z}_{1}^{1s}\), \(\tilde {z}_{2}^{0c}\), and \(\tilde {z}_{2}^{1s}\) coefficients) and 4 (dominated by a sectorial structure of degree 2) present unusual time changes in the 1960s and around 1990, in connection with the above observation made about the original \(\tilde {v}_{n}^{m}\) and \(\tilde {h}_{n}^{m}\) series. This illustrates that the structure of the inducing field is more complex than the ring current, as discussed for instance by Olsen (1998) for periods shorter than 30 days, or by Fujii and Schultz (2002) for periods ranging from 5 to 106.7 days.

Figure 11
figure 11

Maps of the Z component at the Earth’s surface for the first eigen modes. Maps (centered on the Greenwich meridian) of the Z component at the Earth’s surface for the first four normalized eigen modes, \({\mathcal {Z}}_{1}\) to \({\mathcal {Z}}_{4}\), obtained from the \(\tilde {z}_{n}^{m}\) coefficient series in dipole coordinates. They represent 70.6% (top left), 9.4% (top right), 5.2% (bottom left), and 4.5% (bottom right) of the \(\tilde {z}_{n}^{m}\) variances. For each mode, the color scale is saturated at the maximum value of the mode.

The UMSL throughout 1930 to 2010

We finally focus on the evolution, in space and time, of the UMSL as defined by Equation 10. Two main sources contribute to the UMSL: instrumental noise and unmodeled regional signals. Since the number of available observatories becomes limited for retrieving large-scale external coefficients from the very early twentieth century, we present our UMSL estimate from 1930 onward. Examples of the final residuals and of the obtained UMSL series are presented in Figures 3,4 and 5.

In several observatories (see for instance the series for Alibag or Chambon-la-Forêt in Figures 3 and 4), periodic signals in the Z component tend to decrease towards recent epochs, which coincides with a decrease of the UMSL in Z. A link to instrumental effects is suspected, such as the introduction of proton magnetometers for absolute measurements after 1960 (Turner et al. 2007) and of fluxgate sensors for vector variometry, as suggested by Le Mouël et al. (2004b). Indeed, measurement improvements have been recorded by Chulliat and Telali (2007) at this particular time at Alibag and Chambon-la-Forêt. This is also illustrated by Figure four in Gillet et al. (2010) and Figure fourteen in Gillet et al. (2009), where an important decrease is observed for the scatter in the difference between two Z series at nearby observatories (the difference removing the effect of external fields that are coherent at the two sites).

The three components of the UMSL are presented in Figure 12 as a function of geomagnetic colatitude in 1996 (low solar activity) and 2001 (high solar activity). The most important pattern is the higher UMSL values found in (geomagnetic) polar regions, with a sharp transition between moderate and large values around 55° geomagnetic latitude. We notice only a weak dependence with longitude. The global maps of the UMSL in 1996 and 2001, presented in Figure 13, illustrate with more detail the geographical distribution of the final residual level.

Figure 12
figure 12

Latitudinal dependence of the UMSL. Amplitude of the UMSL in 1996 (a, solar minimum) and 2001 (b, solar maximum) for all three components of the magnetic field, as a function of geomagnetic latitude.

Figure 13
figure 13

Global distribution of the UMSL at observatories. From left to right: X, Y and Z components, in geographic frame; at solar minimum (1996, top) and at solar maximum (2001, bottom).

We encapsulate in Figure 14 the UMSL obtained for the three components at all observatories throughout the selected time span, separately for sites with geomagnetic latitudes lower and higher than 60°. We see in the general trend a clear signature of the solar cycle at high latitudes. It is also there, to a lesser extent, at lower latitudes. On average, the UMSL from the low to mid-latitude sites is similar for all three geomagnetic components (a couple of nT towards the most recent epochs), even though perturbations may be larger at some periods and locations for Z. Depending on the observatory, it ranges between 1 and 10 (6) nT for Z (X and Y) data. We notice a slight decrease of the median UMSL in Z after 1960, while the UMSL in X and Y is rather stationary.

Figure 14
figure 14

Time evolution of the UMSL for all observatories. From left to right: X, Y, and Z components for all observatories (grey dots), with median values over all observatories (black), and 11-year running means of the median (red); (top) low to mid and (bottom) high geomagnetic latitudes.

At latitudes higher than 60°, a significant reduction in the UMSL after 1960 is found for Z data (on average, from 10 to 6 nT). In the series of the horizontal components, the median UMSL is more stationary (around 7 and 4 nT for X and Y, respectively). The decrease observed on average before 1950 is essentially due to the fact that we only have access to a restricted group of sites from this period, which does not include the observatories that display the higher UMSL values later on. Note that some sites display a UMSL as large as 20 nT in X at periods of maximum solar activity, a value much larger than found for Y and Z at recent epochs.

Conclusions

We present in this study a sequential analysis of OMM series, where we isolated three contributions:

  1. 1.

    A local, stochastic regression for the main field, constructed in the framework of Gaussian processes. We use an a priori function based on a −4 slope temporal power spectral density, observed at periods longer than 5 years and extrapolated for shorter periods. This method appears robust, in the sense that it is only weakly sensitive to the a priori data errors that enter the definition of the BLUE in Equation 1.

  2. 2.

    A local analysis of annual and semiannual signals. We retrieve their modulation over long time scales, previously put forward by Le Mouël et al. (2004a,b). We suspect that part of these slow variations in amplitude is instrumental in origin, especially in the Z component. The geographic distribution of the amplitude of periodic signals is essentially zonal in geomagnetic coordinates, presenting a sharp gradient in the auroral zone. We see only a little longitudinal dependence, possibly due to regional induction effects.

  3. 3.

    A global spherical harmonic analysis of the remaining large length-scale fields originating from external and induced sources. In order to avoid the ambiguity between internal and external sources, we fit horizontal and vertical components separately, before reconstructing the internal and external coefficients. We then describe the main patterns of the inducing field in Z, showing a non-negligible role of components that are not carried by \({P_{1}^{0}}\). We observe an unusual morphology of this field in the 1960s and around 1990, periods of particularly intense solar activity.

We end up with the final residual series, from which we define the UMSL. At recent epochs, it is mainly a function of the geomagnetic colatitude, and its magnitude is rather similar for all three components (at least at mid- to low latitudes). This reflects the standards required to enter the INTERMAGNET network. If some UMSL variations denote changes in the quality of the instrumentation, as clearly seen in the Z component around 1960, part of the UMSL fluctuations is related to the solar activity. This highlights the difficulty of modeling globally the complex magnetic fields at high latitudes from monthly values. A strategy similar to that presented here could be followed for hourly values if one were to also model Sq variations (see Stening and Winch 2013). This would not only have implications for the analysis of annual and semiannual signals; indeed, it would also require accounting for periodic fields at higher frequencies (see Love and Rigler 2014).

The UMSL could enter the data error covariance matrix in future field modeling studies if one wishes to invert jointly the internal, induced, and external signals at periods shorter than a few years. Indeed, modeling internal and external coefficients through the \(\tilde {h}_{n}^{m}\) and \(\tilde {v}_{n}^{m}\) associated, respectively, with horizontal and vertical data, it is possible to constrain the mantle conductivity through C responses (Olsen 1999). Furthermore, it seems possible that some core signals may reach the Earth’s surface at short periods (Velímský and Finlay 2011), making the joint inversion of core and induced signals relevant. Since the signature of core flows are relatively weak at short periods, accurate estimates of error covariances are required to assess the reliability of the reconstructed motions. In this situation, it may be of importance to consider the UMSL time changes at every site, since accurate error statistics are required to best extract the information contained in OMM series. For this purpose, one may consider the strategy put forward by Haines (1993) when building data error covariance matrices for differentiated series.

Endnote

a Note that the SDE (3) differs from that defined by Equation (12) in (Gillet et al. 2013). Both SDEs define processes with the same auto-covariance function (2), but their Equation (12), contrary to our Equation (3), corresponds to a non-stationary process. This does not affect their results that have been computed from the covariance function and not the SDE.

References

  • Bloxham, J, Jackson A (1992) Time-dependent mapping of the magnetic field at the core-mantle boundary. J Geophys Res 97: 19,537–63.

    Article  Google Scholar 

  • Chapman, S, Bartels J (1940) Geomagnetism, Vol. 2. Clarendon Press, Oxford.

    Google Scholar 

  • Chulliat, A, Telali K (2007) World monthly means database project. Publs Inst Geophys Pol Acad ScC-99(398): 537–552.

    Google Scholar 

  • Cliver, E, Kamide Y, Ling A (2000) Mountains versus valleys: semiannual variation of geomagnetic activity. J Geophys Res 105(A2): 2413–2424.

    Article  Google Scholar 

  • Currie, RG (1968) Geomagnetic spectrum of internal origin and lower mantle conductivity. J Geophys Res 73: 2779–2786.

    Article  Google Scholar 

  • De Santis, A, Barraclough D, Tozzi R (2003) Spatial and temporal spectra of the geomagnetic field and their scaling properties. Phys Earth Planet Int 135(2): 125–134.

    Article  Google Scholar 

  • Elias, AG, Silbergleit VM, de Gonzalez ALC (2011) Long-term variation of the semi-annual amplitude in the aa index. J Atmos Sol Terr Phys 73: 1492–1499.

    Article  Google Scholar 

  • Finlay, CC, Jackson A, Gillet N, Olsen N (2012) Core surface magnetic field evolution 2000–2010. Geophys J Int 189: 761–781. doi:10.1111/j.1365-246X.2012.05395.x.

    Article  Google Scholar 

  • Fujii, I, Schultz A (2002) The 3d electromagnetic response of the earth to ring current and auroral oval excitation. Geophys J Int 151(3): 689–709.

    Article  Google Scholar 

  • Gillet, N, Pais MA, Jault D (2009) Ensemble inversion of time-dependent core flow models. Geochem Geophys Geosyst 10(6): Q06,004. doi:10.1029/2008GC002,290.

    Article  Google Scholar 

  • Gillet, N, Lesur V, Olsen N (2010) Geomagnetic core field secular variation models. Space Sci Rev 155: 129–145. doi:10.1007/s11,214-009-9586-6.

    Article  Google Scholar 

  • Gillet, N, Jault D, Finlay CC, Olsen N (2013) Stochastic modelling of the earth’s magnetic field: inversion for covariances over the observatory era. Geochem Geophys Geosyst 14(4): 766–786. doi:10.1029/2012GC004355.

    Google Scholar 

  • Haines, GV (1993) Modelling geomagnetic secular variation by main-field differences. Geophys J Int 114(3): 490–500. doi:10.1111/j.1365-246X.1993.tb06982.x.

    Article  Google Scholar 

  • Hulot, G, Olsen N, Sabaka TJ (2007) The present field In: In: Kono M, Schubert G (eds) Treatise in Geophysics, Geomagnetism, 33–75.. Elsevier, New York.

    Google Scholar 

  • Jackson, A, Jonkers ART, Walker MR (2000) Four centuries of geomagnetic secular variation from historical records. Phil Trans R Soc Lond A 358: 957–990.

    Article  Google Scholar 

  • Le Mouël, J, Blanter E, Chulliat A, Shnirman M (2004a) On the semiannual and annual variations of geomagnetic activity and components. Ann Geophys 22: 3583–3588.

    Article  Google Scholar 

  • Le Mouël, J, Blanter E, Shnirman M (2004b) The six-month line in geomagnetic long series. Ann Geophys 22: 985–992.

    Article  Google Scholar 

  • Lesur, V, Wardinski I, Asari S, Minchev B, Mandea M (2010) Modelling the Earth’s core magnetic field under flow constraints. Earth Planets Space 62: 503–516.

    Article  Google Scholar 

  • Liu, L, Zhao B, Wan W, Ning B, Zhang ML, He M (2009) Seasonal variations of the ionospheric electron densities retrieved from constellation observing system for meteorology, ionosphere, and climate mission radio occultation measurements. J Geophys Res: Space Phys 114(A02302). doi:10.1029/2008JA013819.

  • Love, JJ, Chulliat A (2013) An international network of magnetic observatories. Eos, Trans Am Geophys Union 94(42): 373–374.

    Article  Google Scholar 

  • Love, JJ, Rigler EJ (2014) The magnetic tides of Honolulu. Geophys J Int 197(3): 1335–1353.

    Article  Google Scholar 

  • Lühr, H, Maus S (2010) Solar cycle dependence of quiet-time magnetospheric currents and a model of their near-Earth magnetic field. Earth Planets Space 62: 843–848.

    Article  Google Scholar 

  • Lyatsky, W, Tan A (2003) Latitudinal effect in semiannual variation of geomagnetic activity. J Geophys Res: Space Phys 108(A5): 1978–2012.

    Article  Google Scholar 

  • Lyatsky, W, Newell P, Hamza A (2001) Solar illumination as cause of the equinoctial preference for geomagnetic activity. Geophys Res Lett 28(12): 2353–2356.

    Article  Google Scholar 

  • Malin, S, Isikara AM (1976) Annual variation of the geomagnetic field. Geophys J R Astron Soc 47(3): 445–457.

    Article  Google Scholar 

  • Malin, S, Winch D (1996) Annual variation of the geomagnetic field. Geophys J Int 124(1): 170–174.

    Article  Google Scholar 

  • Matzka, J, Chulliat A, Mandea M, Finlay C, Qamili E (2010) Geomagnetic observations for main field studies: from ground to space. Space Sci Rev 155: 29–64. doi:10.1007/s11214-010-9693-4.

    Article  Google Scholar 

  • McLeod, MG (1996) Spatial and temporal power spectra of the geomagnetic field. J Geophys Res 101: 2745–2763.

    Article  Google Scholar 

  • Menvielle, M, Berthelier A (1991) The k-derived planetary indices: description and availability. Rev Geophys 29(3): 415–432.

    Article  Google Scholar 

  • Mishima, T, Owada T, Moriyama T, Ishida N, Takahashi K, Nagamachi S, Yoshitake Y, Minamoto Y, Muromatsu F, Toyodome S (2013) Relevance of magnetic properties of soil in the magnetic observatories to geomagnetic observation. Earth Planets Space 65(4): 337–342.

    Article  Google Scholar 

  • Olsen, N (1998) The electrical conductivity of the mantle beneath Europe derived from c-responses from 3 to 720 hr. Geophys J Int 133(2): 298–308.

    Article  Google Scholar 

  • Olsen, N (1999) Long-period (30 days–1 year) electromagnetic sounding and the electrical conductivity of the lower mantle beneath Europe. Geophys J Int 138(1): 179–187.

    Article  Google Scholar 

  • Olsen, N (2002) A model of the geomagnetic field and its secular variation for epoch 2000 estimated from ørsted data. Geophys J Int 149(2): 454–462.

    Article  Google Scholar 

  • Olsen, N, Lowes F, Sabaka TJ (2005a) Ionospheric and induced field leakage in geomagnetic field models, and derivation of candidate models for DGRF 1995 and DGRF 2000. Earth Planets Space 57: 1191–1196.

    Article  Google Scholar 

  • Olsen, N, Sabaka TJ, Lowes F (2005b) New parameterization of external and induced fields in geomagnetic field modeling, and a candidate model for igrf 2005. Earth Planets Space 57(12): 1141–1149.

    Article  Google Scholar 

  • Olsen, N, Glassmeier KH, Jia X (2010a) Separation of the magnetic field into external and internal parts. Space Sci Rev 152(1-4): 135–157.

    Article  Google Scholar 

  • Olsen, N, Mandea M, Sabaka T (2010b) Tøffner-Clausen L. The CHAOS-3 geomagnetic field model and candidates for the 11th generation IGRF. Earth Planets Space 62: 719–727.

    Google Scholar 

  • Olsen, N, Luhr H, Finlay CC, Sabaka TJ, Michaelis I, Rauberg J (2014) Tøffner-Clausen, L. The CHAOS-4 geomagnetic field model. Geophys J Int 197(2): 815–827.

    Google Scholar 

  • Rasmussen, CE, Williams CKI (2006) Gaussian processes for machine learning. The MIT Press, Cambridge.

    Google Scholar 

  • Russell, C, McPherron R (1973) Semiannual variation of geomagnetic activity. J Geophys Res 78(1): 92–108.

    Article  Google Scholar 

  • Sabaka, T, Langel R, Baldwin R, Conrad J (1997) The geomagnetic field 1900-1995, including the large-scale field from magnetospheric sources, and the NASA candidate models for the 1995 revision of the IGRF. J Geomag Geoelectr 49: 157–206.

    Article  Google Scholar 

  • Sabaka, TJ, Olsen N, Langel RA (2002) A comprehensive model of the quiet-time, near-earth magnetic field: phase 3. Geophys J Int 151(1): 32–68. doi:10.1046/j.1365-246X.2002.01774.x.

    Article  Google Scholar 

  • Sabaka, TJ, Olsen N, Purucker ME (2004) Extending comprehensive models of the Earth’s magnetic field with Oersted and CHAMP data. Geophys J Int 159: 521–547.

    Article  Google Scholar 

  • Stening, R, Winch D (2013) The ionospheric sq current system obtained by spherical harmonic analysis. J Geophys Res: Space Phys 118(3): 1288–1297.

    Article  Google Scholar 

  • Turner, GM, Rasson JL, Reeves CV (2007) Observation and measurement techniques In: Kono M, Schubert G (eds) Treatise in Geophysics, 93–146, Geomagnetism.

  • Velímský, J, Finlay CC (2011) Effect of a metallic core on transient geomagnetic induction. Geochem Geophys Geosyst 12(5): Q05,011. doi:10.1029/2011GC003557.

    Article  Google Scholar 

  • Wardinski, I, Holme R (2006) A time-dependent model of the Earth’s magnetic field and its secular variation for the period 1980–2000. J Geophys Res 111. doi:10.1029/2006JB004,401.

  • Wardinski, I, Holme R (2011) Signal from noise in geomagnetic field modelling: denoising data for secular variation studies. Geophys J Int 185: 653–662.

    Article  Google Scholar 

  • Wardinski, I, Lesur V (2012) An extended version of the C 3FM geomagnetic field model – application of a continuous frozen-flux constraint. Geophys J Int 189: 1409–1429. doi:10.1111/j.1365-246X.2012.05384.x.

    Article  Google Scholar 

  • Wardinski, I, Mandea M (2006) Annual and semi-annual variations of the geomagnetic field components analysed by the multi-taper method. Earth Planes Space 58(6): 785.

    Article  Google Scholar 

  • Yaglom, AM (1962) An introduction to the theory of stationary random functions. Prentice-Hall, Upper Saddle River.

    Google Scholar 

Download references

Acknowledgements

We thank the national institutions that support ground magnetic observatories and INTERMAGNET for promoting high standards of practice and making the data available. Jiaming Ou’s two visits to France were financed by the Institute of Geology and Geophysics, Chinese Academy of Sciences. Discussions with Dominique Jault, Nils Olsen, Chris Finlay, and Jurgen Matzka have helped us throughout the evolution of the present study. We thank the two anonymous referees, whose comments helped improve the quality of the manuscript. This work was supported in part by the French Centre National d’Etudes Spatiales (CNES) for the preparation of the Swarm mission of ESA. This work has been partially supported by the French ‘Agence Nationale de la Recherche’ under grant ANR-11-BS56-011. This work is also supported by the National Basic Research Program of China (2014CB845903 and 2012CB825604) and the National Natural Science Foundation of China (41174122, 41031066, 41104093).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Nicolas Gillet.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

JO carried out all the cleaning of geomagnetic series from hourly values. He performed numerically the several steps of the sequential analysis and proposed to consider separately periodic signals. NG supervised the work of JO at ISTerre. He proposed to use the Gaussian processes framework and the global modeling of residuals in terms of horizontal/vertical components to estimate external fields coefficients. NG also carried out the principal components analysis on the external field coefficients and led the writing of the paper. AD supervised the work of JO at IGGCAS, introduced him to geomagnetic observatory records, and made geomagnetic data available for this study. All authors read and approved the final manuscript.

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ou, J., Gillet, N. & Du, A. Geomagnetic observatory monthly series, 1930 to 2010: empirical analysis and unmodeled signal estimation. Earth Planet Sp 67, 13 (2015). https://doi.org/10.1186/s40623-014-0173-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40623-014-0173-z

Keywords