Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth's core surface

1.1 On the Resolution of Core Motions

Some consensus has emerged about the possibility of estimating core surface flows u from the secular variation (defined as the time variation of the magnetic field and henceforth abbreviated to SV) and about their geometry for the most recent epochs. This picture is obtained by inverting for u using the radial component of induction equation at the core surface,

(1)

(B_r is the radial magnetic field, ∇_H· is the horizontal divergence operator). Here magnetic diffusion has been neglected. While stressing the significance of magnetic diffusion as a source for SV, Holme and Olsen [2006] concluded that it does not contribute a larger fraction of the high-degree SV than of the low-degree SV. A recent study based on geodynamo modeling supports this statement [Aubert, 2014a]. Meanwhile, the SV associated with the unresolved part of the radial magnetic field has been diagnosed as the main source of uncertainty in single-epoch analyses of equation 1 [Eymin and Hulot, 2005; Pais and Jault, 2008; Baerenzung et al., 2014]. Consequently, our ability to reconstruct core flows is reduced.

There is nevertheless growing evidence of the equatorial symmetry (ES) of large-scale core surface flows (but see Whaler and Beggan [2015] for an alternative view). Gillet et al. [2011] find an increase of the percentage of ES flow as better SV data are available. Baerenzung et al. [2014] observe from satellite SV models that without imposing any topological constraints, more than 80% of the kinetic energy of core surface motions for the single epoch 2005.0 is stored into their ES component [see also Wardinski et al., 2008]. Furthermore, Aubert et al. [2013] have recently shown that the eccentric gyre proposed by Pais and Jault [2008] can also emerge from three-dimensional geodynamo simulations as a columnar structure that persists over centuries. With this in mind, flows have been continued within the core assuming quasi-geostrophy [Pais and Jault, 2008; Aubert, 2013]. Our knowledge of time variations of the flow is less advanced, partly because time series of satellite data are relatively short.

Early calculations of time-dependent flow models were largely conducted in order to estimate the time changes of core angular momentum [e.g., Jackson, 1997; Pais and Hulot, 2000]. Domination of the steady component of the large length-scale flow in the kinetic energy budget was suggested by Bloxham [1988]. Conversely, emergence of large-scale vortices on 50 year time scales [Amit and Olson, 2006] as well as sudden local accelerations over a few months have been reported [Olsen and Mandea, 2008]. Interpretations of core flow variations have, however, been hampered because most previous studies assumed that the magnetic field entering the forward problem was perfectly known and relied on ad hoc regularizations of the core surface velocity. A first attempt at the estimation of model uncertainties (on which we try to improve here) gave confidence in the detection of fast torsional waves [Gillet et al., 2010] (see also Asari and Lesur [2011] for a description of these waves from 2000 to 2010). Here, we attempt to go further using more consistent forms of prior information concerning the core flow.

1.2 Deterministic Versus Stochastic Modeling

One attractive solution to prescribing prior information is to use a dynamical model, i.e., resorting to geomagnetic data assimilation [Fournier et al., 2010]. Relying on the dynamics of a 3-D numerical geodynamo simulation, Fournier et al. [2011] and Aubert and Fournier [2011] made inferences about the core surface flow. Dynamics entered their linear estimation in the form of covariances between the velocity and the magnetic field, the interior and the surface, all at a single epoch. Aubert [2013] remarked thereafter that linear estimation may have limited applicability because it works well only if the initial guess is not too far from the final solution. Hence, he resorted to a classical frozen-flux inversion of the core surface flow and used the first and second statistical moments obtained from geodynamo simulations to build prior information about the core surface flow and to determine the error covariance matrix, including contributions from model uncertainties that arise as a consequence of truncations. Along the same lines, Aubert [2014a, 2014b] has recently presented images of core flows using statistics provided by the magnetic field model COV-OBS [Gillet et al., 2013], which is also our source of data here. Remarkably, these calculations do not require ad hoc penalization of small-scale flows, which is an awkward feature of classical core surface flow estimation. On the other hand, the important issue of temporal correlations in unresolved scales is not addressed. Furthermore, the present generation of three-dimensional geodynamo simulations probably lacks part of the short time scale variability observed in geomagnetic series [Aubert, 2014a].

As an alternative, we propose here that a stochastic framework can be employed, following a strategy similar to that used during the derivation of the COV-OBS magnetic field model [Gillet et al., 2013]. The radial magnetic field at the core surface B_r(t) was then considered to be a realization of a stationary stochastic (i.e., random) process ξ(t). Knowledge about the spectral density of B_r(t) (in the frequency interval ) was translated into prior information about the time correlation properties of the process ξ sampled by B_r via a correlation function of the form:

(2)

