# The flow at the Earth's core-mantle boundary under weak prior constraints

Julien Baerenzung, Matthias Holschneider, Vincent Lesur

► To cite this version:

Julien Baerenzung, Matthias Holschneider, Vincent Lesur. The flow at the Earth's core-mantle boundary under weak prior constraints. Journal of Geophysical Research: Solid Earth, 2016, 121, pp.1343-1364. 10.1002/2015JB012464. insu-03581670

HAL Id: insu-03581670

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

Submitted on 20 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1002/2015JB012464

**Key Points:**

* Bayesian inversion for the flow at the Earth's core-mantle boundary

* General prior constraint

**Correspondence to:**

J. Baerenzung,
baerenzung@gmx.de

**Citation:**

Baerenzung, J., M. Holschneider, and V. Lesur (2016), The flow at the Earth's core-mantle boundary under weak prior constraints, *J. Geophys. Res. Solid Earth*, **121**, 1343–1364, doi:10.1002/2015JB012464.

Received 21 AUG 2015

Accepted 12 FEB 2016

Accepted article online 15 FEB 2016

Published online 11 MAR 2016

# The flow at the Earth's core-mantle boundary under weak prior constraints

Julien Baerenzung¹,², Matthias Holschneider¹, and Vincent Lesur³

¹Institute for Mathematics, University of Potsdam, Potsdam, Germany, ²Earth's Magnetic Field, GeoForschungsZentrum Potsdam, Potsdam, Germany, ³Institut de Physique du Globe de Paris, Paris, France

**Abstract** Prior information in ill-posed inverse problem is of critical importance because it is conditioning the posterior solution and its associated variability. The problem of determining the flow evolving at the Earth's core-mantle boundary through magnetic field models derived from satellite or observatory data is no exception to the rule. This study aims to estimate what information can be extracted on the velocity field at the core-mantle boundary, when the frozen flux equation is inverted under very weakly informative, but realistic, prior constraints. Instead of imposing a converging spectrum to the flow, we simply assume that its poloidal and toroidal energy spectra are characterized by power laws. The parameters of the spectra, namely, their magnitudes, and slopes are unknown. The connection between the velocity field, its spectra parameters, and the magnetic field model is established through the Bayesian formulation of the problem. Working in two steps, we determined the time-averaged spectra of the flow within the 2001–2009.5 period, as well as the flow itself and its associated uncertainties in 2005.0. According to the spectra we obtained, we can conclude that the large-scale approximation of the velocity field is not an appropriate assumption within the time window we considered. For the flow itself, we show that although it is dominated by its equatorial symmetric component, it is very unlikely to be perfectly symmetric. We also demonstrate that its geostrophic state is questioned in different locations of the outer core.

## 1. Introduction

The geomagnetic field is sustained by the dynamo process acting in the Earth's outer core. The strong degree of turbulence the magnetohydrodynamic flow is exhibiting, the fast rotation of the Earth, and the intense convection [see Kono and Roberts, 2002; Christensen et al., 2010] make it a dynamical system extremely difficult to model and study. In addition, and in contrast with the Earth's atmosphere, directly measuring the velocity field associated with the fluid motions in the outer core is technically excluded. Nevertheless, thanks to satellite and observatory measurements of the magnetic field, probing the outer core becomes feasible. Indeed, assuming the mantle to be insulating, core field models deriving from satellite or observatory data can be estimated everywhere outside the core. At the core-mantle boundary, the radial component of the magnetic field, which is the only component that can continuously cross the jump in electrical conductivity between the mantle and the outer core [see Roberts and Scott, 1965; Jault and Le Mouël, 1991], is advected by the two dimensional fluid motions (the flow is only two dimensional because it cannot penetrate the solid boundary formed by the mantle). Yet this process can be described by the frozen flux (FF) approximation, an equation where diffusion effects are neglected [see Backus et al., 1996]. Therefore, by inverting this equation, the flow at the core-mantle boundary (CMB) can be recovered. Unfortunately, the operation is not straightforward, the problem being ill posed for two reasons. First, only one equation is connecting the two components of the velocity field to the radial component of the magnetic field, leading to an underdetermined problem. The second reason arises from the nature of the available core field. Indeed, when measured at the Earth's surface, the magnetic field is composed of different contributions, in particular, the external field, the lithospheric or crustal field, and the core field. Whereas recent models derived from satellite or observatory data are able to extract the external contribution from the total field, the core and lithospheric field cannot be separated. The lithospheric field being dominant above degree $l=13$ when expanded in spherical harmonics (SH), the core magnetic field and secular variation are undetermined at scales smaller than $l=13$. Yet any velocity field scales are likely to produce large-scale secular variation when interacting nonlinearly with the magnetic field. So once again the number of independent equations (one per secular variation spectral component) is lower than the velocity scales populating the flow in the outer core, leading to an underdetermined problem.
Because the inverse problem is ill posed, it is necessary to impose prior velocity constraints to solve it. To build such constraints, different options can be considered. One possibility is to rely on physical assumptions on the outer core dynamics to analytically construct some constraints. Such an approach led to various formulations, including tangential-geostrophy, tangential-magnetostrophy, general, steady, columnar, helical, or purely toroidal flow (see Holme [2007] and Finlay et al. [2010] for reviews of the underlying assumptions leading to such constraints). However, these constraints being only sufficient to reduce the degree of nonuniqueness of the velocity field, the solution of the inverse problem still needs to be regularized [see Chulliat and Hulot, 2000], an operation which is usually achieved by a priori imposing a certain power law behavior for the total kinetic energy spectrum of the flow [see Holme, 2007]. Another possibility to built prior velocity constraints is to rely on the statistical properties of three-dimensional simulations of the geodynamo, as proposed by Aubert and Fournier [2011], Fournier et al. [2011], and Aubert [2013, 2014]. One of the main advantages of this method is that the full statistical picture of the outer core velocity field and magnetic field is available. Therefore, through their cross correlations, a knowledge of the velocity and the magnetic field at the core-mantle boundary is sufficient to reconstruct them in the entire outer core. Nevertheless, the method has some drawback, the principle one being certainly that up to now, the degree of turbulence reachable in such simulations remains very low. As a consequence, the regimes which are obtained may not allow for a proper reproduction of the geodynamo's dynamical behavior, and the statistical properties of the simulated flow may not be consistent with the real ones.

Yet to obtain realistic posterior statistics, prior assumptions have to reflect the effective knowledge available for the velocity field at the CMB. Since it can also be argued that any of the physical assumptions previously mentioned could describe the dynamics of the fluid motions at the CMB (or the entire outer core), but since these constraints are not all consistent with one another, imposing them in a strong manner is certainly not the best option if one wants to evaluate the reliability of certain features the flow is exhibiting. As an example, if one a priori forces the velocity field to be equatorial symmetric, it will be impossible to a posteriori assess how likely deviations from this symmetric state are. In a similar manner, using unrealistic spectra parameters for the regularization process will undoubtedly lead to erroneous posterior statistics.

This is why we propose to consider that the flow evolving at the CMB is completely general and only characterized by independent poloidal and toroidal power law energy spectra. Although it is believed that the velocity field has dominant large scales [see Holme, 2007] so that it presents a decaying energy spectrum, this hypothesis has only been reinforced, as in Celaya and Wahr [1996], but it has never been proven. Therefore, we assume that the poloidal and toroidal energy spectra are characterized by power laws, with slopes and magnitudes a priori unknown. Whereas incorporating such information in a traditional least squares approach would be impossible, the flexibility of Bayesian methods allows for it. Because of the complexity of the system to be solved, we will proceed in two steps. The energy spectra of the velocity field at the CMB will be first determined, to be then used as a prior information to construct the posterior distribution of the flow.

The article is organized as follows. In section 2 the statistical parametrization of every quantity entering the inverse problem is detailed, and the connection between these quantities is established in a Bayesian framework. In section 3, we first evaluate, with synthetic data, the formalism we developed as well as the numerical approach we proposed to solve the inverse problem. We then apply our method to determine the time-averaged poloidal and toroidal energy spectra of the flow at the CMB for the period 2001–2009.5. These spectra are then used as a prior information to study the flow itself and its possible variations for the epoch 2005. Finally, we give some concluding remarks in section 4.

## 2. Mathematical Formulation of the Inverse Problem

### 2.1. The Magnetic Field and the Flow at the Core-Mantle Boundary

In the Earth's outer core, the dynamical evolution of the magnetic field **B** is prescribed by the induction equation:

$$ \partial_t \mathbf{B} = \nabla \times (\mathbf{u} \times \mathbf{B}) - \nabla \times (\eta \nabla \times \mathbf{B}), \quad (1) $$

