Bayesian inversion for the filtered flow at the Earth’s Core Mantle Boundary

J. Baerenzung¹, M. Holschneider¹, V. Lesur² ¹ Institute for Mathematics, University of Potsdam, Potsdam, Germany
² Helmholtz Center Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany

Abstract

The inverse problem which consists of determining the flow at the Earth’s Core Mantle Boundary according to an outer core magnetic field and secular variation model, has been investigated through a Bayesian formalism. To circumvent the issue arising from the truncated nature of the available fields, we combined two modelization methods. In the first step, we applied a filter on the magnetic field to isolate its large scales by reducing the energy contained in its small scales, we then derived the dynamical equation, referred as filtered frozen flux equation, describing the spatio-temporal evolution of the filtered part of the field. In the second step, we proposed a statistical parametrization of the filtered magnetic field in order to account for both its remaining unresolved scales and its large scale uncertainties. These two modelization techniques were then included in the Bayesian formulation of the inverse problem. To explore the complex posterior distribution of the velocity field resulting from this development, we numerically implemented an algorithm based on Markov Chain Monte Carlo methods. After evaluating our approach on synthetic data and comparing it to previously introduced methods, we applied it on real data for the single epoch $2005.0$ .

I Introduction

The core magnetic field (MF) of the Earth is sustained by the dynamo action taking place in its outer core. Here, the variations in chemical composition and in temperature of the liquid metal allows convection to develop. Since the fluid is electrically conducting, it interacts non-linearly with the magnetic field. While the flow is advecting the MF, this latter is constraining the fluid motions through the Lorentz forces. The evolution of the fluid velocity field (VF) and the magnetic field are therefore entirely connected to one another through energy exchanges between them.

Studying such a system is difficult in many aspects. From a numerical point of view, simulating directly the dynamic of the outer core is a challenging task. Because of the strong regime of turbulence that the magnetohydrodynamic (MHD) flow exhibits, the separation between the smallest and the largest scales of the system is extremely broad. Yet, to properly describe the evolution of both the VF and the MF, all these scales should be considered in simulations since they interact non-linearly together. But this is, with the actual computation power available, impossible.

From an experimental point of view, observing directly the evolution of the outer core is also impossible. Nevertheless, measurements of the MF from ground observatories or satellites may allow to indirectly infer some dynamical properties of the flow and the MF in the core. In particular, because of the low conductivity of the mantle (Velímský (2010)), a knowledge of the core MF at the Earth’s surface is sufficient to evaluate it at the level of the core-mantle boundary (CMB). At this very specific location, the MF is coupled to the outer core VF through the frozen flux equation (a simplified version of the induction equation introduced by Roberts and Scott (1965), and in which the diffusion effects are neglected). By inverting this equation it is therefore possible to evaluate the VF at the CMB. Unfortunately, the problem is ill-posed for different reasons. First, the two components of the velocity field are connected to the radial component of the secular variation (SV) through a unique equation. Then, the available secular variation given by MF models such as GRIMM 2 of Lesur et al. (2010) is only resolved at large scales whereas any scale of the velocity field can contribute to this resolved SV. Finally, as it is the case for the SV, only the large scale core MF can be determined at the Earth surface, yet interactions between the unknown small scale MF with the VF can generate large scale SV (see Eymin and Hulot (2005)).

As shown by Backus (1968), to resolve the non-uniqueness of the velocity field in this inverse problem, constraints have to be imposed to the flow behavior. Different formulations have been proposed over the past decades. This includes tangential-geostrophy, tangential-magnetostrophy, columnar flow, helical flow, or purely toroidal flow (see Holme (2007) and Finlay et al. (2010) for an exhaustive review of the different constraints usually applied and their physical implications). Nevertheless, these physical constraints are not sufficient to provide a unique flow solution (see Chulliat and Hulot (2000) and references therein), and additional regularization assumptions have to be introduced in the problem (see Holme (2007)).

Although identified for more than twenty years (Hulot et al. (1992)), it is only recently that the effects of the unknown small scale MF on the large scale SV have been modeled in the inversion of the Frozen Flux approximation. In particular two methods have given promising results, namely the ensemble approach of Gillet et al. (2009), and the iterative method of Pais and Jault (2008). The philosophy of these two approaches is quite distinct from one another. Whereas in the method of Gillet et al. (2009) an ensemble of magnetic field containing small scales is generated and directly used to evaluate an ensemble of velocity field, in the method of Pais and Jault (2008) the effects of the unknown magnetic field is transposed into a modeling error which is iteratively estimated. In this study, we propose a development of these approaches in the context of Bayesian modeling. The article is therefore constructed in the following manner.

In section II, after describing the principle of derivation of the Frozen Flux equation, we recall the issue raised in the inversion of this equation by the truncated nature of the available MF. We then describe how to isolate the large scales of the MF and present the approximated dynamical equation, referred as Filtered Frozen Flux (FFF) equation, which determines the evolution of these large scales. At the end of the section we present how to formulate the inverse problem in a Bayesian framework when variations of the MF around the prescribed one are allowed to occur. In section III tests on synthetic data are performed. In the first one, we evaluate the improvement brought by using the FFF equation in the inverse problem. In the second test, the Bayesian formalism developed in section II is considered to recover the velocity field from a set of artificially generated data, and the results are compared to the flow obtained with three other approaches. The methodology we developed is then used to determine the velocity field and its underlying uncertainties for the epoch $2005.0$ . Finally we present our conclusions in section IV.

II Governing equations

II.1 The Frozen-Flux approximation

In this section we introduce briefly the hypothesis necessary to derive the Frozen Flux approximation. For a more detailed description see Holme (2007).

Outside the core, the geodynamo’s MF $\bf B$ is irrotational, it can therefore be expressed through a potential $\phi$ according to the relation:

{\bf B}=-{\bf\nabla}\phi\ .

(1)

As mentioned previously, the low conductivity of the mantle (Velímský (2010)) allows to evaluate the MF and therefore its radial component $B_{r}=-\partial_{r}\phi$ at the level of the core mantle boundary. There its evolution is prescribed by the induction equation:

\partial_{t}B_{r}=-{\bf\nabla_{{}_{H}}}({\bf u}{B}_{r})+\eta\left(\Delta{\bf B}\right)\cdot{\bf e_{r}}\ ,

(2)

with ${\bf\nabla_{H}}$ the horizontal divergence operator, ${\bf u}$ the two-dimensional velocity field, $\eta$ the magnetic diffusivity, and ${\bf e_{r}}$ a radial unitary vector. Because for the Earth, on short period of time the dissipation effects are dominated by advection effects (see Holme (2007)) and since our study is limited to single epoch inversion, the fluid can be considered as a perfect conductor. Roberts and Scott (1965) showed that under this assumption, known as the Frozen Flux (FF) approximation, the induction equation can be simplified as:

\partial_{t}{B}_{r}=-{\bf\nabla_{{}_{H}}}({\bf u}{B}_{r})\ .

(3)

II.2 The unresolved scale issue

To perform a consistent inversion of the Frozen-Flux equation at the CMB, in addition to imposing a certain behavior to the flow and to regularizing it (see Holme (2007); Finlay et al. (2010)), every scale composing the secular variation and the magnetic field have to be known. Unfortunately, in the actual models describing the spatio-temporal evolution of the Earth’s magnetic field derived from satellite and observatory data, the resolution of both the SV and the MF is limited. For the GRIMM 2 model of Lesur et al. (2010), for example, the fields do not excess the degree $13$ when expanded in spherical harmonics. So if $g_{l,m}$ corresponds to the spherical harmonics (SH) coefficient at degree $l$ and order $m$ associated with the scalar potential $\phi$ , such as:

\phi=R\sum_{l=0}^{l=+\infty}\sum_{m=-l}^{m=+l}\left(\frac{c}{r}\right)^{l+1}g_{l,m}Y_{l,m}\ ,

(4)

where $R$ is the core radius, and $Y_{l,m}$ is the Schmidt semi-normalized SH, then according to equation (1), the available MF $B_{r}^{<}$ , and its unknown part $B_{r}^{>}$ are respectively given by:

	$\displaystyle B_{r}^{<}$	$\displaystyle=$	$\displaystyle-\sum_{l=1}^{l=lc}(l+1)\sum_{m=-l}^{m=+l}g_{l,m}^{<}Y_{l,m}$		(5)
	$\displaystyle B_{r}^{>}$	$\displaystyle=$	$\displaystyle-\sum_{l=lc+1}^{l=+\infty}(l+1)\sum_{m=-l}^{m=+l}g_{l,m}^{>}Y_{l,m}\ ,$		(6)

with $lc$ the cut-off scale (which is equal to $13$ for the MF provided by the GRIMM 2 model). Note that the total radial component of the MF corresponds to the sum of these two quantities.

To account for the truncated nature of the available MF and SV, the Frozen Flux approximation has to be rewritten as:

	$\displaystyle\partial_{t}{B}_{r}^{<}$	$\displaystyle=$	$\displaystyle-\left({\bf\nabla_{{}_{H}}}({\bf u}{B}_{r})\right)^{<}$		(7)
		$\displaystyle=$	$\displaystyle-\left({\bf\nabla_{{}_{H}}}\left({\bf u}{B}_{r}^{<}\right)\right)^{<}-\left({\bf\nabla_{{}_{H}}}({\bf u}{B}_{r}^{>})\right)^{<}\ ,$		(8)

where the advection term, on the right hand side of equation (8), is decomposed into two parts, one depending on the resolved magnetic field and the other function of the undetermined field. Hulot et al. (1992) were the first to highlight the issue raised by the unknown part of the MF in the inversion of the FF equation. Nevertheless, because at this epoch the uncertainties on SV measurements were large, the contribution of the term $\left({\bf\nabla_{{}_{H}}}({\bf u}{B}_{r}^{>})\right)^{<}$ could be neglected in the inverse problem. The recent increase in quality of both the measurements and models describing the evolution of the core MF (see Hulot et al. (2002)) invalidate this latter statement, and Eymin and Hulot (2005) showed that the unresolved part of the MF could not be neglected anymore.

II.3 Parametrization of the unresolved magnetic field

To model the effects of the unresolved part of the magnetic field on the large scale secular variation, assumptions on this unknown field behavior have to be made. In particular one can prescribe to it a certain energy spectrum. This operation can be performed, for instance, by extrapolating the spectrum associated with the resolved scales which is defined by:

E_{B^{<}}(l)=(l+1)\sum_{m=-l}^{m=l}\left(g_{l,m}^{<}\right)^{2}\ .

(9)

The resulting spectrum $E_{B^{>}}$ can then be used to statistically model the unknown MF $B_{r}^{>}$ . Following the development of Hulot et al. (1992) where the field is assumed to be isotropically distributed, the covariance of the coefficients $g_{l,m}^{>}$ is directly proportional to the extrapolated spectrum through the relation:

E[g^{>}_{l,m},g^{>}_{l^{\prime},m^{\prime}}]=\frac{E_{B^{>}}(l)}{(l+1)(2l+1)}\delta(l-l^{\prime})\delta(m-m^{\prime})\ ,