where τ₀ is a characteristic time. The process ξ with correlation function 2 is the solution of the stochastic differential equation [Yaglom, 1962, equations (2.153) and (2.155')]

(3)

where ζ(t) is the Brownian motion (or Wiener) process. Note that equation 3 corrects the equation (12) of Gillet et al. [2013], which is wrong and corresponds to a process that is not stationary (the mistake in writing this equation did not affect the results presented by Gillet et al. [2013] because the expression for their correlation function (2) was correct).

Here we incorporate preexisting knowledge about the time evolution of the core surface velocity field u via a time correlation function and consider the time series of the flow coefficients to be realizations of a stochastic process. Our specific choice of process for u, which is defined by the stochastic differential equation

(4)

follows from our selection of 3 as the process for B_r, assuming the same differentiability properties for u and ∂B_r/∂t. This choice of prior stochastic process implies that the flow coefficient series are continuous, albeit nondifferentiable (their increments during a time interval of length τ are defined but are not proportional to τ as τ is decreased within the range of time scales that we consider). Realizations of the process 4 have a simple exponential time correlation function of the form

(5)

where τ_u is the characteristic time scale for flow changes. The crucial point here is that the choice of prior for u impacts the calculation of SV model errors arising due to unresolved small scales.

1.3 Ensemble Calculation of SV Model Errors

SV model errors arise when estimating core surface flows because the radial magnetic field that enters the right-hand side of 1 is an uncertain parameter [Jackson, 1995]. These SV model errors should be carefully distinguished from the SV observation errors that are provided by COV-OBS in our case. The foremost example of this difficulty is our ignorance of the small-scale core surface magnetic field B_r with harmonic degree above 14. Workers have traditionally circumvented this problem by searching for flows u presenting rapidly converging spectra [e.g., Hulot et al., 1992], for which the interaction with the small-scale magnetic field is artificially reduced. This strategy was tempting as long as only the very largest scales of the secular variation were known (up to degree about 8 from ground observations). However, magnetic field models obtained from satellite data now display resolved SV up to degree 14 [e.g., Olsen et al., 2010; Lesur et al., 2010; Finlay et al., 2012]. This observation has made it necessary to relax the large-scale flow hypothesis, placing the SV model error problem foremost. Eymin and Hulot [2005] stressed this point and dealt with the problem in the framework of regularized single-epoch core surface flow inversions, tuning the trade-off parameter such that the SV residuals have similar amplitude to the estimated SV modeling errors. Pais and Jault [2008] later recursively estimated a diagonal approximation (i.e., for a single epoch) of the SV model error covariance matrix, with elements being a function of the harmonic degree only, adding this to the observation error covariance matrix to obtain the required covariance matrix [see also Aubert, 2013]. Alternatively, Baerenzung et al. [2014] represented the subgrid processes in the induction equation 1 as a function of the large-scale fields. However, none of these procedures account for time correlations of model uncertainties, which turns out to be essential if one wishes to resolve flow time variations. Consideration of such temporal correlations is a crucial issue because of the high coherence, on decadal time scales, of both the flow and the unresolved magnetic field that enter the induction equation 1. Neglecting the temporal correlation of SV model errors should actually result in an underweighting of the information in geomagnetic data concerning the rapid flow changes and may lead to a biased solution for the core flow.

Here we adopt a new ensemble estimation of SV model errors, taking due account of their time correlation. Gillet et al. [2009] made a first attempt at an ensemble estimation of time-dependent surface core flows from a magnetic field model and its covariance. They generated an ensemble of core surface magnetic fields (including its unresolved component at small length scales), substituted these for B_r in equation 1, and estimated an ensemble of flow solutions from SV data associated with observation errors only. The most probable solution was then calculated as the ensemble average of flows . Considering each realization of the unresolved field as perfectly known when inverting equation 1, SV model errors were estimated a posteriori once and for all. This assumption keeps the inverse problem linear, thus avoiding an iterative estimate of SV model errors as in Pais and Jault [2008]. However, there was a major drawback in the approach of Gillet et al. [2009]: substituting for u in equation 1 systematically underpredicted the SV monitored in ground observatories even though each flow member u^p of the ensemble adequately predicted the observed SV. The cause of this deficiency was that the SV error covariance used during the estimation of the flow contained only the contribution from observation errors. The correct way to calculate ensemble statistics is to generate an ensemble of forward models and to build in a consistent fashion the error covariance matrix by adding both the SV model error and observational error covariances [Evensen, 2009]. In agreement with this statement, Baerenzung et al. [2014] remark in their synthetic tests that omitting SV model error covariances, as in the study of Gillet et al. [2009], causes a degradation in performance compared to the iterative scheme initiated by Pais and Jault [2008]. Note that SV model error covariance matrix has to be calculated recursively because it involves the quantity u that is estimated. Furthermore, in order to properly account for time correlation in the SV model error, we must estimate the flow over the entire time interval in a single batch. For the SV observation error covariance matrix, we rely on the COV-OBS magnetic field model.

Section 2 is devoted to methodology and constitutes the heart of the present work. We develop a time-dependent ensemble approach to calculate a covariance matrix for the SV model errors, following the conclusions of the above discussion. We describe how this transforms the kinematic inversion of core flows into a nonlinear problem. We furthermore extend the ensemble method to the calculation of the covariance matrix for the uncertainties in the model parameters. It is at this stage that prior knowledge about time-correlations is incorporated. In section 3 we analyze the predictions of the resulting flow for geophysical observations, namely, length-of-day and SV changes at observatories. Then, the geometry and the time dependence of the flow are discussed (section 4), together with their uncertainties. This section ends with considerations about torsional waves, electromagnetic core-mantle coupling, and a focus on the dynamics in the equatorial region. We conclude this paper (section 5) with discussions of possible methodological improvements and perspectives concerning core physics.

2.1 Notations

We expand all quantities (B_r, ∂B_r/∂t, and the poloidal and toroidal scalars S and T describing the surface flow u=∇×(Tr)+∇_H(rS), where r is the position vector) in spherical harmonics:

(6)

((θ,φ) are spherical coordinates, Y_ℓm are Schmidt quasi-normalized functions, and is the complex conjugate of x). The expansions of the magnetic field and its time derivative are truncated, respectively, at degrees L_b (for B_r) and L_y (for ∂B_r/∂t). Contrary to Jackson [1997] and Gillet et al. [2009], we do not expand in time in terms of a B-spline basis and consider instead the spherical harmonic coefficients of the time variable fields as discrete time parameter sets regularly sampling every δt = (t_N−t₁)/(N − 1) = 1 year (N = 71). At each epoch t_n, we store B_r, ∂B_r/∂t, and core flow spherical harmonic coefficients in vectors b(t_n), y(t_n), and x(t_n), respectively. From these we build vectors , , and , which are linked through the forward problem:

(7)

with the SV error vector and the operator

(8)

H(b(t_n)) results from the transformation of the snapshot equation 1 at epoch t_n in matrix form.

2.2 Formulation of the Quasi-Geostrophic Topological Constraint

We assume quasi-geostrophy and incompressibility in the outer core volume, which results in the columnar flow constraint at the core surface [Amit and Olson, 2004; Amit and Pais, 2013],

(9)

together with the equatorial symmetry constraint [Pais and Jault, 2008, equation (27)]. In contrast with our previous attempts at calculating core flows [Pais and Jault, 2008; Gillet et al., 2009], we do not impose no penetration across the cylindrical surface tangent to the solid inner core. We note below that our method makes it easy to incorporate this constraint at a later stage if desired.

The above hypotheses yield a set of linear constraints that can be formally written as

(10)

Following Jackson [1997], we use the matrix G obtained from the QR decomposition of Q in order to project the vector of unknowns onto a reduced basis:

(11)

The solution for which we invert is composed of vectors whose size, L_x(L_x+1)/2 for L_x even, is reduced by a factor about 4 compared to the size 2L_x(L_x+2) of the original unknown vectors.

2.3 Ensemble Core Flow Estimation

We use an ensemble approach, following the example of previous studies in oceanic and atmospheric dynamics [Evensen, 2003], to recursively estimate stationary second-order statistics for the flow coefficients and the model errors. We derive an ensemble of P solutions (typically here P = 20) to the forward problem 7 under the constraint 10 from the ensemble of replications drawn from COV-OBS. Preexisting knowledge about w is represented by the covariance matrix C_w, while SV error covariances are described by C_e. Then, for each replication , the least squares solution W^p minimizes the cost function

(12)

Here is the forward operator of equation 8 rotated into the reduced basis using 11. Rather than define and fix the two matrices C_e and C_w prior to the inversion, here we update these at each iteration using the current ensemble of flow solutions . Thus, the minimization of the functional 12 becomes a nonlinear problem, and we calculate the ensemble of models recursively, with at each iteration k:

(13)

We keep the same ensemble of replications (B^p,Y^p) from one iteration to the next.

Concerning the construction of , we estimate the a priori covariances between the coefficients w_i(t_n) and at iteration k + 1 from time and ensemble averages of the variances of the coefficients of the flow models , as calculated at iteration k. In addition, for the required temporal correlation function, we simply adopt the exponential function 5. Combining these gives

(14)

Here the notation means ensemble averaging:

(15)

The choice of a common value for all coefficients of the characteristic flow time scale τ_u is made for the sake of simplicity. Spatial cross covariances (except those associated with the QG assumption) are ignored. To summarize, equation 14 prescribes the elements of the flow model covariance matrix that is used as prior information at iteration k + 1. In order to initialize the calculation, is built from a flat core-mantle boundary (CMB) spatial power spectrum (∀m∈[−ℓ,ℓ],E(t_ℓm)²=10²/ℓ(ℓ + 1)) and the time autocorrelation function 5. The final flow solutions are found to be insensitive to small changes of .

Use of equation 14 means that we build the a priori probability distribution of the flow from the expected values obtained from the existing ensemble of flow models. This approximation is made to reduce the numerical cost; indeed, one should in principle also account for the posterior uncertainties associated with each flow realization [see Baerenzung et al., 2014]. Tests on problems of small dimensions demonstrate that with this simplification we only slightly underestimate the posterior uncertainties on the flow model and the related SV model errors.

We have no prior knowledge of the covariance matrices that enter the cost function 12 to be minimized. The iterative process yields (i) the flow model, (ii) the SV model errors (as described in section 2.4), and (iii) the flow covariances. We seek to exhibit possible solutions within this framework. On the other hand, the inversion process is inherently nonlinear, and we do not claim unicity of the flow solution.

2.4 Accounting for Time-Correlated SV Model Errors

The SV model error covariance matrix is also estimated recursively using an ensemble approach. Recall that the error vector E in equation 7 contains both the SV observation errors E^o and the SV model errors E^m that arise from unresolved interactions between the core flow and the magnetic field. SV model errors E^m result from the unresolved flow δW interacting with the entire (resolved or not) magnetic field, plus the resolved flow interacting with the unresolved field δB. We omit model errors at the first iteration. We calculate the resolved flow and magnetic field as the ensemble averages of the flow solutions, (calculated at iteration k), and of the COV-OBS field models, , respectively.

An ensemble of Q = 40 realizations of δB is obtained from the product of normally distributed random vectors with the Choleski decomposition of the covariance matrix . For ℓ≤14, we simply take the COV-OBS error covariance matrix as C_b. Note that these are potentially underestimated at degrees ℓ ∼ 14 due to the signature of unmodeled lithospheric field at large length scales [see Jackson, 1990; Thébault and Vervelidou, 2015]. For ℓ > 14, C_b is constructed from the correlation function 2 proposed in Gillet et al. [2013], using correlation times and variances extrapolated with power laws from those obtained in satellite field models for ℓ≤12. We truncate the unresolved field at L_b=30 since increasing L_b further does not significantly affect the results. The ensemble of Q realizations is calculated independently for the P replications B^p that are used in 13. The dispersion within the ensemble of flow solutions defines the unresolved flow, i.e., . It should be accounted for when estimating the covariances of SV model errors because these depend on the flow (the reason that our problem is nonlinear, see Pais and Jault [2008]). We thus estimate PQ = 800 realizations of the SV errors from all possible pairs (δW^p,δB^q):

(16)

from which we build the matrix

(17)

using the ensemble average . We have found it important that the realizations of δB that enter equation 16 are independent of the P realizations used in 13. This yields an unbiased estimate of δY (i.e., in practice ) and thus of . Otherwise, SV model errors tend to contribute constructively to the observed SV (they are positively correlated with the SV data), and variances in 17 tend to be dramatically underestimated.

The matrix is positive semidefinite but unfortunately most of its eigenvectors correspond to null eigenvalues because the size PQ of the ensemble is smaller than the number of rows R_y=NL_y(L_y+2). We have

(18)

where V is any nonzero column vector. As a result, is not a valid covariance matrix and is not invertible. It would require PQ ≫ R_y in order to have a converged estimate of all cross covariances in ; their accuracy is limited by the size of the ensemble. To avoid considering spurious cross covariances that would bias the measure of the SV model errors, we have to modify . First, spatial correlations, in the spherical harmonic domain, of the flow model uncertainties are neglected in order to avoid overestimates of these quantities, and we consider only temporal correlations. Tests on a smaller problem show that ignoring these spatial correlations does not significantly affect the time changes of the output flow models. Second, we guard against artificial correlations between SV model errors at distant epochs by using a covariance localization approach [Gaspari and Cohn, 1999]. This consists in taking the Hadamard product (element-by-element multiplication) of and L,

(19)

where L is a covariance matrix constructed from a correlation function ρ_loc. The matrix is a valid covariance matrix because the Hadamard product of a positive definite matrix (such as L) and of a positive semidefinite matrix with all its diagonal elements strictly positive (such as ) is positive definite (Schur's theorem; see Horn [1990]). This heuristic approach has been extensively used in data assimilation studies of the dynamics in the atmosphere and the ocean, to filter the background covariance matrix as a function of distance [Hamill et al., 2001]. Instead, we use it here to filter the SV model error covariance matrix as a function of time separation.

As in oceanic applications, L is constructed from a correlation function ρ_gc defined by equation (4.10) in Gaspari and Cohn [1999], which involves a cutoff period τ_loc. As illustrated in Figure 1, this allows time differentiability properties compatible with the SV spectrum (see the smooth behavior of ρ_gc at zero lag). Furthermore, being defined on a compact support prevents overestimation of cross covariances at large lag. The cutoff period should be slightly larger than empirical estimates of the correlation length [Oke et al., 2007]. We find through the analysis of that the SV model error typically have a decay time between 10 and 20 years (depending on the SV coefficients). We therefore adopt a cutoff period τ_loc=30 years. Figure 1 illustrates the impact of the localization process on the ensemble estimate of time correlations. The time correlation functions at zero lag is less sharp than the Laplacian function obtained for a stochastic process such as that defined by equation 4: this corresponds to a case intermediate between the two tutorial examples illustrated in Appendix A. Note that the most extreme localization would involve keeping only the diagonal elements of , i.e., ρ_loc(τ) = δ(τ), which we also test for comparison purposes.

Finally, the combined SV data error covariance matrix is then . Covariances for E^o are derived from the error covariance matrix C_spl for the COV-OBS spline model coefficients [Gillet et al., 2013]: if D(t) is the operator relating the SV coefficients at epoch t to the spline coefficients, then the matrix is built from N × N blocs that contain the covariances for the observation errors between all SV spherical harmonic coefficients at epochs t_n and .

Note that Gillet et al. [2009, 2011] assumed . Our procedure to obtain C_e also differs from the iterative method of Pais and Jault [2008] inasmuch as they needed only one realization (B¹,Y¹) to estimate how the diagonal elements of C_e vary with the degree. In practice less than 10 iterations suffice to obtain a converged ensemble of solutions and thus converged covariance matrices C_e and C_w. For all frequencies and length scales, the rate of change of the average solution, from one iteration to the next, is smaller than the dispersion within the ensemble of solutions (about 5% of the kinetic energy), and we observe no drift of the solution through the iterative process. The variances that enter C_w in equation 14, averaged per degree ℓ and over time, change by about 4% from one iteration to the next (similarly, with no drift).

3.1 Importance of Time Correlations of Flow Coefficients and SV Model Errors

Table 1. Statistics of the Derived Flow Models for Several Values of τ_u, and Several Localization Functions ρ_loca

Case	ρloc(τ)	τu (y)	χ2					rγ
A10	ρgc(τ)	10	0.213	0.86	1.02	1.72	1.56	0.87	0.77
A100	ρgc(τ)	100	0.406	0.81	0.86	1.34	1.03	0.86	0.80
A300	ρgc(τ)	300	0.664	0.75	0.74	1.19	0.86	0.85	0.69
D10	δ(τ)	10	0.045	0.93	0.75	1.91	0.88	0.92	0.79
D100	δ(τ)	100	0.110	0.88	0.44	1.49	0.37	0.89	0.75

• ^a ρ_gc stands for the localization function from Gaspari and Cohn [1999], used with a cutoff frequency τ_loc=30 years. r_γ is the correlation coefficient between the LOD data and the LOD prediction from the ensemble average solution (similarly, is defined for the series filtered at periods between 4 and 9.5 years). Other symbols are defined in the text (section 3.1).

We find that ignoring temporal covariances of model errors (cases “D” with ρ_loc(τ) = δ(τ)) leads to a loss of information. Table (1) shows that this omission leads us to either overpredict LOD changes on decadal time scales (see D₁₀) or underpredict LOD changes on inter annual time scales (see D₁₀₀).

We also find that the amplitude of LOD predictions is a very useful diagnostic in assessing the appropriate flow correlation time τ_u. For τ_u=10 years (or less), the predicted LOD fluctuations are significantly more intense than the actual LOD fluctuations (see A₁₀). On the other hand, it appears impossible to account for rapid SV fluctuations with τ_u=300 years, or more (see A₃₀₀). We thus consider the solution obtained for τ_u=100 years as providing an acceptable compromise, producing reasonable amplitudes of the SV and LOD changes at both decadal and interannual periods. Henceforth, we focus on our preferred flow ensemble A₁₀₀, and discuss it in detail.

3.2 Fit to SV Observations

In this section we analyze how predictions from our ensemble of flow solutions fit the SV data. We consider both the misfit to COV-OBS Gauss coefficients and comparisons to annual differences of observatory annual means. The former are analyzed in terms of time-averaged SV spatial power spectra, presented in Figure 2. Overall, the SV spectra of the a posteriori and a priori errors calculated through our ensemble scheme superimpose. The latter is similar to that found by Pais and Jault [2008] for a snapshot problem. We obtain a posteriori errors slightly larger (smaller) than the prior errors at long (short) length scales, which translates into a misfit less than unity (see Table 1).

We focus on degrees 1 and 2 of the SV that show large posterior errors compared to the prior errors. Analyzing SV coefficients individually, we find anomalously enhanced misfits for . SV predictions for these coefficients are biased toward zero as illustrated by Figure 3 for the axial dipole. Extrema of (e.g., from 1960 to 1980) are particularly poorly reproduced (see the red curves in the Figure 3). The difficulty to reproduce large values for some low-degree coefficients (for instance and in Figure 3) explains why SV predictions from our flows have slightly yet systematically lower amplitude than SV measured at observatories, especially when the magnitude of the SV is high (e.g., dY/dt at Niemegk or dZ/dt at Kakioka and M'Bour, in Figure 4). The biases are less pronounced for recent epochs (cf. Figures 3 and 4), due to the improved accuracy of SV data. In contrast with the slow changes, the high-frequency fluctuations are very well reproduced, once we account for time-correlated SV model errors.

Note that these results are independent of the imposition of the equatorial symmetry constraint. If we omit the time correlation of SV model errors (see the blue curve in Figure 3) or if we artificially decrease these errors for , we are able to better fit these coefficients, which indicates that the topological constraints are not the origin of the difficulty. Furthermore, looking at the spatial distribution of residuals at the core-mantle boundary (not shown), we find no evidence for any particular region displaying preferentially large SV misfit. Despite the improvements we have introduced for the calculation of the SV model errors, a bias in the distribution of SV residuals persists. To some extent, this is unavoidable if the prior distributions of flow coefficients are all centered around 0. The optimization scheme then preferentially selects models with SV predictions biased towards zero, within the specified error bars (see the tutorial example in Appendix A). Unfortunately, the relatively large amplitude of the SV model errors for the low-degree coefficients exacerbates this effect. In section 5 we discuss possible approaches for handling this issue.

In Figure 4 we compare the SV predicted by our ensemble of flow models with annual differences of observatory annual means in Kakioka, Niemegk, and M'Bour. We observe a larger dispersion of the ensemble at earlier epochs: as noted by Pais and Jault [2008], the more knowledge we have about the observed SV (at recent epochs), the smaller SV model errors become. Nevertheless, within this dispersion, the interannual changes are rather coherent between realizations. Interannual changes can be recovered, despite the relatively large SV model errors, due to the consideration of time correlations in C_e (see Table 1). Of course, the bias mentioned above in the spherical harmonic domain also maps into the predictions at observatory locations (see Figure 4) albeit to a smaller extent than in the study of Gillet et al. [2009].

4.1 Resolution of the Calculated Core Flows

The spatial structure in our ensemble of flow models is predominantly steady throughout the studied interval of 1940–2010. By way of illustration, the correlation coefficient of core surface velocity maps calculated for pairs of epochs 10 years apart remains very close to 0.93 during [1945,2005]. Assuming that the correlation function has the form 5, this corresponds to a correlation time of at least 140 years. Remarkably, the time-averaged flow also includes significant small-scale constituents, as can be seen from the time-averaged spatial power spectrum of the flow:

(23)

Similarly, we define to be the spectrum for the time-dependent component of the flow, after the time-averaged part has been removed. In Figure 5 we display the ensemble average of the two power spectra, and , together with the power spectra for the ensemble averages, and . We find the flow coefficients are resolved within the ensemble up to degree about 10, above which departs gradually from . There is likely some useful information on the time average flow up to degree 13. A similar scale of spatial resolution is obtained for the time-dependent component of the flow when comparing with , even though large length-scales seem slightly less resolved for the flow fluctuations than for its time average. The energy of the time-variable part of the flow peaks at degree 10, above which there are larger uncertainties on the flow coefficients. Thus, the time-variable flow is predominantly small scale.

In order to obtain a more detailed measure of the flow resolution as a function of harmonic degree and frequency f, we introduce the quantity

(24)

where the complex coefficients (τ_ℓm,σ_ℓm) are the Fourier transforms of the flow coefficients time series t_lm(t),s_lm(t). This quantity is 1 (respectively 0) when the spread in the flow harmonic coefficients reaches 100% (respectively 0%) of the model variances. A low value of is a necessary but not sufficient condition for the true core flow to be captured.

Distributions of are shown in Figure 6 for the two intervals 1940–1975 and 1975–2010. We observe a clear improvement at recent epochs, with well-resolved flows at both shorter periods and smaller length scales. Indeed, for the most recent time interval, flow variations with a period of about 6 years are well resolved up to degree 6, with some information provided up to degree 12, whereas slower decadal variations are well constrained until degree 12, with some information provided up to degree 16. Though resolution is poorer for older epochs, some interannual changes are constrained even during the first 35 years. Again, we wish to emphasize that our ability to reliably retrieve interannual flows is a direct consequence of accounting for time-correlated SV model errors (see Appendix A).

4.2 Planetary Gyre and and Midlatitude Eddies

Under the incompressible QG hypothesis (see section 2.2), the flow in the whole volume can be represented through a stream function ψ(s,φ) [Jault and Finlay, 2015; Canet et al., 2014, equations (14) and (15)]:

(25)

(1_s,1_φ,1_z) are the unit vectors in cylindrical polar coordinates. is the half height of a geostrophic cylinder with c the outer core radius. We give in Appendix B the relation between coefficients describing the stream function and the toroidal and poloidal coefficients. The fluid flow in the equatorial plane, the first term on the right-hand side of 25, is parallel to the isolines of ψ, but its intensity is not proportional to the density of curves.

In Figure 7 (top), we present maps at the CMB and in the equatorial plane of the time average flow showing a planetary-scale gyre similar to that described in earlier studies [Pais and Jault, 2008; Gillet et al., 2009; Aubert, 2013]. We find that the gyre possesses a detailed structure, in agreement with the findings of Amit and Pais [2013] for flows based on the incompressible QG hypothesis. The most conspicuous features within the gyre are two anticyclones centered at (45^∘E, 60^∘N, S) and (60^∘W, 45^∘N, S) at the core surface. Large velocities are observed where the flow is in the azimuthal direction. In particular, the westward flow in the Atlantic hemisphere is split into two branches, one within ±10^∘ latitude at the equator and the other at latitudes from 30 to 45^∘. These are separated by a region of lower azimuthal velocity.

In addition to the time-averaged gyre structure, we also observe time-dependent features at decadal time scales. There is a general increase in the westward solid-body rotation from 1940 onward (see Figure 8, top). Assuming that the system of solid Earth (crust, mantle, and core) is closed, we find that our ensemble of geostrophic flows, of the form u_G(θ,t) at the core surface, account rather well for the observed decadal changes in the length of day (Figure 8, middle). A correlation coefficient r_γ close to 0.9 is found between the data and the predictions at decadal periods (see Table 1).

We also observe transient nonzonal circulations. The flow perturbation in 2005 (see Figure 7, bottom right) enhances anticyclonic eddies centered near (50^∘E, 40^∘N, S), (40^∘W, 65^∘N, S) and (90^∘W, 45^∘N, S), and cyclonic eddies around (60^∘E, 60^∘N, S), (60^∘W, 45^∘N, S), and (150^∘W, 45^∘N, S). These are reminiscent of the two main time-dependent structures isolated by Pais et al. [2014] on time scales about 70 years and longer. However, the most energetic time variable flows are nonzonal azimuthal jets located around 30^∘ latitude and in the equatorial belt (see Figure 7, bottom left). If the former appears related to the gyre, the latter are difficult to describe using equatorial projections of the stream function. The dynamics in the equatorial region is further discussed in section 4.5.

4.3 Taylor's State in the Earth's Core and Excitation of Torsional Waves

Nonzonal flows account for the majority of the energy of time variable flows, as can be inferred from the power spectral densities:

(26)

of respectively zonal and nonzonal motions ( stands for the power spectral density of the series t_ℓm(t), with similar definition for ). Nevertheless, the ratio of the zonal to nonzonal kinetic energies shows distinctive spectral bands centered on 6–8 years and around 3 years where zonal motions are relatively more intense (see Figure 9). We shall not attach much importance to the 3 year spectral band which is not adequately resolved (see Figure 6), but the well-resolved peak at 6–8 years provides an indication that there may be torsional waves present in the derived core flows.

In Figure 10 we show a time cylindrical radius map of the zonal velocity from 1940 to 2010, for the period band between 4 and 9.5 years where we find enhanced zonal energy in Figure 9. The contribution of these flows to LOD changes is displayed in Figure 8 (bottom), which shows the predictions from individual members of the flow ensemble and from the ensemble mean. The contribution of external fluid envelopes (mainly the atmosphere and to a lesser extent the ocean) to angular momentum changes is known to dominate for periods up to about 3 years [Gross et al., 2004]. However, the interannual variability of atmospheric angular momentum [Paek and Huang, 2012] appears to be too small to account for a prominent signal observed with a period of around 6 years [Abarca del Rio et al., 2000; Chao et al., 2014]. The zonal core motions that we have isolated readily provide a explanation for this signal. We corroborate the results of Gillet et al. [2010], which were limited to the time interval 1955–1975 (note that there is a delay compared to their Figure 2b, due to the fact that they used a noncausal filter and omitted to shift the time axis). The time variability of the observed and predicted LOD in the period range [4–9.5] years [see also Chao et al., 2014] is in apparent conflict with the finding of Holme and De Viron [2013] who decomposed iteratively the LOD data (from 1962 to 2012) into a decadally varying signal and a 5.9 year oscillation of almost constant amplitude (compare their Figure 2 with our Figure 8, bottom). In our opinion, the relatively small amplitude of the oscillation (in comparison with that of decadal changes) makes it difficult to decide whether it is long standing or heavily damped. Figure 8 (bottom) displays all the LOD changes produced by geostrophic flows in the frequency range [4–9.5] years. They need not all be attributed to the propagation of torsional waves. In any case, the remarkable agreement between our predictions and the geodetic data encourages us in the interpretation of the flow model down to periods about 4 years.

Geostrophic motions appear very clear over 1995–2010, particularly as the torsional wave approaches the equator, at latitudes below 40^∘, with a node of the waveform at about 10^∘ latitude (see Figure 10, bottom). The amplitude of the motions in this region is significantly larger than the spread in the flow ensemble (even at earlier epochs), yet the better resolution of the field model at recent epochs may have increased the sensitivity in the relatively small (in latitudinal extent) equatorial area. We confirm the slower propagation inferred by Gillet et al. [2010] as the wave gets closer to the equator and find no evidence for reflection at the equator.

Now the theory of “magnetostrophic dynamos” [see, e.g., Roberts and Wu, 2015], which has been developed to account for the Earth's magnetic field, gives us a tool to interpret the ratio as a function of frequency. We note above that this ratio remains small and does not vary much for periods larger than 8 years (see Figure 9). Taylor [1963] demonstrated that in the absence of inertia and viscosity we have

(27)

with Σ(s) the geostrophic cylinders (see Roberts and Aurnou [2012] or Jault and Finlay [2015] for modern accounts of Taylor's theory). Differentiating in time 27 and substituting ∂B/∂t with its expression from the induction equation

(28)

Taylor [1963] obtained a linear relationship between the geostrophic zonal flow u_G and nongeostrophic motions u_NG, which depends on the magnetic field inside the core (see his equation (4.5); η is the magnetic diffusivity),

(29)

with

(30)

The quantity C(s) has the dimension of a velocity and is proportional to the r.m.s. value of the cylindrical radial field B_s averaged over geostrophic cylinder (ρ is the density of the outer core, and μ is the free space magnetic permeability). Both and are quadratic functions of the magnetic field, depends linearly on the nonzonal velocity, and stems from the diffusion term η∇²B. In a quasi-geostrophy framework, we can identify u_NG with the nonzonal motions u_NZ. Equation 29 then determines a linear relationship between u_G and u_NZ that may explain the uniform ratio between their energies at long periods. There is perhaps no need to seek another mechanism for decadal variations of the geostrophic velocity and for decadal signals in the length of a day.

Restoring the inertial acceleration ∂u_G/∂t [Roberts and Aurnou, 2012, equation (24b)], Taylor's relationship is transformed into

(31)

The homogeneous part (on the left-hand side) corresponds to the torsional waves equation of Braginsky [1970], where C can now be interpreted as the torsional wave velocity. Comparing 29 and 31 makes clear that Taylor's differential equation 29 is valid on time scales that are long compared to the period of the torsional waves. The nonzonal velocities on the right-hand side of 31 appear as a possible source term for the torsional waves provided they have the appropriate time scale, as mentioned also by Teed et al. [2014] who searched for torsional waves in numerical simulations of the geodynamo.

Over the best resolved era (the last 15 years), the outward propagation of geostrophic motions, which can be interpreted as a torsional wave in most of the outer core, present a node 10^∘ away from the equator where u_G remains of small amplitude (Figure 10, bottom). The mechanism responsible for the special behavior close to the equator remains unclear. Latitudes below 10^∘ involve geostrophic motions only up to 50 km cylindrical depth from the equator. These geostrophic motions, which become important after 1995, carry a tiny fraction (of the order of 1%) of the outer core angular momentum. We may lack resolution to detect them further back in time, and their existence at older epochs cannot be ruled out solely on the basis of the good fit to LOD variations from 1950 onward.

4.4 On Electromagnetic Core-Mantle Coupling

Putting aside the fast torsional waves governed by 31, we read equation 29 as an indication that u_G and u_NZ have similar time scales. This equation, however, does not directly constrain the solid body rotation part of the core flow, for which ∂(u_G/s)/∂s = 0. The time evolution of this rotation is instead governed by the coupling mechanisms acting at the outer core boundaries.

The torque Γ^o(t) acting between the core and the mantle is linearly related to the time derivative of the observed LOD. In a series of papers, Holme [1998a, 1998b] found that core surface flow models exist that explain geomagnetic data between 1900 and 1980 and can also account for decadal changes in Γ^o(t) through electromagnetic coupling if the conductance of the mantle

(32)

is 10⁸ S or greater (δ is the thickness of the conducting layer; for the figures given below, the conductivity of the mantle, as a function of radius, is chosen as κ_m(r) = 3000(c/r)³⁰ S/m). We test here whether this result holds for our ensemble of flow models, while accounting for their uncertainties. Assuming that these are correct, we can follow a direct approach rather than the inverse approach of Holme [1998a]. We associate each member u^p of the ensemble of core flows with a time series Γ^p(t) of the electromagnetic torque acting at the CMB (see Appendix C). We find that the mantle conductance G has to be about 7 × 10⁸ S to make the typical amplitude of Γ^p(t) match the torque values inferred from LOD changes at decadal periods (about 10¹⁸ N m peak to peak), with a correlation coefficient of 0.36 between observations and the ensemble average predictions. It seems difficult to reconcile the time variation of our flow model with an electromagnetic explanation for the core-mantle coupling if the mantle conductance is less than 3 × 10⁸ S. The difference with the results of Holme [1998b] comes from the longer characteristic times of our large-scale flows, which are responsible for most of the electromagnetic torque (when ignoring time correlations of the SV model errors, in case D₁₀, we find G ∼ 4×10⁸ S for the magnitude of the predicted torque to match that of Γ^o). We also calculate the torque at interannual periods. The correlation coefficient between the observed torque and the ensemble average predictions is 0.56. A large conductance G ∼ 3×10⁹ S is required to account for the amplitude of the observed torque, at these periods.

The phase of the predicted series fits relatively well that of the geodetic series over both period ranges. However, the conductances required to make electromagnetic coupling a viable mechanism at, respectively, 6 and 25 year periods are strongly at variance. This does not come as a surprise since the solid body rotation accounts for a significant portion of both the electromagnetic torque and the LOD changes. Assuming that fluctuations of were linearly related to the variation of angular momentum at periods from 6 to 60 years, they could not also be linearly related also to the evolution of its time derivative (and hence of Γ^o) over the same period range. The high value of G that is required for electromagnetic core-mantle coupling at either period would certainly hinder the propagation of torsional waves across the core [Dumberry and Mound, 2008; Gillet et al., 2010].

Finally, we have also calculated the cross correlation between observed and predicted LOD changes from our ensemble of core flows (Figure 11). We find that the correlation is maximum for a delay between observed and predicted values of τ = 0.26 ± 0.29 year. The obtained value is the ensemble average of time lags, and the error is estimated as the standard deviation within the ensemble of time lags. The positive value of the delay means observations are, in average, ahead of predictions. Since diffusion in the mantle may cause a negative delay, we estimate that the lag most probably lays between −0.3 and 0 year (within 2 standard deviations) ; its estimation may help constrain the electrical conductivity of the lowermost mantle.

4.5 Dynamics of the Equatorial Region

We have observed in sections 4.2 and 4.3 a remarkable morphology of the azimuthal velocity in the equatorial belt at various time scales: the time-averaged velocity, the decadal flow changes, and the interannual geostrophic motions all show a minimum in amplitude at about 10^∘ latitude (see Figures 7 and 10, bottom). Decadal, nonzonal (i.e., nonaxisymmetric), azimuthal jets in the equatorial region do not seem directly related to u_φ at higher latitudes, even though they are linked to u_θ at latitudes around 10^∘ through mass conservation (see for instance around longitudes 60^∘W, 70^∘E, or 130^∘E in Figure 7, bottom left). These low-latitude features and their time evolution are consistently replicated within the flow ensemble (see Figure 12, top).

At interannual periods, nonzonal azimuthal motions have a minimum in amplitude at 10^∘ latitude, as illustrated with the latitudinal profiles of the r.m.s. of u_φ in Figure 13. Longitudinal profiles of u_φ at the equator present localized structures, with particularly intense flows around longitude 85^∘W at recent epochs (see Figure 12, bottom). Their peak-to-peak amplitude is much larger than that found for the geostrophic waves discussed in section 4.3. Even though the largest of these jets below 10^∘ can be traced back only to the mid-1990s, be they either axisymmetric (zonal) or else nonzonal (see Figures 10 and 12, bottom), we cannot exclude that such active low-latitude dynamics was present at earlier times, on account of the improved resolution of magnetic field models constructed from satellite observations. To our knowledge, the intriguing velocity minimum at 10^∘ latitude has not been reported from geodynamo simulations. In our opinion, further theoretical and numerical studies of this particular region are called for.

5.1 Core Flow Time Changes

We have presented an attempt to consistently reconstruct time changes of core motions. These are primarily associated with disturbances of the westward eccentric gyre identified by Pais and Jault [2008]. We obtain a weaker temporal variability compared to previous studies [e.g., Amit and Olson, 2006]. Indeed, we find the gyre to be largely steady over 1940–2010.

Our model does not perfectly account for the low-frequency changes of the dipole term and of a few other low-degree coefficients. This seems to be a result of the relatively large SV model errors, associated with our flow, affecting the slow variations of these coefficients. This indicates that it may be necessary to model core motions as a perturbation centered on a nonzero background flow (the equivalent of the climatic mean for the oceans dynamics). We must acknowledge, however, that there may be contributions to these slow, large-scale, field changes that are neglected in our quasi-geostrophic frozen-flux model. For example, one consequence of quasi-geostrophy is that the longitudinal average of the meridional flow is zero. A contrasting view has been offered by Buffett [2014] who calculated waves at the top of the core assuming that it is stably stratified. He found that the zonal flow is coupled to an axisymmetric meridional flow that causes the dipole magnetic field to fluctuate with a period of about 60 years. On the other hand, our quasi-geostrophic flows are able to reproduce much of the dipole decay observed in recent decades, when the impact of SV modeling errors is less pronounced, through nonzonal meridional flows acting on the longitudinally asymmetric core field.

5.2 Limits on the Time Resolution of Core Flow Models

Information on the time variability of SV model errors amounts to a time-dependent constraint on core flow calculations. Our investigations show that consistently accounting for time-correlated SV model errors has very encouraging consequences regarding the rapid flow changes that may be inferred from satellite magnetic data. We found that the dispersion within our ensemble of flow solutions at subdecadal periods was significantly reduced at recent epochs, indicating that a more detailed picture of the rapid core flow changes is emerging from the availability of continuous satellite data.

Our confidence in the estimated flow changes is supported by their ability to reproduce LOD variations at both decadal and interannual periods. Consistent modeling of interannual and decadal LOD changes was possible here for the first time thanks to our modeling of time variable SV model errors. In previous core flow reconstructions investigating this issue [Wardinski, 2004; Gillet et al., 2010], there was no constraint on the flow changes between epochs (i.e., ) and no time-correlated SV model error was considered (ρ_loc=δ). As a result, decadal LOD changes were found to be significantly overestimated.

We identify the approximately 6 years periodic LOD signal with torsional waves in the outer core. Our analysis is compatible with a modulation of their amplitude inside the time interval covered by our study, with particularly intense velocities at latitudes below 40^∘ after 1995. These torsional waves may be excited by the interaction of nonzonal motions with the magnetic field inside the core. Both the geostrophic flow and the nonzonal azimuthal velocity present a minimum in amplitude at 10^∘ latitude. The understanding of the interaction between zonal and nonzonal motions in the equatorial region requires further theoretical and numerical investigations.

Torsional waves are not the main source of the intense and rapidly changing SV patterns recently identified near the equator from the analysis of satellite measurements during the past 15 years [Olsen and Mandea, 2007; Chulliat and Maus, 2014]. Rather, enhanced nonzonal flow at low latitudes, particularly around 90^∘W as well as 60^∘ and 120^∘E, appear to be responsible.

5.3 Possible Additions to the Procedure

Determining the posterior core flow probability density given the information contained in the SV data, as attempted here, not only provides an estimate of the uncertainty on the core motions but also provides an opportunity to construct flow models subject to extra constraints. Indeed, this can be formalized as an optimization problem of the form

(33)

with given by the posterior covariance matrix for the flow model errors.

A first instance of an extra constraint is at s = c, which is the boundary condition for torsional waves under the insulating mantle hypothesis [see Schaeffer et al., 2012; Jault and Finlay, 2015]. Similarly, one could produce flow solutions satisfying the nonpenetration condition u_s=0 at the cylindrical surface tangent to the inner core.

We have estimated the modeling error E^m from the unresolved magnetic field δB. Alternatively, one could also consider estimating directly E^m, using “augmented state” data assimilation schemes [see, for instance, Evensen, 2003]. In this framework, time covariances of SV model errors can be accounted for by advecting E^m with a stochastic equation that is consistent with their correlation functions, as presented in Figure 1.

The stochastic approach employed throughout this study can be modified with the incorporation of deterministic elements involving either 3-D geodynamo simulations [Aubert, 2014a] or QG models (such as that of Canet et al. [2009]). For instance, the information about time covariances governed by a stochastic process could be combined with spatial covariances inferred from geodynamo simulations. Another improvement may be to estimate simultaneously the velocity and magnetic fields (as in Lesur et al. [2015]) using sequential data assimilation.

The flow model presented here, together with its associated covariance matrix, is readily available from the address http://isterre.fr/recherche/equipes/geodynamo/themes-de-recherche/article/analyse-de-donnees-geomagnetiques?lang=fr.