where $\eta$ is the magnetic diffusivity and **u** is the velocity of the electrically conducting fluid. As mentioned previously, because of the low electrical conductivity of the mantle [see Velímský, 2010; Jault, 2015], the connection between magnetic field models derived from satellite or observatory data and the magnetic field sustained by dynamo action in the outer core can be formalized at the level of the core-mantle boundary.
At this location, since magnetic field models such as GRIMM-3 of Lesur and Wardinski [2012] only provide the spatiotemporal evolution of the large-scale part of the field (the small-scale field being obscured by litho-spheric field at the Earth's surface), and since the diffusion effects on these large scales are negligible in comparison to the advection effects [see Holme, 2007], the radial component of the induction equation can be simplified into the frozen flux (FF) equation:

$$ \partial_t B_r = - \nabla_H (\mathbf{u} B_r), \qquad (2) $$

where $\nabla_H$ is the horizontal divergence operator.

In spectral space, the radial component of the magnetic field as well as the secular variation $\partial_t B_r$ can be expanded in spherical harmonics such as

$$ B_r = - \sum_{l=1}^{l=+\infty} (l+1) \sum_{m=-l}^{m=l} g_{l,m}(t) Y_{l,m}, \qquad (3) $$

$$ \partial_t B_r = - \sum_{l=1}^{l=+\infty} (l+1) \sum_{m=-l}^{m=l} \partial_t (g_{l,m}(t)) Y_{l,m}, \qquad (4) $$

with $Y_{l,m}$ the Schmidt seminormalized spherical harmonics (SH) of degree $l$ and order $m$, and $g_{l,m}$ the spherical harmonics coefficients of the magnetic field outside the core. According to this decomposition, two vectors $b$ and $\gamma$, respectively, containing the spherical harmonics coefficients associated with $B_r$ and $\partial_t B_r$ are built. In a similar manner, vectors $b_0$ and $\gamma_0$ are defined and filled with the spherical harmonics coefficients associated with the GRIMM-3 magnetic field and secular variation.

The velocity field $\mathbf{u}$ at the CMB is decomposed into a poloidal scalar field $\phi$ and a toroidal scalar field $\psi$ such as

$$ \mathbf{u} = \mathbf{e}_r \times \nabla_H \psi + \nabla_H \phi. \qquad (5) $$

In spectral space, the poloidal and toroidal fields, respectively, derive from the coefficients $\phi_{l,m}$ and $\psi_{l,m}$ through the formulation:

$$ \phi = \sum_{l=1}^{l=+\infty} \sum_{m=-l}^{m=l} \phi_{l,m} \hat{Y}_{l,m}, \qquad (6) $$

$$ \psi = \sum_{l=1}^{l=+\infty} \sum_{m=-l}^{m=l} \psi_{l,m} \hat{Y}_{l,m}, \qquad (7) $$

where $\hat{Y}_{l,m}$ are orthonormalized spherical harmonics. The coefficients $\phi_{l,m}$ and $\psi_{l,m}$ are then inserted in a vector referred to $\mathbf{u}$.

The frozen flux approximation can now be expressed for the vectors $\gamma$, $b$, and $\mathbf{u}$ as follows:

$$ \gamma = -A_b u = -A_u b, \qquad (8) $$

where the linear operators $A_b$ and $A_u$ allow us to calculate the SH coefficients associated with the advection term $\nabla_H(\mathbf{u}B_r)$ when they are, respectively, applied to $\mathbf{u}$ and $b$.

## 2.2. Usual Formulation of the Inverse Problem

Models of the Earth's secular variation and magnetic field can be used to estimate the flow at the CMB through the inversion of the FF equation. In a Bayesian formalism [see Jackson, 1995; Baerenzung et al., 2014], this operation consists of determining the maximum of the posterior distribution $p(\mathbf{u}|\gamma_0)$ defined as follows:

$$ p(\mathbf{u}|\gamma_0) = \frac{1}{p(\gamma_0)} \int p(\gamma_0|\mathbf{u}, b)p(\mathbf{u})p(b)db, \qquad (9) $$

where $p(\mathbf{u})$ and $p(b)$ are, respectively, the prior distributions of the velocity field and the magnetic field, $\gamma_0$ is the secular variation model with its associated distribution $p(\gamma_0)$, and $p(\gamma_0|\mathbf{u}, b)$ is the likelihood distribution.
In this method $p(u)$ is prescribed and usually chosen as a Gaussian distribution with a mean $E[u]=0$ and a covariance $E[uu^T]=\Sigma_u$. The derivation of $\Sigma_u$ is commonly based on two assumptions: the flow at the CMB is isotropically distributed, and its energy spectrum $E_u(l)$ is characterized by a certain power law, which, for orthonormalized spherical harmonics gives the following:

$$E_u(l) = l(l+1) \sum_{m=-l}^{m=l} (\phi_{l,m}^2 + \psi_{l,m}^2), \qquad (10) \\ \sim A^2 l^{-p},$$

where $A$ is the magnitude of the spectrum and the exponent $P$ determines the strength of the penalty on the high-degree flow components [Holme, 2007]. Whereas the parameter $P$ is generally prescribed, *Pais and Hulot* [2000] proposed to estimate $A$ a posteriori. This latter operation is performed by evaluating the flow maximizing $p(u)|_{\gamma_0}$ with different values for $A$ and keeping the one which is consistent with both the secular variation uncertainties and the prior belief on the flow energy.

### 2.3. General Prior Information on the Flow Behavior

In reality, if the spectrum of the flow at the CMB is characterized by a certain power law, neither its magnitude nor its slope are a priori known. Furthermore, the poloidal and toroidal parts of the flow certainly present distinct spectral behavior. To account for these particular aspects of the velocity field in its prior distribution, we assume the following:

1. The poloidal and toroidal part of the total velocity field are uncorrelated and isotropically distributed over the CMB.

2. The poloidal and toroidal spectra are individually characterized by power laws with possibly different spherical harmonics ranges, where the slopes of the spectra may vary from one range to the other.

3. The velocity field is normally distributed over the CMB with a 0 mean and a covariance depending on the above statements.

Under these latter hypothesis, the poloidal and toroidal spectra of the flow can be written as follows:

$$E_{\phi}(l) = C_{\phi_i} A_{\phi}^2 l^{-P_{\phi_i}} \quad \text{for } l \in \Delta_{\phi_i}, \qquad (11)$$

$$E_{\psi}(l) = C_{\psi_i} A_{\psi}^2 l^{-P_{\psi_i}} \quad \text{for } l \in \Delta_{\psi_i}, \qquad (12)$$

where $A_{\phi}$ and $A_{\psi}$ are, respectively, the magnitude of the poloidal and toroidal energy spectra and $P_{\phi_i}$ and $P_{\psi_i}$ are the slopes of the spectra within the spherical harmonics ranges $\Delta_{\phi_i}$ and $\Delta_{\psi_i}$. The number of ranges characterizing the poloidal and toroidal spectra, respectively, $N_{\phi}$ and $N_{\psi}$, may differ from one case to the other. Finally, the constant $C_{\phi_i}$ and $C_{\psi_i}$ which allow the spectra to be continuous between two different ranges are given by

$$C_{\phi_i} = \prod_{a=2}^{a=i} \exp \left( \log \left( I_{\phi_{a-1}} \right) \left( P_{\phi_a} - P_{\phi_{a-1}} \right) \right), \qquad (13)$$

$$C_{\psi_i} = \prod_{a=2}^{a=j} \exp \left( \log \left( I_{\psi_{a-1}} \right) \left( P_{\psi_a} - P_{\psi_{a-1}} \right) \right), \qquad (14)$$

$I_{\phi_a}$ and $I_{\psi_a}$ being the spherical harmonic degrees where transitions in slopes occur.

With the assumptions on the flow we detailed at the beginning of this section and with the energy spectra given in equations (11) and (12), the covariances of the poloidal and toroidal fields can be expressed as follows:

$$E[\psi_{l,m} \psi_{l',m'}] = \frac{E_{\psi}(l)}{l(l+1)(2l+1)} \delta_{ll'} \delta_{mm'} \qquad (15)$$

$$E[\phi_{l,m} \phi_{l',m'}] = \frac{E_{\phi}(l)}{l(l+1)(2l+1)} \delta_{ll'} \delta_{mm'} \qquad (16)$$

$$E[\phi_{l,m} \psi_{l',m'}] = 0 \quad \forall l, l', m, m'. \qquad (17)$$
These correlation functions depend directly on the poloidal and toroidal spectra of the flow; they are therefore connected to the spectra parameters referred to the following:

$$ \mathcal{M} = \{ A_{\phi}, P_{\phi_i}, \Delta_{\phi_i}, A_{\psi}, P_{\psi_j}, \Delta_{\psi_j} \}. \qquad (18) $$

Since the covariance matrix of the flow can only be constructed once the parameters $\mathcal{M}$ are known, we express it as $\Sigma_{u|M}$. As a consequence, the prior distribution of the velocity field is also conditioned by $\mathcal{M}$ and reads as follows:

$$ p(u|\mathcal{M}) = \mathcal{N}(0, \Sigma_{u|\mathcal{M}}). \qquad (19) $$

## 2.4. Joint Posterior Distribution of the Flow, the Magnetic Field, and the Spectra Parameters

The different quantities entering the inverse problem are connected together through the joint distribution $p(u, b, \mathcal{M}, \gamma_0)$. According to Bayes [1763], and neglecting a priori the correlations between the magnetic field and the velocity field, $p(u, b, \mathcal{M}, \gamma_0)$ can be split as follows:

$$ p(u, b, \mathcal{M}, \gamma_0) = p(u, b, \mathcal{M}|\gamma_0)p(\gamma_0) \qquad (20) $$

$$ = p(\gamma_0|u, b, \mathcal{M})p(b)p(u|\mathcal{M})p(\mathcal{M}), \qquad (21) $$

where $p(u, b, \mathcal{M}|\gamma_0)$ is the joint posterior distribution of the flow, the magnetic field, and the spectra parameters. The remaining distributions entering equations (20) and (21) are detailed in the following.

### 2.4.1. Likelihood Distribution

The uncertainties on the secular variation model being assumed to be normally distributed with a 0 mean and a covariance $\Sigma_{\gamma}$, the likelihood distribution $p(\gamma_0|u, b, \mathcal{M})$ reads as follows:

$$ p(\gamma_0|u, b, \mathcal{M}) = \frac{\exp\left[-\frac{1}{2}(\gamma_0 + A_b u)^T \Sigma_{\gamma}^{-1} (\gamma_0 + A_b u)\right]}{(2\pi)^{\frac{d}{2}}|\Sigma_{\gamma}|^{\frac{1}{2}}}, \qquad (22) $$

where $d$ is the dimension of the secular variation vector. The construction of the secular variation covariance matrix is based on an estimated spectrum of the secular variation uncertainties at the level of the Earth's surface. This spectrum is chosen to be flat (as would be the spectrum of spatially uncorrelated measurement errors), with a magnitude of $A_{\gamma}^2 = 0.01$ (nT/yr)$^2$, a value approximately corresponding to the energy of the secular variation at SH degree $l = 14$, where the noise level of high-latitude data starts to generate significant artifacts. At the CMB the covariance matrix is therefore diagonal with the following structure:

$$ E \left[ (\gamma_0(l, m) - \gamma(l, m))(\gamma_0(l', m') - \gamma(l', m')) \right] = \frac{A_{\gamma}^2 (R_s/R)^{(2l+4)}}{(l+1)(2l+1)} \delta_{ll'} \delta_{mm'}, \qquad (23) $$

with $R_s = 6371.2$ km the radius at the Earth's surface, $R = 3485$ km the core radius, and $\gamma$ the real secular variation.

### 2.4.2. Conditional Prior Distribution of the Velocity Field Given the Spectra Parameters

As already expressed in section 2.3, we recall that the conditional prior distribution of the velocity field given the spectra parameters $\mathcal{M}$, $p(u|\mathcal{M})$, is assumed to be normally distributed with a 0 mean and a covariance $\Sigma_{u|\mathcal{M}}$ such as

$$ p(u|\mathcal{M}) = \frac{\exp\left[-\frac{1}{2}u^T \Sigma_{u|\mathcal{M}}^{-1} u\right]}{(2\pi)^{\frac{d}{2}}|\Sigma_{u|\mathcal{M}}|^{\frac{1}{2}}}, \qquad (24) $$

with $d_u$ the dimension of the velocity field vector.

### 2.4.3. Prior Distribution of the Spectra Parameters

To reflect the complete lack of information we possess on $\mathcal{M} = \{A_{\phi}, P_{\phi_i}, \Delta_{\phi_i}, A_{\psi}, P_{\psi_j}, \Delta_{\psi_j}\}$, we suppose that all these parameters are independent from one another and uniformly distributed over infinite or semi-infinite ranges such as the prior distribution $p(\mathcal{M})$ reads as follows:

$$ p(\mathcal{M}) = p(A_{\psi})p(A_{\phi}) \prod_{i=1}^{i=N_{\phi}} (p(P_{\phi_i})p(\Delta_{\phi_i})) \prod_{j=N_{\psi}-1}^{j=N_{\psi}} (p(P_{\psi_j})p(\Delta_{\psi_j})), \qquad (25) $$
with

$$p(P_{\phi_i}) = p(P_{\psi_j}) = \mathcal{U}(-\infty, \infty), \qquad (26)$$

$$p(A_{\psi_j}) = p(A_{\phi_i}) = \mathcal{U}(0, \infty), \qquad (27)$$

$$p(\Delta_{\phi_i}) = p(\Delta_{\psi_j}) = \mathcal{U}(1, \infty), \qquad (28)$$

and where the notation $\mathcal{U}(a, b)$ stands for uniformly distributed between a and b.

#### 2.4.4. Prior Distribution of the Magnetic Field

The magnetic field is assumed to be normally distributed with a mean $b_0$ corresponding to the field provided by models fitting satellite or observatory data (in our case the GRIMM-3 model of Lesur and Wardinski [2012]), and a covariance $\Sigma_b$ such as

$$p(b) = \frac{\exp\left[-\frac{1}{2}(b-b_0)^T \Sigma_b^{-1} (b-b_0)\right]}{(2\pi)^{\frac{d_b}{2}} |\Sigma_b|^{\frac{1}{2}}}, \qquad (29)$$

where $d_b$ is the dimension of the vector b.

To build the covariance matrix $\Sigma_b$, one has to analyze the possible variations of the magnetic field around its mean $b_0$. These variations mainly arise from two phenomena. At large scales ($1 \le l \le 13$), magnetic field models such as GRIMM-3 predict a biased core field since this latter is contaminated by the lithospheric field. The model also contains uncertainties due to measurement errors. However, the fact that the spectrum associated with the difference between the GRIMM3 model and the CHAOS-3 model of Olsen et al. [2010] is generally below the lithospheric spectrum (except for the SH degree $l = 1$, an issue which will not be considered here) indicates that uncertainties due to measurement errors can be neglected. At small scales ($l > 13$), because at the Earth's surface the intensity of the lithospheric field is larger than the core one, this latter field is completely obscured.

To account for this uncertain nature of the available magnetic field in the inverse problem, we split it into two parts as $b = b^< + b^>$, where $b^<$ and $b^>$, respectively, correspond to the large- and small-scale parts of the total field. Each part of the field is then treated statistically in a distinct manner. Whereas the average value of $b^<$ is given by $E[b^<] = b_0$, the expectation of $b^>$ is assumed to be $E[b^>] = 0$. The covariance matrices associated with each field are built accordingly to two different spectra.

For the large-scale field, the theoretical spectrum of the lithospheric field proposed by Thébault and Vervelidou [2015] is considered. The latter expressed in Tesla² reads as follows:

$$E_L(l) = \frac{1}{2} (l+1) (\mu_0 |M| F_l(\epsilon))^2 l^{\gamma} C_l, \qquad (30)$$

with $|M| = 0.45$ A m⁻¹ the mean apparent induced magnetization, $\epsilon = 27$ km the equivalent magnetized layer thickness, and the constants $\mu_0 = 4\pi 10^{-7}$ m kg s⁻² A⁻² and $\gamma = 1.4$. The two functions $F_l(\epsilon)$ and $C_l$ are given by

$$C_l = \frac{l(20l^3 + 8l^2 - 13l + 3)}{3(2l + 3)(2l + 1)(2l - 1)} \\ F_l(\epsilon) = \frac{1 - (1 - \epsilon/R_s)^{(l-1)}}{l-1}.$$

From this spectrum the covariance matrix $\Sigma_{b^<}$ of the large-scale field $b^<$ is defined by

$$\Sigma_{b^<} = E \left[ (b^< - b_0) (b^< - b_0)^T \right] \\ = \frac{E_L(l)(R_s/R)^{(2l+4)}}{(l+1)(2l+1)}, \qquad (31)$$
For the small-scale field, an extrapolation of the resolved scale spectrum $E_{b_0}$ is used to build $\Sigma_{b^>}$, the covariance matrix of $b^>$. To extrapolate $E_{b_0}$, the analytical formulation of the magnetic field spectrum at the CMB from *Buffett and Christensen* [2007] is employed. It reads the following:

$$E_{b^>} (l) = C_1 \chi^l, \qquad (32)$$

where $\chi = 0.99$ and where the constant $C_1$ is determined by fitting $E_{b_0}$ between the SH degrees $2 < l \le 13$. The covariance matrix of $b^>$ is therefore given by

$$\Sigma_{b^>} = E[b^>b^{>T}] = \frac{E_{b^>} (l)}{(l+1)(2l+1)}. \qquad (33)$$

The total covariance of the magnetic field $\Sigma_b$ is finally obtained by combining $\Sigma_{b<}$ and $\Sigma_{b^>}$.

## 2.5. Posterior Distribution of the Flow Spectra Parameters

Before evaluating the flow at the CMB, we wish to determine the optimal parametrization its poloidal and toroidal energy spectra are exhibiting. Therefore, we have to evaluate the maximum of the posterior distribution $p(M|\gamma_0)$, defined as follows:

$$\begin{aligned} p(M|\gamma_0) &= \frac{1}{p(\gamma_0)} \iint p(\gamma_0|u, b, M)p(b)p(u|M)p(M)dbdu \\ &= p(\gamma_0|M) \frac{p(M)}{p(\gamma_0)}. \end{aligned} \qquad (34)$$

The analytical calculation of this integral is extremely complicated, and the resulting distribution would be very hard to handle numerically. Therefore, we propose to approximate the marginalized likelihood distribution $p(\gamma_0|M)$ by a normal distribution, where the mean $\tilde{\gamma}_0$ and the covariance $\Sigma_{\gamma_0}$ are given by

$$\tilde{\gamma}_0 = \iiint \gamma_0 p(\gamma_0|u, b, M)p(b)p(u|M)dbdudy_0 \qquad (35)$$

$$\Sigma_{\gamma_0} = \iint \int (\gamma_0 - \tilde{\gamma}_0) (\gamma_0 - \tilde{\gamma}_0)^T \times p(\gamma_0|u, b, M)p(b)p(u|M)dbdudy_0. \qquad (36)$$

Introducing the third-order tensor $A_{ijk}$ such as

$$A_{ijk} u_j b_k = (A_b u)_i = (A_u b)_i \qquad (37)$$

the solutions of integrals (35) and (36) are the following, respectively,

$$\tilde{\gamma}_{0i} = 0, \qquad (38)$$

$$\begin{aligned} (\Sigma_{\gamma_0})_{ij} = (\Sigma_\gamma)_{ij} &+ (A_{b_0} \Sigma_{u_{1,M}} A_{b_0}^T)_{ij} \\ &+ A_{imn} (\Sigma_{u_{1,M}})_{mr} (\Sigma_b)_{ns} A_{jrs}. \end{aligned} \qquad (39)$$

The posterior distribution $p(M|\gamma_0)$ is then approximated by

$$p(M|\gamma_0) \sim \frac{\exp\left[-\frac{1}{2}\gamma_0^T \Sigma_{\gamma_0}^{-1} \gamma_0\right]}{(2\pi)^{\frac{d}{2}} |\Sigma_{\gamma_0}|^{\frac{1}{2}}} \times \frac{p(M)}{p(\gamma_0)}. \qquad (40)$$

Note that if the magnetic field was exactly known, and given by $b_0$, the last term on the right-hand side of equation (39) would vanish and relation (40) would be exact.
### 3. Evaluation and Results

#### 3.1. Numerical Method

To evaluate the flow at the CMB and its posterior variability, the first step consists in determining the optimal parameters (magnitudes and slopes), its poloidal, and toroidal energy spectra are exhibiting. These parameters are referred to $\hat{\mathcal{M}}$ and obtained with a gradient method approximating the maximum of distribution $p(\mathcal{M}|\gamma_0)$ given in equation (40).

Once the optimal spectra parameters are determined, they are used as a prior information for the estimation of the flow at the CMB. The distribution to be considered is therefore the posterior distribution of the velocity field $u$, conditioned by the secular variation model $\gamma_0$, and the spectra parameters $\hat{\mathcal{M}}$, namely, $p(u|\gamma_0, \hat{\mathcal{M}})$. This probability density function can be obtained by marginalizing with respect to the magnetic field $b$, the following distribution:

$$p(u, b|\gamma_0, \hat{\mathcal{M}}) = p(\gamma_0|u, b, \hat{\mathcal{M}})p(b)p(u|\hat{\mathcal{M}})/p(\gamma_0). \quad (41)$$

The result, presented by Baerenzung et al. [2014] without the $\hat{\mathcal{M}}$ dependency, reads as follows:

$$p(u|\gamma_0, \hat{\mathcal{M}}) = \frac{(2\pi)^{-\frac{3d}{4}}}{|\Sigma_{\hat{\gamma}}^{\frac{d}{4}}|} \exp \left[ -\frac{1}{2} (\gamma_0 + A_u b_0)^T \Sigma_{\hat{\gamma}}^{-1} (\gamma_0 + A_u b_0) \right] \times \frac{1}{|\Sigma_u|^{\frac{d_u}{2}}} \exp \left[ -\frac{1}{2} u^T \Sigma_{u|\hat{\gamma}}^{-1} u \right] \times \frac{1}{p(\gamma_0)} \quad (42)$$

$$\Sigma_{\hat{\gamma}} = \Sigma_{\gamma} + A_u \Sigma_b A_u^T. \quad (43)$$

To sample this distribution, Baerenzung et al. [2014] proposed to use a Markov chain Monte Carlo algorithm. Unfortunately, this method is numerically costly and necessitate a fine tuning of the velocity field random walk. Here we propose another approach, much more efficient numerically, and requiring no tuning; it is referred to the Gibbs sampling algorithm [see Gamerman and Lopes, 2006]. Instead of sampling $p(u|\gamma_0, \hat{\mathcal{M}})$, this latter method is used to draw pairs of magnetic field and velocity field from the joint posterior distribution $p(u, b|\gamma_0, \hat{\mathcal{M}}) = p(\gamma_0|u, b, \hat{\mathcal{M}})p(b)p(u|\hat{\mathcal{M}})/p(\gamma_0)$. To do so, the algorithm proceeds iteratively such that at each step $n$, a velocity field $u^n$ is randomly drawn from the conditional probability distribution:

$$p(u^n|b^{n-1}, \gamma_0, \hat{\mathcal{M}}) = \mathcal{N}(\hat{u}^n, \Sigma_{u^n}) \quad (44)$$

with

$$\hat{u}^n = - (A_{b^{n-1}}^T \Sigma_{\gamma}^{-1} A_{b^{n-1}} + \Sigma_{u|\hat{\gamma}}^{-1})^{-1} A_{b^{n-1}}^T \Sigma_{\gamma}^{-1} \gamma_0 \quad (45)$$

$$\Sigma_{u^n} = (A_{b^{n-1}}^T \Sigma_{\gamma}^{-1} A_{b^{n-1}} + \Sigma_{u|\hat{\gamma}}^{-1})^{-1} \quad (46)$$

and a magnetic field $b^n$ is randomly drawn from the following:

$$p(b^n|u^n, \gamma_0, \hat{\mathcal{M}}) = \mathcal{N}(\hat{b}^n, \Sigma_{b^n}) \quad (47)$$

with

$$\hat{b}^n = (A_{u^n}^T \Sigma_{\gamma}^{-1} A_{u^n} + \Sigma_b^{-1})^{-1} (\Sigma_b^{-1} b_0 - A_{u^n}^T \Sigma_{\gamma}^{-1} \gamma_0) \quad (48)$$

$$\Sigma_{b^n} = (A_{u^n}^T \Sigma_{\gamma}^{-1} A_{u^n} + \Sigma_b^{-1})^{-1}. \quad (49)$$

Note that the Gaussian nature of distributions (44) and (47) allows efficient sampling them numerically.
<table><caption>**Table 1.** Flow Spectra Parameters *M* Imposed and Estimated in the Evaluation Phaseª</caption><thead><tr><th>Case</th><th>Δ<sub>φ</sub></th><th>A<sub>φ</sub></th><th>P<sub>φ</sub></th><th>Δ<sub>ψ</sub></th><th>A<sub>ψ</sub></th><th>P<sub>ψ</sub></th></tr></thead><tbody><tr><td>a</td><td>[1 26]</td><td>5</td><td>1</td><td>[1 26]</td><td>20</td><td>2</td></tr><tr><td>b</td><td>[1 26]</td><td>5.02</td><td>0.99</td><td>[1 26]</td><td>19.91</td><td>1.95</td></tr><tr><td>c</td><td>[1 26]</td><td>2.98</td><td>0.62</td><td>[1 26]</td><td>16.8</td><td>1.78</td></tr><tr><td>d</td><td>[1 26]</td><td>34.57</td><td>3</td><td>[1 26]</td><td>34.57</td><td>3</td></tr></tbody></table>

Calculating the expectation of some function $f(u)$ using the posterior ensemble is similar to evaluating $E[f(u)]$ according to the distribution $p(u|\gamma_0, \hat{M})$ since

$$ 
\begin{aligned}
E[f(u)] &= \iint f(u)p(u, b|\gamma_0, \hat{M})dbdu \\
&= \int f(u) \left( \int p(u, b|\gamma_0, \hat{M})db \right) du \\
&= \int f(u)p(u|\gamma_0, \hat{M})du.
\end{aligned} 
\quad (50) $$

Numerically, all computations are performed on a discrete CMB as in Baerenzung et al. [2014], and the evaluation of the different mathematical operators is performed through a pseudospectral code.

For every simulation, the prescribed magnetic field $b_0$ and the secular variation $\gamma_0$ are defined up to SH degree $l = 13$, the magnetic field is extrapolated up to SH degree $l=26$ accordingly to the formulation of Buffett and Christensen [2007], and the velocity field is expanded up to SH degree $l=26$.

## 3.2. Reliability of the Spectra Parameters Evaluation

In this section, we test, with artificial data, the ability of the approximated distribution $p(M|\gamma_0)$ given in equation (40), and the algorithm which determines its maximum, to provide the spectra parameters of the flow at the CMB. To do so, the maximum of the posterior distribution $p(M|\gamma_0)$ is evaluated numerically for an ensemble of $N_e = 5000$ artificially generated velocity fields, magnetic fields and secular variations.

For each realization, a velocity field is randomly drawn from $p(u|M) = N'(0, \Sigma_{u|M})$ where the spectra parameters *M* are prescribed and given in Table 1. Note that in this evaluation phase both poloidal and toroidal energy spectra are characterized by unique slopes, therefore the indices of the spherical harmonics ranges $\Delta_{\phi}$ and $\Delta_{\psi}$, and slopes $P_{\phi}$ and $P_{\psi}$ can be dropped, such that these spectra parameters become, respectively, $\Delta_{\phi}, \Delta_{\psi}, P_{\phi},$ and $P_{\psi}$.

The large-scale magnetic field $b_0$ is taken from the GRIMM-3 model at the epoch 2002.0. It is contaminated by a lithospheric field $b_L$ randomly drawn from the distribution $p(b_L) = N'(0, \Sigma_{b_L})$ discussed in section 2.4.4. To generate small-scale magnetic fields, the energy spectrum associated with $b_0$ is extrapolated up to SH degree $l = 26$. From this spectrum, the covariance matrix $\Sigma_{b^>}$ is built as shown by equation (33). The ensemble of small-scale magnetic fields are then randomly drawn from the distribution $p(b^>) = N'(0, \Sigma_{b^>})$.

For each pair of artificially generated magnetic field and velocity field, the associated secular variation is calculated and contaminated by a Gaussian noise with a 0 mean and a covariance $\Sigma_{\gamma}$. The maximum of $p(M|\gamma_0)$ is then numerically approximated for every artificial secular variation and large-scale magnetic field and referred to $\hat{M} = \{\hat{A}_{\phi}, \hat{P}_{\phi}, \hat{A}_{\psi}, \hat{P}_{\psi}\}$.

The results of the numerical simulations are shown in Figure 1 (top row) through the two-dimensional histograms $p(\hat{A}_{\phi}, \hat{P}_{\phi})$ (top left) and $p(\hat{A}_{\psi}, \hat{P}_{\psi})$ (top right). As it can be observed, the maximum of $p(M|\gamma_0)$ does not always point toward the exact values of the spectra parameters chosen for the experiment and given in Table 1. Instead, they spread around these values. Nevertheless, when averaged over the entire ensemble

<footnote>ªCase a: parameters prescribed to randomly generate velocity fields in both tests of section 3.2; case b: average values of the recovered parameters; case c: estimated parameters *M* associated with the pair of magnetic field and secular variation selected in section 3.3; case d: same as case c, but with imposed spectral slopes and estimated total magnitude.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Histograms (color maps) and marginalized distributions (side plot) of the (left column) poloidal and (right column) toroidal spectra magnitudes and slopes estimated with (top row) synthetic secular variations and (middle and bottom rows) truncated velocity fields. The truncation scales are set to spherical harmonic degree $l=10$ (Figure 1, middle) and $l=5$ (Figure 1, bottom) for the poloidal field and to $l=9$ (Figure 1, middle) and $l=5$ (Figure 1, bottom) for the toroidal field. We recall that the spectra parameters imposed in the experiment are $\{A_{\phi}, P_{\phi}, A_{\psi}, P_{\psi}\} = \{5, 1, 20, 2\}$.</figcaption></figure>

of solutions we found, the mean values of the estimated parameters $\{E[\hat{A}_{\phi}], E[\hat{P}_{\phi}], E[\hat{A}_{\psi}], E[\hat{P}_{\psi}]\} = \{5.02, 0.99, 19.91, 1.95\}$ become very close to the prescribed ones $\{A_{\phi}, P_{\phi}, A_{\psi}, P_{\psi}\} = \{5, 1, 20, 2\}$, indicating that the parameter estimator is not biased.

Two mechanisms may explain why the estimated parameters can strongly deviate from the imposed ones. First, the posterior variance of the parameters is large, and the deviations are coming from the combination between data uncertainties and parameters prior variability. Another possibility is directly linked to the design of the experiment. Indeed, when a velocity field is randomly drawn from the distribution $p(u|M)$, the energy associated with any of its component will usually not be equal to the prior variance of the component. As a consequence, the slopes and magnitudes the poloidal and toroidal energy spectra exhibit may differ from the ones chosen a priori. If, in addition, only a portion of the velocity field scales contributes enough to the large-scale secular variation to be detected, the variations in the observed flow spectra will even be larger.
In other words, the slope and the magnitude a poloidal or a toroidal energy spectrum exhibit may also depend on the size of the spectral range which is considered. To measure the importance of these effects, we performed a second test. There, $N_e=20,000$ velocity fields are randomly drawn from $p(u|M)$ according to the parameters $M$ given in Table 1. Instead of conditioning the distribution of $M$ to the secular variation, we directly connect it to the randomly generated velocity fields. Furthermore, we take into account the dependency between the spectral parameters and the amount of velocity field scales considered. Therefore, for every realization we evaluate numerically the maximum $\hat{M}_{l_{\phi}}^{l_{\psi}} = \{\hat{A}_{\phi}^{l_{\phi}}, \hat{P}_{\phi}^{l_{\phi}}, \hat{A}_{\psi}^{l_{\psi}}, \hat{P}_{\psi}^{l_{\psi}}\}$ of distribution:

$$p(M_{l_{\phi}}^{l_{\psi}}|u) = p(u|M_{l_{\phi}}^{l_{\psi}})p(M_{l_{\psi}}^{l_{\phi}})/p(u), \qquad (51)$$

where $M_{l_{\phi}}^{l_{\psi}}$ corresponds to the spectra parameters of the poloidal and toroidal fields, respectively, truncated at spherical harmonic degree $l_{\phi}$ and $l_{\psi}$. Estimating $\hat{M}_{l_{\phi}}^{l_{\psi}}$ is similar to evaluating the parameters of power laws (slopes and magnitudes) which fit best the poloidal and the toroidal energy spectra between, respectively, $l=[1, l_{\phi}]$ and $l=[1, l_{\psi}]$.

The ensembles of $N_e=20,000$ estimated parameters indicate that whereas the averaged values of the parameters do not depend on $l_{\phi}$ and $l_{\psi}$ and are in any case very close to the imposed ones; the higher the truncation degree, the lower the variance of the parameters, and conversely. By comparing the statistical properties of the recovered parameters for every value of $l_{\phi}$ and $l_{\psi}$ with the statistical properties of the parameters estimated through the artificial secular variation, we could observe a close match when the truncation degrees were set to $l_{\phi}=10$ and $l_{\psi}=9$. Figure 1 (middle row), which displays the two-dimensional histograms $p(\hat{A}_{\phi}^{10}, \hat{P}_{\phi}^{10})$ (middle left) and $p(\hat{A}_{\psi}^{9}, \hat{P}_{\psi}^{9})$ (middle right), illustrates this result. To highlight the difference between optimal and suboptimal solutions, the two-dimensional histograms $p(\hat{A}_{\phi}^{5}, \hat{P}_{\phi}^{5})$ and $p(\hat{A}_{\psi}^{5}, \hat{P}_{\psi}^{5})$, where $l_{\phi}=l_{\psi}=5$, are also shown in Figures 1 (bottom left) and 1 (bottom right), respectively.

The combination of the two numerical tests we presented in this section shows that both the approximated posterior distribution $p(M|\gamma_0)$ and the numerical approach we employed to determine its maximum are well suited to evaluate the spectral characteristics of the flow at the Earth's CMB. Nevertheless, restrictions on the method abilities have to be considered. In particular, the range of validity of the recovered parameters is not infinite and is certainly conditioned by the real nature of the poloidal and toroidal fields themselves as well as the level of uncertainties on the secular variation model.

### 3.3. Sampling of the Posterior Distribution of the Velocity Field According to Optimal and Imposed Flow Spectra Parameters

In usual inversions of the frozen flux equation, prior constraints are imposed. Here contrary to this approach, we propose a formalism to estimate parameters of the flow energy spectra in order to construct a realistic prior distribution of the velocity field. To highlight the impact of imposed or derived prior constraints on the posterior statistics of the flow, the two approaches are now compared on synthetic data.

To do so, a pair of magnetic field and secular variation is taken from the previous test. The optimal spectra parameters $\hat{M}$ associated with these fields, the ones we numerically approximated and which are given in Table 1, case c, are considered to construct the prior covariance of the velocity field following the development given in equations (11)–(17). According to the associated prior distribution, a first ensemble of 50,000 pairs of velocity and magnetic fields $\{u_{ep}, b_{ep}\}$ is randomly drawn from the posterior distribution $p(u, b|\gamma_0, \hat{M})$ using the Gibbs sampling algorithm presented in section 3.1, equations (44)–(49). This ensemble is referred to the ep ensemble. A second ensemble $\{u_{ip}, b_{ip}\}$ of the same size is drawn from the posterior distribution $p(u, b|\gamma_0, M_{ip})$. In this case, the slopes belonging to the spectra parameters $M_{ip}$ are imposed and set to $P_{\phi} = P_{\psi} = 3$, and the magnitudes $A_{\phi}$ and $A_{\psi}$ are forced to be equal. Only the total amplitude $A = \sqrt{A_{\phi}^2 + A_{\psi}^2}$ is evaluated such that $p(M_{ip}|\gamma_0)$ is maximized. The resulting spectra parameters are given in Table 1, case d. Note that except for the estimation of the magnitude $A$, such a prior parametrization of the velocity field has already been employed in the inversion of the frozen flux equation, as in Eymin and Hulot [2005]. This ensemble is referred to the ip ensemble.

Before comparing the two subensembles $\{u_{ep}\}$ and $\{u_{ip}\}$ of, respectively, $\{u_{ep}, b_{ep}\}$ and $\{u_{ip}, b_{ip}\}$ with $u_0$, the velocity field used to construct the artificial secular variation, we check that both $\{u_{ep}\}$ and $\{u_{ip}\}$ are able to explain the secular variation at the right level of uncertainty. In usual least squares inversions, the balance
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spectrum of the artificial secular variation $\gamma_0$ at the Earth's surface (black lines), together with the imposed spectrum of secular variation errors (blue lines). The gray lines correspond to some spectra of fit residuals to the artificial secular variation $\epsilon_{\gamma}^i = \gamma_0 + \nabla_H(u^i b^i)$, where pairs of magnetic field and velocity field $(u^i, b^i)$ are taken from the ensembles (top) $\{u_{ep}, b_{ep}\}$ deriving from estimated prior constraints and (bottom) $\{u_{ip}, b_{ip}\}$ deriving from imposed prior constraints. The red lines correspond to the average spectra (over the entire ensembles) of the fit residuals to the artificial secular variation.</figcaption></figure>

between fit to the observed secular variation and prior belief on the flow energy, is controlled by a scalar factor $\lambda$ which is a priori unknown. In general, inversions are performed with different values of $\lambda$ until an optimal balance is reached [see *Pais and Hulot*, 2000]. In the Bayesian method we propose that the quantity which can be associated with the balance factor of the least squares approach is the magnitude of the total kinetic energy spectrum $A = \sqrt{A_{\phi}^2 + A_{\psi}^2}$ with $\lambda \sim 1/A^2$. Yet, in our case this quantity is determined before estimating any velocity field. For $A$ to be consistent with the synthetic secular variation $\gamma_0$, the secular variation $\gamma^i = -\nabla_H(u^i b^i)$ which is reconstructed accordingly to any given element $(u^i, b^i)$ of the ensembles $\{u_{ep}, b_{ep}\}$ or $\{u_{ip}, b_{ip}\}$ has to fit $\gamma_0$ at the right level of uncertainty. This means that the residual to the artificial secular variation, $\epsilon_{\gamma}^i = \gamma_0 - \gamma^i$, has to present similar statistical properties to the imposed measurement errors. We recall that the latter were assumed to present a flat spectrum at the Earth's surface with a squared magnitude $A_{\gamma}^2=0.01 \text{ nT}^2/\text{yr}^2$, as illustrated with a blue line in both Figures 2 (top) and 2 (bottom). Looking at the same figure, one can observe that the spectra displayed with gray lines are spreading around this imposed error spectrum. They correspond to the spectra of the residuals $\epsilon_{\gamma}^i$ associated with some pairs of velocity and magnetic field $(u^i, b^i)$ taken from the ensembles $\{u_{ep}, b_{ep}\}$ (Figure 2, top) and $\{u_{ip}, b_{ip}\}$ (Figure 2, bottom). So qualitatively, both ensembles seem to produce pairs of magnetic field and velocity field consistent with the synthetic data. When the spectra of the residuals are averaged over the entire ensembles, the result, shown with red lines in Figure 2 and referred to the global fit residual to the secular variation, become very close to the imposed error spectrum. To quantify the degree of consistency of our ensembles with the synthetic data, we define the average normalized misfit as

$$ \chi_{\gamma}^2 = \frac{1}{d} E \left[ (\gamma_0 + \nabla_H(ub))^T \Sigma_{\gamma}^{-1} (\gamma_0 + \nabla_H(ub)) \right]_{\text{post}} \quad (52) $$

with $d$, the dimension of the secular variation vector, $\Sigma_{\gamma}$, the covariance of the synthetic secular variation, and where the expectation $E[\dots]_{\text{post}}$ is performed over the joint posterior distributions of the velocity and the magnetic field. The values of this misfit are of $\chi_{\gamma}^2 = 0.9997$ and $\chi_{\gamma}^2 = 1.06$ for, respectively, the ensembles
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (left column) Toroidal and (right column) poloidal spectra, associated with the artificial velocity field $u_0$ and with the ensemble averaged velocity fields (top row) $\hat{u}_{ep}$ and (bottom row) $\hat{u}_{ip}$. Thick black lines: spectra of the artificial field $u_0$. Thin black lines: prior spectra. Thick grey lines: spectra of the ensemble averaged velocity fields $\hat{u}_{ep}$ (Figure 3, top) and $\hat{u}_{ip}$ (Figure 3, bottom). Grey dashed lines: rescaled scale-by-scale variances estimated with the velocity field posterior ensembles. Crosses: spectra of the error fields $\epsilon_{ep} = \hat{u}_{ep} - u_0$ (Figure 3, top) and $\epsilon_{ip} = \hat{u}_{ip} - u_0$ (Figure 3, bottom).</figcaption></figure>

$\{u_{ep}, b_{ep}\}$ and $\{u_{ip}, b_{ip}\}$. The fact that these quantities are very close to 1 indicates that both velocity field models are able to accurately explain the synthetic secular variation. However, as we will see in the following, the two flow solutions are not equivalent.

The velocity fields of both ensembles $\{u_{ep}, b_{ep}\}$ and $\{u_{ip}, b_{ip}\}$ are now examined in spectral space and compared to $u_0$, the artificial velocity field of this test. The average velocity field of $\{u_{ep}, b_{ep}\}$ is referred to $\hat{u}_{ep}$, whereas the average velocity field of $\{u_{ip}, b_{ip}\}$ is referred to $\hat{u}_{ip}$. Figure 3 displays the spectra of different quantities which derive from the toroidal (left) and poloidal (right) parts, of the velocity field $u_0$ and the subensemble $\{u_{ep}\}$ (top), and of the velocity field $u_0$ and the subensemble $\{u_{ip}\}$ (bottom). The thick black lines correspond to the spectra of the artificial toroidal and poloidal fields, whereas the spectra of the estimated fields $\hat{u}_{ep}$ and $\hat{u}_{ip}$ are shown with thick grey lines in Figures 3 (top row) and 3 (bottom row), respectively. The first observation one can make is that whereas at large scales, $\hat{u}_{ep}$ and $\hat{u}_{ip}$ present comparable levels of energies to the ones $u_0$ exhibit, at small scales, the energy spectra of $\hat{u}_{ep}$ and $\hat{u}_{ip}$ are decaying much faster than the ones associated with $u_0$. Comparing the two graphics on the right side of Figure 3 seems to indicate that the poloidal spectrum of $\hat{u}_{ip}$ fits much better the poloidal spectrum of $u_0$ than the poloidal spectrum of $\hat{u}_{ep}$. However, when computing the error fields $\epsilon_{ep} = \hat{u}_{ep} - u_0$ and $\epsilon_{ip} = \hat{u}_{ip} - u_0$, one can observe through their spectra (shown with crosses in Figure 3) that the energy they contain in any scale is of the same order for both average velocity fields. Even for the poloidal fields, where the energy associated with the degrees $l=2, 3, 4$ is lower for $\epsilon_{ip}$ than for $\epsilon_{ep}$ but larger for degrees $l=1, 5, 6, 7$.

From the subensembles $\{u_{ep}\}$ and $\{u_{ip}\}$ predictions on the uncertainties associated with $\hat{u}_{ep}$ and $\hat{u}_{ip}$ should be available. Here these predictions are computed in terms of scale-by-scale posterior variances. After being rescaled by the factor $l(l+1)$, they are plotted in Figure 3 with grey dashed lines. From the comparison between this quantity and the rescaled scale-by-scale prior variances (shown with thin black lines in Figure 3 and also referred to prior spectra) two spherical harmonics ranges can be identified. One where data have an influence on the posterior distribution of the flow, up to SH degree $l=8$ to $l=13$, and one where the prior distribution completely dominates posterior statistics, typically above SH degree $l=13$. For the part where data still have an influence on the estimated flow, one can see that the spectrum associated with the error field $\epsilon_{ep} = \hat{u}_{ep} - u_0$,
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparisons between the artificial velocity field $u_0$ and the ensemble solution $\hat{u}_{ep}$ deriving from estimated prior spectra parameters. Black arrows: (top row) artificial velocity field $u_0$ and (bottom row) ensemble averaged velocity field $\hat{u}_{ep}$, through their (left column) toroidal and (right column) poloidal components. The color maps in Figure 4 (top row) correspond to the absolute value of the difference in orientation between the artificial $u_0$ and estimated $\hat{u}_{ep}$ fields, whereas in Figure 4 (bottom row), they are associated with the 90% confidence interval on flow orientation according to the posterior velocity field sample {$u_{ep}$}.</figcaption></figure>

spreads around the scale-by-scale posterior variances, meaning that the prediction on flow uncertainties is consistent with the effective error. This is not the case anymore for the prediction associated with the {$u_{ip}$} subensemble, where the prior distribution is imposed. Here the predicted errors are clearly underestimating the effective errors with even a toroidal scale-by-scale posterior variance lower than the toroidal error spectrum at any SH degree. To quantify the departure or consistency of the flow predictions with the synthetic velocity field, one can evaluate the following normalized velocity field misfit:

$$ \chi_u^2 = \frac{1}{d_u} E \left[ (\hat{u} - u_0)^T \Sigma_{u|y_0:M}^{-1} (\hat{u} - u_0) \right]_{\text{post}} \quad (53) $$

where $d_u$ is the size of the velocity field vector and $\Sigma_{u|y_0:M} = E \left[ (u - \hat{u}) (u - \hat{u})^T \right]_{\text{post}}$ is the posterior covariance of the velocity field. Note that to evaluate $\chi_u^2$, only components of the poloidal and toroidal fields up to SH degree $l=13$ are taken into account since above this limit data have no more influence on the posterior velocity fields. The values $\chi_u^2$ is exhibiting confirm the visual diagnosis we established with the different spectra in Figure 3. Indeed, for the velocity field ensemble deriving from the estimated prior distribution, the value of $\chi_u^2 = 1.07$ is close to one, indicating that the prediction of the uncertainties associated with $\hat{u}_{ep}$ is consistent with the effective error $\epsilon_{ep} = \hat{u}_{ep} - u_0$. The value of $\chi_u^2 = 3.05$ associated with the velocity field generated accordingly to the imposed prior distribution indicates that the error $\epsilon_{ep} = \hat{u}_{ep} - u_0$ is not consistent with the posterior covariance of the velocity field. The error predictions are therefore strongly underestimated in this case.

Comparison between {$u_{ep}$}, {$u_{ip}$}, and $u_0$ is now performed in physical space. Since above SH degree $l=13$ data have almost no influence on the estimated flows at the CMB, the posterior distributions of the velocity field being completely dominated by the prior distributions, the different velocity fields ($\{u_{ep}\}$, $\{u_{ip}\}$, and $u_0$) are truncated at degree $l_c$. Here we present an evaluation of the ability of the method to correctly estimate the flow orientation, and if not to predict its possible deviations. Note that the results we will present are similar for the flow magnitude. Figure 4 is associated with the velocity field subensemble {$u_{ep}$} depending on the estimated spectra parameters. It displays the toroidal (left) and poloidal (right) parts of the artificial velocity field $u_0$ (top) and estimated velocity field $\hat{u}_{ep}$ (bottom) with black arrows. The color maps in Figure 4 (top row) correspond to the absolute value of the angle (in degree) between $u_0$ and $\hat{u}_{ep}$; it is therefore a measure of the error on the orientation of $\hat{u}_{ep}$. In Figure 4 (bottom row), the color maps are associated with the 90% confidence interval (according to the velocity ensemble) on the flow orientation. These latter maps must be read as follows. At a given location, the ensemble predicts that the orientation of the real flow is given by the orientation of the average velocity field vector $\hat{u}_{ep}$, plus or minus the value associated with the color at
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as Figure 4, but comparing the artificial velocity field $u_0$ and the ensemble solution {$u_{ip}$} deriving from imposed prior spectra parameters.</figcaption></figure>

this place, and this, with a 90% confidence. Globally, most of the locations where the ensemble is predicting little uncertainties in flow orientation (violet and blue colors) correspond to locations where the direction of $\hat{u}_{ep}$ remains close to the direction of $u_0$. Conversely, locations where the orientation of $\hat{u}_{ep}$ strongly deviates from the effective flow orientation (red patches) are usually belonging to areas of large orientation uncertainty.

The same analysis in physical space is now performed for the subensemble {$u_{ip}$} deriving from imposed prior spectra parameters. When looking at Figure 5 which displays the same quantities as Figure 4 but for {$u_{ip}$}, identical conclusions to the ones we just mentioned can be drawn. However, the comparison between Figures 5 and 4 provides much more interesting information.

Reflecting what we already observed in spectral space, the poloidal and toroidal parts of the average fields $\hat{u}_{ep}$ and $\hat{u}_{ip}$ are very close to one another. The orientation errors associated with these fields are also sharing similar locations and present comparable levels of intensity. Nevertheless, the orientation error is lower for the poloidal and toroidal parts of $\hat{u}_{ep}$ than for the poloidal and toroidal parts of $\hat{u}_{ip}$ on, respectively, 50.9% and 54.8% of the surface of the core-mantle boundary.

The most striking differences between Figures 4 and 5 appear in the color maps associated with the 90% confidence interval on the flow orientation. Although predictions of reliable or uncertain flow orientation are located in similar places, the intensity of these predictions strongly varies from one case to the other. According to Figures 4 and 5 the velocity field $\hat{u}_{ip}$ presents weaker posterior variability in orientation than $\hat{u}_{ep}$. However, when providing the 90% confidence interval on flow orientation, one has to expect that on approximately 90% of the surface of the CMB, the real velocity field orientation belong to the predicted interval. This is the case for predictions provided by the {$u_{ep}$} subensemble, where the real poloidal and toroidal field orientations (the one associated with $u_0$) are lying within the predicted interval on, respectively, 92.6% and 88.5% of the total core surface. For predictions arising from the {$u_{ip}$} subensemble, these proportions are strongly dropping to, respectively, 77.7% and 70.4% of the area of the CMB, indicating that this ensemble clearly underestimates the appropriate posterior variability of the velocity field.

From this synthetic test, two main conclusions can be drawn. First, the method we developed to estimate the spectra parameters in order to build prior constraints for the flow allows us to generate velocity fields which are consistent with the observed secular variation. Second, we demonstrated that although the spectra parameters chosen to construct the prior distribution of the flow do not strongly affect the average velocity field of the posterior ensemble (at least with the parameters we chose), taking appropriate values for these parameters is mandatory to obtain accurate estimations of flow uncertainties. Yet the parameters evaluation we propose in this study is very well suited to obtain realistic posterior statistics of the flow evolving at the Earth's core-mantle boundary.
### 3.4. The Poloidal and Toroidal Energy Spectra at the Core-Mantle Boundary

The method tested in the two previous sections is now applied to the GRIMM-3 secular variation and magnetic field model for the period 2001–2009.5. Instead of calculating the spectra parameters at a single epoch, we propose to evaluate their time-averaged values within the 2001–2009.5 period. To do so, the magnetic field and secular variation time series are sampled every 0.1 years. The likelihood distribution of each secular variation $\gamma_i$ taken from the sample is given by

$$p(\gamma_i|\mathcal{M}) = \frac{\exp\left[-\frac{1}{2}\gamma_i^T\Sigma_{\gamma_i|\mathcal{M}}^{-1}\gamma_i\right]}{(2\pi)^{\frac{d}{2}}|\Sigma_{\gamma_i|\mathcal{M}}|^{\frac{1}{2}}}. \qquad (54)$$

Neglecting the temporal dependency of both the secular variation and the magnetic field, the posterior distribution of the spectra parameters given the ensemble of $N_s = 86$ secular variations $\gamma_i$ and magnetic fields $b_i$ becomes the following:

$$p(\mathcal{M}|\Gamma = \gamma_1, \dots, \gamma_i, \dots, \gamma_{N_s}) = \left[ \prod_{n=1}^{n=N_s} \frac{p(\gamma_n|\mathcal{M})}{p(\gamma_n)} \right] \times p(\mathcal{M}). \qquad (55)$$

The maximum $\hat{\mathcal{M}}$ of this distribution is then numerically approximated for different cases. In the simulations we performed, up to two spherical harmonics ranges are used to characterize the toroidal and poloidal energy spectra. This choice was dictated by recent studies as in Aubert [2013], Baerenzung et al. [2014], and Gillet et al. [2015] showing that the flow at the core-mantle boundary was exhibiting a very energetic large-scale eccentric gyre. Yet the spectral signature of this structure may differ from the one associated with the remaining part of the flow. Therefore, by modeling the spectra with two ranges, if they exist, the two different spectral regimes will be distinguishable. Numerically, only the slopes and the magnitudes of the toroidal and poloidal energy spectra are evaluated, the separation of the spectra in different SH ranges being prescribed a priori. To determine at which scales the optimal changes in spectral slopes should occur, a large set of possible configurations are compared to one another through the maximum posterior probability density they exhibit. This operation is relevant since we a priori assess equal chances to any possible decomposition of the poloidal and toroidal energy spectra. Table 2 summarizes the different parameters we employed in the simulations as well as the results we obtained.

In this parametric study we proceeded gradually. In the first simulation (referred to I), the poloidal and toroidal spectra are both assumed to be characterized by unique slopes. The result we obtained shows that the toroidal field strongly dominates the poloidal one in the entire spectral range we modeled (i.e., between SH degree $l=1$ and $l=26$). In addition, whereas the poloidal spectrum is increasing, the toroidal one is decaying with a slope of the order of $P_\psi \sim 1$. This first result contradicts results of Celaya and Wahr [1996] who deduced that the total energy spectrum of the flow at the CMB should decay as $l^{-2}$ or faster. In the second set of simulations (simulations II to VIII), the poloidal spectrum is described by a single slope whereas the toroidal one is decomposed into two ranges. The degree where the transition in slope of the toroidal energy spectrum occurs is varied between $l=2$ and $l=8$. For each case the magnitudes of both the poloidal and the toroidal spectra are evaluated as well as their different slopes. From Table 2, one can observe that shifting the toroidal transition scale toward small scales coincides with a diminution of the spectral slope within the spherical harmonic range $\Delta_{\psi_1}$, whereas the slope associated with the spectrum at small scales remains always close to 0. This phenomenon may indicate that the toroidal field can be separated into two parts, one very energetic at large scales (probably the large-scale eccentric gyre), and another one presenting an almost flat spectrum at small scales. The computation of $\log(p(\mathcal{M}|\Gamma))$ for each simulation, indicates that the degree where the transition in slope is most likely to happen is $l=3$, corresponding to simulation III. In the final ensemble of simulations (simulations IX to XX), both the poloidal and toroidal energy spectra are assumed to present two different ranges. Whereas the change in slope in the toroidal case is always prescribed around SH degree $l=3$, values between $l=4$ and $l=9$ are tested for the poloidal case. It can first be observed that splitting the toroidal spectrum into the two ranges $1 \le l \le 3$ and $3 \le l \le 26$ always present a higher posterior probability density than any other decomposition. Then, the most probable scale at which the slope of the poloidal spectrum is changing is SH degree $l=8$.
<table><caption>**Table 2.** Optimal Spectra Parameters $\hat{M}$ for Different Prior Decompositions of the Time-Averaged Toroidal and Poloidal Energy Spectra Within the 2001–2009.5 Period<sup>a</sup></caption><thead><tr><th>Simulation</th><th>Range Index i</th><th>Δ<sub>φ<sub>i</sub></sub></th><th>A<sub>φ</sub></th><th>P<sub>φ<sub>i</sub></sub></th><th>Δ<sub>ψ<sub>i</sub></sub></th><th>A<sub>ψ</sub></th><th>P<sub>ψ<sub>i</sub></sub></th><th>log (p(<i>M</i>|γ<sub>0</sub>)) + c</th></tr></thead><tbody><tr><td>I</td><td>1</td><td>[1 26]</td><td>1.24</td><td>-0.76</td><td>[1 26]</td><td>20.2</td><td>0.92</td><td>0</td></tr><tr><td rowspan="2">II</td><td>1</td><td>[1 26]</td><td>1.36</td><td>-0.65</td><td>[1 2]</td><td>26.2</td><td>2.19</td><td rowspan="2">21.3</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[2 26]</td><td>15.37</td><td>0.65</td></tr><tr><td rowspan="2">III</td><td>1</td><td>[1 26]</td><td>1.61</td><td>-0.44</td><td>[1 3]</td><td>34.05</td><td>2.78</td><td rowspan="2">89.5</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[3 26]</td><td>7.45</td><td>0.01</td></tr><tr><td rowspan="2">IV</td><td>1</td><td>[1 26]</td><td>1.66</td><td>-0.42</td><td>[1 4]</td><td>27.45</td><td>1.95</td><td rowspan="2">63.2</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[4 26]</td><td>6.92</td><td>-0.04</td></tr><tr><td rowspan="2">V</td><td>1</td><td>[1 26]</td><td>1.68</td><td>-0.37</td><td>[1 5]</td><td>25.67</td><td>1.69</td><td rowspan="2">57</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[5 26]</td><td>5.5</td><td>-0.23</td></tr><tr><td rowspan="2">VI</td><td>1</td><td>[1 26]</td><td>1.68</td><td>-0.4</td><td>[1 6]</td><td>23.2</td><td>1.39</td><td rowspan="2">31.6</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[6 26]</td><td>5.55</td><td>-0.2</td></tr><tr><td rowspan="2">VII</td><td>1</td><td>[1 26]</td><td>1.61</td><td>-0.44</td><td>[1 7]</td><td>22.23</td><td>1.23</td><td rowspan="2">21.5</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[7 26]</td><td>5.66</td><td>-0.17</td></tr><tr><td rowspan="2">VIII</td><td>1</td><td>[1 26]</td><td>1.65</td><td>-0.43</td><td>[1 8]</td><td>21.73</td><td>1.17</td><td rowspan="2">19.9</td></tr><tr><td>2</td><td>-</td><td>-</td><td>-</td><td>[2 26]</td><td>4.74</td><td>-0.29</td></tr><tr><td rowspan="2">IX</td><td>1</td><td>[1 4]</td><td>0.28</td><td>-3.31</td><td>[1 3]</td><td>33.75</td><td>2.77</td><td rowspan="2">107.3</td></tr><tr><td>2</td><td>[4 26]</td><td>3.28</td><td>0.23</td><td>[3 26]</td><td>7.26</td><td>-0.03</td></tr><tr><td rowspan="2">X</td><td>1</td><td>[1 5]</td><td>0.3</td><td>-2.94</td><td>[1 3]</td><td>33.98</td><td>2.8</td><td rowspan="2">117.9</td></tr><tr><td>2</td><td>[5 26]</td><td>5.99</td><td>0.77</td><td>[3 26]</td><td>7.1</td><td>-0.05</td></tr><tr><td rowspan="2">XI</td><td>1</td><td>[1 6]</td><td>0.37</td><td>-2.54</td><td>[1 3]</td><td>34.15</td><td>2.82</td><td rowspan="2">129</td></tr><tr><td>2</td><td>[6 26]</td><td>13.78</td><td>1.51</td><td>[3 26]</td><td>6.96</td><td>-0.07</td></tr><tr><td rowspan="2">XII</td><td>1</td><td>[1 7]</td><td>0.43</td><td>-2.23</td><td>[1 2]</td><td>25.86</td><td>2.22</td><td rowspan="2">54.5</td></tr><tr><td>2</td><td>[7 26]</td><td>26</td><td>1.98</td><td>[2 26]</td><td>14.46</td><td>0.54</td></tr><tr><td rowspan="2">XIII</td><td>1</td><td>[1 7]</td><td>0.48</td><td>-2.14</td><td>[1 3]</td><td>34.64</td><td>2.87</td><td rowspan="2">136.1</td></tr><tr><td>2</td><td>[7 26]</td><td>51.51</td><td>2.68</td><td>[3 26]</td><td>6.76</td><td>-0.1</td></tr><tr><td rowspan="2">XIV</td><td>1</td><td>[1 7]</td><td>0.5</td><td>-2.07</td><td>[1 4]</td><td>27.42</td><td>1.98</td><td rowspan="2">107.7</td></tr><tr><td>2</td><td>[7 26]</td><td>49.98</td><td>2.65</td><td>[4 26]</td><td>6.24</td><td>-0.16</td></tr><tr><td rowspan="2">XV</td><td>1</td><td>[1 8]</td><td>0.52</td><td>-1.96</td><td>[1 2]</td><td>25.87</td><td>2.22</td><td rowspan="2">54.8</td></tr><tr><td>2</td><td>[8 26]</td><td>192.21</td><td>3.73</td><td>[2 26]</td><td>14.37</td><td>0.53</td></tr><tr><td rowspan="2">XVI</td><td><b>1</b></td><td><b>[1 8]</b></td><td><b>0.62</b></td><td><b>-1.77</b></td><td><b>[1 3]</b></td><td><b>34.67</b></td><td><b>2.87</b></td><td rowspan="2"><b>136.7</b></td></tr><tr><td><b>2</b></td><td><b>[8 26]</b></td><td><b>427.59</b></td><td><b>4.52</b></td><td><b>[3 26]</b></td><td><b>6.76</b></td><td><b>-0.1</b></td></tr><tr><td rowspan="2">XVII</td><td>1</td><td>[1 8]</td><td>0.64</td><td>-1.72</td><td>[1 4]</td><td>27.46</td><td>1.97</td><td rowspan="2">108.6</td></tr><tr><td>2</td><td>[8 26]</td><td>387.53</td><td>4.43</td><td>[4 26]</td><td>6.24</td><td>-0.16</td></tr><tr><td rowspan="2">XVIII</td><td>1</td><td>[1 9]</td><td>0.74</td><td>-1.53</td><td>[1 2]</td><td>25.82</td><td>2.22</td><td rowspan="2">50.5</td></tr><tr><td>2</td><td>[9 26]</td><td>17471.7</td><td>7.63</td><td>[2 26]</td><td>14.25</td><td>0.51</td></tr><tr><td rowspan="2">XIX</td><td>1</td><td>[1 9]</td><td>0.77</td><td>-1.45</td><td>[1 3]</td><td>35.23</td><td>2.91</td><td rowspan="2">136</td></tr><tr><td>2</td><td>[9 26]</td><td>174982</td><td>9.77</td><td>[3 26]</td><td>6.61</td><td>-0.13</td></tr><tr><td rowspan="2">XX</td><td>1</td><td>[1 9]</td><td>0.84</td><td>-1.38</td><td>[1 4]</td><td>27.89</td><td>2.02</td><td rowspan="2">110.4</td></tr><tr><td>2</td><td>[9 26]</td><td>193315</td><td>9.85</td><td>[4 26]</td><td>5.97</td><td>-0.21</td></tr></tbody></table>

<footnote><sup>a</sup>A corresponds to the magnitudes and P to the slopes of the poloidal (with φ subscripts) and toroidal (with ψ subscripts) energy spectra. The last column provides the posterior probability density of the spectra parameters (the larger the density, the most likely the decomposition), and the constant c is given by – log (p(<i>M</i>|γ<sub>0</sub>)) of simulation I. Note that for poloidal and toroidal fields expanded in Schmidt seminormalized spherical harmonics, the magnitudes given in this table have to be divided by a factor of $\sqrt{4\pi}$ (the slopes remain similar). Bold data are associated with the optimal decomposition.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (top) Toroidal and (bottom) poloidal, energy spectra (black lines), prior spectra (grey lines), and rescaled scale-by-scale variances (circles), associated with the velocity field $\hat{u}$ at the epoch 2005.0.</figcaption></figure>

To provide confidence intervals on the spectra parameters of simulation XVI (the one presenting the highest posterior probability density), we numerically evaluated the standard deviation:

$$ \sigma_M = \left[ \int (M - \hat{M})^2 \rho(M) dM \right]^{\frac{1}{2}} \qquad (56) $$

with a Monte Carlo integration algorithm. Combining this uncertainty estimation to the results of simulations I to XX, the following conclusions on the spectral behavior of the poloidal and toroidal flows can be drawn:

1. The toroidal energy spectrum has a magnitude of 34.7 ± 2.7 and can be characterized by two ranges, one between 1 ≤ l ≤ 3 where the spectrum is decaying with a slope of 2.9 ± 0.2 and one between 3 ≤ l ≤ 26 where the spectrum is almost flat with a slope of -0.1 ± 0.1.

2. The poloidal energy spectrum which is characterized by a magnitude of 0.6 ± 0.1 is first increasing with a slope of 1.8 ± 0.2 to reach its maximum at SH degree l = 8 where it starts to decay with a strong slope of 4.5 ± 1.1.

### 3.5. The Velocity Field at the Core-Mantle Boundary and Its Uncertainties in 2005.0

The optimal spectra parameters found in the previous section, and referred to $\hat{M}$, are now used as a prior information for the estimation of the flow at the CMB at the epoch 2005.0. For this study we generated an ensemble of N = 42, 500 (u, b) fields. The average velocity field of this ensemble is referred to $\hat{u}$.

In spectral space, the posterior variances of the poloidal and toroidal fields are computed for each degree l. These latter are then rescaled by the factor l(l + 1) and plotted in Figure 6 (circles) together with the spectra of $\hat{u}$ (black lines) and with the prior spectra of the velocity field (gray lines). In comparison with the results of Baerenzung et al. [2014], where the same analysis was performed, the posterior variance of both poloidal and toroidal fields is much more important in this study. This is certainly due to the choice of the prior information on the flow, which in Baerenzung et al. [2014] was imposing a much stronger slope to the kinetic energy spectrum than in the present case. Note that the strength of the penalty on the small scales imposed by the prior
<figure><img src="image_8.png" /></figure>

**Figure 7.** (top) Toroidal field (black arrows) and (bottom) poloidal field (black arrows) for the epoch 2005.0 together with the 90% confidence interval on their orientation (color maps). The blue and violet contours, respectively, correspond to the boundary of patches within which the tangential-geostrophic constraint and the surface incompressible quasi-geostrophic constraint are violated with a 90% confidence. The white circles in the Northern and Southern Hemispheres correspond to the intersections between the tangent cylinder and the CMB.

velocity field spectrum is not the only possible factor reducing the posterior variability of the solution; indeed, forcing the flow to fulfill some physical constraint will also affect the posterior variance of the velocity field. This effect can be observed in the flow obtained by Gillet *et al.* [2015], where although their ensemble of flow solution suggests a flat prior spectrum, the quasi geostrophic constraint applied to the velocity field seems to drastically reduce the variability of the solution. Going back to Figure 6, one can observe on its upper part that only a few scales of the toroidal field (SH degree $l=1, 2, 4, 6$) present a lower or equivalent variance level than the energy they exhibit. For the poloidal field (Figure 6, bottom) the situation is even more dramatic, since the energy contained in any scale is always lower than the variance. In the best case (SH degree $l=7$ and $l=8$), these two quantities are of the same order of magnitude. Another observation one can make with Figure 6 is that above SH degree $l=13$, the posterior distribution of the flow is completely dominated by prior properties. This is why the velocity field in physical space will only be expanded up to SH degree $l=13$.

In physical space, Figure 7 displays, with black arrows, the average toroidal (top) and poloidal (bottom) fields, together with the 90% confidence interval on their orientation in degree (color maps). The color transition from violet to red is associated with, respectively, reliable to uncertain flow orientations. The locations of reliable or uncertain flow patterns strongly differ between the poloidal and the toroidal fields. Indeed, the most robust feature of the toroidal field is the large planetary scale eccentric gyre already identified in previous studies with various prior constraints [see Pais and Jault, 2008; Gillet *et al.*, 2009; Aubert, 2013], whereas only part of the poloidal field in the Eastern Hemisphere can be accurately estimated. In agreement with the results of Baerenzung *et al.* [2014], the flow at the level of the CMB below the Pacific Ocean, the northern and southern part of the Atlantic and Europe cannot be precisely determined through the inversion of the FF equation. This means, in particular, that the counterclockwise gyre below the tip of Africa or the clockwise gyre east of Québec, recurrently reported [see Finlay *et al.*, 2010; Aubert, 2013], and which are also present in the flow we recovered, are not necessary to explain the observed secular variation. The toroidal field we obtained also contradicts many prior physical assumptions usually employed in the inverse problem. In particular, the extremely robust parts of the flow crossing the equator below India and Brazil violate the condition imposed
by quasi-geostrophy, or tangential-geostrophy, which are forcing the polar component of the velocity field to vanish at the equator. Another feature contradicting the geostrophic assumption is the clear flow penetration of the cylinder tangent to the inner core (in both the Southern and Northern hemispheres) whose intersections with the CMB are represented with white lines in Figure 7. Finally, the property of equatorial symmetry of the flow requested by quasi-geostrophy is also not fulfilled by the velocity field we calculated, the intensity of the flow in the Southern Hemisphere below Africa and below the Atlantic, being larger than in the Northern Hemisphere. Note that this feature is usually observed when equatorial symmetry is not a priori imposed to the velocity field as in *Amit and Pais* [2013].

Recently, Aubert [2013]; Aubert et al. [2013] has shown that the eccentric nature of the planetary scale gyre was due to a difference in buoyancy release between the Eastern Hemisphere and the remaining part of the inner core boundary (ICB). In their 3-D simulations of the geodynamo, they observed that imposing a longitudinal anomalous mass flux on the ICB presenting a maximum beneath Indonesia, allowed them to recover the particular localizations of the geomagnetic secular variation. This heterogeneous forcing of the flow leads to fast upwelling movements of fluid particles above locations where it is applied. A signature of this mechanism can be observed in the flow we estimated. Indeed, upwelling movements of the fluid in the outer core are converted into poloidal field when meeting the core-mantle boundary. Yet the poloidal field we obtained in eastern regions is much more energetic, with some velocities reaching 7 km yr⁻¹, than in the remaining part of the CMB, where the poloidal velocities do not exceed 2 km yr⁻¹. The presence of such intense and reliable poloidal fields indicates that if there were to be a stratified layer at the top of outer core, as proposed by Pozzo et al. [2012] and Gubbins and Davies [2013], the latter would have only weak effects on large-scale upwelling and downwelling motions in the outer core and would only impact small-scale structures as predicted by Vidal and Schaeffer [2015].

Thanks to the ensemble of velocity fields we possess, we can estimate that the reliability of particular features û is exhibiting. In the following we measure the posterior probability of different properties characterizing some particular behaviors of the flow.

*Equatorial Symmetry.* The toroidal part of û clearly presents strong equatorial symmetries. The range of [0.72; 0.84] corresponding to the 90% confidence interval on the ratio between the energy of the equatorial symmetric and total toroidal field confirm this aspect of the flow. This particular behavior of the velocity field is less pronounced for its poloidal part. Indeed, the 90% confidence interval for the ratio between the energy of the equatorial symmetric and total poloidal field is of [0.47; 0.67].

*Northern and Southern Intensity of the Westward Drift.* We define the westward drift intensity as the root mean square azimuthal velocity lying between longitudes ±60° and colatitudes ±30°, its northern part $u_N$ and southern part $u_S$ being delimited by the equator. For 94.2% of the velocity fields belonging to the posterior ensemble, $u_S$ is greater than $u_N$. To quantify the departure from symmetry of the westward drift, the posterior distribution of the ratio $u_N/u_S$ has been computed. The fact that this distribution is characterized by a maximum at $u_N/u_S = 0.8$, and a 90% confidence interval of [0.63; 1.00], confirms that the westward drift is very likely to be more intense in the Southern Hemisphere than in the northern one.

*Geostrophy.* To evaluate if the geostrophy assumption is valid everywhere in the outer core, we studied the posterior distribution of the following two quantities: $∇_H(u \cos \theta)$ associated with the tangential-geostrophy (TG) constraint and $∇_H(u_{es} \cos^2 \theta)$, where $u_{es}$ is the equatorial symmetric part of the velocity field, associated with the surface incompressible quasi-geostrophic (SIQG) constraint. For the TG or the SIQG constraints to be fulfilled, either $∇_H(u \cos \theta)$ or $∇_H(u_{es} \cos^2 \theta)$ has to vanish at the core-mantle boundary. Contours of patches in which the 0 value does not belong to the 90% confidence intervals of $∇_H(u \cos \theta)$ and $∇_H(u_{es} \cos^2 \theta)$ are shown in Figure 7 (bottom). Blue and violet contours, respectively, correspond to locations in which the TG and SIQG constraints are violated. These locations mostly coincide with areas of intense upwelling or downwelling motions of the outer core fluid. This may indicate that the difference in buoyancy release between the Eastern Hemisphere and the remaining part of the inner core boundary could break the geostrophic state of the outer core.

## 4. Conclusion

In this study, we proposed a consistent formalism to determine both the poloidal and toroidal energy spectra of the velocity field evolving at the Earth's core-mantle boundary. We showed that using these spectra
as a prior information in the inversion of the frozen flux approximation, together with employing an appropriate sampling algorithm, allowed recovering realistic posterior statistics of the flow. Applying our method to real data allowed us to extract different information on the velocity field at the core-mantle boundary. First, the large-scale approximation as it is usually considered cannot be applied in a Bayesian framework to determine the flow between 2001.0 and 2009.5. Indeed, the toroidal field, which dominantly contributes to the total velocity field, exhibits an energy spectrum characterized by two ranges: one between spherical harmonic degrees $l = 1$ and $l = 3$ where the spectrum is strongly decaying, in agreement with the large-scale approximation, but also one above degree $l = 3$, contradicting the assumption. The large-scale assumption is also violated by the poloidal field which presents an energy spectrum increasing between spherical harmonics degrees $l = 1$ and $l = 8$. Then, we could estimate the reliability of certain features of the flow. Whereas the toroidal part of the velocity field is dominantly equatorial symmetric, it is very unlikely to be perfectly symmetric, as reflected by the difference in intensity between the northern and southern part of the westward drift. Finally, we could demonstrate that if the flow was to be geostrophic in the outer core, this property is probably breaking in locations where intense upwelling and downwelling motions of fluid particles occur.

From the results of this study, improvements on the methodology could be considered. In particular, the strong degree of equatorial symmetry of the flow suggests that its symmetric and antisymmetric parts could be treated distinctively a priori. Accounting for the temporal dependency of the velocity field in the inverse problem is also clearly lacking in this study. If well characterized, time correlations would certainly improve the quality of the flow solution by reducing its posterior variability. The method we propose, because of its numerical efficiency and mathematical consistency, could provide a solid basis to tackle the time-dependant problem.

## Acknowledgments

This work has been supported by the German Research Foundation (DFG) within the Priority Program SPP1488 "Planetary Magnetism." The data for this paper are available by contacting the corresponding author at baerenzu@uni-potsdam.de. This is IPGP contribution 3720.

## References

* Amit, H., and M. A. Pais (2013), Differences between tangential geostrophy and columnar flow, *Geophys. J. Int.*, **194**, 145–157, doi:10.1093/gji/ggt077.

* Aubert, J. (2013), Flow throughout the Earth's core inverted from geomagnetic observations and numerical dynamo models, *Geophys. J. Int.*, **192**, 537–556, doi:10.1093/gji/ggs051.

* Aubert, J. (2014), Earth's core internal dynamics 1840–2010 imaged by inverse geodynamo modelling, *Geophys. J. Int.*, **197**, 1321–1334, doi:10.1093/gji/ggu064.

* Aubert, J., and A. Fournier (2011), Inferring internal properties of Earth's core dynamics and their evolution from surface observations and a numerical geodynamo model, *Nonlinear Processes Geophys.*, **18**, 657–674, doi:10.5194/npg-18-657-2011.

* Aubert, J., C. Christopher, and A. Fournier (2013), Bottom-up control of geomagnetic secular variation by the Earth's inner core, *Nature*, **502**, 219–223, doi:10.1038/nature12574.

* Backus, G., R. Parker, and C. Constable (1996), *Foundations of Geomagnetism*, Cambridge Univ. Press, San Diego, Calif.

* Baerenzung, J., M. Holschneider, and V. Lesur (2014), Bayesian inversion for the filtered flow at the Earth's core mantle boundary, *J. Atmos. Sci.*, **68**, 2757–2770, doi:10.1175/2010JAS3445.1.

* Bayes, T. (1763), An essay towards solving a problem in the doctrine of chances, *Philos. Trans. R. Soc. London*, **53**, 370–418.

* Buffett, B. A., and U. R. Christensen (2007), Magnetic and viscous coupling at the core-mantle boundary: Inferences from observations of the Earth's nutations, *Geophys. J. Int.*, **171**, 145–152, doi:10.1111/j.1365-246X.2007.03543.x.

* Celaya, M., and J. Wahr (1996), Aliasing and noise in core-surface flow inversions, *Geophys. J. Int.*, **126**, 447–469, doi:10.1111/j.1365-246X.1996.tb05302.x.

* Christensen, U. R., J. Aubert, and G. Hulot (2010), Conditions for Earth-like geodynamo models, *Earth Planet. Sci. Lett.*, **296**, 487–496, doi:10.1016/j.epsl.2010.06.009.

* Chulliat, A., and G. Hulot (2000), Local computation of the geostrophic pressure at the top of the core, *Phys. Earth Planet. Inter.*, **117**, 309–328, doi:10.1016/S0031-9201(99)00104-1.

* Eymin, C., and G. Hulot (2005), On core surface flows inferred from satellite magnetic data, *Phys. Earth Planet. Inter.*, **152**, 200–220, doi:10.1016/j.pepi.2005.06.009.

* Finlay, C. C., M. Dumberry, A. Chulliat, and M. A. Pais (2010), Short timescale core dynamics: Theory and observations, *Space Sci. Rev.*, **155**, 177–218, doi:10.1007/s11214-010-9691-6.

* Fournier, A., J. Aubert, and E. Thébault (2011), Inference on core surface flow from observations and 3-D dynamo modelling, *Geophys. J. Int.*, **186**, 118–136, doi:10.1111/j.1365-246X.2011.05032.x.

* Gamerman, D., and H. F. Lopes (2006), *Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference*, Chapman and Hall/CRC Texts in Statistical Science Series, Taylor and Francis.

* Gillet, N., M. A. Pais, and D. Jault (2009), Ensemble inversion of time-dependent core flow models, *Geochem. Geophys. Geosyst.*, **10**, Q06004, doi:10.1029/2008GC002290.

* Gillet, N., D. Jault, and C. C. Finlay (2015), Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth's core surface, *J. Geophys. Res. Solid Earth*, **120**, 3991–4013, doi:10.1002/2014JB011786.

* Gubbins, D., and C. J. Davies (2013), The stratified layer at the core-mantle boundary caused by barodiffusion of oxygen, sulphur and silicon, *Phys. Earth Planet. Inter.*, **215**, 21–28, doi:10.1016/j.pepi.2012.11.001.

* Holme, R. (2007), Large-scale flow in the core, in *Treatise on Geophysics*, edited by G. Schubert, pp. 107–130, Elsevier, Amsterdam, doi:10.1016/B978-044452748-6.00127-9.

* Jackson, A. (1995), An approach to estimation problems containing uncertain parameters, *Phys. Earth Planet. Inter.*, **90**, 145–156, doi:10.1016/0031-9201(95)05080-U.

* Jault, D., and J. L. Le Mouelé (1991), Physical properties at the top of the core and core surface motions, *Phys. Earth Planet. Inter.*, **68**, 76–84, doi:10.1016/0031-9201(91)90009-7.
Jault, D. (2015), Illuminating the electrical conductivity of the lowermost mantle from below, *Geophys. J. Int.*, **202**, 482–496, doi:10.1093/gji/ggv152.

Kono, M., and P. H. Roberts (2002), Recent geodynamo simulations and observations of the geomagnetic field, *Rev. Geophys.*, **40**, 1–14–53, doi:10.1029/2000RG000102.

Lesur, V., and I. Wardinski (2012), Evaluating secular acceleration in geomagnetic field model GRIMM-3, Abstract GP24A-08 presented at 2012 Fall Meeting AGU, San Francisco, Calif.

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

Pais, A., and G. Hulot (2000), Length of day decade variations, torsional oscillations and inner core superrotation: Evidence from recovered core surface zonal flows, *Phys. Earth Planet. Inter.*, **118**, 291–316, doi:10.1016/S0031-9201(99)00161-2.

Pais, M. A., and D. Jault (2008), Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field, *Geophys. J. Int.*, **173**, 421–443, doi:10.1111/j.1365-246X.2008.03741.x.

Pozzo, M., C. Davies, D. Gubbins, and D. Alfé (2012), Thermal and electrical conductivity of iron at Earth's core conditions, *Nature*, **485**, 355–358, doi:10.1038/nature11031.

Roberts, P. H., and S. Scott (1965), On analysis of the secular variation, *J. Geomagn. Geoelec.*, **17**(2), 137–151, doi:10.5636/jgg.17.137.

Thébault, E., and F. Vervelidou (2015), A statistical spatial power spectrum of the Earth's lithospheric magnetic field, *Geophys. J. Int.*, **201**, 605–620, doi:10.1093/gji/ggu463.

Velimský, J. (2010), Electrical conductivity in the lower mantle: Constraints from CHAMP satellite data by time-domain EM induction modelling, *Phys. Earth Planet. Inter.*, **180**, 111–117, doi:10.1016/j.pepi.2010.02.007.

Vidal, J., and N. Schaeffer (2015), Quasi-geostrophic modes in the Earth's fluid core with an outer stably stratified layer, *Geophys. J. Int.*, **202**(3), 2182–2193, doi:10.1093/gji/ggv282.