(10)

where $E[...]$ corresponds to the mathematical expectation.

The issue with such a modelization is that no universal spectrum of the geodynamo’s MF at the core mantle boundary is available. Nevertheless, different formulations have been proposed over the past decades, but unfortunately most of them strongly differ from one another. An illustration of this statement is given in figure 1 where the spectrum of the MF prescribed by the GRIMM 2 model for the epoch $2005.0$ (between SH degree $1\leq l\leq 13$ ) together with three different extrapolations (thin lines) are plotted. In this example, the laws derived by Buffett and Christensen (2007), Roberts et al. (2003) and Voorhies (2004) were used to extrapolate the resolved scale spectrum. They respectively read:

$\displaystyle E_{B}^{B}(l)$	$\displaystyle=$	$\displaystyle C_{1}\chi^{l}$	(11)
$\displaystyle E_{B}^{R}(l)$	$\displaystyle=$	$\displaystyle C_{2}e^{-Sl}$	(12)
$\displaystyle E_{B}^{V}(l)$	$\displaystyle=$	$\displaystyle C_{3}\frac{l+1/2}{l(l+1)}\,$	(13)

with $\chi=0.99$ and $S=0.055$ . The constants $C_{1}$ , $C_{2}$ and $C_{3}$ were determined by fitting the resolved magnetic field spectrum between the degree $2<l\leq 13$ .

Refer to caption — Figure 1: Energy spectrum of the GRIMM 2 magnetic field at the CMB for the epoch $2005.0$ (between SH degree $1\leq l\leq 13$ ) and its filtered (thick lines) and non filtered (thin lines) extrapolations. The laws used to extrapolate the large scale spectrum were taken from Roberts et al. (2003) (plus symbols), Voorhies (2004) (circles) and Buffett and Christensen (2007) (triangles).

As it is confirmed in figure 1 the three extrapolated spectra are presenting distinct behaviors. Furthermore, the law proposed by Buffett and Christensen (2007) predicts a strong concentration of energy at small scales. If this latter formulation was to be the closest one from reality, the unresolved MF would have to be modeled at very high degree (up to $l\sim 500$ ), which is nowadays numerically impossible. We propose therefore to reduce the impact of the small scale magnetic field on the large scale secular variation by applying a filter on the magnetic field and by determining the dynamical equation governing the evolution of the resulting filtered field.

II.4 Filtering of the fields

As it is the case for the Earth’s dynamo, the degree of turbulence in astrophysical or geophysical systems is most of the time extremely high. This implies that the fields are populated by a great variety of scales. Capturing all these scales in observations or in numerical simulations is usually impossible. Therefore, focusing on the large scales, and their dynamics, is probably the best solution to study such systems. This approach has been widely followed in the context of hydrodynamic (HD) and magnetohydrodynamic (MHD) turbulence (see Sagaut (2006); Lesieur (2008); Baerenzung et al. (2008a, b); Fabre and Balarac (2011)). To extract the large scales from a given field, one has to filter it. In HD or MHD turbulence three different filters are generally used, the top-hat or the Gaussian filter for studies in Cartesian space, or the sharp cut-off filter in spectral space. While the top-hat or Gaussian filter are continuous filters, the cut-off filter truncate the field above a certain predefined scale. In spherical coordinates the top-hat filter is usually preferred (Sun and Su (2007)) since until recently no equivalent of the Gaussian filter existed for such geometries (see however Bülow (2004) ). In our study, we decided to consider the isotropic filter developed by Bülow (2004). The principle of derivation of this filter is detailed in appendix A. Applying this filter to a scalar field ${\xi}$ leads to:

\overline{\xi}=G\star\xi\ ,

(14)

where $\overline{\xi}$ is the filtered field, the $\star$ symbol denotes the convolution product over the sphere between convolution kernel $G$ and the scalar field $\xi$ . According to the convolution theorem (see Driscoll and Healy (1994)), in spectral space equation (14) becomes:

\overline{\xi}_{l,m}=\sqrt{\frac{4\pi}{2l+1}}\xi_{l,m}G_{l}^{0}\ ,

(15)

where $\overline{\xi}_{l,m}$ and $\xi_{l,m}$ are respectively the filtered and the total SH coefficients associated with $\xi$ , and $G_{l}^{0}$ is convolution kernel expressed in spectral space at degree $l$ . Since for the filter developed by Bülow (2004), $G$ is zonal, only its coefficients at order $m=0$ are different from zero. As shown in appendix A equation (101), equation (15) can be expressed as:

\overline{\xi}_{l,m}=\xi_{l,m}\exp{\left(-\frac{l(l+1)\overline{\Delta}^{2}}{24R^{2}}\right)}\ .

(16)

where $\overline{\Delta}$ corresponds to the filter width. With this formulation one can directly connect the spectrum of the scalar field $E_{\xi}$ to the spectrum of its filtered part $E_{\overline{\xi}}$ as:

E_{\overline{\xi}}=E_{\xi}\exp{\left(-\frac{l(l+1)\overline{\Delta}^{2}}{12R^{2}}\right)}\ .

(17)

We then applied this filter to the three different extrapolated spectra presented in the previous section (see equations (11)-(13)). In figure 1 the non filtered (thin lines), and the filtered (thick lines) spectra are plotted. The width of the filter has been set to $\overline{\Delta}=500$ km in order to preserve the large scale fields and to suppress the small scale ones. One can observe that most of the energy that was contained in the small scales of the initial fields has now vanished in the filtered fields. Furthermore, the variations between the different extrapolated filtered spectra are now much less pronounced than in the non filtered case.

II.5 The Filtered Frozen-Flux (FFF) approximation

Applying the spatial filter presented in the preceding section to the velocity and the magnetic field allows to decompose them into an averaged part $\overline{{\bf u}}$ and $\overline{B}_{r}$ , and a fluctuating part ${\bf u}^{\prime}$ and $B_{r}^{\prime}$ such as:

	$\displaystyle{\bf u}$	$\displaystyle=$	$\displaystyle\overline{{\bf u}}+{\bf u}^{\prime}$		(18)
	$\displaystyle B_{r}$	$\displaystyle=$	$\displaystyle\overline{B}_{r}+B_{r}^{\prime}\ .$		(19)

Since the spatial filter considered here is homogeneous and time invariant, it commutes with all the differential operators encountered in the induction equation. Therefore applying this filter to the Frozen Flux approximation (equation (3)) leads to:

\partial_{t}\overline{B}_{r}=-{\bf\nabla_{{}_{H}}}(\overline{{\bf u}{B}}_{r})\ .

(20)

Following the decomposition introduced by Leonard (1974), the part of the subgrid stress tensor which allows interactions between the tangential components of the velocity field and the radial component of the MF reads:

{\bf\tau}=\overline{{\bf u}{B}}_{r}-\overline{{\bf u}}\overline{B}_{r}\ .

(21)

Now equation (20) can be rewritten as:

\partial_{t}\overline{B}_{r}=-{\bf\nabla_{{}_{H}}}(\overline{{\bf u}}\overline{B}_{r})-{\bf\nabla_{{}_{H}}}{\bf\tau}\ .

(22)

The subgrid stress tensor $\tau$ , through the averaged product $\overline{{\bf u}{B}}_{r}=\overline{(\overline{\bf u}+u^{\prime})(\overline{B}_{r}+B_{r}^{\prime})}$ , incorporates interactions between averaged and fluctuating part of the fields, thus equation (22) is not closed in the sense that it does not only contains large scale quantities. To close this equation, expressions of the fluctuating quantities depending on the averaged ones have to be derived.

As mentioned in appendix A, the filtered velocity and magnetic field are both solutions of the diffusion equations:

$\displaystyle\partial_{s}{\bf u}({\bf x},s)$	$\displaystyle=$	$\displaystyle D_{u}\Delta_{{}_{H}}{\bf u}({\bf x},s)$	(23)
$\displaystyle\partial_{s}B_{r}({\bf x},s)$	$\displaystyle=$	$\displaystyle D_{{}_{B}}\Delta_{{}_{H}}B_{r}({\bf x},s)$	(24)
$\displaystyle\partial_{s}\left({\bf u}({\bf x},s)B_{r}({\bf x},s)\right)$	$\displaystyle=$	$\displaystyle D_{{}_{B}}\Delta_{{}_{H}}\left({\bf u}({\bf x},s)B_{r}({\bf x},s)\right)$	(25)

where the filtered and total fields are respectively $\overline{\bf u}={\bf u}({\bf x},s)$ and $\overline{B}_{r}=B_{r}({\bf x},s)$ , and ${\bf u}={\bf u}({\bf x},0)$ and $B_{r}=B_{r}({\bf x},0)$ , and the product between $s$ and the two diffusion coefficients $D_{B}$ and $D_{u}$ determines the width of the filter. Note that since $D_{B}$ is the diffusion coefficient applied to each term of the Frozen Flux equation it has to be applied to the MF but not necessarily to the VF. To link fluctuating to filtered fields, Taylor expansions at the first order in $s$ of ${\bf u}$ , $B_{r}$ , and $\overline{{\bf u}{B}}_{r}$ , using relations (23) to (25) are performed, leading to:

$\displaystyle{\bf u}$	$\displaystyle\sim$	$\displaystyle\overline{\bf u}-sD_{u}\Delta_{{}_{H}}\overline{\bf u}$	(26)
$\displaystyle B_{r}$	$\displaystyle\sim$	$\displaystyle\overline{B}_{r}-sD_{{}_{B}}\Delta_{{}_{H}}\overline{B}_{r}$	(27)
$\displaystyle\overline{{\bf u}B}_{r}$	$\displaystyle\sim$	$\displaystyle{\bf u}B_{r}+sD_{{}_{B}}\Delta_{{}_{H}}({\bf u}B_{r}).$	(28)

By analogy to the usual Gaussian filters used in Large Eddy Simulations, one can define the following characteristic lengths:

	$\displaystyle{\overline{\Delta}}^{2}_{u}$	$\displaystyle=$	$\displaystyle 24sD_{u}$		(29)
	$\displaystyle{\overline{\Delta}}^{2}_{{}_{B}}$	$\displaystyle=$	$\displaystyle 24sD_{{}_{B}}$		(30)

Letting ${\overline{\Delta}}^{2}_{u}={\overline{\Delta}}^{2}_{{}_{B}}={\overline{\Delta}}^{2}$ , injecting equations (26) and (27) into equation (28), and keeping only the first order terms, one gets:

\overline{{\bf u}B}_{r}\sim\overline{{\bf u}}\overline{B}_{r}-\frac{{{\overline{\Delta}}^{2}}}{24}\left(\overline{{\bf u}}\Delta_{{}_{H}}\overline{B}_{r}+\overline{B}_{r}\Delta_{{}_{H}}\overline{{\bf u}}-\Delta_{{}_{H}}(\overline{{\bf u}}\overline{B}_{r})\right)\ ,

(31)

which can be reduced to:

\overline{{\bf u}B}_{r}\sim\overline{{\bf u}}\overline{B}_{r}+\frac{{{\overline{\Delta}}^{2}}}{12}\left(\left({\bf\nabla_{{}_{H}}}\overline{B}_{r}\right){\bf\nabla_{{}_{H}}}\right)\overline{\bf u}\ .

(32)

So the subgrid stress tensor $\tau$ in its approximated closed form reads:

{\bf\tau}=\frac{{{\overline{\Delta}}^{2}}}{12}\left(\left({\bf\nabla_{{}_{H}}}\overline{B}_{r}\right){\bf\nabla_{{}_{H}}}\right)\overline{\bf u}\ .

(33)

The Filtered Frozen Flux equation can now be written for the filtered secular variation $\partial_{t}\overline{B}_{r}^{<}$ truncated at degree $l_{c}=13$ as:

\partial_{t}\overline{B}_{r}^{<}=-\left({\bf\nabla_{{}_{H}}}(\overline{{\bf u}}\overline{B}_{r})\right)^{<}-\left({\bf\nabla_{{}_{H}}}{\bf\tau}\right)^{<}\ ,

(34)

with $\overline{B}_{r}=\overline{B}_{r}^{<}+\overline{B}_{r}^{>}$ .

In our study, the extrapolated filtered magnetic field $\overline{B}_{r}^{>}$ is extended up to the spherical harmonic degree $l=30$ . Above this scale, the filtered magnetic field (for $\overline{\Delta}=500$ km) is extremely weak (see figure 1) and its influence on the large scale secular variation is neglected.

II.6 Bayesian formulation of the inverse problem

The problem of determining the velocity field at the CMB knowing the exact MF and the SV together with its uncertainties is an ill-posed inverse problem: in a discrete approximation, the number of unknown is twice as large as the number of equations. One of the most common methods to tackle this problem consists in minimizing an energy functional composed of two main terms; a quantity measuring the discrepancies between the model and the data, balanced, through a regularization parameter, with a quantity expressing a prior knowledge on the expected solution. This method has been widely used to evaluate the flow at the CMB and for a review of the different parametrization employed see Holme (2007).

Another option consists in formulating the problem in a Bayesian framework. The solution becomes then the full posterior distribution of the velocity field given the secular variation. This method allows the estimation of a model for the flow together with the quantification of its uncertainties.

When variations around the prescribed MF are allowed to occur, the inverse problem becomes more complicated. For models of the Earth’s core MF derived from satellite or observatory data, the nature of these variations is diverse. At large scales they can be due to a leakage of the external and lithospheric field into the core field, whereas at small scales, the entire MF is undetermined because of the dominance of the lithospheric field at the Earth’s surface. Recent models such as GRIMM 2 (Lesur et al. (2010)) are able to separate at large scales (between SH degree $1\leq l\leq 13$ ) the external from the core field, but not the lithospheric field from the core field. As a consequence, the large scale lithospheric field becomes a source of uncertainty on the geodynamo’s MF.

To parametrize the unresolved part of the MF in the inverse problem different approaches have been recently developed. This includes the ensemble method of Gillet et al. (2009) and the iterative algorithm of Pais and Jault (2008). In this study we propose to extend these methods to the context of Bayesian modeling, following the development of Jackson (1995). Furthermore, in addition to parametrizing the unresolved MF, we also consider the uncertainties on the large scale MF due to the lithospheric field.

II.6.1 CMB velocity distribution

From now on, in order to simplify the notations, the filtered MF and VF as well as the filtered and truncated SV will be written as:

$\displaystyle\overline{B}_{r}$	$\displaystyle=$	$\displaystyle b$	(35)
$\displaystyle\overline{\bf u}$	$\displaystyle=$	$\displaystyle u$	(36)
$\displaystyle\dot{\overline{B}}_{r}^{<}$	$\displaystyle=$	$\displaystyle\gamma$	(37)

In this section, the distribution we want to characterize is the posterior distribution of the velocity field given the secular variation $p\left(u|\gamma\right)$ . But since we want to account for the unknown small scale magnetic field together with the uncertainties on the large scale field, this distribution cannot be expressed directly. Nevertheless, it can be obtained by marginalizing the joint posterior distribution of the velocity field and the magnetic field as following:

p\left(u|\gamma\right)=\int{p(u,b|\gamma)\mathrm{d}b}\ .

(38)

According to Bayes (1763) the distribution on the right hand side of relation (38) can be decomposed into:

p(u,b|\gamma)=\frac{p(\gamma|u,b)p\left(u,b\right)}{p(\gamma)}\ ,

(39)

with $p\left(u,b\right)$ the joint prior distribution of the VF and the MF, $p(\gamma)$ the distribution of the SV, which is constant with respect to both $u$ and $b$ , and finally $p(\gamma|u,b)$ the likelihood distribution. Because $u$ and $b$ are a priori assumed to be independent random variables their joint distribution can be split into two distributions such as:

p\left(u,b\right)=p\left(u\right)p\left(b\right)\ .

(40)

The prior distribution of the VF $p(u)$ is assumed to be Gaussian with the following general form:

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

(41)

with $d$ the dimension of the VF vector, which in our case is twice as large as the dimension of the MF and SV, and $\Sigma_{u}$ the velocity covariance matrix chosen to enforce the spatial smoothness of the flow. For this purpose we impose that:

	$\displaystyle u^{T}\Sigma_{u}^{-1}u$	$\displaystyle=$	$\displaystyle\int_{\Omega}\left(\|{\bf\nabla_{{}_{H}}}\left({\bf\nabla_{{}_{H}}}\cdot u\right)\|^{2}+\|{\bf\nabla_{{}_{H}}}\left({\bf r}\times{\bf\nabla_{{}_{H}}}\cdot u\right)\|^{2}\right)\mathrm{d}\omega$		(42)
		$\displaystyle+$	$\displaystyle\alpha^{2}\int_{\Omega}\|u\|^{2}\mathrm{d}\omega\ .$		(42)

The right hand-side of equation (42) is composed of two different norms on $u$ , namely the Bloxham’s ”strong norm” (Bloxham (1988); Jackson et al. (1993)) for the first one and the standard $L^{2}$ -norm for the second one. The domain of integration $\Omega$ is the surface of the sphere describing the CMB, and the balance factor $\alpha^{2}$ is chosen to be small enough not to modify the correlation length induced by the Bloxham norm. Furthermore, the covariance is rescaled such as the averaged standard deviation of the velocity field intensity is $20$ km/yr. This implies that at any location of the core-mantle boundary, the probability for the flow intensity to excess the value of $50$ km/yr, an upper limit calculated by Finlay and Amit (2011), is of the order of $0.01$ .

As for the velocity field, the magnetic field $b$ is also assumed to be normally distributed, but with a mean $b_{0}$ corresponding to the resolved MF, and a covariance $\Sigma_{b}$ . Its prior distribution can therefore be expressed as:

p(b)=\frac{\exp\left[-\frac{1}{2}\left(b-b_{0}\right)^{T}\Sigma_{b}^{-1}\left(b-b_{0}\right)\right]}{(2\pi)^{\frac{d}{4}}|\Sigma_{b}|^{\frac{1}{2}}}\ .

(43)

In this equation the two parts of the MF $b=b^{<}+b^{>}$ have different behaviors. The truncated field $b^{<}$ is taken from the GRIMM 2 model with $E[b^{<}]=b_{0}$ whereas the unknown MF $b^{>}$ is characterized by the averaged value $E[b^{>}]=0$ . Since in the GRIMM 2 model, the core field and the litospheric field are overlapping at large scales, this latter field can be viewed as a source of uncertainties on $b^{<}$ . We propose therefore to use the theoretical spectrum of the litospheric field given by Thebault and Vervelidou (2013) to build the covariance matrix $\Sigma_{b^{<}}$ for the resolved scales magnetic field. This spectrum reads:

E_{B}^{L}(l)=\left(l+1\right)\left(\mu_{0}|M|F_{l}^{a}\left(\epsilon\right)\right)^{2}l^{-\delta}C_{l}

(44)

with $|M|=0.4225A.m^{-1}$ the averaged crust magnetization, $\epsilon=27$ km the equivalent magnetized layer thickness, and the constants $\mu_{0}=4\pi 10^{-7}$ , $\delta=1.28$ , and $a=6371.2$ km. The two functions $F_{l}^{a}\left(\epsilon\right)$ and $C_{l}$ are given by:

	$\displaystyle C_{l}$	$\displaystyle=$	$\displaystyle\frac{l(l+1)(160l^{5}+264l^{4}-192l^{3}-130l^{2}+96l-9)}{6(2l+3)^{2}(2l+1)^{2}(2l-1)^{2}}$
	$\displaystyle F_{l}^{a}(\epsilon)$	$\displaystyle=$	$\displaystyle\frac{1-(1-\epsilon/a)^{(l-1)}}{l-1}\ .$

Letting the truncated MF $b^{<}$ being linked to its spectral counterpart $g_{l,m}^{<}$ through the relation $b^{<}=FPg_{l,m}^{<}$ , where the operator $P$ projects the coefficients in physical space, and the operator $F$ filters the field, and assuming that the litospheric field is isotropically distributed, the covariance of the truncated MF is:

$\displaystyle\Sigma_{b^{<}}\!$	$\displaystyle\!\!=\!\!$	$\displaystyle\!\!E\left[\left(b^{<}-b_{0}\right)\left(b^{<}-b_{0}\right)^{T}\right]$	(45)
	$\displaystyle\!\!=\!\!$	$\displaystyle\!\!E\left[PF\left(g_{l,m}^{<}-g_{l,m}^{0}\right)\left(g_{l,m}^{<}-g_{l,m}^{0}\right)^{T}F^{T}P^{T}\right]$	(46)
	$\displaystyle\!\!=\!\!$	$\displaystyle\!\!PF\frac{E_{B}^{L}(l)}{(l+1)(2l+1)}F^{T}P^{T}\ ,$	(47)

where $g_{l,m}^{0}$ are the spherical harmonics coefficients associated with $b_{0}$ .

The extrapolated SH coefficients of the MF $g_{l,m}^{>}$ are assumed to individually have a $0$ mean and a variance depending on the extrapolated spectrum $E_{B^{>}}^{B}(l)$ from Buffett and Christensen (2007) as shown in equation (10). Therefore, in physical space, the MF also has a $0$ mean at any spatial location and a covariance given by:

\Sigma_{b^{>}}=PF\frac{E_{B^{>}}^{B}(l)}{(l+1)(2l+1)}F^{T}P^{T}\ .

(48)

By combining $\Sigma_{b^{<}}$ with $\Sigma_{b^{>}}$ , one gets the total covariance $\Sigma_{b}$ for the MF.

The last distribution to characterize is the likelihood distribution which measures the discrepancies between the model (the Filtered Frozen Flux equation) and the data (the secular variation). It reads:

p(\gamma|b,u)=\frac{\exp\left[-\frac{1}{2}\left(\gamma+A_{u}b\right)^{T}\Sigma_{\gamma}^{-1}\left(\gamma+A_{u}b\right)\right]}{(2\pi)^{\frac{d}{4}}|\Sigma_{\gamma}|^{\frac{1}{2}}}\ ,

(49)

where the operator $A_{u}$ , when applied to $b$ , allows to calculate the non linear term of the filtered Frozen Flux equation $\left({\bf\nabla_{{}_{H}}}\left(ub+\tau\right)\right)^{<}$ , and the covariance matrix $\Sigma_{\gamma}$ is the SV posterior covariance of the GRIMM 2 model.

All the distributions entering the velocity field posterior distribution being detailed, this latter can be evaluated. The integral given in equation (38) has already been calculated by Jackson (1995), so we only present the result which reads:

	$\displaystyle p(\gamma\|u)$	$\displaystyle=$	$\displaystyle\int p(\gamma\|b,u)p(b)\mathrm{d}b$		(50)
		$\displaystyle\sim$	$\displaystyle\frac{(2\pi)^{\frac{d}{4}}}{\|N\|^{\frac{1}{2}}}\exp\left[-\frac{1}{2}\left(c-r^{T}N^{-1}r\right)\right]$		(51)

with:

$\displaystyle N$	$\displaystyle=$	$\displaystyle A_{u}^{T}\Sigma_{\gamma}^{-1}A_{u}+\Sigma_{b}^{-1}$	(52)
$\displaystyle r$	$\displaystyle=$	$\displaystyle\Sigma_{b}^{-1}b_{0}-A_{u}^{T}\Sigma_{\gamma}^{-1}\gamma$	(53)
$\displaystyle c$	$\displaystyle=$	$\displaystyle\gamma^{T}\Sigma_{\gamma}^{-1}\gamma+b_{0}^{T}\Sigma_{b}^{-1}b_{0}$	(54)

By multiplying the density (51) with the prior distribution (41) one gets the posterior distribution:

	$\displaystyle p(u\|\gamma)$	$\displaystyle=$	$\displaystyle\frac{p(\gamma\|u)p(u)}{p(\gamma)}$
		$\displaystyle\sim$	$\displaystyle\frac{1}{\|N\|^{\frac{1}{2}}}\exp\left[-\frac{1}{2}\left(c-r^{T}N^{-1}r+u^{T}\Sigma_{u}^{-1}u\right)\right]$

Since this expression does not allow to easily apprehend the effects arising from the modelization of the MF variations on the posterior distribution of the velocity field, we decided to rewrite it into a more intuitive form. It reads:

$\displaystyle p(u\|\gamma)$	$\displaystyle=$	$\displaystyle\frac{(2\pi)^{-\frac{3d}{4}}}{\|\Sigma_{\tilde{\gamma}}\|^{\frac{1}{2}}}\exp\left[-\frac{1}{2}\left(\gamma+A_{u}b_{0}\right)^{T}\Sigma_{\tilde{\gamma}}^{-1}\left(\gamma+A_{u}b_{0}\right)\right]$	(56)
	$\displaystyle\times$	$\displaystyle\frac{1}{\|\Sigma_{u}\|^{\frac{1}{2}}}\exp\left[-\frac{1}{2}u^{T}\Sigma_{u}^{-1}u\right]\times\frac{1}{p(\gamma)}$	(56)
$\displaystyle\Sigma_{\tilde{\gamma}}$	$\displaystyle=$	$\displaystyle\Sigma_{\gamma}+A_{u}\Sigma_{b}A_{u}^{T}\ .$	(57)

This formulation is very similar to the posterior distribution of the velocity field in the case where the magnetic field is exactly know. Indeed this latter distribution can be obtained by simply replacing the covariance matrix $\Sigma_{\tilde{\gamma}}$ by $\Sigma_{\gamma}$ . One can therefore observe that accounting for the small scale magnetic field and the lithospheric field when formulating the inverse problem in a Bayesian framework leads to an increase of the secular variation uncertainties through the quantity $A_{u}\Sigma_{b}A_{u}^{T}$ . Because of the dependency of this latter term on the velocity field $u$ , the maximum of the posterior distribution cannot be analytically calculated as already mentioned by Jackson (1995). Nevertheless, it is numerically possible to extract the main statistical characteristics of this posterior distribution using a Markov Chain Monte Carlo method (see Mosegaard and Rygaard-Hjalsted (1999); Rygaard-Hjalsted et al. (2000)). The algorithm we chose to explore the posterior distribution and the results we obtained are presented in the next section.

III Numerical method and results

The SV, MF and VF are expressed in physical space, therefore the surface of the CMB has been discretized by recursively dividing an initial icosahedron (see figure 2). The grid construction and its properties, as well as the approximation of the differential operators are explicited in appendix B.

Given that the grid is composed of $N$ nodes, the vectors $b$ , ${\gamma}$ , and ${u}$ are given by:

$\displaystyle b^{T}$	$\displaystyle=$	$\displaystyle\left(b_{0}...b_{i}...b_{N-1}\right)\ ,$
$\displaystyle{\gamma}^{T}$	$\displaystyle=$	$\displaystyle\left(\gamma_{0}...\gamma_{i}...\gamma_{N-1}\right)\ ,$
$\displaystyle u^{T}$	$\displaystyle=$	$\displaystyle\left(u_{0}...,u_{i},...,u_{N-1}\right)\ .$

III.1 Evaluation of the Filtered Frozen flux model

To evaluate the FFF model, an inversion of this equation was performed using GRIMM 2 as well as artificially generated data. The large scale magnetic field was taken from the GRIMM 2 model at the epoch $2001.0$ , whereas the small scales were randomly generated for SH degree lying between $14$ and $160$ according to the exponential law (11). The MF was then filtered (with $\overline{\Delta}=80$ km) such as its smallest scales, which cannot be properly represented on the grid, exhibited a low energy level. This MF is referred as $b_{0}$ . To create a SV associated with this MF, we drew randomly a velocity field and used it to advect the MF with the FF equation (3). The velocity field was decomposed into a poloidal and a toroidal part such as:

{\bf{u}}={\bf{\nabla}}\Phi+{\bf{e_{r}}}\times{\bf\nabla}\psi\ .

(58)

In spectral space, the spherical harmonics coefficients for the poloidal and toroidal field respectively read $\Phi_{l,m}$ and $\psi_{l,m}$ . In order to promote interactions between small and large scales, $\Phi_{l,m}$ and $\psi_{l,m}$ were extended up to spherical harmonic degree $80$ with the following statistical properties:

E[\Phi_{l,m}]=E[\psi_{l,m}]=0\quad\forall l,m

(59)

$\displaystyle E[\Phi_{l,m}\Phi_{l^{\prime},m^{\prime}}]$	$\displaystyle=$	$\displaystyle Cl^{-14/3}\|m\|^{-11/3}\delta_{ll^{\prime}}\delta_{mm^{\prime}}$	(60)
$\displaystyle E[\psi_{l,m}\psi_{l^{\prime},m^{\prime}}]$	$\displaystyle=$	$\displaystyle Cl^{-14/3}\|m\|^{-11/3}\delta_{ll^{\prime}}\delta_{mm^{\prime}}$	(61)
$\displaystyle E[\Phi_{l,m}\psi_{l^{\prime},m^{\prime}}]$	$\displaystyle=$	$\displaystyle 0\quad\forall l,l^{\prime},m,m^{\prime}\ ,$	(62)

where $C$ is a normalization constant. The choice of a $l^{-14/3}$ power law imposed onto the poloidal and toroidal field correlation allowed to generate a flow with similar statistical properties than a two-dimensional turbulent flow (see Sukoriansky et al. (2002)).

To directly simulate the advection of the MF, a fourth order Runge-Kutta scheme had been implemented. The integration time was taken to be $0.05$ year, and the computation was performed on a grid refined $7$ times (with $163842$ nodes).

In figure 3, the spectrum of the resulting secular variation is plotted. As it can be observed, the energy of the SV is maximum at scales lying between $l=60$ and $l=100$ indicating that the MF is strongly affected by the velocity field at these scales.

To mimic the situation encountered when using models fitting satellite data where both the MF and SV cannot be entirely taken into account, the artificial fields were truncated at degree 40. The resulting MF and SV are respectively referred as $\gamma$ and $b_{0}^{<}$ . These fields were then used as an input for the inversion of the FFF equation. Two different filter widths were tested, $\overline{\Delta}=500$ km, and $\overline{\Delta}=0$ km. In the latter case the equation reduces then to the usual FF approximation. Since the statistical properties of the velocity field were exactly known, they were directly injected in the prior information $p({u})$ . No measurement errors on the secular variation and the magnetic field had been generated, therefore the likelihood and the MF prior distributions respectively read:

	$\displaystyle p\left(\gamma\|u,b\right)$	$\displaystyle=$	$\displaystyle\delta\left(\gamma+{\bf\nabla_{{}_{H}}}(bu+\tau)\right)\ ,$		(63)
	$\displaystyle p(b)$	$\displaystyle=$	$\displaystyle\delta(b-b_{0}^{<})\ .$		(64)

Multiplying together these two distributions and marginalizing the result with respect to $b$ leads to:

p\left(\gamma|u\right)=\delta\left(\gamma+{\bf\nabla_{{}_{H}}}\left(b_{0}^{<}u+\tau\left(b_{0}^{<}\right)\right)\right)\ ,

(65)

with $\tau(b_{0}^{<})$ the subgrid stress tensor evaluated with $b_{0}^{<}$ . The posterior distribution of the velocity field is then proportional to:

p\left(u\right|\gamma)\sim p\left(\gamma|u\right)p(u)\ .

(66)

The Bayesian formulation of the problem being described, the discrete velocity field that maximizes the posterior distribution could be determined for the two cases. To do this, a particular solution of the equation $\gamma+{\bf\nabla_{{}_{H}}}\left(b_{0}^{<}u+\tau\left(b_{0}^{<}\right)\right)=0$ was calculated. In addition, the null space of the operator $A_{b_{0}^{<}}$ defined such as $A_{b_{0}^{<}}u=\nabla_{{}_{H}}\left(b_{0}^{<}u+\tau\left(b_{0}^{<}\right)\right)$ was parametrized. The final solution corresponded then to the sum of the particular solution with the null-space one which minimized the prior information on the VF. The grid used to realize the computation was composed of $10242$ nodes (approximately $16$ times less than the one taken to advect the MF). For the results to be comparable, the three different velocity fields were truncated at SH degree $l=40$ . Furthermore the artificially generated field as well as the field obtained by inverting the FF approximation ( $\overline{\Delta}=0$ ), were filtered with $\overline{\Delta}=500$ km.

On the left side of figure 4 the poloidal and toroidal spectra respectively defined by:

	$\displaystyle E_{\Phi}$	$\displaystyle=$	$\displaystyle l(l+1)\sum_{m=-l}^{m=l}(\Phi_{l,m})^{2}$		(67)
	$\displaystyle E_{\psi}$	$\displaystyle=$	$\displaystyle l(l+1)\sum_{m=-l}^{m=l}(\psi_{l,m})^{2}\ ,$		(68)

are plotted for the three velocity fields. The behavior of the spectra associated with the artificial field (full line) is correctly reproduced at large scale ( $1\leq l\leq 25$ ) by both the FF (crosses) and the FFF (circles) flow models, whereas at SH degree close to the cutoff $l_{c}=40$ , the discrepancies between the exact spectra and the spectra of the inverted fields become larger. Nevertheless, the spectra of the difference between exact and inverted fields, displayed on the right side of figure 4, show that the use of the FFF equation allows to reproduce more accurately the artificial field at almost any scale.

To quantify the spatial error of the two inverted fields, the following quantities were computed:

$\displaystyle\mathcal{E}_{FF}$	$\displaystyle=$	$\displaystyle\frac{\int\|u-u_{{}_{FF}}\|^{2}\mathrm{d}\Omega}{\int\mathrm{d}\Omega}=4,92\;\text{km}^{2}.\text{yr}^{-2}$	(69)
$\displaystyle\mathcal{E}_{FFF}$	$\displaystyle=$	$\displaystyle\frac{\int\|u-u_{{}_{FFF}}\|^{2}\mathrm{d}\Omega}{\int\mathrm{d}\Omega}=2,83\;\text{km}^{2}.\text{yr}^{-2}$	(70)
$\displaystyle E_{tot}$	$\displaystyle=$	$\displaystyle\frac{\int\|u\|^{2}\mathrm{d}\Omega}{\int\mathrm{d}\Omega}=130\;\text{km}^{2}.\text{yr}^{-2}\ ,$	(71)

where the integration domain is the surface of the CMB, $u$ is the exact filtered velocity field, $u_{FF}$ and $u_{FFF}$ are respectively the filtered velocity field obtained by inverting the FF and the FFF equations. The result of these computations shows that although the energy associated to the error fields is weak in comparison to the total energy of the exact flow, performing an inversion of FFF approximation reduces the global error on the VF.

III.2 Sampling of the velocity posterior distribution

In this part we present a method to sample the posterior distribution $p(u|\gamma)$ given in equation (56). Since this distribution exhibits a complex form with respect to the velocity field $u$ , directly drawing sample from it is impossible. Nevertheless, by building an appropriate Markov Chain on the VF one can map the distribution (for a complete description of Markov Chain Monte Carlo methods see Gamerman and Lopes (2006)). For this study we chose an algorithm of the Metropolis-Hastings type to construct the chain. The principle of the method is the following:

(a)

in the initial step, a VF $u^{n}$ with $n=0$ is generated. No particular property have to be imposed on this field, but choosing a field which is as close as possible to the one maximizing the target distribution will allow the chain to converge faster.
(b)

from $u^{n}$ a field $u^{n+1}$ is constructed according to some arbitrary transition kernel $q(u^{n+1},u^{n})$ .

(c)

the next step consists in accepting or rejecting the move from $u^{n}$ to $u^{n+1}$ . Therefore, an acceptance probability $\alpha(u^{n+1},u^{n})$ is defined. In the case of the Metropolis-Hastings algorithm this probability is expressed as:

\alpha(u^{n+1},u^{n})=\min\left\{1,\frac{p(u^{n+1}|\gamma)q(u^{n+1},u^{n})}{p(u^{n}|\gamma)q(u^{n},u^{n+1})}\right\}\ .

(72)

(d)

the process returns then to step (b) with $u^{n+1}$ if the move from $u^{n}$ to $u^{n+1}$ is accepted, and with $u^{n}$ otherwise.

The ensemble is then assumed to be representative of the posterior distribution, once its averaged field as converged towards a fixed vector.

Note that this algorithm has already been employed by Mosegaard and Rygaard-Hjalsted (1999); Rygaard-Hjalsted et al. (2000) for sampling the posterior distribution given in equation (II.6.1), but with a different parametrization of the magnetic field and the velocity field. In particular, only the uncertainties on the large scale MF were considered in these two studies.

As already mentioned previously, analytically calculating the maximum of the posterior distribution $p(u|\gamma)$ is not feasible. We therefore decided to approximate it by taking the averaged VF of the ensemble generated by the Markov Chain, such as:

\operatorname*{arg\,max}_{u}p(u|\gamma)\sim\int u\;p(u|\gamma)\mathrm{d}u\ .

(73)

III.2.1 Evaluation and comparison of the method with artificial data

To evaluate our method, we created an artificial set of data, and performed the inversion of the FFF equation using these data. We also compared our results to the ones obtained with alternative approaches. The construction of the synthetic fields was performed as following:

Artificial magnetic field

The large scale MF (between spherical harmonic degree $1\leq l\leq 13$ ) was taken from the GRIMM 2 model at the epoch $2004.0$ . Its spectrum (the thick black line in the top of figure 5) was then extrapolated up to degree $l=30$ according to the formulation (11) of Buffett and Christensen (2007). From this small scale spectrum, and under the assumption of isotropy and $0$ mean of the field, a MF was randomly generated and added to the GRIMM 2 large scale field. In figure 5 (top) the spectrum of the small scale MF is represented by the plus symbols.
Artificial velocity field

The coefficients, in spectral space, of the poloidal ( $\Phi$ ) and toroidal ( $\psi$ ) fields were assumed to be isotropically distributed with a $0$ mean and a covariance $\tilde{\Sigma}_{u}$ derived from the following power law spectrum:

$E_{\Phi}(l)=E_{\psi}(l)=A^{2}l^{-5/3}\ ,$ (74)

where the value of the amplitude $A$ was chosen such as the averaged velocity intensity at the CMB was of $17$ km.yr^-1. A velocity field extending up to degree $l=26$ was then randomly drawn accordingly to these statistical properties.
Artificial secular variation

To generate an artificial large scale SV (extending up to degree $l=13$ ), the MF was advected by the VF through the non-linear term of the Frozen Flux approximation (3). The energy spectrum of the resulting SV is presented with the thick line in figure 5 (bottom).

The next step of this evaluation was to recover the velocity field according to the magnetic field and the secular variation. But for this inverse problem to be more realistic, uncertainties were added to the data. A large scale lithospheric field ( $1\leq l\leq 13$ ) was randomly generated accordingly to the theoretical spectrum $E_{B}^{L}$ of Thebault and Vervelidou (2013) (see equation 44) and added to the artificial MF. This contaminated field was then truncated at degree $l=13$ . It is referred as $b_{0}$ and its energy spectrum is plotted on top of figure 5 with circles. The uncertainties on the secular variation were built by randomly drawing a field from a Gaussian distribution with a $0$ mean and a covariance given by the posterior covariance matrix of the GRIMM 2 secular variation $\Sigma_{\gamma}$ . The resulting field was then superimposed on the artificial field. The spectrum associated with the total SV is shown with circles on the bottom of figure 5.

Except for the prior distribution of the velocity field $p(u)$ , all the distributions derived in section II.6.1 are consistent to describe the posterior distribution of the velocity field in this test. We recall that they read:

	$\displaystyle p(\gamma\|u,b)$	$\displaystyle=$	$\displaystyle\mathcal{N}\left(-\left({\bf\nabla_{{}_{H}}}\left(ub+\tau\right)\right)^{<},\Sigma_{\gamma}\right)\ ,$		(75)
	$\displaystyle p(b)$	$\displaystyle=$	$\displaystyle\mathcal{N}(b_{0},\Sigma_{b})\ ,$		(76)

where $\mathcal{N}(x_{0},\Sigma_{x})$ corresponds to the normal distribution centered in $x$ and with covariance $\Sigma_{x}$ , and ${\Sigma}_{b}$ is the covariance of the MF in which are included the uncertainties due to the lithospheric field and the modelization of the unresolved MF. For this evaluation phase, the prior distribution of the VF was given by:

p(u)=\mathcal{N}(0,\tilde{\Sigma}_{u})\ ,

(77)

with $\tilde{\Sigma}_{u}$ the covariance of the VF derived from the power laws (74).

The estimation of the flow with respect to the artificial MF and SV was then realized with the four algorithms presented below.

MCMC method

In this approach the full posterior distribution of the velocity field $p(u|\gamma)$ is explored with the Metropolis-Hastings algorithm described in the beginning of the section. The transition kernel $q(u^{n+1},u^{n})$ entering the algorithm (step b) is derived from the prior distribution of the VF as following:

q(u^{n+1},u^{n})=\frac{exp\left[{-\frac{1}{2\lambda^{2}}(u^{n+1}-u^{n})^{T}\tilde{\Sigma}_{u}^{-1}(u^{n+1}-u^{n})}\right]}{(2\pi)^{\frac{d}{2}}|\tilde{\Sigma}_{u}|^{\frac{1}{2}}}\ ,

(78)

where the factor $\lambda$ allows to rescale the covariance $\tilde{\Sigma}_{u}$ in order to limit the distance of the move from one field to the other. Since this kernel is symmetric with respect to $u^{n+1}$ and $u^{n}$ , the acceptance probability $\alpha(u^{n+1},u^{n})$ of equation (72) reduces to:

\alpha(u^{n+1},u^{n})=\min\left\{1,\frac{p(u^{n+1}|\gamma)}{p(u^{n}|\gamma)}\right\}\ .

(79)

The initial field was set to $0$ and an ensemble of $N_{m}=230000$ VF $u_{i}$ was generated. We then approximated the flow maximizing the posterior distribution through the averaging operation:

u_{0}=\frac{1}{N_{m}}\sum_{i=1}^{i=N_{m}}u_{i}\ .

(80)

Least square method

In this method, the magnetic field $b$ is assumed to be exactly known such as:

p(b)=\delta(b-b_{0})\ .

(81)

As a consequence, the posterior distribution of the velocity field becomes proportional to:

	$\displaystyle p(u\|\gamma)$	$\displaystyle\sim$	$\displaystyle\exp\left[-\frac{1}{2}\left(\gamma+A_{b_{0}}u\right)^{T}\Sigma_{{\gamma}}^{-1}\left(\gamma+A_{b_{0}}u\right)\right]$		(82)
		$\displaystyle\times$	$\displaystyle\exp\left[-\frac{1}{2}u^{T}\tilde{\Sigma}_{u}^{-1}u\right]\ ,$		(82)

where $A_{b_{0}}$ is the operator allowing to evaluate the non linear term of the FFF equation when it is applied to $u$ . The maximum of this distribution can be calculated analytically and corresponds to the usual least square solution:

u_{0}=-(A_{b_{0}}^{T}\Sigma_{{\gamma}}^{-1}A_{b_{0}}+\tilde{\Sigma}_{u}^{-1})^{-1}A_{b_{0}}^{T}\Sigma_{{\gamma}}^{-1}\gamma\ .

(83)

Ensemble method

This approach was developed by Gillet et al. (2009) and consists in generating an ensemble of magnetic field and calculating their associated velocity field. In this test, the ensemble of MF was randomly drawn from the distribution $p(b)$ given in equation (76). In total we generated $N_{e}=100$ magnetic fields ${b_{i}}$ (with $1..i..N_{e}$ ) extending up to degree $l=30$ . Each velocity field $u_{i}$ were then determined by the relation:

u_{i}=-(A_{b_{i}}^{T}\Sigma_{{\gamma}}^{-1}A_{b_{i}}+\tilde{\Sigma}_{u}^{-1})^{-1}A_{b_{i}}^{T}\Sigma_{{\gamma}}^{-1}\gamma\ .

(84)

where the operator $A_{b}$ now depends on each realization $b_{i}$ of the magnetic field. The flow solution of the complete inverse problem is given by:

u_{0}=\frac{1}{N_{e}}\sum_{i=1}^{i=N_{e}}u_{i}\ .

(85)

Iterative method

In this approach, Lesur et al. (2013) proposed to determine the CMB velocity field through the following iterative process:

u^{n+1}=-(A_{b_{0}}^{T}\left(\Sigma_{\tilde{\gamma}}^{n}\right)^{-1}A_{b_{0}}+\tilde{\Sigma}_{u}^{-1})^{-1}A_{b_{0}}^{T}\left(\Sigma_{\tilde{\gamma}}^{n}\right)^{-1}\gamma

(86)

where $n$ is the index of the iteration, and the covariance $\Sigma_{\tilde{\gamma}}^{n}$ is given by:

\Sigma_{\tilde{\gamma}}^{n}=\Sigma_{\gamma}+\left(A_{u^{n}}\right)\Sigma_{b}\left(A_{u^{n}}\right)^{T}

(87)

with $A_{u^{n}}$ the operator depending on the velocity field $u^{n}$ , and which allows to evaluate the non linear term of the FFF equation when it is applied to $b$ . In the numerical implementation of this method, the initial field $u^{0}$ was taken as the one solution of the least square approach.

Note that this algorithm provides an estimation of the maximum of the posterior distribution given in equation (56) if the quantity $A_{u}\Sigma_{b}A_{u}$ is assumed to vary slowly with respect to $u$ .

In the previous section we mentioned the algorithm of Pais and Jault (2008) which also allows to take into account the variations of the MF in the inverse problem. We decided not to implement it in this evaluation since it is an approximation of the method proposed by Lesur et al. (2013).

To compare the different approaches, the artificial velocity field $u$ is decomposed into a toroidal ( $\psi$ ) and a poloidal ( $\Phi$ ) field, with $\psi_{l,m}$ and $\Phi_{l,m}$ their respective spectral counterpart. The same operation is performed on the velocity fields $u_{0}$ given by the four inversion methods. Their toroidal and poloidal parts are then referred as $\psi_{0}$ and $\Phi_{0}$ in physical space, and $\psi_{l,m}^{0}$ and $\Phi_{l,m}^{0}$ in spectral space.

To measure the accuracy, scale by scale, of the various fields, we define the poloidal and toroidal error fields as:

	$\displaystyle\epsilon_{l,m}^{\Phi}$	$\displaystyle=$	$\displaystyle\Phi_{l,m}-\Phi_{l,m}^{0}$		(88)
	$\displaystyle\epsilon_{l,m}^{\psi}$	$\displaystyle=$	$\displaystyle\psi_{l,m}-\psi_{l,m}^{0}\ ,$		(89)

and their associated energy spectra:

	$\displaystyle E_{\epsilon_{\Phi}}$	$\displaystyle=$	$\displaystyle l(l+1)\sum_{m=-l}^{m=l}(\epsilon_{l,m}^{\Phi})^{2}$		(90)
	$\displaystyle E_{\epsilon_{\psi}}$	$\displaystyle=$	$\displaystyle l(l+1)\sum_{m=-l}^{m=l}(\epsilon_{l,m}^{\psi})^{2}\ .$		(91)

In figure 6 the spectra of the error fields (symbols) are plotted together with the poloidal and toroidal spectra of the artificial velocity field (solid lines). The first observation one can make is that above the SH degree $l=10$ , the energy associated with the poloidal and toroidal error is of the same order than the energy of the artificial fields. This implies that above this degree, the estimation of the flow is not reliable whatever the method employed to determine it. At large scale, the performance of the different approaches varies strongly. Whereas the flow obtained with the least square method (triangles) is the one which deviates the most from the artificial flow, the velocity fields evaluated with the iterative method (circles) and the MCMC algorithm (plus symbols) are the ones presenting the lowest error intensities. In between, in terms of accuracy, lies the ensemble method (crosses). It can be noted that at degree $l=5$ and between $l=5$ and $l=6$ for respectively the ensemble and the least square method, the energy associated with the toroidal error field is larger than the energy of the artificial field itself. This is not the case anymore when using the iterative or MCMC approaches.

To evaluate the different models in physical space, the coefficients associated with the poloidal and toroidal error fields (equations 88-89) were truncated at degree $l=10$ , since at smallest scales the uncertainties are maximal, and projected in real space. Figure 7 shows the intensity of these fields at the level of America for the four approaches. Although the locations of the errors are similar, their intensities differs strongly from one flow to the other. A computation of the averaged energy associated with the poloidal and toroidal error fields (see table LABEL:tableError) shows that globally the best approximation of the artificial velocity field is provided by the MCMC algorithm, followed in order by the iterative, ensemble and least square methods. Note that the differences between the MCMC and iterative approach are very low, and since the computation time required to sample the full posterior distribution is much larger than the one to approximate its maximum using the iterative method, if one wants to determine the flow without its underlying uncertainties, using the algorithm of Lesur et al. (2013) is certainly more appropriate.

Table 1: Averaged energy of the error in km

{}^{2}.

yr^-2 associated with the different inverted velocity fields. These values have to be compared with the averaged energies of the artificial poloidal and toroidal fields which are respectively of

48.07

{}^{2}.

yr^-2 and

45.28

{}^{2}.

yr^-2.

	MCMC	Iterative	Ensemble	Least square
Poloidal Field	$2.20$	$2.42$	$3.04$	$3.72$
Toroidal Field	$6.45$	$6.72$	$8.69$	$10.1$

Since the MCMC algorithm provides an information on the flow uncertainties, we extracted the standard deviation of the velocity intensity $\sigma{|u|}$ , from the variance $\sigma_{u}^{2}$ as following:

\sigma_{u}^{2}=\int diag\left[\left(u-u_{0}\right)\left(u-u_{0}\right)^{T}\right]p(u|\gamma)\mathrm{d}u

(92)

where the term $diag$ means that only the diagonal elements of the matrix lying within the brackets are kept. At each node of the discrete CMB, the VF $u$ is composed of a polar and an azimuthal component, so to get the variance associated with the velocity intensity $\sigma_{|u|}^{2}$ , the variance of each component has to be summed up. In figure 8 the quantity $2\sigma_{|u|}$ , corresponding to the $95\%$ confidence interval on the flow intensity, is displayed at the level of $America$ . This picture presents a pessimistic view of the uncertainties to be expected when evaluating the flow maximizing the posterior distribution. Because the probability for the real flow to lie within the tails of the posterior distribution is very low, the predicted error will globally be larger than the effective one. Nevertheless, locations where the differences between the solution flow and the real one are important corresponds to areas of high posterior variance.

Finally, we wanted to measure the impact of the velocity field prior information on the solution of the inverse problem. We therefore realized a simulation with the iterative algorithm of Lesur et al. (2013), in which only the covariance matrix of the velocity field $\tilde{\Sigma}_{u}$ of equation (87) has been modified. Instead of being the covariance imposed to generate the artificial velocity field, this latter was replaced by the covariance ${\Sigma_{u}}$ defined in equation (42) and derived from the Bloxham’s ”strong norm”. We then calculated, as previously, the averaged energy of the error field using the velocity field we obtained. It is of $3.45$ km ${}^{2}.$ yr^-2 and $9.23$ km ${}^{2}.$ yr^-2 for respectively the poloidal and toroidal part of the error. This shows us that the choice of the prior information on the velocity field influences in a significant manner the results of the inverse problem. Nevertheless, it can be noticed that the intensity of the error in this case remains lower than the intensity of the error for the flow resulting from the least square approach, showing the importance of modeling the possible variations of the MF.

III.2.2 Application of the MCMC algorithm for the epoch 2005.0

In the simulation we realized, the MF and the SV as well as the covariance for the SV, are given up to SH degree $l=13$ by the GRIMM 2 model of Lesur et al. (2010) for the epoch $2005.0$ . For the magnetic field, the covariance associated with its large scales ( $0<l\leq 13$ ) is derived from the theoretical spectrum of the litospheric field of Thebault and Vervelidou (2013) as shown in equations (44)-(47), whereas the extrapolation of the MF spectrum proposed by Buffett and Christensen (2007) is used to build the covariance matrix of the small scale MF ( $13<l\leq 30$ ) as presented in equation (48). The velocity field is extended up to SH degree $l=26$ , and its prior distribution is given by the relation (41). The width of the filter has been set to $\overline{\Delta}=500$ km, and the initial field of the Markov chain $u_{0}$ , is the velocity field solution of the least square approach:

u_{0}=-\left(A_{B}^{T}\Sigma_{\gamma}^{-1}A_{B}+\Sigma_{u}^{-1}\right)^{-1}A_{B}^{T}\Sigma_{\gamma}^{-1}\gamma

(93)

with $A_{B}u=\left({\bf\nabla_{{}_{H}}}\left(ub+\tau\right)\right)^{<}$ .

As for the synthetic test of section III.2.1, the transition kernel $q(u^{n+1},u^{n})$ was derived from the prior distribution of the VF as shown in equation (78).

The Metropolis-Hastings algorithm was then numerically simulated on a discrete CMB refined 4 times. An ensemble of $130000$ VF $u$ mapping the posterior distribution have been generated. The velocity field $\hat{u}$ maximizing the posterior, is approximated by taking the average field of the ensemble we created.

In figure 9, the vector field $\hat{u}$ and its intensity are displayed in different locations of the core mantle boundary. Many features of the flow we obtained have already been reported in previous studies. In particular, the eccentric and planetary scale anticyclonic gyre observed by Pais and Jault (2008) and Gillet et al. (2009) is also present in the flow we calculated. This observation reinforces the hypothesis that the fluid motions in the outer core can be well-described by the compressible quasi-geostrophic assumption, a constraint a priori applied in both studies. Nevertheless, according to our results, deviations from quasi-geostrophy (QG) have also to be expected. Indeed, under this hypothesis, the flow is forced to be symmetric with respect to the equator outside the tangential cylinder, a condition which is not fulfilled everywhere in our case. Although the symmetric part of the velocity field is dominant in our simulation ( $82\%$ of the energy of the total VF is concentrated in its symmetric components) certain patterns, such as the flow crossing the equator below India or the larger intensity of the westward drift in the southern hemisphere are violating this property. The flow we obtained also exhibits a much smoother spatial behavior than the ones presented by Pais and Jault (2008); Gillet et al. (2009). This is certainly due to the choice we made to characterize a priori the velocity field. Through the Bloxham’s ”strong norm”, we imposed a very steep spectrum (in $l^{-5}$ ) to both the poloidal and toroidal field, as a consequence the intermediate scales of the VF were probably over-damped in our simulation.

Possessing the full posterior distribution of the VF allows to extract many useful information on the flow. In particular the uncertainties on the velocity field intensity and orientation (with respect of course to the prescribed modelization) can be evaluated at any spatial location on the grid. To compute the standard deviation of the velocity field intensity $\sigma_{|u|}$ we followed the protocol given in section III.2.1, whereas the standard deviation of the velocity orientation $\sigma_{\Gamma}$ , is derived from the formula:

\sigma_{\Gamma}^{2}=\int diag\left[\Gamma\Gamma^{T}\right]p(u|\gamma)\mathrm{d}u\ ,

(94)

where $\Gamma$ is the angle, in degree, between the velocity field $u$ and $\hat{u}$ . The quantities $\sigma_{|u|}$ and $\sigma_{\Gamma}$ , are displayed in figure 10 (on the left for the $\sigma_{|u|}$ and on the right for the $\sigma_{\Gamma}$ ). This figure first shows that a strong (resp. weak) uncertainty on the intensity coincides with a strong (resp. weak) uncertainty on the orientation of the flow. Secondly one can notice that the uncertainties are not homogeneously distributed on the surface of the CMB. While the planetary scale eccentric gyre seems to be very robust, the VF at the level of the Pacific ocean is much more uncertain. It is known that in this latter area the magnetic field activity is moderate (Hulot et al. (2002)). As a consequence the secular variation is low, and its associated uncertainties, due to the inaccuracy of the measurements and to the interactions between the unresolved MF and the VF, may be larger than the signal itself. It is therefore very difficult to evaluate the velocity field in this region. Nevertheless, this maps shows that there is a large part east of Australia and around the longitude of New-Zealand where the VF can be accurately estimated. Another particular feature can be noticed at the level of the North-Atlantic ocean, where the uncertainties on both the flow intensity and orientation are very large. As already questioned by Finlay et al. (2010), the robustness of the clockwise gyre usually observed in this area by models assuming Tangential-Geostrophy seems to be very weak according to our results.

Since we are using truncated MF and SV it may be interesting to investigate the spectral properties of the posterior VF. The poloidal $\Phi$ and toroidal $\psi$ part of the VF $u$ are expanded in spherical harmonics. The resulting fields are respectively referred as $\Phi_{l,m}$ and $\psi_{l,m}$ . The same operation is performed on $\hat{u}={\bf{\nabla}}\hat{\Phi}+{\bf{e_{r}}}\times{\bf\nabla}\hat{\psi}$ , with $\hat{\Phi}_{l,m}$ and $\hat{\psi}_{l,m}$ its poloidal and toroidal SH coefficients. In figure 11 are plotted the spectra associated with $\hat{\Phi}_{l,m}$ and $\hat{\psi}_{l,m}$ as well as the variance on these coefficient summed up over the order $m$ and rescaled by the factor $l(l+1)$ . We can observe that the largest scales of the flow are dominated by the toroidal field. A computation of the poloidal and toroidal energy shows that more than $88\%$ of the total energy is of toroidal nature. The other information which can be extracted from figure 11 is that, as for the synthetic test we realized previously, above the SH degree $l=10$ , the intensity of the flow uncertainties becomes larger than the intensity of the flow itself. As a consequence the evaluation of the small scale velocity field cannot be considered as reliable.

It has to be emphasized that all the results we obtained are conditioned by the choice of the prior information imposed to the flow, and should not be considered as absolute. In particular the globally low level of uncertainties on the flow solution is certainly a consequence of the strong regularization we employed.

IV Conclusions

In this study we have presented a new method to determine the velocity field at the Earth’s core mantle boundary according to an outer core magnetic field and secular variation model. We showed that using an appropriate dynamical equation to prescribe the large scale magnetic field evolution in the inverse problem, permitted to reduce the modeling errors arising from the truncated nature of the available fields. We also demonstrated that the Bayesian formalism we developed to account for the large scale uncertainties on the magnetic field and to model the unresolved small scale MF, allowed to properly describe the inverse problem as soon as the information introduced a priori were accurate. Through the evaluation of our method and the comparison with other approaches, we could indirectly confirm that the unresolved part of the magnetic field contributed significantly to the observed secular variation, and that its modelization was necessary to obtain a more accurate description of the flow at the core mantle boundary.

When we applied our method to real data, provided by the GRIMM 2 model for the epoch $2005.0$ , we could recover many features of the flow already observed in previous studies where a different prior information on the velocity field had been considered. In particular we could retrieve the planetary-scale eccentric gyre characteristic of flow evaluated under the compressible quasi-geostrophy assumption (see Pais and Jault (2008); Gillet et al. (2009)). Nevertheless, according to our simulation, the flow is crossing the equator below India and the intensity of the westward drift is the larger in the southern hemisphere, indicating that the equatorial symmetry imposed by the quasi-geostrophy hypothesis is broken. Through this observation one can conclude that deviations from quasi-geostrophy should be allowed to occur when this latter constraint is imposed in the inverse problem. Another specificity of the velocity field we obtained is its very smooth spatial behavior. This property is certainly induced by the prior distribution of the velocity field we chose, since this latter imposes a very steep spectrum to both the poloidal and toroidal part of the velocity field.

Finally, thanks to the ensemble of velocity field we generated to map the posterior distribution, we could evaluate the uncertainties associated with the flow solution of the inverse problem. According to our results, whereas on the one hand, the robustness of the flow is questioned in many area, and particularly in almost the entire Pacific ocean, and in the northern part of the Atlantic ocean, on the other hand, the planetary-scale eccentric gyre seems to be a very robust structure. From the evaluation of the uncertainties in spectral space, we concluded that the flow at the CMB could only be accurately estimated at large scales (between spherical harmonics degree $1$ and $10$ ).

Appendix A Spherical diffusion

In Cartesian space, applying a Gaussian filter to a scalar or a vector field, or letting the field evolve through a diffusion process can be interpreted as being similar operations. Indeed, the kernel of the Gaussian filter is:

G({\bf x}-{\bf x^{\prime}})=\left(\frac{\gamma}{\pi\overline{\Delta}^{2}}\right)^{\frac{d}{2}}\exp\left({-\frac{\gamma|{\bf x}-{\bf x^{\prime}}|^{2}}{\overline{\Delta}^{2}}}\right)\ ,

(95)

where $d$ corresponds to the spatial dimension, $\overline{\Delta}$ is the filter width, and $\gamma$ is a constant usually set to $6$ . The solution of the diffusion equation:

\partial_{t}\xi({\bf x},t)=D\Delta\xi({\bf x},t)\ ,

(96)

where $D$ is a diffusion coefficient and $t$ the time, is the convolution between the scalar field $\xi({\bf x^{\prime}},t=0)$ , and the Green function:

G({\bf x}-{\bf x^{\prime}},t)=\frac{1}{(4\pi Dt)^{\frac{d}{2}}}\exp\left(-\frac{|{\bf x}-{\bf x^{\prime}}|^{2}}{4Dt}\right)\ .

(97)

So if one sets $Dt=\frac{\overline{\Delta}^{2}}{4\gamma}$ , diffusing a scalar field through equation (98) is equivalent to filtering that field with the convolution kernel expressed in (95). Based on this observation, Bülow (2004) derived a Gaussian-like filter on the surface of a sphere of radius $R$ , by determining the convolution kernel of the spherical diffusion equation:

\partial_{t}\xi({\bf x},t)=D\Delta_{{}_{H}}\xi({\bf x},t)\ .

(98)

which reads:

G=\sum_{l\epsilon\mathbb{N}}\sqrt{\frac{2l+1}{4\pi}}Y_{l0}\exp\left(-\frac{l(l+1)Dt}{R^{2}}\right)\ ,

(99)

where $Y_{l0}$ is the spherical harmonic of degree $l$ and order $m=0$ . So in spectral space, the filtering operation of a scalar field $\xi$ expanded in spherical harmonics:

\xi=\sum_{l\epsilon\mathbb{N}}\sum_{m=-l}^{m=l}\xi_{l,m}Y_{lm}

(100)

simply reduces to the operation:

\overline{\xi}_{l,m}=\xi_{l,m}\exp\left(-\frac{l(l+1)\overline{\Delta}^{2}}{24R^{2}}\right)\ .

(101)

Appendix B Discretization of the Core Mantle Boundary

B.1 Construction of the grid

The grid describing the Core Mantle Boundary is obtained by recursively subdividing an initial icosahedron as explained in Baumgardner and Frederickson (1985); Stuhne and Peltier (1999). For each grid refinement procedure, a node is added in the middle of the geodesic arc linking every two neighboring points. The refinement degree $rd$ of the grid corresponds to the number of time this procedure has been applied. Therefore $rd=0$ corresponds to the icosahedron itself which possess $N_{p}=12$ nodes, and $N_{c}=20$ spherical triangle cells. As the refinement degree increases, the number of grid points and cells increases as:

	$\displaystyle N_{p}$	$\displaystyle=$	$\displaystyle 2+10\times 4^{rd}$		(102)
	$\displaystyle N_{c}$	$\displaystyle=$	$\displaystyle 20\times 4^{rd}\ .$		(103)

To approximate differential operators, a Voronoi-based finite volume method is chosen. This approach has been widely used in advection-diffusion problems such as in Heikes and Randall (1995); Lazarov et al. (1996); Satoh et al. (2008) and has proven to be efficient to tackle these kind of problems. Since in finite volume methods, differential operators are converted to surface integrals, control volumes (or Voronoi cells) surrounding each grid points have to be defined.

As shown in figure 12 a), each grid point is surrounded by $6$ (or $5$ when the point corresponds to a generator of the initial icosahedron) nodes in its direct neighborhood. From this cluster of node, one can build an ensemble of spherical triangle cells $P_{0}P_{i}P_{i+1}$ all connected together by the common vertex $P_{0}$ . Taking the gravity center $G_{i}$ of each cell allow then to draw a hexagonal (or pentagonal) volume control around the central node $P_{0}$ .

Tomita et al. (2001); Du et al. (2003b) have shown that moving the grid points in a manner that they coincide with the gravity center of their control volume, allows to improve the accuracy of the differential operators to the second-order. The procedure employed here is the Constrained Centroidal Voronoi Tessellations(CCVT) developed by Du et al. (2003a). The principle of this procedure as we implemented it, known as the Lloyd’s method, is the following:

(1)

starting with an initial distribution of nodes on the sphere, taken here as the different subdivision of the icosahedron, as explained previously in this part.
(2)

building the Voronoi cells associated with each grid point.
(3)

moving each node to the gravity center of the cell it is belonging to.
(4)

returning to Step 2 until some convergence criteria is reached. In our case, we imposed that the averaged geodesic distance between the node and the gravity center of the Voronoi cell has to be smaller than $10^{-10}$ .

B.2 Approximation of the differential operators

The horizontal divergence and curl applied to a vector field ${\bf u}({\bf x})$ , and the gradient applied to a scalar field $\Phi$ can be expressed in their integral form as:

$\displaystyle{\bf\nabla_{{}_{H}}}\cdot{\bf u}({\bf x})$	$\displaystyle=\lim_{\Omega\to 0}\frac{1}{\Omega}\int_{\delta_{\Omega}}{\bf u}({\bf x})\cdot{\bf t}\;\mathrm{d}\delta_{\Omega}$	(104)
$\displaystyle{\bf\nabla_{{}_{H}}}\times{\bf u}({\bf x})$	$\displaystyle=\lim_{\Omega\to 0}\frac{1}{\Omega}\int_{\delta_{\Omega}}{\bf u}({\bf x})\cdot{\bf n}\;\mathrm{d}\delta_{\Omega}$	(105)
$\displaystyle{\bf\nabla_{{}_{H}}}\Phi$	$\displaystyle=\lim_{\Omega\to 0}\frac{1}{\Omega}\int_{\delta_{\Omega}}\Phi{\bf n}\;\mathrm{d}\delta_{\Omega}$	(106)

where the control surface as an area $\Omega$ , delimited by the contour $\delta_{\Omega}$ , and the unit tangential and normal vector to the contour are respectively ${\bf t}$ and ${\bf n}$ .

On the discrete sphere, vector and scalar fields are known on the nodes, therefore, to differentiate them, one need to approximate them on the contour of the Voronoi cells. First the different quantities are evaluated on the corners of the cells:

	$\displaystyle{\bf u}(G_{i})$	$\displaystyle=$	$\displaystyle\frac{\alpha_{i}{\bf u}(P_{0})+\beta_{i}{\bf u}(P_{i})+\gamma_{i}{\bf u}(P_{i+1})}{\alpha_{i}+\beta_{i}+\gamma_{i}}$		(107)
	$\displaystyle\Phi(G_{i})$	$\displaystyle=$	$\displaystyle\frac{\alpha_{i}\Phi(P_{0})+\beta_{i}\Phi(P_{i})+\gamma_{i}\Phi(P_{i+1})}{\alpha_{i}+\beta_{i}+\gamma_{i}}\ ,$		(108)

where the areas $\alpha_{i}$ , $\beta_{i}$ , $\gamma_{i}$ are shown in figure 12 b).

Then, following the notation of figure 12 a) the discrete approximation of the different differential operators becomes:

$\displaystyle{\bf\nabla_{{}_{H}}}\cdot{\bf u}({\bf x})$	$\displaystyle\sim$	$\displaystyle\frac{1}{A(P_{0})}\sum_{i=1}^{N_{G}}\widetilde{G_{i}G}_{i+1}\frac{\left({\bf u}(G_{i})+{\bf u}(G_{i+1}\right)}{2}\cdot{\bf n}_{i}$
$\displaystyle{\bf\nabla_{{}_{H}}}\times{\bf u}({\bf x})$	$\displaystyle\sim$	$\displaystyle\frac{1}{A(P_{0})}\sum_{i=1}^{N_{G}}\widetilde{G_{i}G}_{i+1}\frac{\left({\bf u}(G_{i})+{\bf u}(G_{i+1}\right)}{2}\cdot{\bf t}_{i}$
$\displaystyle{\bf\nabla_{{}_{H}}}\Phi$	$\displaystyle\sim$	$\displaystyle\sum_{i=1}^{N_{G}}\widetilde{G_{i}G}_{i+1}\frac{1}{A(P_{0})}\frac{\left(\Phi(G_{i})+\Phi(G_{i+1}\right)}{2}{\bf n}_{i}$
		$\displaystyle-\frac{\Phi(P_{0})}{A(P_{0})}\sum_{i=1}^{N_{G}}\widetilde{G_{i}G}_{i+1}{\bf n}_{i}$

where $A(P_{0})$ corresponds to the area of the control volume, $\widetilde{G_{i}G}_{i+1}$ is the geodesic length between the points $G_{i}$ and $G_{i+1}$ , and $N_{G}$ is the number of vertices of the control volume. Note that when the subscript $i+1=N_{G}+1$ then $i+1=1$ .

References

Velímský (2010) J. Velímský, Physics of the Earth and Planetary Interiors 180, 111 (2010).
Roberts and Scott (1965) P. H. Roberts and S. Scott, Journal of Geomagnetism and Geoelectricity 17, 137 (1965).
Lesur et al. (2010) V. Lesur, I. Wardinski, M. Hamoudi, and M. Rother, Earth, Planets, and Space 62, 765 (2010).
Eymin and Hulot (2005) C. Eymin and G. Hulot, Physics of the Earth and Planetary Interiors 152, 200 (2005).
Backus (1968) G. E. Backus, Royal Society of London Philosophical Transactions Series A 263, 239 (1968).
Holme (2007) R. Holme, in Treatise on Geophysics, edited by E. in Chief:Gerald Schubert (Elsevier, Amsterdam, 2007), pp. 107 – 130, ISBN 978-0-444-52748-6.
Finlay et al. (2010) C. C. Finlay, M. Dumberry, A. Chulliat, and M. A. Pais, Space Science Reviews 155, 177 (2010).
Chulliat and Hulot (2000) A. Chulliat and G. Hulot, Physics of the Earth and Planetary Interiors 117, 309 (2000).
Hulot et al. (1992) G. Hulot, J. L. L. Mouël, and J. Wahr, Geophysical Journal International 108, 224 (1992).
Gillet et al. (2009) N. Gillet, M. A. Pais, and D. Jault, Geochemistry, Geophysics, Geosystems 10, Q06004 (2009).
Pais and Jault (2008) M. A. Pais and D. Jault, Geophysical Journal International 173, 421 (2008).
Hulot et al. (2002) G. Hulot, C. Eymin, B. Langlais, M. Mandea, and N. Olsen, Nature 416, 620 (2002).
Buffett and Christensen (2007) B. A. Buffett and U. R. Christensen, Geophysical Journal International 171, 145 (2007).
Roberts et al. (2003) P. H. Roberts, C. Jones, and A. Calderwood, Earth’s Core and Lower Mantle pp. 100–129 (2003).
Voorhies (2004) C. V. Voorhies, Journal of Geophysical Research 107, 5034 (2004).
Sagaut (2006) P. Sagaut, Large Eddy Simulation for Incompressible Flows: An Introduction, Scientific Computation (Springer, 2006), ISBN 9783540263449.
Lesieur (2008) M. Lesieur, Turbulence in Fluids, Fluid Mechanics and Its Applications (Springer, Berlin, 2008), ISBN 1402064357, 9781402064357.
Baerenzung et al. (2008a) J. Baerenzung, H. Politano, Y. Ponty, and A. Pouquet, Physical Review E 77, 046303 (2008a), eprint 0707.0642.
Baerenzung et al. (2008b) J. Baerenzung, H. Politano, Y. Ponty, and A. Pouquet, Physical Review E 78, 026310 (2008b), eprint 0803.4499.
Fabre and Balarac (2011) Y. Fabre and G. Balarac, Physics of Fluids 23, 115103 (2011).
Sun and Su (2007) O. S. Sun and L. K. Su, in APS Division of Fluid Dynamics Meeting Abstracts (2007), p. K6.
Bülow (2004) T. Bülow, Transactions on Pattern Analysis and Machine Intelligence 26 (2004).
Driscoll and Healy (1994) J. R. Driscoll and D. M. Healy, Advances in Applied Mathematics 15, 202 (1994).
Leonard (1974) A. Leonard, in Turbulent Diffusion in Environmental Pollution (1974), pp. 237–248.
Jackson (1995) A. Jackson, Physics of the Earth and Planetary Interiors 90, 145 (1995).
Bayes (1763) T. Bayes, Phil. Trans. of the Royal Soc. of London pp. 370–418 (1763).
Bloxham (1988) J. Bloxham, Geophysical Research Letters 15, 585 (1988).
Jackson et al. (1993) A. Jackson, J. Bloxham, and D. Gubbins, Washington DC American Geophysical Union Geophysical Monograph Series 72, 97 (1993).
Finlay and Amit (2011) C. C. Finlay and H. Amit, Geophysical Journal International 186, 175–192 (2011).
Thebault and Vervelidou (2013) E. Thebault and F. Vervelidou, submitted to Physics of the Earth and Planetary Interiors (2013).
Mosegaard and Rygaard-Hjalsted (1999) K. Mosegaard and C. Rygaard-Hjalsted, Inverse Problems 15, 573 (1999).
Rygaard-Hjalsted et al. (2000) C. Rygaard-Hjalsted, K. Mosegaard, and N. Olsen, in Methods and Applications of Inversion, edited by P. Hansen, B. Jacobsen, and K. Mosegaard (Springer Berlin Heidelberg, 2000), vol. 92 of Lecture Notes in Earth Sciences, pp. 255–275, ISBN 978-3-540-65916-7, URL http://dx.doi.org/10.1007/BFb0010296.
Sukoriansky et al. (2002) S. Sukoriansky, B. Galperin, and N. Dikovskaya, Physical Review Letters 89, 124501 (2002).
Gamerman and Lopes (2006) D. Gamerman and H. F. Lopes, Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference, Chapman and Hall/CRC Texts in Statistical Science Series (Taylor & Francis, 2006), ISBN 9781584885870.
Lesur et al. (2013) V. Lesur, I. Wardinski, and K. Whaler, in IAGA Scientific Assembly (2013).
Baumgardner and Frederickson (1985) J. R. Baumgardner and P. O. Frederickson, SIAM Journal on Numerical Analysis 22, 1107 (1985).
Stuhne and Peltier (1999) G. R. Stuhne and W. R. Peltier, Journal of Computational Physics 148, 23 (1999).
Heikes and Randall (1995) R. Heikes and D. A. Randall, Monthly Weather Review 123, 1862 (1995).
Lazarov et al. (1996) R. Lazarov, P. Mishev, and P. Vassilevski, Journal of Numerical Analysis 33, 31 (1996).
Satoh et al. (2008) M. Satoh, T. Matsuno, H. Tomita, H. Miura, T. Nasuno, and S. Iga, Journal of Computational Physics 227, 3486 (2008).
Tomita et al. (2001) H. Tomita, M. Tsugawa, M. Satoh, and K. Goto, Journal of Computational Physics 174, 579 (2001).
Du et al. (2003b) Q. Du, M. Gunzburger, and L. Ju, Computational Methods in Applied Mechanics and Engineering 192, 3933 (2003b).
Du et al. (2003a) Q. Du, M. Gunzburger, and L. Ju, Journal of Scientific Computation 24, 1488 (2003a).

	$\displaystyle p(u\|\gamma)$	$\displaystyle=$	$\displaystyle\frac{p(\gamma\|u)p(u)}{p(\gamma)}$
		$\displaystyle\sim$	$\displaystyle\frac{1}{\|N\|^{\frac{1}{2}}}\exp\left[-\frac{1}{2}\left(c-r^{T}N^{-1}r+u^{T}\Sigma_{u}^{-1}u\right)\right]$