Fluid Dynamics of Earth's core: geodynamo, inner core dynamics, core formation

TL;DR

This chapter provides an overview of Earth's core dynamics, including the geodynamo, inner core behavior, and core formation, based on lectures covering magnetohydrodynamics, magnetic field modeling, and core energetics.

Contribution

It offers a comprehensive introduction to Earth's core processes, integrating historical context, theoretical models, and current understanding of geodynamo and inner core dynamics.

Findings

Magnetic field generated by convection-driven geodynamo.

Inner core dynamics influenced by magnetic field effects.

Insights into Earth's core formation processes.

Abstract

This chapter is build from three 1.5 hours lectures given in Udine in april 2018 on various aspects of Earth's core dynamics. The chapter starts with a short historical note on the discovery of Earth's magnetic field and core (section 1). We then turn to an introduction of magnetohydrodynamics (section 2), introducing and discussing the induction equation and the form and effects of the Lorentz force. Section 3 is devoted to the description of Earth's magnetic field, introducing its spherical harmonics description and showing how it can be used to demonstrate the internal origin of the geomagnetic origin. We then move to an introduction of the convection-driven model of the geodynamo (section 4), discussing our current understanding of the dynamics of Earth's core, obtaining heuristically the Ekman dependency of the critical Rayleigh number for natural rotating convection, and introducing the equations and non-dimensional parameters used to model a convectively driven dynamo. The following section deals with the energetics of the geodynamo (section 5). The final two section deal with the dynamics of the inner core, focusing on the effect of the magnetic field (section 6), and with the formation of the core (section 7). Given the wide scope of this chapter and the limited time available, this introduction to Earth's core dynamics is by no means intended to be comprehensive. For more informations, the interested reader may refer to Jones (2011), Olson (2013), or Christensen and Wicht (2015) on the geomagnetic field and the geodynamo, to Sumita and Bergman (2015), Deguen (2012) and Lasbleis and Deguen (2015) on the dynamics of the inner core, and to Rubie et al. (2015) on core formation.

Tables

Table 1. Table 1: Representative values of some properties of the outer core.

radius of the core-mantle boundary	$r_{\mathrm{cmb}}$	3480 km
radius of the inner core boundary	$r_{\mathrm{icb}}$	1221 km
acceleration of gravity at CMB	$g_{\mathrm{cmb}}$	$10.68$ m.s-2
kinematic viscosity	$\nu$	${10}^{-6}$ m2.s-1
thermal diffusivity	$\kappa$	$5\times {10}^{-6}$ m2.s-1
compositional diffusivity	${\kappa }_{\chi }$	${10}^{-9\pm 2}$ m2.s-1
magnetic diffusivity	$\eta$	$1$ m2.s-1
thermal expansion coefficient	$\alpha$	${10}^{-5}$ K-1
compositional expansion coefficient	$\beta$	$1$ wt%-1

Equations

$$\nabla \cdot \mathbf{E}$$

$$\nabla \times \mathbf{E}$$

$$\nabla \cdot \mathbf{B}$$

$$\nabla \times \mathbf{B}$$

$$\frac{\partial {\rho }_c}{\partial t}+\nabla \cdot \mathbf{j}=0,$$

$$\mathbf{j}={\rho }_c\mathbf{u}+\sigma \left(\mathbf{E}+\mathbf{u}\times \mathbf{B}\right),$$

$$\frac{{ϵ}_0}{\sigma }\frac{D{\rho }_c}{Dt}+{\rho }_c=-{ϵ}_0\nabla \cdot \left(\mathbf{u}\times \mathbf{B}\right),$$

$${\rho }_c=-{ϵ}_0\nabla \cdot \left(\mathbf{u}\times \mathbf{B}\right).$$

$$\frac{|{\rho }_c\mathbf{u}|}{|\sigma \mathbf{u}\times \mathbf{B}|}\sim \frac{|{ϵ}_0\nabla \cdot \left(\mathbf{u}\times \mathbf{B}\right)\mathbf{u}|}{|\sigma \mathbf{u}\times \mathbf{B}|}\sim \frac{{ϵ}_0/\sigma }{L/U}.$$

$$\mathbf{j}=\sigma \left(\mathbf{E}+\mathbf{u}\times \mathbf{B}\right).$$

$$\frac{\left|\frac{1}{c^2}\frac{\partial \mathbf{E}}{\partial t}\right|}{\left|\nabla \times \mathbf{B}\right|}\sim \frac{L}{c^{2}T}\frac{E}{B}.$$

$$\frac{\left|\frac{1}{c^2}\frac{\partial \mathbf{E}}{\partial t}\right|}{\left|\nabla \times \mathbf{B}\right|}\sim {\left(\frac{L/T}{c}\right)}^2.$$

$$\left|\frac{1}{c^2}\frac{\partial \mathbf{E}}{\partial t}\right|\ll \left|\nabla \times \mathbf{B}\right|.$$

$$\nabla \cdot \mathbf{E}$$

$$\nabla \times \mathbf{E}$$

$$\nabla \cdot \mathbf{B}$$

$$\nabla \times \mathbf{B}$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{u}\times \mathbf{B})+\eta {\nabla }^2\mathbf{B},$$

$$\eta =\frac{1}{{\mu }_0\sigma }.$$

$$\frac{\partial \mathbf{B}}{\partial t}+\underset{\text{Advection}}{\underbrace{\left(\mathbf{u}\cdot \nabla \right)\mathbf{B}}}=\underset{\text{Stretching}}{\underbrace{\left(\mathbf{B}\cdot \nabla \right)\mathbf{u}}}+\eta {\nabla }^2\mathbf{B},$$

$$\frac{D\mathbf{B}}{Dt}=\left(\mathbf{B}\cdot \nabla \right)\mathbf{u}+\eta {\nabla }^2\mathbf{B}.$$

$$\frac{|\nabla \times (\mathbf{u}\times \mathbf{B})|}{|\eta {\nabla }^2\mathbf{B}|}\sim \frac{UB/L}{\eta B/L^2}=\frac{UL}{\eta }=R_m,$$

$$\frac{\partial \mathbf{B}}{\partial t}=\eta {\nabla }^2\mathbf{B}.$$

$${\tau }_{\eta }=\frac{L^2}{\eta }={\mu }_0\sigma L^2.$$

$$\frac{D\mathbf{B}}{Dt}=\frac{\partial \mathbf{B}}{\partial t}+\underset{\text{Advection}}{\underbrace{\left(\mathbf{u}\cdot \nabla \right)\mathbf{B}}}=\underset{\text{Stretching}}{\underbrace{\left(\mathbf{B}\cdot \nabla \right)\mathbf{u}}}.$$

$$B_{0}\left[\begin{array}[]{c}\partial_{z}u_{x}\\ \partial_{z}u_{z}\end{array}\right]$$

$$B_x$$

$$B_z$$

$$B_x$$

$$B_z$$

$$\mathbf{\delta }(t+dt)$$

$$=\mathbf{\delta }(t)+\mathbf{\delta }\cdot \nabla \mathbf{u}dt+\mathcal{O}({\mathbf{\delta }}^2).$$

$$\frac{D\mathbf{\delta }}{Dt}=\left(\mathbf{\delta }\cdot \nabla \right)\mathbf{u}.$$

$$\frac{d}{dt}{\int }_S\mathbf{B}\cdot d\mathbf{S}={\int }_S\frac{\partial \mathbf{B}}{\partial t}\cdot d\mathbf{S}+{\oint }_{\mathcal{C}}\mathbf{B}\cdot \mathbf{u}\times d\mathbf{l},$$

$${\int }_S\frac{\partial \mathbf{B}}{\partial t}\cdot d\mathbf{S}={\int }_S\nabla \times (\mathbf{u}\times \mathbf{B})\cdot d\mathbf{S},$$

Full text

Given the wide scope of this chapter and the limited time available, this introduction to Earth’s core dynamics is by no means intended to be comprehensive. For more informations, the interested reader may refer to Jones (2011), Olson (2013), or Christensen and Wicht (2015) on the geomagnetic field and the geodynamo, to Sumita and Bergman (2015), Deguen (2012) and Lasbleis and Deguen (2015) on the dynamics of the inner core, and to Rubie et al. (2015) on core formation.

Chapter 1 Introduction

0.1 The birth of geomagnetism

Magnetism, the study of magnetic fields, may be argued to start with the discovery and description of natural magnets, or lodestones. Magnets have been known since at least $\sim 600$ BC in Greece, and $\sim 300$ BC in ancient China. Aristote (384-322 BC) attributes to Thales ($\sim 600$ BC) the observation that “magnets exert a force on iron”. The birth of geomagnetism requires the realisation that there is an ambient magnetic field at Earth’s scale. The primary observation is that a magnet allowed to rotate freely in a horizontal plane points toward the north (as defined by the direction of the polar star). This seems to have been known since at least the I${}^{\text{st}}$ century in China; the compass was used in navigation since at least the X${}^{\text{th}}$ century in China, and since the XII${}^{\text{th}}$ century in Europe.

In the western world, the first scientific treatise describing the properties of magnets is a letter (Epistola) written by Pierre de Maricourt (also known as Petrus Peregrinus) in 1269 to a friend of him, Sygerus de Foucaucourt. Pierre de Maricourt wrote his epistola during a military campaign led by Charles d’Anjou, while the french army was besieging the town of Lucera, in southern Italia. In his letter, Pierre de Maricourt describes some of the most important properties of natural magnets: (i) a magnet has two poles (North and South); (ii) the South pole of a magnet attracts the North pole of another magnet, two identical poles repel each others; (iii) the two poles of a magnet cannot be isolated: breaking a magnet into two pieces gives two magnets, each with two poles (in modern terms, this means that there is no magnetic monopole); (iv) a magnet free to rotate points toward the North pole. In a second part of his letter, Pierre de Maricourt describes the design of a perpetual motion machine using the properties of magnets.

The exact geometry of Earth’s magnetic field turns out to be more complex than what early observations suggested:

Compasses do not point exactly toward the geographic north. The magnetic declination – the angle between the direction of the true geographic north and the direction given by a compass – has been known since at least the end of XIth’s century in China. Christopher Colombus was the first to observe that the magnetic declination varies spatially: while sailing westward from the old world, he observed that the direction given by his compasses changed with longitude along his path. The declination was positive in Europe (i.e. a compass points slightly to the west of the true north), decreased gradually along Columbus’ path, reached zero somewhere in the middle of the Atlantic ocean, and became negative on the western side of the Atlantic. 2. 2.

Furthermore, the actual direction of the magnetic field is tilted from the horizontal plane: when allowed to rotate around an horizontal axis, a magnetised needle points toward a direction which makes a well defined angle with the horizontal, called the magnetic inclination. This has been noticed by Georg Hartmann in 1544 and precisely measured by Robert Norman, who found a declination of 71o50’ in London in 1581.

Perhaps the first scientifically grounded theory of the origin of the geomagnetic field is Gilbert’s theory published in 1600 in his De Magnete, which states that Earth is a magnet. Gilbert’s claim was based on the similarity between the latitudinal variations of the inclination of Earth’s magnetic field and the variations of the magnetic field orientation he measured on a model of Earth, his terella, obtained by carving a lodestone into the shape of a sphere (figure 1). The geometry of Earth’s magnetic field indeed strongly resembles that of a magnet, but it has soon been clear that this theory is at best incomplete. Perhaps the strongest argument against Gilbert’s theory comes from the observation that the geomagnetic field varies with time. Measurements of the declination in London across the XVIIth century have shown variations of several degrees in a few decades; variations of the declination and inclination on a decadal timescale have been well documented by the mid XVIIIth century. Other evidences of fast magnetic field variations include the westward drift of a line of zero declination inferred to be in the middle of the Atlantic ocean in 1701 by Halley (figure 2), which had reached the South American continent about 60 years after the publication of Halley’s map. Modern measurements and inferences have confirmed that the geomagnetic field varies on timescales spanning a very wide range, including monthly to decadal variations. Such fast time variations are hardly compatible with Gilbert’s theory.

0.2 The discovery of the core

The basic structure of the Earth – a solid, rocky mantle surrounding an iron-alloy core which is molten except for a solid inner core at its center – has not been known before the early XXth century. In the XVIIIth and early XIXth centuries, it was widely believed that most of the Earth is molten, except for a thin crust on which we live. A qualitative argument in favour of this view is the observation of molten lava erupting from volcanoes, which is straightforwardly explained if the interior of the Earth is indeed molten. A more quantitative argument is given by the observation that the temperature increases with depth in mines, at a rate such that the melting temperature of most crustal rocks would be reached at a depth lower than 50 km.111The temperature actually does not reach the melting temperature because: (i) the melting temperature increases with pressure; and (ii) the high temperature gradient observed in mines is limited to depth of $\sim 30$ km or less, the temperature gradient becoming much less pronounced at deeper depth because of convective motions in the mantle.

This view was challenged by scientists including Ampère, Poisson, or Kelvin, Ampère arguing that such a thin crust cannot withstand the deviatoric stresses induced by the lunar tides. The fact that the interior of the Earth is solid to a large extent has been demonstrated in the middle of the XIXth century by Lord Kelvin and George Darwin (the son of Charles), by considering the deformation of the interior of the Earth in response to tidal forcing. The oceanic tides that we see at Earth’s surface are the relative motion of the oceans with respect to the crust. We would not see them if the Earth was fully molten since it would deform in the same way as the oceans. The fact that we do see oceanic tides means that the interior of the Earth resists deformation, i.e. has a finite rigidity. Calculations based on available data led Kelvin to state that “the Earth is as rigid as steel”.

The first quantitative model of the interior of the Earth including a rocky mantle and an iron-rich core was proposed in 1896 by Emil Wiechert based on the following constraints and assumptions. The mean density of the Earth (5.6, obtained from the value $g$ at Earth’s surface) is significantly higher than the density of rocks near the surface ($\sim 3$). Wiechert believed that compression alone cannot increase the density of rocks to such extent, and argued that a change of chemical composition at some depth is required. Since only metals were known to have densities higher than 5.6, it was natural for Wiechert to assume the deepest layer (the core) of his model to be metallic. The estimated moment of inertia of the Earth provides further constraints on the repartition of mass within the Earth. By considering a two-layer model with constant density in each layer, the constraints coming from the mean density and moment of inertia of the Earth enable to determine the size of the core and densities of the core and mantle. Wiechert obtained a core radius of about 4970 km, and mantle and core densities of $3.2$ and $8.2$, respectively. His estimate of the core density being close to that of iron at low pressure ($7.8$), Wiechert proposed that the core is composed of iron, and attributes the difference of density to compression.

The end of XIXth century and the beginning of the XXth century witnessed the emergence of seismology as a new and growing branch of geophysics. Milne, Wiechert, Oldham and others realised that the study of seismic waves222Compression waves, called P-waves, and shear waves, called S-waves. propagating within the Earth offers a way to probe the interior of the Earth. By studying the time of arrival of P-waves as a function of the distance from Earthquakes’ epicenters, Oldham (1906) discovered that no P-wave were observed at epicentral distances between $\sim 105^{\mathrm{o}}$ and $\sim 140^{\mathrm{o}}$. He explained this shadow zone by the presence of a discontinuity (which he interpreted as the boundary between the core and mantle of Wiechert’s model), at which the wave propagation velocity decreases, thus refracting the waves downward. Gutenberg (who was a student of Wiechert) estimated the radius of the core to be about $3470\pm 20$ km, which is consistent with the modern estimate of 3480 km. At that time the core was believed to be solid, mostly because of Kelvin and Darwin’s estimates of the rigidity of the Earth. The molten state of the core has been established by Jeffrey in 1926 by comparing the mean rigidity of the Earth to the rigidity of the mantle obtained from estimates of the speed of propagation of seismic waves in the mantle. The mean rigidity of the Earth is significantly smaller than the rigidity of the mantle, which implies that the rigidity of the core must be low.

The final major piece of the basic structure of the Earth has been discovered by Inge Lehmann, a danish seismologist. In spite of what has been said in the previous paragraph, P-waves are sometimes observed in the core shadow zone. These anomalous arrivals, called P’, have been attributed by Lehmann (1936) to P-waves reflected at a discontinuity inside the core at which the seisimic waves velocity increases. She thus discovered the inner core, a sphere 1221 km in radius (Inge Lehmann’s estimate was 1400 km), which was later shown to be solid (Dziewonski and Gilbert, 1971).

By comparing seismological estimates of the density and elastic moduli in the core to the results of high pressure experiments, later work established that the liquid and solid parts of the core are both made predominantly of iron, alloyed with lighter elements (Si, O, S, C, H, …) (e.g. Birch, 1952, 1964; Jephcoat and Olson, 1987; Alfè et al., 2000; Badro et al., 2007). These light elements are more abundant in the liquid outer core (about 10 wt.%) than in the solid inner core (less than 5 wt.%), which is consistent with the idea that the inner core is crystallising from the outer core (Birch, 1940; Jacobs, 1953) since the impurities in an alloy are usually partitioned into the melt upon crystallisation.

Chapter 2 A short introduction to magnetohydrodynamics

Magnetohydrodynamics concerns the interactions between the flow of an electrically conducting fluid and a magnetic field. It thus combines hydrodynamics, which is modelled with the Navier-Stokes and continuity equations, and electromagnetism, which is modelled with Maxwell’s equations, supplemented with Ohm’s law and the conservation of electric charges. Interactions between the flow and magnetic field are twofold: (i) the motion of an electrically conducting material can act as a source term for the evolution of the magnetic field, which can be quantified with the induction equation, a prognostic equation for the evolution of the magnetic field; and (ii) the magnetic field generates a force on the electrically conducting material, the Lorentz force. The goal of this section is to introduce the induction equation, and discuss the form and effects of the Lorentz force.

1 Classical electromagnetism

Maxwell’s equations are

[TABLE]

where $\mathbf{E}$ is the electric field, $\mathbf{B}$ is the magnetic field, $\rho_{c}$ is the electric charge density, $\mathbf{j}$ is the electric current density, $\mu_{0}$ is the permeability of free space, and $\epsilon_{0}$ is the permittivity of free space. Note that $\mu_{0}\epsilon_{0}=1/c^{2}$, where $c$ is the speed of light. Maxwell’s equations are supplemented by an equation expressing that electric charges are conserved, which we write

[TABLE]

and by Ohm’s law, which can be written in a moving electrically conducting material as

[TABLE]

where $\sigma$ is the electrical conductivity, and $\mathbf{u}$ the velocity of the conducting material.

1.1 The charge density in an electrically conducting material

Replacing $\mathbf{j}$ in equation (5) by its expression from Ohm’s law and then using Gauss law (equation (1)), we obtain the following equation for $\rho_{c}$,

[TABLE]

which means that the charge density follows time variations of $\epsilon_{0}\boldsymbol{\nabla}\cdot\left(\mathbf{u}\times\mathbf{B}\right)$ with a response time, or time lag, $\epsilon_{0}/\sigma$. In most liquid metals, $\epsilon_{0}/\sigma$ is between $10^{-16}$ and $10^{-18}$ s and the charge density responds almost instantaneously to changes in $\mathbf{u}$ and $\mathbf{B}$. The charge density must therefore be very close to

[TABLE]

This estimate can be used to discuss the importance of the term $\rho_{c}\mathbf{u}$ in Ohm’s law. Denoting by $L$ and $U$ typical length and velocity scales of a given problem, comparing $\rho_{c}\mathbf{u}$ and $\sigma\mathbf{u}\times\mathbf{B}$ shows that

[TABLE]

The term $\rho_{c}\mathbf{u}$ can therefore be neglected if the macroscopic timescale $L/U$ is large compared to the charge response time $\epsilon_{0}/\sigma$. Since, again, $\epsilon_{0}/\sigma\sim 10^{-18}-10^{-16}$ s, $\rho_{c}\mathbf{u}$ can be neglected in any situation typically encountered in MHD problems, including planetary core dynamics. Ohm’s law can therefore be written as

[TABLE]

1.2 The non-relativistic limit

When applied to the geodynamo problem (as well as most MHD problems), Ampère’s law (equation (4)) can be simplified as follows in the non-relativistic limit. Consider a problem with $\mathbf{B}$ and $\mathbf{E}$ having magnitude $B$ and $E$ varying on typical lengthscale $L$ and timescale $T$. The ratio of the displacement current $\mu_{0}\epsilon_{0}\partial_{t}\mathbf{E}$ to the curl of $\mathbf{B}$ is on the order of

[TABLE]

Using Faraday’s law of induction (equation (2)), which implies that $E/L\sim B/T$, we obtain

[TABLE]

The ratio $L/T$ is an estimate of typical velocities, which in the non-relativistic limit is $\ll c$. This implies that

[TABLE]

This will always be assumed here, and we will use the non-relativistic Maxwell’s equations:

[TABLE]

2 The induction equation

A prognostic equation for $\mathbf{B}$ can be obtained from Ohm’s law and Maxwell’s equations. Taking the curl of Ohm’s law and using Faraday and Ampère’s laws, one can obtain the so-called induction equation, which writes

[TABLE]

where $\eta$, the magnetic diffusivity, is defined as

[TABLE]

An alternative – and sometimes more useful – form of the induction equation can be obtained by noting that $\boldsymbol{\nabla}\times(\mathbf{u}\times\mathbf{B})=\left(\mathbf{B}\cdot\boldsymbol{\nabla}\right)\mathbf{u}-\left(\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{B}$, which allows to transform equation (18) into

[TABLE]

or

[TABLE]

This shows that the Lagrangian derivative of $\mathbf{B}$ depends on a competition between two terms: a stretching term $\left(\mathbf{B}\cdot\boldsymbol{\nabla}\right)\mathbf{u}$, and a diffusion term $\eta\nabla^{2}\mathbf{B}$. Diffusion will always tend to smooth out spatial variations of $\mathbf{B}$. In contrast, the stretching term can increase the magnitude of $\mathbf{B}$.

To estimate the relative importance of these two terms, we denote by $B$ and $U$ the magnitudes of the magnetic and velocity fields, which are both assumed to vary over the same length scale $L$. Forming the ratio of the stretching and diffusion terms, we obtain

[TABLE]

which defines the magnetic Reynolds number $R_{m}$. $R_{m}$ is a measure of the relative importance of stretching and advection of the magnetic field to its diffusion (it could also have been called a “magnetic Péclet number”).

To understand the induction equation, we will first consider the limits of small and large $R_{m}$, before considering a more general case.

2.1 The $R_{m}\ll 1$ limit

If $R_{m}\ll 1$, the advection and stretching terms are both negligible compared to the diffusion term, and the induction equation becomes a simple diffusion equation:

[TABLE]

In this limit, a magnetic field varying over a length scale $L$ would smooth out by diffusion in a timescale

[TABLE]

In the Earth’s core ($L\sim 3000$ km, $\eta\sim 1$ m2.s*-1*), the diffusion timescale is $\tau_{\eta}\sim 300\,000$ years. Properly taking into account the spherical geometry of the core yields a smaller timescale, $\tau_{\eta}=L^{2}/(\pi^{2}\eta)\sim 30\,000$ years, which implies that a spatially varying magnetic field can perdure only during a few tenths of thousand years. Since this is small compared to the age of the Earth (4.56 Gy), this implies that the geomagnetic field cannot be of primordial origin.

In contrast, the magnetic field would diffuse very rapidly in any reasonable size laboratory experiment: in a one meter size experiment using liquid sodium (which has a very low magnetic diffusivity, $\eta\sim 0.1$ m2), the diffusion time is $\sim 10$ s.

2.2 The $R_{m}\gg 1$ limit

At $R_{m}\gg 1$, the induction equation becomes:

[TABLE]

Let us consider the evolution of an initially uniform magnetic field $\mathbf{B}=B_{0}\mathbf{e}_{z}$ in response to a flow with velocity $\mathbf{u}=(u_{x},u_{z})$. At time $t=0$, at which $\mathbf{B}=B_{0}\mathbf{e}_{z}$, the term $\left(B\cdot\boldsymbol{\nabla}\right)\mathbf{u}$ writes

[TABLE]

From this one can see that:

shearing the magnetic field ($\mathbf{u}=u_{x}(z)\mathbf{e}_{x}$) produces a magnetic field component perpendicular to the initial direction of the magnetic field. Solving the induction equation with a simple shear flow ($u_{x}=\dot{\gamma}z$, where $\dot{\gamma}$ is the shear rate) gives

[TABLE]

The magnetic lines are tilted by the shear, and tends to be aligned with the shear direction. Note also that $|\mathbf{B}|$ increases as $|\mathbf{B}|=B_{0}\sqrt{1+(\dot{\gamma}t)^{2}}$. There is therefore a net production – and not just a redistribution – of magnetic energy. 2. 2.

stretching (resp. compressing) the magnetic field increases (resp. decreases) its magnitude. Solving the induction equation with $\mathbf{u}=u_{z}(z)\mathbf{e}_{z}$ gives

[TABLE]

A constant $\partial u_{z}/\partial z$ results in exponential growth (if it is positive) or decrease (if negative) of $B_{z}$.

In the limit of infinite Rm, one can obtain two useful theorems:

Helmholtz theorem:

*magnetic lines are material lines.333 Proof: consider a small vector $\boldsymbol{\delta}$ having material end-points, advected by the flow. At time $t+dt$ this small vector will be equal to

$\displaystyle\boldsymbol{\delta}(t+dt)$ $\displaystyle=-\mathbf{u}(\mathbf{x},t)dt+\boldsymbol{\delta}(t)+\mathbf{u}(\mathbf{x}+\boldsymbol{\delta},t)dt,$

(31)

$\displaystyle=\boldsymbol{\delta}(t)+\boldsymbol{\delta}\cdot\boldsymbol{\nabla}\mathbf{u}\,dt+\mathcal{O}(\boldsymbol{\delta}^{2}).$

(32)

Taking the $dt\rightarrow 0$, $\boldsymbol{\delta}\rightarrow 0$ limit gives

$\frac{D\boldsymbol{\delta}}{Dt}=\left(\boldsymbol{\delta}\cdot\boldsymbol{\nabla}\right)\mathbf{u}.$

(33)

$\boldsymbol{\delta}$ therefore evolves according to the same equation as $\mathbf{B}$ (equation (25)). *

Alfven’s frozen flux theorem:

a magnetic tube is material and its magnetic flux is conserved: denoting by $S$ any cross-section of a magnetic tube, then the magnetic flux $\int_{S}\mathbf{B}\cdot d\mathbf{S}$ (which is constant along the tube because $\boldsymbol{\nabla}\cdot\mathbf{B}=0$) does not vary with time.444Proof: The fact that a magnetic tube is material follows directly from Helmholtz theorem. To show that its magnetic flux does not vary with time, write its time derivative as

$\frac{d}{dt}\int_{S}\mathbf{B}\cdot d\mathbf{S}=\int_{S}\frac{\partial\mathbf{B}}{\partial t}\cdot d\mathbf{S}+\oint_{\mathcal{C}}\mathbf{B}\cdot\mathbf{u}\times d\mathbf{l},$

(34)

use the diffusion-free induction equation to write the first term on the RHS as

$\int_{S}\frac{\partial\mathbf{B}}{\partial t}\cdot d\mathbf{S}=\int_{S}\boldsymbol{\nabla}\times(\mathbf{u}\times\mathbf{B})\cdot d\mathbf{S},$

(35)

and the identity $(\mathbf{A}\times\mathbf{B})\cdot\mathbf{C}=-\mathbf{B}\cdot(\mathbf{A}\times\mathbf{C})$ plus Stokes’ theorem to write the second term as

$\displaystyle\oint_{\mathcal{C}}\mathbf{B}\cdot\mathbf{u}\times d\mathbf{l}$ $\displaystyle=-\oint_{\mathcal{C}}(\mathbf{u}\times\mathbf{B})\cdot d\mathbf{l},$

(36)

$\displaystyle=-\int_{S}\boldsymbol{\nabla}\times(\mathbf{u}\times\mathbf{B})\cdot d\mathbf{S},$

(37)

which is equal to the opposite of the first term on the RHS of equation (34).

2.3 An example of a kinematic solution of the full induction equation

Let us now consider one simple example of solution of the full induction equation. We consider a velocity field of the form $\mathbf{u}=U\sin(2\pi z/L)\mathbf{e}_{x}$ (i.e. a pure shear flow with a shear rate varying periodically with $z$). With a velocity field of this form, the induction equation writes

[TABLE]

The magnetic field is assumed to be initially uniform, $\mathbf{B}(t=0)=B_{0}\mathbf{e}_{z}$. Solving these equations (by looking for a solution of the form $B_{x}=f(t)\cos\left(2\pi z/L\right)$) gives

[TABLE]

where

[TABLE]

The magnitude of $B_{x}$ increases linearly with $t$ when $t\ll\tau$, and saturates at a value $R_{m}B_{0}$ when vertical diffusion balances magnetic field production, which happens at $t\sim\tau$. Compared to the infinite $R_{m}$ limit, there is still a tendency to align $\mathbf{B}$ with the shear, but this is now mitigated by diffusion.

3 The Lorentz force

Experiments show that a charged particle ($q$) moving with velocity $\mathbf{u}$ in an electric field $\mathbf{E}$ and magnetic field $\mathbf{B}$ experiences a force

[TABLE]

called the Lorentz force. This can be generalised to the case of a charged continuous medium: electric and magnetic fields produce a volumetric force given by

[TABLE]

3.1 The non-relativistic limit

In the non-relativistic limit, one can first use Faraday’s law (without the current displacement) to write

[TABLE]

Using the estimate of $\rho_{c}$ for electrical conductors (equation (8)), the relative importance of the two terms of the Lorentz force can be estimated as

[TABLE]

where we have used the scaling relation $E/B\sim L/T\sim U$ obtained from Faraday’s law. In a non-relativistic problem, we can therefore safely neglect the $\rho_{c}\mathbf{E}$ term and write the Lorentz force as

[TABLE]

3.2 Magnetic pressure and tension

A bit of algebra on $(\boldsymbol{\nabla}\times\mathbf{B})\times\mathbf{B}$ allows to decompose the Lorentz force into two terms:

[TABLE]

The first term acts as a pressure term : a force equal to the gradient of the magnetic pressure $B^{2}/(2\mu_{0})$. The second term is called magnetic tension, for reasons which will be explained below. This decomposition is often useful because in many cases the gradient of the magnetic pressure can be balanced to a large extent by pressure gradients, so that only the magnetic tension has a strong dynamic effect (see section 11 for example). Only the magnetic tension can produce vorticity.

We restrict here our analysis to 2D for the sake of simplicity, but the following reasoning can be generalised to 3D. Let us consider a curvilinear Frenet coordinate system $(\mathbf{t},\mathbf{n})$ attached to a given magnetic field line (figure 3). We denote by $s$ the coordinate along this field line. In this coordinate system, we can write

[TABLE]

where $\mathcal{K}=\boldsymbol{\nabla}\cdot\mathbf{n}$ is the (signed) curvature of the field line, the derivative of $\mathbf{t}$ according to $s$ being equal to $-\mathcal{K}\mathbf{n}$ . The operator $\boldsymbol{\nabla}_{\mathbf{t}}(...)=\boldsymbol{\nabla}(...)-\mathbf{n}(\mathbf{n}\cdot\boldsymbol{\nabla})(...)$ is the component of the gradient parallel to the field line.

This equation bears a strong resemblance with the expression of the stress jump induced by interfacial tension across the interface between two immiscible fluids, which is equal to

[TABLE]

where $\gamma$ is the interfacial tension. The first term is the Laplace pressure term, which expresses the fact that interfacial tension induces a pressure jump across a curved interface. This pressure jump is proportional to the interfacial tension and to the curvature of the interface. The second term is a stress tangential to the interface associated with gradients of interfacial tension, which is responsible of Marangoni effects.

Equation (51) has a similar mathematical form with $B^{2}/(2\mu_{0})$ taking the role of interfacial tension. The analog of the Laplace pressure implies that deforming a magnetic line induces a volumetric force proportional to $B^{2}$ normal to the line, and directed toward the center of curvature of the field lines. This force thus acts against deformation of the magnetic lines and tends to straighten them. The analog of the Marangoni tension is a force acting parallel to the magnetic lines toward regions of higher magnetic field intensity. It produces a tension parallel to the fields lines if the magnetic field intensity varies along field lines.

Putting back the magnetic pression and magnetic tension terms together, the Lorentz force writes

[TABLE]

where $\boldsymbol{\nabla}_{\mathbf{n}}$ is the component of the gradient normal to the field line. Both terms acts perpendicularly to the field lines (consistently with the $\frac{1}{\mu_{0}}(\boldsymbol{\nabla}\times\mathbf{B})\times\mathbf{B}$ expression of the Lorentz force).

3.3 Alfvèn waves

The relationship between the Lorentz force and the magnetic field lines curvature suggests that a magnetised fluid can carry waves: if magnetic lines are deformed from an initial equilibrium state by a flow with velocity $\mathbf{u}$, deformation of the lines will produce a restoring force (the Lorentz force) normal to the line and of magnitude proportional to the curvature of the lines (figure 4). Only a velocity field perpendicular to the magnetic lines would curve the lines and thus produce a restoring force. Taking into account that *a velocity perturbation propagating as a plane wave in an incompressible fluid must be a transverse wave555 To show this, consider a plane wave propagating in the $z-$direction, with fluid velocity

$\mathbf{u}=\mathcal{R}\{\mathbf{u}_{0}\exp[i(k_{z}z-\omega t)]\},$

(56)

where $\mathcal{R}$ indicates the real part. If the fluid is incompressible, $\boldsymbol{\nabla}\cdot\mathbf{u}=0$, and hence

$\frac{\partial u_{z}}{\partial z}=-\frac{\partial u_{x}}{\partial x}-\frac{\partial u_{y}}{\partial y},$

(57)

which is equal to 0 since the plane wave assumption implies that the velocity field is a function of $z$ and $t$ only. This implies that the component $u_{z}$ of the velocity field must be spatially uniform. In other words, only the component of the velocity field perpendicular to the propagation direction can oscillate spatially: the wave must be transverse. , we can thus expect velocity perturbations perpendicular to the field lines to propagate in the form of a transverse wave travelling in the direction of magnetic lines. These are called Alfvèn waves. The celerity of these waves can be obtained from dimensional analysis: the propagation velocity $v_{A}$ can be a function of the magnetic field intensity $B_{0}$ [T=kg s-2* A*-1*], of the magnetic permeability $\mu_{0}$ [m kg s*-2* A*-2*], of the density $\rho$ [kg m*-3*] of the fluid, and of the wave number $k$ [m*-1*] of the wave, which are the only parameters of the problem. There is only one way to build a group of parameters with the dimension of velocity from this list of parameters (this is proved with the Vaschy-Buckingham theorem), which is

[TABLE]

An important point is that it is not possible to build a velocity involving the wave number $k$, which implies that Alfvèn waves must be non-dispersive.

A classical way of deriving the dispersion equation of these waves is to work directly from the Navier-Stokes (including the Lorentz force) and induction equations.666Consider small perturbations of the velocity and magnetic fields, linearise the Navier-Stokes and induction equations, take the curl of these two equations, and combine them after taking the time derivative of the curled Navier-Stokes equation (vorticity equation).

As an alternative, we will here obtain the wave equation from an analysis based on magnetic tension. We consider a conducting, incompressible fluid of uniform density permeated by a magnetic field, and neglect magnetic diffusion in the induction equation (thus taking the infinite Rm limit). Noting that $\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{u}=0$ (since we are considering plane transverse waves), the Navier-Stokes and induction equations then write

[TABLE]

We consider a uniform background magnetic field $\mathbf{B}_{0}=B_{0}\mathbf{e}_{z}$, which is slightly perturbed by a small velocity field perturbation. By “slightly pertubed”, we mean here that the curvature $\mathcal{K}$ of the magnetic lines is assumed to remain small. The restoring Lorentz force being perpendicular to $\mathbf{B}$, we expect oscillations of the fluid perpendicular to $\mathbf{B}$, in the $x$ direction, and, since only transverse waves can be carried by an incompressible fluid, propagation parallel to the magnetic fields lines, in the $z$ direction. A transverse wave propagating along a field line will shear the magnetic line and produce by induction a magnetic field in the $x$ direction. We thus consider small velocity and magnetic fields perturbations of the form $\mathbf{u}=u_{x}(z,t)\mathbf{e}_{x}$ and $\mathbb{b}=b_{x}(z,t)\mathbf{e}_{x}$, small curvature $\mathcal{K}(z)$, and linearize equations (59) and (60) to obtain

[TABLE]

We denote by $\delta_{x}$ the $x$-displacement of a magnetic field line. Since in the infinite Rm limit magnetic field lines are material lines, the horizontal displacement of the field lines is linked to the velocity field by

[TABLE]

With $\mathbf{n}\simeq\mathbf{e}_{x}-\frac{\partial\delta_{x}}{\partial z}\mathbf{e}_{z}$, the curvature is

[TABLE]

The $x$-component of Navier-Stokes then writes

[TABLE]

which is a non-dispersive wave equation, with celerity

[TABLE]

The period $T$ of theses waves is related to their wavelength $\lambda$ by

[TABLE]

Chapter 3 The geometry of Earth’s magnetic field

4 Spherical harmonics decomposition

Figure 5a shows a map of the intensity of the radial component $B_{r}$ of the geomagnetic field at the surface of the Earth. As already discussed in the introduction, and first noted by Gilbert in 1600, Earth’s magnetic field strongly resembles the magnetic field which would be produced by a magnet inside the Earth and aligned with the rotation axis: the field is strongly dipolar, close to axisymmetric, with negative $B_{r}$ in the northern hemisphere and positive $B_{r}$ in the southern hemisphere. However, closer inspection of figure 5a reveals deviations from a NS-oriented dipolar field and smaller scales details. For example, the position of the magnetic equator (the line at the surface where $\mathbf{B}$ is horizontal) wanders quite far from the geographic equator. Figure 5a also shows patches of stronger field intensity under North America, Siberia, and the south Indian ocean.

The fact that the mantle, crust, and lower atmosphere of the Earth can be considered as current-free regions (they are good electric insulators) greatly simplifies the mathematical description of the geomagnetic field: Ampère’s law (equation (17)) with $\mathbf{j}=0$ writes $\boldsymbol{\nabla}\times\mathbf{B}=0$, which implies that $\mathbf{B}$ can be written as the gradient of a scalar field, the geomagnetic potential $V$:

[TABLE]

Since $\mathbf{B}$ is a divergence free vector field and $\boldsymbol{\nabla}\cdot\boldsymbol{\nabla}(...)=\Delta(...)$, the geomagnetic potential obeys Laplace equation,

[TABLE]

The general solution of this equation can be written as

[TABLE]

where $(r,\theta,\phi)$ are the usual spherical coordinate systems (radius, colatitude, longitude) as defined on figure 6, and $Y_{lm}(\theta,\phi)$ are the spherical harmonics, $l$ and $m$ being the degree and order, which form a complete set of orthogonal functions on the sphere.

4.1 Demonstrating the internal origin of the geomagnetic field

Written this way, the geomagnetic potential can be seen as the sum of terms of internal origins, which decrease with increasing $r$ (the $A_{lm}$ terms), and terms of external origins, which increase with increasing $r$ (the $B_{lm}$ terms). By determining the $A_{lm}$ and $B_{lm}$ coefficients, it is therefore possible to determine whether the magnetic field observed at the surface of the Earth is predominantly of internal or external origin. In his General Theory of Terrestrial Magnetism published in 1839, Gauss introduced the geomagnetic potential and its spherical harmonics expansion. Considering only terms of internal origin, Gauss determined the coefficients of the expansion up to degree $l=4$ using least-square inversion from intensity, inclination, and declination maps available at that time. From the excellent agreement between the observations and his spherical harmonics expansion, he concluded that the geomagnetic field must be of internal origin.

The internal and external terms only differ by the way they vary with $r$, and it may not be obvious at first sight how it is possible to differentiate between them from measurements all made at the same radius. But the actual data are measurements of the magnetic field $\mathbf{B}=-\boldsymbol{\nabla}V$, and not $V$. The latitudinal variations of the direction of $\mathbf{B}$ at a given $r$ happens to be sensitive to the origin of the magnetic field. To get a feeling of how this works, let us consider the $l=1$ terms only, and write the potential as

[TABLE]

From this we obtain the corresponding magnetic field components:

[TABLE]

and calculate the inclination $I$ (the angle between $\mathbf{B}$ and the horizontal plane) as

[TABLE]

where $\lambda=\pi/2-\theta$ is the latitude. If the geomagnetic field is either of internal origin only ($B_{10}=0$), or of external origin only ($A_{10}=0$), then

[TABLE]

Figure 7 shows the inclination $I$ as a function of $\lambda$ for the IGRF 2010 geomagnetic model (grey diamonds), and for the data used by Gauss (blue diamonds). The data points are quite dispersed around the $\tan I=2\tan\lambda$ curve, but this is mostly due to the fact that the geomagnetic dipole is actually tilted from the rotation axis, and also because Earth’s magnetic field includes non-negligible higher degree terms.

From now on we will thus keep only the contributions of internal origin. The spherical harmonics can be written as $Y_{lm}=\mathrm{e}^{im\phi}P_{lm}(\cos\theta)$, where $P_{lm}$ are the associated Legendre polynomials. Keeping only the terms in $r^{-l-1}$ in equation (70) and using real coefficients, one can write the magnetic potential as

[TABLE]

where $g_{l}^{m}$ and $h_{l}^{m}$ are the Gauss coefficients, and $r_{\Earth}$ is the radius of the Earth.

It is useful to look at the spectrum of magnetic energy as a function of the spherical harmonic degree $l$. Using the fact that spherical harmonics form an orthogonal basis, one can show that the magnetic energy at a radius $r$ corresponding to all components of degree $l$ is given by

[TABLE]

The resulting spectrum is called the Lowes-Mauersberger spatial power spectrum. Figure 5c shows this spectrum at Earth’s surface in blue circles. The energy of the dipolar components ($l=1$) is more than an order of magnitude larger than the $l=2$ component. $R_{l}$ then decreases with $l$ up to $l$ equal 13 or 14 where it reaches a plateau. At $l>13$, the magnetic energy is dominated by contributions of the crustal magnetisation, which obscures the field from deeper origin. Only the $l\leq 13$ spherical harmonics components can be attributed to the magnetic field originating from the core.

4.2 The magnetic field at the core-mantle boundary

One very interesting property of a curl-free magnetic field is that it is possible to extrapolate to other radii the magnetic potential once we know the coefficients of the spherical harmonics decomposition. Since the mantle of the Earth can be considered to be current-free, it is in particular possible to extrapolate the surface magnetic field down to the base of the mantle, though this can be done only for harmonic degrees $l\leq 13$ since the higher degrees components of the core field are obscured by the crustal field. The potential for $l\leq 13$ at the core-mantle boundary can be obtained by writing $V$ at $r=r_{\mathrm{cmb}}$ from equation (78). From this one can obtain the associated magnetic field map and spectrum.

Figure 5b shows the radial component of the field at the core-mantle boundary (CMB). Though the field is still quite dipolar, it exhibits much more smaller scale variations than the surface field (figure 5a). This is confirmed by inspection of the energy spectrum (figure 5c, orange circles): the spectrum is about flat for $l\geq 2$, with the energy of the dipolar components still about an order of magnitude larger than the higher $l$ components. The fact the dipolar components still stand out at the CMB is an important result. Since the spherical components decrease with $r$ as $r^{-l-1}$, any locally produced magnetic field would appear dipolar when seen from far enough. The fact that the field at the CMB is still dominated by the $l=1$ terms suggests that the dipolar nature of the geomagnetic field is a robust feature of the geodynamo.

4.3 The field within the core: poloidal-toroidal decomposition

In planetary cores, $\mathbf{B}$ cannot be considered to be curl-free anymore ($\mathbf{j}\neq 0$), and thus cannot be written as the gradient of a potential. However, the fact that $\mathbf{B}$ is a divergence-free vector field allows to write it as

[TABLE]

where $T$ and $P$ are the poloidal and toroidal potentials. Note that the toroidal part has no radial component, which means that the magnetic field reconstructed at the CMB only corresponds to the poloidal part. We have no direct constraints on the toroidal part of the magnetic field in the core.

5 Core flow inversion

The magnetic field at the CMB evolves with time in a measurable way. These variations can be interpreted with the help of the induction equation, remembering that at the CMB only the radial component of the magnetic field is known. Taking the dot product of the induction equation with $\mathbf{e}_{r}$ gives

[TABLE]

where $\boldsymbol{\nabla}_{h}\cdot(...)$ is the horizontal part of the divergence operator, and $\mathbf{u}_{h}=(0,u_{\theta},u_{\phi})$ is the horizontal part of the velocity field. The diffusion term is usually neglected, which can be justified a posteriori on the basis that the magnetic Reynolds number based on the smallest spatial scale considered ($\sim 1000$ km) and estimated velocity is $\sim 500$. Dropping the diffusion term, equation (81) writes

[TABLE]

Knowing $B_{r}$ and its time derivative, one can in principle invert equation (82) to obtain the horizontal velocity field just below the CMB. This happens to be a severely ill-posed inverse problem, the most obvious reason being that we are trying to estimate a two-components vector field ($\mathbf{u}_{h}$) from a scalar field ($B_{r}$). In other words, we have only one equation for two unknowns. Inverting equation (82) thus requires a second equation for $\mathbf{u}_{h}$. Various assumptions have been made (e.g. steady flow, toroidal flow, tangentially geostrophic flow, columnar flow, quasi-geostrophic flow, …), and the choice does impact the resulting flow pattern (see Holme (2015) for a review). Robust features of the inverted flows include a strong westward flow under the Atlantic, a much weaker flow under the Pacific, and some degree of symmetry about the equatorial plane. The RMS flow velocity is around $12$ to $14$ km per year, or about $4\times 10^{-4}$ m.s*-1*.

Chapter 4 Basics of planetary core dynamics

6 The geodynamo hypothesis

As discussed in section 1, Gilbert’s claim that Earth is a magnet has been dismissed by the observation of fast temporal changes of the geomagnetic field. In addition, we now know that permanent magnets (i.e. ferromagnetic or ferrimagnetic materials) lose their permanent magnetic properties (by becoming paramagnetic) above a critical temperature called the Curie temperature. Magnetite and iron have Curie temperatures of 858 K and 1043 K respectively. The temperature in the Earth exceeds these temperatures at depths larger than around $100-150$ km (Jaupart and Mareschal, 2010), which confines ferromagnetism to rather shallow depths. Magnetisation of crustal material can in fact be quite strong, but cannot explain the large scale part of Earth’ s magnetic field.

Leaving aside ad-hoc theories, we are thus left with the MHD equations introduced in section 2 to understand the origin of Earth’s magnetic field. We have already shown (section 2) that the geomagnetic field cannot be of primordial origin, and must be sustained in some way against ohmic dissipation: solving the induction equation with no velocity field indeed shows that spatial variations of the magnetic field would be smoothed out by diffusion on a timescale of $\sim 30\,000$ years, while we know from paleomagnetism that the geomagnetic field has been sustained for (at least) the last $\sim 3.5$ Gy.

6.1 Self-exciting dynamos

The current theory for the origin of the Sun and Earth’s magnetic fields has been proposed in 1919 by Sir Joseph Larmor in a meeting of the British Association for the Advancement of Science. In the short report of his presentation (Larmor, 1919), he wrote:

“In the case of the Sun, surface phenomena point to the existence of a residual internal circulation mainly in meridian planes. Such internal motion induces an electric field acting on the moving matter; and if any conducting path around the solar axis happens to be open, an electric current will flow round it, which may in turn increase the inducing magnetic field. In this way it is possible for the internal cyclic motion to act after the manner of the cycle of a self-exciting dynamo, and maintain a significant magnetic field from insignificant beginnings, at the expense of some of the energy of the internal circulation.

[…]

The very extraordinary feature of the Earth magnetic field is its great and rapid changes, comparable with its whole amount. […] [a self-exciting dynamo] would account for magnetic change, sudden or gradual, on the Earth merely by change of internal conducting channels: though, on the other hand, it would require fluidity and residual circulation in deep-seated regions.”

Larmor’s proposition was inspired by the self-exciting dynamos developed in the second half of 19th century, which were the first efficient electric generators. The underlying mechanism can be understood as a positive feedback loop involving Faraday’s law of induction and Ampère’s law. Consider an electrically conducting material moving into an arbitrarily small seed magnetic field. Faraday’s law of induction tells that time variations of a magnetic field in an electrically conducting material produces electric currents. Equivalently, the motion of the conducting material into the seed magnetic field produces electric currents. According to Ampère’s law, these electric currents would themselves produce a magnetic field. If the orientation of the electric current is such that the induced magnetic field reinforces the initial seed magnetic field, then the intensity of the magnetic field can grow.

The word dynamo was coined by Faraday to name the electric generator he invented, now called the Faraday disk. Faraday’s disk consists in a copper disk rotating within a magnetic field produced by a horseshoe magnet. The rotating motion of the copper disk produces a difference of electric potential between its centre and periphery by virtue of Faraday’s law of induction, and hence an electrical current if the center and periphery of the disk are linked trough an electric circuit. Self-exciting dynamos are based on a similar principle except that the permanent magnet is replaced by electromagnets fed by the induced electric current. This creates a positive feedback loop which increases the intensity of the magnetic field, thus increasing the currents produced by induction. The concept has been formulated by Anyos Jedlik around 1856; practical designs of working self-exciting dynamos have been patented by Varley in 1866, and presented by Wheatstone and Siemens in 1867. These self-exciting dynamos were able to produce a much higher power output than permanent magnet dynamos, which opened the way to the industrial use of electricity. A simple conceptual model of a self-exciting dynamo was proposed by Bullard and Gellman (1954) as a toy model of the geodynamo (figure 8).

Compared to self-exciting dynamos such as Bullard’s, the concept of geodynamo faces several additional difficulties. To quote Bullard and Gellman (1954): “A central problem […] is to determine whether there exist motions of a simply connected, symmetrical fluid body which is homogeneous and isotropic that will cause it to act as a self-exciting dynamo […]. We call such dynamos ’homogeneous’ to distinguish them from the dynamos of the electrical engineer, which are multiply connected and of low symmetry.” In industrial dynamos, the electric currents produced by induction are fed into electric wires which are arranged in such a way that the induced magnetic field indeed reinforces the seed field. In Earth’s core, no such wiring exists and the electric currents path is determined from Ohm’s law by the direction of the electromotive force (and thus ultimately by the geometry of the velocity and magnetic fields). The motion of molten iron in the core has to be such that, on average, the induced electric currents can maintain the geomagnetic field.

As formulated by Bullard and Gellman (1954), the question is not (yet) whether a free-flowing fluid can maintain a magnetic field in spite of the expected feedback of the Lorentz force onto the flow. At first no such feedback was considered: early proponents of the geodynamo hypothesis first tried to find kinematic dynamos, where only the induction equation is solved, the velocity field being chosen in the hope that it could produce dynamo action. This is a difficult problem, even if the velocity field is prescribed, but working kinematic dynamos have indeed been found. Classical exemples include the dynamos found by Herzenberg (1958), G.O. Roberts (1972), and Ponomarenko (1973).

The first laboratory demonstration of homogeneous dynamo action was given by Lowes and Wilkinson (1963) with an apparatus inspired by Herzenberg’s dynamo. Lowes & Wilkinson’s apparatus consists in two rotating solid cylinders made of a high magnetic permeability iron alloy (“Perminvar”) with their axes at right angles, embedded in a solid block of the same material (figure 8). The electrical contact between the cylinders and the block is ensured by a thin layer of liquid mercury. Denoting by $\omega$ the rotation rate of the cylinders, $a$ their radii, and $\lambda=1/(\mu\sigma)$ the magnetic diffusivity of the “Perminvar” alloy, Lowes and Wilkinson (1963) found that their apparatus operates as a self-exciting dynamo when the non dimensional number $\omega a^{2}/\lambda$ exceeds a critical value. This number is a magnetic Reynolds number (with velocity scale $\omega a$), so the fact that it must exceed a critical value for dynamo action is consistent with our analysis of the induction equation (section 2). In spite of being made of two materials (Perminvar plus a thin layer of mercury), this is effectively a homogeneous dynamo in the sense that the electric currents are not forced into specific paths.

A next critical step for establishing the viability of the geodynamo hypothesis has been to obtain dynamos where the velocity field is not imposed a priori, and where the feedback of the magnetic field on the flow is taken into account (via the Lorentz force). Fluid dynamos obtained numerically by self-consistently solving the coupled Navier-Stokes, induction, and heat transfer equations have started to appear in the eighties in the contexts of the solar dynamo (Gilman and Miller, 1981; Glatzmaier, 1984, 1985a, 1985b) and geodynamo (Zhang and Busse, 1988; Glatzmaier and Roberts, 1995; Kageyama et al., 1995). In spite of a number of shortcomings, these numerical dynamos driven by convective motions have been successful in reproducing some of the more salient features of Earth’s magnetic field (more on this in section 8).

Experimental liquid homogeneous self-exciting dynamos have been obtained more recently, in two apparatus developed in Riga (Latvia) (Gailitis et al., 2001) and Karlsruhe (Germany) (Stieglitz and Müller, 2001). In both experiments liquid sodium was forced by propellers into a system of pipes, which effectively imposed the velocity field. Liquid sodium is the best electrically conducting liquid available in large quantities (it has a magnetic diffusivity $\eta=0.1$ m2.s*-1*, about ten times lower than molten iron), and has been systematically used as working fluid in dynamo experiments. In both apparatus the imposed flow is inspired by known kinematic dynamos: the Riga experiment is based upon the Ponomarenko dynamo, and the Karlsruhe experiment is based upon the G.O. Roberts dynamo. The Riga experiment used 2 m3 of sodium and required a power input of 200 kW; the Karlsruhe experiment used $1.6$ m3 of sodium and required a power input of 630 kW.

Obtaining experimentally a dynamo with a much less constrained flow proved to be even more arduous. To date the only successful free-flowing experiment is the VKS experiment developed in Cadarache (France) (Monchaux et al., 2007; Berhanu et al., 2007). Though numerical dynamo calculations are now well established, liquid metal experiments are still very valuable tools for understanding planetary core dynamics (e.g. Nataf and Gagnière, 2008; Cabanes et al., 2014). On certain aspects liquid metal experiments are dynamically closer to Earth’s core conditions than numerical simulations. In particular, liquid sodium has a magnetic Prandtl number (the ratio of kinematic viscosity to the magnetic diffusivity) of about $10^{-6}$, similar to that of Earth’s core, while working numerical dynamos have been so far limited to magnetic Prandtl number values above $\sim 0.1$.

6.2 What drives the geodynamo?

In parallel to the question of the feasibility of a liquid homogeneous dynamo, a central question has been the source of power of the geodynamo. What drives the motion and provides the energy that is being lost by ohmic dissipation?

In industrial devices or laboratory experiences, the power is provided by an external mean and therefore is not, theoretically speaking, an issue (though it can be a practical issue, since driving a liquid dynamo in the laboratory does require a quite large power input).

In Earth’s core, the possible sources of motion and power fall into two broad categories: (i) natural convection, either of thermal origin (Bullard, 1949, 1950; Verhoogen, 1961) or compositional origin (Braginsky, 1963; Gubbins, 1977; Loper, 1978; O’Rourke and Stevenson, 2016); and (ii) stirring produced by astronomical forcing, i.e. motions forced in the core by either tidal forcing or changes of the mantle rotation vector (in direction – precession or nutation – or intensity – libration) (Bondi and Lyttleton, 1948; Bullard, 1949; Malkus, 1963, 1968). The second class of models is the subject of the chapter by Le Reun and Le Bars; we will focus here on natural convection.

Core convection, whether it be thermal or compositional, is controlled by the rate at which the core loses heat to the mantle. From a thermal point of view, the core is the slave of the mantle: The heat flux from the core to the mantle is dictated by the efficiency of heat transport by mantle convection, and core convection is thus tied to mantle convection. Cooling of the core can drive convective motions is several ways:

Cooling the core from above can potentially produce thermal convection, if the imposed heat flux is larger than the flux which can be conducted along an adiabat in the core (see section 7). 2. 2.

The secular cooling of the core is responsible for its slow solidification and the formation of the solid inner core. In spite of the core being colder at the CMB than at deeper depth, the core started to solidify at its center and the inner core is now growing outward. The reason for this is that the solidification temperature of the core material increases with pressure faster than the actual temperature (Jacobs, 1953). If the core started hot and fully molten and then gradually cooled down, the temperature first reached the solidification temperature at the center of the core, which allowed the inner core to nucleate. Further cooling results in the slow growth of the inner core.

One key point here is that the core is made of an iron-rich alloy rather than pure iron. Since the density of the core is lower than that of pure iron, we know that the main impurities are “light elements” (i.e. lighter than iron), likely a mixture of mainly oxygen, sulfur, and silicium. Upon solidification, these elements are partitioned preferentially into the liquid outer core: solidification of the inner core thus results in an outward flux of light elements at the inner core boundary (ICB), which can drive compositional convection. In addition, the latent heat of solidification helps thermal convection by slowing down the cooling of the inner core boundary. 3. 3.

Finally, cooling of the core may also be at the origin of a compositional flux across the core-mantle boundary.

Light elements like oxygen, silicium, or magnesium, are present both in the core and mantle (composed predominently of silicates and oxydes of iron and magnesium). Since the core and mantle materials in contact at the core-mantle boundary should be very close to thermodynamic equilibrium, the relative abundance of these elements in the core and silicates at the CMB is set by the partitioning coefficient of these elements between silicates and the core alloy. This in general depends on temperature: cooling of the core thus modifies this chemical equilibrium and results in a flux of elements between mantle and core. If light elements are transferred from the mantle to the core, then a stably stratified layer may form below the CMB. If in the other way, the flux of light elements would drive compositional convection in the core. The actual flux of element is limited by the rate at which convection in the mantle provides “fresh” material which can react with the molten iron of the core; again, mantle convection controls to a large part the buoyancy flux available to drive core convection.

Another possible mechanism is exsolution of light elements from the core. It has been proposed that MgO (O’Rourke and Stevenson, 2016; Badro et al., 2016) and SiO2 (Hirose et al., 2017) can precipitate (or exsolve) from the core alloy if initially abundant enough. The saturation concentration of these species happens to decrease with decreasing temperature, so a gradual cooling of the core results in the progressive removal of magnesium, oxygen, and silicon from the core. If exsolution happens predominantly in the vicinity of the core-mantle boundary, then removing MgO and SiO2 from the core induces a buoyancy flux (by leaving behind an iron-rich, dense melt) which can drive convection in the core.

7 Rotating convection

7.1 Governing equations and non-dimensional parameters

The liquid outer core is modelled as a spherical shell of outer and inner radii $r_{\mathrm{cmb}}$ and $r_{\mathrm{icb}}$. The shell thickness is denoted by $H=r_{\mathrm{cmb}}-r_{\mathrm{icb}}$. In what follows, we will use either a spherical coordinate system $(r,\theta,\phi)$ or a cylindrical coordinate system $(s,\phi,z)$ (figure 6). The core is rotating at a rate $\Omega$, and the gravity field is $\mathbf{g}=-g(r)\mathbf{e}_{r}$. We assume that there is a positive temperature difference $\Delta T$ between the inner and outer boundaries. We denote by $\nu$ the kinematic viscosity of the liquid core, by $\rho$ its density, by $\kappa$ its thermal diffusivity, by $\alpha$ its thermal expansion coefficient, and by $c_{p}$ is specific heat capacity. Though this is a somewhat dubious assumption, the outer core is often modelled as a Boussineq fluid (which in particular includes the assumption that the fluid is incompressible). Molten iron in the outer core is also typically assumed to behave as a newtonian fluid with temperature-independent kinematic viscosity.

Rotating thermal convection is governed by the Navier-Stokes equation in a rotating frame of reference, mass conservation, and a transport equation for temperature $T$. Denoting by $\mathbf{u}$ the velocity field, by $P$ the pressure, and by $T$ the temperature, this set of equations can be written (under the Boussinesq approximation), as

[TABLE]

These equations can be made dimensionless with a variety of different choices of scales. For exemple, scaling lengths by the outer core thickness $H$, time by $H^{2}/\nu$, velocity by $\nu/H$, temperature by the difference of temperature $\Delta T$ across the shell, and pressure by $\rho\Omega\nu$, gives

[TABLE]

where the Ekman, Rayleigh, and Prandtl numbers are defined as

[TABLE]

The quantity $RaE/Pr$ is sometimes taken as a “modified Rayleigh number” defined as

[TABLE]

If a heat flux $q$ is imposed at the boundaries rather than a temperature difference, $Ra$ and $Ra^{*}$ can be modified to give flux-based Rayleigh numbers by replacing $\Delta T$ by $qH/k$, where $k$ is thermal conductivity.

The Ekman number is a measure of the ratio of the viscous forces and Coriolis acceleration, if the viscous forces are estimated assuming a flow varying spatially on a length scale $\sim H$:

[TABLE]

It is about $10^{-15}$ for the outer core, which means that viscous forces would become of importance only at relatively small length scales (more about this later).

It can also be useful to consider the vorticity equation, obtained by taking the curl of equation (86),

[TABLE]

which can usually be simplified to

[TABLE]

by assuming that the vorticity $\boldsymbol{\zeta}$ is small compared to the rotation rate $\mathbf{\Omega}$. Here the term $2\Omega\frac{\partial\mathbf{u}}{\partial z}$ is the production of vorticity through stretching of planetary vorticity, $\alpha g\boldsymbol{\nabla}T\times\mathbf{e}_{r}$ is the baroclinic production of vorticity, and $\nu\boldsymbol{\nabla}^{2}\boldsymbol{\zeta}$ is the viscous diffusion of vorticity.

7.2 Geostrophy

The very small value of the Ekman number of the core suggests that viscous forces may be neglected, at least as long as only large scale motions are considered. Keeping only the Coriolis acceleration and the pressure gradient (thus neglecting also inertia and buoyancy), the Navier-Stokes equation reduces to

[TABLE]

which is called the geostrophic balance. Taking the curl of equation (96), which gives

[TABLE]

shows that the velocity field is invariant along the rotation axis if this balance holds. This is known as the Taylor-Proudman theorem.

In a rotating container with sloped boundaries, such as planetary cores, the motion is further restricted by the condition that no fluid can cross the boundary. In a spherical container, this condition writes $\mathbf{u}\cdot\mathbf{e}_{r}=0$. In other words, the velocity field at the boundary can only have $\theta$ and $\phi$ components. Of these, only the $\phi$ component is allowed by the Taylor-Proudman constraint. This means that the only possible motions are longitudinal (zonal) motions of the form

[TABLE]

Mass conservation ($\boldsymbol{\nabla}\cdot\mathbf{u}=(r\sin\theta)^{-1}\partial u_{\phi}/\partial\phi=0$) further implies that $\mathbf{u}$ cannot depend on $\phi$. The only possible geostrophic flows in a spherical container thus consist in rigid cylinders rotating about Earth’s axis of rotation, called geostrophic cylinders:

[TABLE]

7.2.1 Torsional oscillations

Though it is somewhat of a side note in the context of a thermal convection-oriented section, it is worth discussing quickly the effect of a magnetic field on geostrophic motions. The addition of a magnetic field permeating the geostrophic cylinders allows the propagation of a subclass of Alfvèn waves called torsional waves, or torsional oscillations (Braginsky, 1970).

The basic mechanism behind these waves can be easily understood using the concept of magnetic tension introduced in section 3. Assume one geostrophic cylinder rotates at a different rate (figure 9a). The differential rotation of this cylinder with respect to its neighbours will bend the magnetic lines as shown on figure 9b. The resulting magnetic tension will produce a restoring force of magnitude inversely proportional to the azimuthal displacement of the cylinder. The net effect on the geostrophic cylinder is a restoring torque. As shown in section 3 Alfvèn waves are longitudinal waves; torsional waves therefore must propagate in the $s$-direction.

A dispersion equation for these waves can be obtained by applying conservation of angular momentum to a geostrophic cylinder of inner radius $s$ and thickness $ds$. By restricting the analysis to small relative motions of the cylinders, one find that the angular velocity $\omega(s)$ obeys the following equation:

[TABLE]

where

[TABLE]

is the r.m.s. value of $B_{s}$ on the geostrophic cylinder of radius $s$ (Braginsky, 1970).

The propagation velocity of these oscillations being proportional to $\bar{B}_{s}$, identifying torsional oscillations in core flow inversions can provide an estimate of the magnetic fields intensity within the core. Gillet et al. (2010) found a 6-year period signal in their core flow inversions which they interpreted as a torsional wave. It takes about 4 years for this wave to propagate through the outer core; the resulting propagation velocity gives $\bar{B}_{s}\sim 2$ mT. For comparison, the radial component of the magnetic field at the CMB peaks at about 1 mT (figure 5).

7.3 Onset of thermal convection

We now turn to the question of the initiation of thermal convection in a rotating planetary core. Linear stability analysis shows that at low Ekman number the critical Rayleigh number $\mathrm{Ra}_{\mathrm{cr}}$ and the horizontal length scale $\ell_{c}$ of the first unstable mode scale as

[TABLE]

The linear stability analysis is challenging (Roberts, 1968; Busse, 1970; Jones et al., 2000; Dormy et al., 2004); to try to understand this result at a lower mathematical cost, we will take a heuristic approach and derive the Ekman number dependency of $\mathrm{Ra}_{\mathrm{cr}}$ from simple physical arguments.

We have seen above that geostrophic motions in a spherical container are restricted to zonal flows of the form $\mathbf{u}=v_{\phi}(s,t)\mathbf{e}_{\phi}$. This obviously cannot carry heat radially; thermal convection therefore necessarily involves deviations from geostrophy. This can result from buoyancy forces, inertia, viscous forces, or Lorentz force. We will not consider here the effect of the Lorentz force since we are interested in the initiation of convection in a planetary core, which here is assumed to predate the generation of the magnetic field. Furthermore, we can safely neglect inertia at the onset of convection since it is quadratic in $\mathbf{u}$ while the Coriolis and viscous forces are linear in $\mathbf{u}$. We therefore cannot rely on inertia to break the geostrophic balance, and we are left with viscous and buoyancy forces.

Let us first consider the vorticity equation (95), which, neglecting the inertia-derived terms, is

[TABLE]

Denoting by $U$ the characteristic scale of convective velocities, $\ell$ the characteristic horizontal length scale of convective motions, and $\delta T$ the characteristic scale of horizontal temperature variations, the first term is on the order of $\Omega U/H$, the baroclinic term on the order of $\alpha g\delta T/\ell$, and the viscous term on the order of $\nu U/\ell^{3}$ (the vorticity $\boldsymbol{\zeta}=\boldsymbol{\nabla}\times\mathbf{u}$ is on the order of $U/\ell$). Since buoyancy is the driver of the flow, we expect the baroclinic vorticity production to be always of importance, and balanced by either the viscous term or the Coriolis term, or both.

A balance between baroclinic production of vorticity and diffusion of vorticity gives

[TABLE]

which is basically a Stokes velocity. Balancing baroclinic production of vorticity and stretching of the planetary vorticity gives

[TABLE]

(a thermal wind balance).

The two above velocity scales are equal when the three terms in equation (101) are of equal importance, which happens when the flow length scale is on the order of $\ell_{c}\sim(\nu H/\Omega)^{1/3}=HE^{1/3}$. At $\ell\gg\ell_{c}$, the viscosity term is small compared to the planetary vorticity stretching term and the relevant velocity scaling is equation (103); at $\ell\ll\ell_{c}$, the viscosity term is large compared to the planetary vorticity stretching term and the relevant velocity scaling is equation (102). Since the velocity scale increases with $\ell$ while $\ell<\ell_{c}$ (equation (102)), and then decreases with increasing $\ell$ (equation (103)), $\ell_{c}$ also happen to be the convection length scale at which the velocity would be maximal.

One key point here is that the $2\Omega\frac{\partial\mathbf{u}}{\partial z}\sim\alpha g\boldsymbol{\nabla}T\times\mathbf{e}_{r}$ balance would yield a velocity field $\mathbf{u}$ with no radial component. This will not carry heat radially. In the absence of inertia, radial convection would therefore requires viscous effects to be important. This may seem counter-intuitive, but viscous effects are actually necessary to initiate a Rayleigh-Bénard-type convection in a rotating sphere ! This would suggest that an horizontal scale smaller than $\ell_{c}$ is necessary for the initiation of convection. Since in addition the velocity scale at $\ell<\ell_{c}$ decreases with decreasing $\ell$, we may expect that the optimal scale for the initiation of convection is on the order of $\ell_{c}$: radial motions are damped by rotation effects at larger scale, while smaller scale motions are slower due to viscous friction.

Let us now consider a liquid parcel of size $\ell$ and temperature $T$ larger than the surrounding temperature $\bar{T}(r)$ (temperature excess $\delta T=T-\bar{T}$). $\bar{T}(r)$ is the background conductive temperature profile. The radial velocity of the parcel is given by

[TABLE]

since, as argued above, viscous effects are necessary to the initiation of convection in a rotating sphere. To see how the velocity of the parcel evolves, take its Lagrangian derivative:

[TABLE]

The first term in the parenthesis can be estimated from the heat equation,

[TABLE]

while the second term is

[TABLE]

Putting everything together, this gives

[TABLE]

This can be rearranged to give an estimate of the growth rate of the vertical velocity:

[TABLE]

With $\frac{d\bar{T}}{dr}\sim-\Delta T/H$, this gives

[TABLE]

A parcel displaced upward (resp. downward) will keep rising (resp. sinking) if $\frac{1}{u_{r}}\frac{Du_{r}}{Dt}>0$, i.e. if the term within the parenthesis in equation (111) is positive. This condition can be recast as a condition for the Rayleigh number, which must exceed a critical value $\mathrm{Ra}_{\mathrm{cr}}$ which is a function of the length scale $\ell$ of the perturbation:

[TABLE]

Perturbations with the largest length scale $\ell$ will have the lowest critical Rayleigh number and will then be favoured. In non-rotating convection, the only limit to $\ell$ is the size $H$ of the convecting layer, so the lowest $\mathrm{Ra}_{\mathrm{cr}}$ will correspond to perturbations with $\ell\sim H$. The initiation of convection then simply requires Ra to exceed some critical value, which depends only on the geometry and boundary conditions (it is for example $27\pi^{4}/4\simeq 657.5$ for Rayleigh-Bénard convection in a plane layer with free-slip boundaries). In a rotating sphere with radial gravity, driving radial motion requires initiating the convection with a horizontal length scale $\lesssim\ell_{c}=HE^{1/3}$. Since the critical Rayleigh number is a decreasing function of $\ell$, we expect the length scale of the fastest growing mode to be $\sim\ell_{c}$, which gives

[TABLE]

as predicted by linear stability analysis.

7.4 Compressibility effect

In the above analysis, we have left aside compressibility effects on the initiation of convection. Rather than the Boussinesq version of the heat transfer equation (85), we now consider its more general form

[TABLE]

where $\dot{\epsilon}:\tau$ is the viscous dissipation ($\tau$ is the stress tensor and $\dot{\epsilon}$ the deformation rate tensor). Equation (114) can be obtained from the entropy balance (for its derivation, see e.g. Ricard, 2007).

7.4.1 Schwarzschild’s stability criterion for thermal convection

To see the effect of the pressure term, let us consider a parcel of fluid displaced upward with negligible viscous dissipation and heat diffusion (thus following an isentropic path). It follows directly from equation (114) that the parcel’s temperature will vary with its pressure as

[TABLE]

which corresponds to adiabatic heating (cooling) due to compression (decompression).777The evolution of the parcel is actually even an isentropic process, since the parcel’s evolution is both adiabatic (no transfer of heat and mass with its surrounding) and reversible (no friction). This is equivalent to a radial gradient, the adiabatic gradient, given by

[TABLE]

The adiabatic temperature difference $\Delta T_{\mathrm{ad}}$ across the outer core is obtained by integrating one of equations (115) or (116) (here assuming that $\alpha/c_{p}$ is constant and that $g(r)=g_{\mathrm{cmb}}r/r_{\mathrm{cmb}}$):

[TABLE]

where $T_{\mathrm{icb}}$ is the temperature at the Inner Core Boundary. In Earth’s core, the adiabatic temperature difference is $\Delta T_{\mathrm{ad}}\sim 1000$ K, so this is a significant effect.

If displaced fast enough for negligible heat diffusion, a parcel of fluid displaced vertically from its initial state will follow an adiabatic path. Its temperature $T$ will then vary with pressure $P$ according to equation (115). If the background temperature gradient $\bar{T}$ is steeper than the adiabatic gradient, i.e. if

[TABLE]

then the temperature of a parcel of fluid displaced upward will becomes larger than the background temperature, with the temperature difference increasing with the upward displacement of the parcel (figure 10a). The parcel will thus become less and less dense than the surrounding, and will keep rising. In contrast, if

[TABLE]

then a parcel of fluid displaced upward will have a temperature lower than the background temperature (figure 10b). It will thus be denser than the surrounding fluid, and its buoyancy will eventually send it back toward its initial position. A temperature profile satisfying the condition (120) is thus stable against thermal convection.

Equation (119) is thus a necessary condition for thermal convection, known as Schwarzschild’s stability criterion888Alternatively, equation (120) is a sufficient condition for stability.. It is not a sufficient condition since it assumes that fluid parcels follow adiabatic and reversible paths when displaced upward or downward. In a real fluid, the parcels would exchange heat with the surrounding by diffusion, and viscous dissipation would be non-zero.

7.4.2 Critical Rayleigh number with compressibility effects

Going back to our heuristic derivation of the critical Rayleigh number, we consider again the evolution of a fluid parcel displaced from its initial position, but now assume that the evolution of its temperature is governed by

[TABLE]

ignoring the viscous dissipation term $\dot{\epsilon}:\tau$ included in equation (114). Viscous dissipation is quadratic in $\mathbf{u}$, and can therefore be neglected when looking at the initiation of convection.

In equation (121), the pressure derivative can be estimated as

[TABLE]

neglecting terms in pressure fluctuations. Equation (121) thus gives

[TABLE]

Using this relation instead of equation (106) and going through the same steps (equations (104) to (113)), the rate of change of the parcel’s velocity is now given by

[TABLE]

instead of equation (110). Denoting by $\Delta T_{\mathrm{>ad}}=\Delta T-\Delta T_{\mathrm{ad}}$ the excess temperature difference above the adiabatic temperature across the outer core (or super-adiabatic temperature difference), we have $\frac{d(\bar{T}-T_{\mathrm{ad}})}{dr}\sim-\Delta T_{\mathrm{>ad}}/H$, and

[TABLE]

where

[TABLE]

is the super-adiabatic Rayleigh number. The instability criterion is thus generalised by simply replacing the classical Rayleigh number by the super-adiabatic Rayleigh number. In a compressible fluid at low Ekman, convection will thus start if

[TABLE]

7.5 Application to Earth’s core

7.5.1 Thermal convection

How much superadiabatic the core needs to be to meet this requirement? Using the definitions of $Ra_{\mathrm{>ad}}$ and $E$, the condition (127) can be recast as a condition for the superadiabatic temperature difference:

[TABLE]

With parameter values from table 1, the required $\Delta T_{\mathrm{>ad}}$ is only on the order of $10^{-7}-10^{-6}$ K: this is tiny, compared for example to the adiabatic temperature variation across the core, which is $\sim 1000$ K.

In practice this means that Schwarzschild’s criterion is a very good indicator of the likelihood of core convection. In a planetary core, the temperature difference across the core is not imposed. What is imposed is the heat flux $Q_{\mathrm{cmb}}$ at the Core mantle Boundary, which is controlled by mantle convection. A requirement for thermal convection is therefore that the heat flux imposed by mantle convection is larger than the heat flux $Q_{\mathrm{cmb}}^{\mathrm{ad}}$ which can be carried by diffusion along an adiabatic temperature profile:

[TABLE]

where $k$ is the thermal conductivity of the core. The value of the thermal conductivity of the core is currently quite debated. Until quite recently, available estimates were in the range $30-60$ W.m*-1*.K*-1* (Stacey and Anderson, 2001; Stacey and Loper, 2007), but much higher values ($>150$ W.m*-1*.K*-1*) have since been proposed by several independent groups (de Koker et al., 2012; Pozzo et al., 2012; Gomi et al., 2013, 2016). Other groups (Seagle et al., 2013; Konôpková et al., 2016) still favour a relatively low value of $k$.

Finally, it is also instructive to estimate the horizontal scale $\ell_{c}$ at which convection is initiated. From equation (100), this is $\ell_{c}\sim HE^{1/3}\sim 2000\ \mathrm{km}\times(10^{-15})^{1/3}\sim 20$ m, a factor $10^{5}$ smaller than the outer core shell thickness!

7.5.2 Compositional convection

Though equations (83)-(85) have been written for thermal convection, the same set of equations can be used to describe compositional convection if $T$ is replaced by a concentration (in light elements) $\chi$, $\alpha$ by a coefficient of compositional expansion $\beta$, and $\kappa$ by a compositional diffusivity $\kappa_{\chi}$. There is no direct compositional analog of the adiabatic gradient, so the condition for compositional convection is simply

[TABLE]

where $\Delta\chi$ is the light element concentration difference across the outer core. Written in terms of $\Delta\chi$, condition (131) is equivalent to the following condition:

[TABLE]

A tiny difference in composition can drive core convection.

8 Convective dynamos

8.1 Governing equations

In its simplest form, the set of governing equations for a convectively driven dynamo consists of the three equations governing rotating convection with the addition of the Lorentz force in the Navier-Stokes equation, plus the induction equation:

[TABLE]

This set of equations can be made dimensionless by using the characteristic scales already used in section 7 in the case of thermal convection for lengths ($H$), time ($H^{2}/\nu$), temperature ($\Delta T$ or $qH/k$) and pressure ($\rho\Omega\nu$). A magnetic field scale can for example be obtained by assuming a balance between the Coriolis acceleration and the Lorentz force. Writing the Lorentz force as $\mathbf{j}\times\mathbf{B}$ and estimating $\mathbf{j}$ from Ohm’s law as $\mathbf{j}\sim\sigma\mathbf{u}\times\mathbf{B}$ gives a Lorentz force $\sim\sigma UB^{2}$. Balancing this with the Coriolis force which is $\sim\rho\Omega U$ gives a magnetic field scale $B_{0}=\sqrt{\rho\Omega/\sigma}$. Using this set of scales, we obtain the following dimensionless set of equations:

[TABLE]

where the magnetic Prandtl number $P_{m}$ is defined as the ratio of the kinematic viscosity $\nu$ to the magnetic diffusivity:

[TABLE]

$P_{m}$ is about $10^{-6}$ in Earth’s core: the magnetic field diffuses much faster than momentum, and we therefore expect the magnetic field to vary over larger length scales that the velocity field.

In addition to the input non-dimensional numbers ($E$, $Ra$, $Pr$, $P_{m}$), it is also often useful to consider output non-dimensional numbers based on measured dynamical quantities such as some averages of the velocity and magnetic field, $\langle u\rangle$ and $\langle B\rangle$. Examples of useful output non-dimensional numbers include the followings:

[TABLE]

In this expression of the Elsasser number, the Lorentz force has been estimated from the expression $\mathbf{j}\times\mathbf{B}$ by taking $\mathbf{j}\sim\sigma\mathbf{u}\times\mathbf{B}$ from Ohm’s law.

These numbers can be estimated for Earth’s core as follows. If we accept that the geomagnetic field is produced by dynamo effect in the core, then the magnetic Reynolds number must be larger that $\mathcal{O}(10)$ (the critical magnetic Reynolds number for dynamo action in a spherical shell is typically $\sim 50$). This implies that $Re$ must be at least of $\mathcal{O}(10^{7})$ and Rossby above $\mathcal{O}(10^{-8})$. The order of magnitude of the velocity corresponding to $Rm\sim 10$ is $10\eta/H\sim 5\times 10^{-6}$ m.s*-1*. If instead we take a velocity scale of $5\times 10^{-4}$ m.s*-1* as obtained from core flow inversions, we obtain $R_{m}\sim 10^{3}$, $Re\sim 10^{9}$, and $Ro\sim 10^{-6}$. With a magnetic field of $\sim 4$ mT (Gillet et al., 2010), the Elsasser number is $\Lambda\sim 10$.

The values of these numbers suggest that in the Navier-Stokes equation the dominant forces would be the Coriolis force, the Lorentz force, and presumably the buoyancy force (which is difficult to estimate from simple arguments, but which is likely non-negligible since it is the source of motion). This corresponds to the so-called MAC balance (Magnetic, Archimedes, Coriolis).

Solving numerically these equations in regimes which are relevant to the geodynamo is difficult; it has in fact not yet been possible to solve them with parameter values approaching that of the Earth. They are two main reasons for this:

The low viscosity of molten iron means that the velocity field likely develops small scale turbulent fluctuations, which would require a high spatial resolution and fine time-stepping to be fully resolved. We have seen in section 7 that thermal convection would initiate at a length scale on the order of $HE^{1/3}$, or about $20$ m for $E\sim 10^{-15}$. This is $10^{5}$ smaller than the outer core thickness. Resolving this scale in 3D numerical simulations would necessitate at least $10^{5}$ grid points in each direction, or $10^{15}$ grid points in total. Today’s most resolved numerical simulations of the geodynamo have a spatial resolution of about 2 km, which allows to reach $E=10^{-7}$ (Schaeffer et al., 2017). 2. 2.

The magnetic Prandtl number is quite small in liquid metals (probably $\sim 10^{-6}$ in Earth’s core), but much larger values are being used in dynamo simulations. Solving the induction equation with a low $P_{m}$ present no intrinsic difficulty: a high magnetic diffusivity usually means smooth variations of the magnetic field and small magnetic fields gradients. What is difficult is to obtain dynamo action at a low $P_{m}$ in numerical simulations. The reason for this is easily understood by noting that the magnetic Reynolds number can be written as the product of the Reynolds number and magnetic Prandtl number: $Rm=Re\times P_{m}$. Since dynamo operation requires reaching $Rm>\mathcal{O}(10)$, doing this with a low $P_{m}$ requires a high value of $Re$. This implies the development of turbulent velocity fluctuations down to small length scales, which are difficult to resolve numerically.

8.2 Successes and challenges

Though numerical dynamos are still quite far from Earth’s conditions in terms of non-dimensional parameters, this does not mean that relevant numerical simulations cannot be done. In the past few decades, the geodynamo modelling community has been quite successful in strengthening the case for a convectively powered geodynamo.

As already discussed in section 6, a major achievement has been to demonstrate that dynamos can indeed be sustained by rotating convection (Zhang and Busse, 1988; Glatzmaier and Roberts, 1995; Kageyama et al., 1995). Since the first few numerical models of the dynamo (which were done at relatively high $E$ and $P_{m}$), the non-dimensional parameter space has been explored toward lower $E$ and $P_{m}$, and higher $Ra$. By performing geodynamo simulations at various values of the governing non-dimensional parameters, scaling laws for typical velocity and magnetic field strength have been obtained. When extrapolated to the core condition, the predictions of these scaling laws are in reasonable agreement with estimates of the core flow (from inversion of secular variation of the CMB magnetic field) and observed field strength (e.g. Starchenko and Jones, 2002; Olson and Christensen, 2006; Christensen and Aubert, 2006; Christensen, 2010). In addition, it has been shown that the magnetic field can be “Earth like” in a well-defined region of the parameter space, which includes the estimated state of Earth’s core (Christensen et al., 2010). “Earth-like” geodynamo simulations also often exhibit polarity reversals.

In spite of the important successes described above, the question of whether numerical simulations do reach a dynamical regime similar to planetary cores has been a constant preoccupation of the geodynamo community. In other words, are geodynamo simulations “Earth-like” for good reasons, or is it just coincidental? Numerical simulations are still far from the Earth in terms of non-dimensional numbers: the most massive numerical simulations of the geodynamo have reached $E=10^{-7}$, $P_{m}=0.1$, with a Rayleigh number $6\times 10^{3}$ times supercritical (Schaeffer et al., 2017). Though quite impressive, this is still far from the Earth (seven or eight orders of magnitude in terms of $E$, five orders of magnitude in $P_{m}$). And most of available geodynamo simulations are significantly further away from Earth’s core parameters values. The main question is whether viscous effects in state of the art simulations play a significant role or not. Christensen and Aubert (2006) have argued, from a large set of numerical dynamos ($E\in 10^{-6}-3\times 10^{-4}$, $P_{m}\in 0.06-10$), that the magnetic field strength and mean velocity follow scaling laws which are independent of diffusivities. This may suggest that an asymptotic regime has been reached, and that further decreasing the viscosity down to Earth’s value would not change the dominant force balance, allowing to use the obtained scaling laws to extrapolate outputs of numerical simulations to planetary cores conditions. However, several authors have questioned whether viscosity effects in current geodynamo simulations are indeed small enough, and whether diffusivity-free scaling laws correctly describe the available set of numerical geodynamos (Stelzer and Jackson, 2013; King and Buffett, 2013; Cheng and Aurnou, 2016). The ratio of magnetic to kinetic energies, which is estimated to be $\sim 10^{3}-10^{4}$ in Earth’s core, is much smaller in numerical models at low ($\lesssim 1$) $P_{m}$: $\mathcal{O}(10)$ at best, and often below 1 (Schaeffer et al., 2017).

Chapter 5 Energetics of the geodynamo

Figure 11 shows schematically the energy flow of a dynamo powered by natural convection (thermal or compositional) or astronomical forcing (tides, precession, nutation, or libration). Natural convection converts the gravitational potential energy of unstable density gradients into kinetic energy; inertial instabilities excited by astronomical forcing converts rotational energy into kinetic energy. If the magnetic Reynolds number is large enough for dynamo action, a fraction of this kinetic energy is converted into magnetic energy through induction. Magnetic energy can be converted back into kinetic energy (through the work of the Lorentz force). Kinetic and magnetic energies are dissipated into heat through viscous and ohmic dissipation.

Dynamo action requires the rate of release of gravitational or rotational energies to be high enough to maintain the magnetic field against viscous and ohmic dissipation. The two energy reservoirs are huge: the gravitational potential energy of the Earth is $\simeq 2\times 10^{32}$ J, and the rotational energy is $\simeq 2\times 10^{29}$ J. Only a fraction of these would be enough to maintain Earth’s magnetic field for the last 4 Gy if the amount of dissipation is in the range 0.1-3 TW as discussed later in this section. However, the fact that enough energy is available in Earth’s system is by no means a sufficient condition for sustaining the geomagnetic field. Only a small fraction of this energy can actually be made available to sustain the geomagnetic field; the goal of the energetics approach of the geodynamo is to estimate this fraction.

An equation for the evolution of the magnetic energy $E_{m}=\int_{\mathcal{V}_{\infty}}\frac{B^{2}}{2\mu_{0}}d{\mathcal{V}}$ can be obtained from Ohm’s law and Maxwell’s equations, or equivalently from the dot product of the induction equation with $\mathbf{B}$. Integration over the whole space $\mathcal{V}_{\infty}$ gives

[TABLE]

where $\mathbf{f}_{L}$ is the Lorentz force, and $\mathcal{V}$ the volume of the core. In this equation, the magnetic energy includes the contribution of the magnetic field outside the core. On the other hand, the integrals of the power of the Lorentz force and of ohmic dissipation can be restricted to the core volume $\mathcal{V}$ under the assumption that the electric conductivity is equal to 0 outside the core. This equation cannot be used alone to predict whether the geodynamo can be dynamically sustained, because the (unknown) velocity field $\mathbf{u}$ appears in the equation. We therefore need a second conservation equation involving $\mathbf{u}$, which is obtained by taking the dot product of Navier-Stokes equation with $\mathbf{u}$, and again integrating over the volume of the core. This gives an equation for the time derivative of the total kinetic energy of the core:

[TABLE]

where $\mathcal{S}$ is the surface area of the core-mantle boundary, $\underline{\tau}$ is the deviatoric stress tensor and $\Phi_{\nu}=\int_{\mathcal{V}}\underline{\tau}:\nabla\mathbf{u}\,d\mathcal{V}$ is the viscous dissipation. There is no power associated with the Coriolis force since $(\rho\Omega\times\mathbf{u})\cdot\mathbf{u}=0$.

Combining the kinetic and magnetic energy equations, we can get rid of the power of the Lorentz force and obtain:

[TABLE]

In this equation:

The term $\int_{\mathcal{V}}\left(\rho\mathbf{g}-\boldsymbol{\nabla}P\right)\cdot\mathbf{u}\,d{\mathcal{V}}$ is the rate of work of non-hydrostatic processes. It is equal to 0 if the core is in a hydrostatic state. In planetary cores, deviations from a hydrostatic state can be due to density, pressure, or gravity field fluctuations, which can arise from either natural convection or astronomical forcing. The integral is positive in the case of natural convection; it can be either positive or negative in the case of astronomical forcing. 2. 2.

The term $\oint_{\mathcal{S}}\mathbf{(}\underline{\tau}\cdot\mathbf{u})\cdot\mathbf{n}\,d{\mathcal{S}}$ is the rate of work of deviatoric stresses (viscous stresses) at the CMB. It is different from zero if the mantle rotates at a different rate than the core (in particular due to precession/nutation/libration) or if the CMB is deformed (tides). 3. 3.

$\Phi$ is the total dissipation due to Joule heating (Ohmic dissipation) and viscous heating (viscous dissipation). It is always positive.

Starting from a state at rest with $E_{k}=E_{m}=0$, the sum of the kinetic and magnetic energy would increase if the sum of the two first terms on the RHS of equation (148), which we denote by $P_{a}$, is positive. Increasing $\mathbf{u}$ and $\mathbf{B}$ will increase the dissipation $\Phi$ up to a point of saturation at which the rate of energy production is on average balanced by dissipation, i.e. $P_{a}\simeq\Phi$. In a statistically steady state, time averaging equation (148) gives

[TABLE]

From an energetic point of view, the geodynamo problem amounts to estimate whether the rate of energy energy release $P_{a}$ is large enough to sustain Earth’s magnetic field against dissipation. To answer this question, we need to estimate: (i) how dissipative is the geodynamo, and (ii) what is the available power $P_{a}$ in the core.

8.3 How dissipative is the geodynamo?

The amount of ohmic dissipation in the core is given by

[TABLE]

To estimate how dissipative is the geodynamo, we thus need to know how $\mathbf{B}$ varies spatially in the core. This is difficult. As discussed in section 3, direct observations of the core magnetic field are restricted to the poloidal field at the core-mantle boundary, up to spherical harmonic degree $l_{\mathrm{max}}=13$, which corresponds to a spatial wavelength of about 1600 km. We have no direct observation of the harmonic components of the poloidal field at $l>13$, and no constraint on the spatial variations of the toroidal part of the magnetic field. Numerical simulations of the geodynamo suggest that ohmic dissipation is dominated by the contributions of smaller scale components of the magnetic field (Roberts et al., 2003), so estimating the geodynamo dissipation from direct observations seems hopeless. Published estimates of $\Phi_{J}$ have been obtained from extrapolations of numerical or experimental dynamos results (Buffett, 2002; Roberts et al., 2003; Christensen and Tilgner, 2004; Christensen et al., 2010; Stelzer and Jackson, 2013), and range from $\sim 0.1$ TW to several TW.

8.4 Estimating the available power

The energetics approach of the geodynamo problem has so far been restricted to the case of a convectively driven dynamo; the case of mechanical forcing has not been treated in a rigorous and usable way, in part because of the difficulty of estimating the $\int_{\mathcal{V}}\left(\rho\mathbf{g}-\boldsymbol{\nabla}P\right)\cdot\mathbf{u}\,d{\mathcal{V}}$ term. We will thus restrict ourselves to the case of a convectively driven dynamo. In this situation the surface integral in equation (148) is equal to 0.

Two different methods have been classically used: either work on the entropy budget rather than on equation (148) (e.g. Gubbins et al., 2004; Labrosse, 2015), or estimate directly the contribution of convective motions on $P_{a}$. These two approaches are equivalent (Lister, 2003) and in both cases require assumptions on the state of the core, which is usually assumed to be continuously well mixed by convective motions. Here we use the later approach, which is perhaps a bit more physically intuitive, and follow the approach of Buffett et al. (1996). We will only sketch the derivation of the model; the details can be found in Buffett et al. (1996).

Denoting by $\psi$ the gravitational potential, with $\mathbf{g}=-\boldsymbol{\nabla}\psi$, one can show by using mass conservation that

[TABLE]

where $\partial\rho/\partial t$ is the rate of change of the density due solely to convective re-arrangement. The next step is to calculate the change of density in the core associated to the convective re-distribution of heat and chemical elements, under the assumption that convection keeps the core well mixed and isentropic. By doing this, one can write $P_{a}$ as

[TABLE]

where $Q_{\mathrm{icb}}^{\mathrm{>ad}}=Q_{\mathrm{icb}}-Q_{\mathrm{icb}}^{\mathrm{ad}}$ and $Q_{\mathrm{cmb}}^{\mathrm{>ad}}=Q_{\mathrm{cmb}}-Q_{\mathrm{cmb}}^{\mathrm{ad}}$ are the super-adiabatic heat flux at the inner core boundary (ICB) and core-mantle boundary (CMB), ${F}_{\mathrm{icb}}=4\pi r_{\mathrm{icb}}^{2}\dot{r}_{\mathrm{icb}}c$ is the flux of light elements concentration (in wt.%.m3.s*-1*) at the ICB ($c$ being the concentration in light elements of the outer core), $\bar{\psi}$ is the mass-averaged value of the gravitational potential in the outer core, and $\psi_{\mathrm{cmb}}$ and $\psi_{\mathrm{icb}}$ are the values of the gravitational potential at the CMB and ICB. The factors $\bar{\psi}-\psi_{\mathrm{cmb}}$ and $\psi_{\mathrm{icb}}-\bar{\psi}$ come from the fact that density perturbations originating from either the CMB or ICB are redistributed over the whole core. The superadiabatic heat flux at the ICB happens to be well approximated by the release of latent heat due to inner core solidification, $Q_{\mathrm{icb}}^{\mathrm{>ad}}\simeq 4\pi r_{\mathrm{icb}}^{2}\rho L\dot{r}_{\mathrm{icb}}$, where $\dot{r}_{\mathrm{icb}}$ is the time-derivative of the inner core radius $r_{\mathrm{icb}}$, and $L$ is the latent heat of iron at core conditions. Since ${F}_{\mathrm{icb}}$ is also proportional to $\dot{r}_{\mathrm{icb}}$, the whole ICB contribution is proportional to the rate of growth of the inner core, which is itself controlled by the rate at which heat is extracted from the core.

The terms in equation (152) are all tied to the core thermal evolution and controlled by the heat flux at the CMB. To use equation (152) in a predictive way, one needs to relate the rate of inner core growth to the heat flux at the CMB, which can be done by considering the energy balance of the core (coming from the first law of thermodynamics), which simply states that the heat flux extracted from the core at the CMB is equal to the sum of sensible heat released by the cooling of the core, latent heat of inner core crystallisation, and compositional energy (associated with the mixing in the core of light elements released by inner core solidification). Doing this allows to write equation (152) as

[TABLE]

where $\epsilon_{\mathrm{icb}}$ and $\epsilon_{\mathrm{cmb}}$ are the so-called dynamo efficiencies.999The CMB efficiency is usually defined in another way, such that $P_{a}=\left(\epsilon_{\mathrm{icb}}+\epsilon_{\mathrm{cmb}}\right)Q_{\mathrm{cmb}}$. This has the advantage of linking $P_{a}$ directly to $Q_{\mathrm{cmb}}$, but then $\epsilon_{\mathrm{cmb}}$ is itself a function of $Q_{\mathrm{cmb}}$ (it includes a factor $Q_{\mathrm{cmb}}^{>\mathrm{ad}}/Q_{\mathrm{cmb}}$). Expressions for these efficiencies can be found in Buffett et al. (1996) or Lister (2003). One important point is that $\epsilon_{\mathrm{icb}}$ is an increasing function of inner core size (figure 12), which means that the contribution of compositional convection increases with inner core size. An approximate expression of $\epsilon_{\mathrm{cmb}}$ is

[TABLE]

Having obtained estimates of $\epsilon_{\mathrm{icb}}$ and $\epsilon_{\mathrm{cmb}}$, it is now possible to estimate the convective power $P_{a}$. One key parameter is the heat flux conducted along the adiabat, $Q_{\mathrm{cmb}}^{\mathrm{ad}}$, which depends mostly on the thermal conductivity of the outer core (see equation (130)). Figures 13a and b show $P_{a}$ as a function of the CMB heat flux and inner core size, for $Q_{\mathrm{cmb}}^{\mathrm{ad}}$ equal to either 5 TW or 15 TW. The 5 TW and 15 TW values are representative of the low ($\sim 30$ W.m*-1*.K*-1*) and high ($\sim 100$ W.m*-1*.K*-1*) thermal conductivity estimates. Figures 13c and d show the evolution with time of inner core size $r_{\mathrm{icb}}$ and convective power $P_{a}$, for $Q_{\mathrm{cmb}}^{\mathrm{ad}}$ equal to either 5 TW or 15 TW, and imposed CMB heat flux $Q_{\mathrm{cmb}}$ equal to 5, 10, or 15 TW, the current CMB heat flux being thought to be in the range 5-15 TW (Lay et al., 2006). On figures 13a and b, the grey line corresponds to the limit above which $P_{a}>0$.

One can see that a positive $P_{a}$ requires the CMB heat flux to be above a threshold value, which is equal to $Q_{\mathrm{cmb}}^{\mathrm{ad}}$ at $r_{\mathrm{icb}}=0$, and decreases with increasing $r_{\mathrm{icb}}$ due to the increasing importance of the ICB buoyancy flux. From equation (153), a dynamo dissipating at a given rate $\Phi_{0}$ requires a CMB heat flux

[TABLE]

The CMB heat flux above which $P_{a}>0$ is given by setting $\Phi_{0}=0$ in equation (155). On figures 13a and b are also shown the $Q_{\mathrm{cmb}}$ values corresponding to $\Phi_{0}=0.1$ TW and 1 TW.

One can see that in the current state of the core, the convective power can easily be of a few terawatts if the CMB heat flux is in the range 5-15 TW as currently believed (Lay et al., 2006). There is therefore no difficulty in powering the geodynamo with the inner core at its current size. However, powering the geodynamo happens to be more problematic when the inner core is smaller, and before its nucleation. With a low core thermal conductivity and $Q_{\mathrm{cmb}}^{\mathrm{ad}}=5$ TW, the convective power $P_{a}$ would be of a few tens of TW. Driving the geodynamo with thermal convection seems therefore possible, though the convective power would be much lower than at present. Whether this would have a significant impact on the large scale part of the geomagnetic field remains an open question. If the high estimates of the core thermal conductivity are correct, the CMB heat flux may well be sub-adiabatic. In this situation thermal convection is not possible, and driving the dynamo requires another source of motion and energy. Possible additional sources of energy include exsolution of light elements from the core (O’Rourke and Stevenson, 2016; Badro et al., 2016; O’Rourke et al., 2017) and astronomical forcing (Andrault et al., 2016). In addition, adding some radioactive heating in the core (possible due to 40K) would help by decreasing the rate of cooling of the core and increasing the age of the inner core (Labrosse, 2015).

Chapter 6 Inner core dynamics

The Earth’s inner core is the deepest layer of our planet: a 1221-km-radius sphere of solid iron-alloy surrounded by molten metal. Its existence was unknown until the first observations of seismic reflexions at the inner core boundary by Inge Lehmann, in one of the shortest-title paper ever: P’ (Lehmann, 1936). The arrivals of P-waves in the core shadow zone, where P-waves are refracted away by the presence of the 3600-km-radius core, have been explained by the existence of a new discontinuity inside the core, the inner core boundary.

Since these first observations, the study of seismic waves travelling through the inner core and normal modes sampling the deepest layers have provided a blurry image of the inner core structure (figure 14). Birch (1940) and Jacobs (1953) have proposed that the inner core is a solid sphere of the same metal constituting the outer core, while the actual solidity has been demonstrated only several years later by Dziewonski and Gilbert (1971). Poupinet et al. (1983) were the first to note the different propagation velocities of waves travelling parallel to the rotation axis and perpendicular to it. This anisotropy of seismic properties has since been extensively studied, demonstrating the existence of a complex structure of the inner core. Among the most robust and surprising features of the inner core structure, we can cite two: a strong anisotropy for the bulk of the inner core, and an uppermost layer of the inner core with a strong hemispherical dichotomy in P-waves velocity, but no detectable anisotropy.

The existence of anisotropy in the inner core is an evidence for crystal orientation within the bulk of the inner core (iron crystals being elastically anisotropic). Such crystal orientation is the main motivation for studying flows in the inner core, as lattice-preferred orientation (LPO) may be deformation-induced (Karato, 2012). Thus, it is thought that the observed structure may be an evidence for flows within the inner core, likely to be combined with initial crystallisation-induced LPO.

We will discuss here the constitutive equations and some of the aspects of inner core dynamics, focusing on the flow forced in the inner core by the magnetic field diffused from the outer core.

9 Constitutive equations

We consider an incompressible fluid in a spherical domain of radius $r_{\text{icb}}$, with a newtonian rheology and uniform viscosity $\eta$. Neglecting inertia and rotation, the equation of continuity and the conservation of momentum are written as

[TABLE]

where $\mathbf{F}$ denotes volume forces, $p$ is the dynamic pressure (including gravity potential) and $\mathbf{u}$ the velocity field.

Among the volume forces, the buoyancy force is written as $\mathbf{F}_{\text{buoyancy}}=\Delta\rho\mathbf{g}$, where $\Delta\rho$ is the density difference compared to a neutral density profile and $\mathbf{g}=g_{\text{icb}}r/r_{\text{icb}}\mathbf{e}_{r}$ the acceleration of gravity. $g_{\text{icb}}$ is the acceleration of gravity at the surface of the inner core.

In this model, the density variations only play a role in the buoyancy term, and they can be related to variations in temperature or concentration in light elements (Si, S, O, …) compared to a reference state. We take as the reference temperature profile the adiabatic temperature profile $T_{\text{ad}}(r)$ anchored at the melting temperature at the radius of the inner core $r_{\text{icb}}$. The deviation of the temperature field compared to this reference profile is the potential temperature $\Theta=T-T_{\text{ad}}$. The potential composition field is defined as $C=c^{s}-c^{s}_{\text{icb}}$, where $c^{s}$ is the concentration of light elements in the inner core and $c^{s}_{\text{icb}}$ its value at the inner core boundary. The density variations are thus $\rho\alpha_{T}\Theta$ or $\rho\alpha_{C}C$ for respectively thermal or compositional stratification, where $\alpha_{T}$ and $\alpha_{C}$ are the thermal and compositional expansion coefficients. As both potential temperature and potential composition are solutions of an advection-diffusion equation and are both set to zero by construction at the inner core boundary, we will consider a general equation for a quantity $\chi$, representing either of these quantities. $\chi$ is solution of

[TABLE]

where $\kappa_{\chi}$ is the diffusivity and $S\chi(t)$ a source term built from the evolution of the reference profile as

[TABLE]

The continuity and momentum equations can be solved using a poloidal-toroidal decomposition of the velocity field $\mathbf{u}=\nabla\times(T\mathbf{r})+\nabla\times\nabla\times(P\mathbf{r})$, where $\mathbf{r}=r\mathbf{e}_{r}$ is the position vector and $T$ and $P$ respectively the toroidal and poloidal components. In the following, we will only consider boundary conditions with a zero vertical vorticity and volume forces without toroidal components. The flow is thus expected to have only non-zero poloidal component, and applying $\mathbf{r}\cdot(\nabla\times\nabla\times)$ to equation (157), we obtain

[TABLE]

where $L^{2}$ is the Laplace horizontal operator defined as

[TABLE]

One could note that for a volume force of the form $\mathbf{F}=F_{r}r\mathbf{e}_{r}$ such as the buoyancy forces, the second term of the right-hand side of the equation simplifies as $\mathbf{r}\cdot(\nabla\times\nabla\times\mathbf{F})=L^{2}F_{r}$. Splitting the volume force term as one term for the buoyancy forces and one term for the other forces, we have $\mathbf{F}=\rho\alpha_{\chi}\chi rg_{\text{icb}}\mathbf{e}_{r}/r_{\text{icb}}+\mathbf{F}_{\text{volume}}$. Expanding the scalar fields $P$ and $\chi$ with horizontal spherical harmonics $Y_{l}^{m}$ satisfying $L^{2}Y_{l}^{m}=-l(l+1)Y_{l}^{m}$ of degree $l$ and order $m$, as $P=\sum_{l,m}P_{l}^{m}Y_{l}^{m}$ and $\chi=\sum_{l,m}\chi_{l}^{m}Y_{l}^{m}$, the equation of interest is eventually written for each $(l,m)$ as

[TABLE]

where $f_{l}^{m}$ is the spherical harmonics decomposition of the poloidal component of the volume force $\mathbf{F}$, and $\text{D}_{l}$ a second-order differential operator defined as

[TABLE]

9.1 Boundary conditions

The inner core boundary is a crystallisation front, where the iron-alloy of the outer core freezes due to the slow secular cooling of our planet. Its exact position is determined by the intersection of the melting temperature profile of the iron-alloy and the temperature profile in the core. Any solid material pushed dynamically further from this intersection would melt, while any liquid pushed inward would freeze. The timescale $\tau_{\phi}$ involved in the freezing or melting of a small topography at the inner core boundary can be estimated from the timescale needed by outer core convection to extract the latent heat released by crystallisation. We note $h$ the topography at the inner core boundary.

From continuity of stress at the ICB, the mechanical boundary conditions are written at $r=r_{\text{icb}}(t)$. We consider that the dynamical topography $h$ is small (compared to the horizontal wavelength) and that the vector normal to the boundary is close to the radial unit vector. The tangential and normal components of the stress tensor are written as

[TABLE]

where $p$ here is the total pressure. The viscosity of the outer core being much smaller than the viscosity of the inner core, we can assume tangential stress free conditions. The boundary conditions are then written as $\tau_{r\theta}(r=r_{\text{icb}})=\tau_{r\phi}(r=r_{\text{icb}})=0$ and continuity of $\tau_{rr}$ across the ICB, which for a small topography amounts to state that the normal stress on the inner core side at $r=r_{\text{icb}}$ is equilibrated by the weight of the topography:

[TABLE]

$p^{\prime}$ being the dynamical pressure on the inner core side of the ICB.

To close the system of equations, we consider the time evolution of the topography. The topography can be formed by deformation of the inner core boundary by the underlying flow $u_{r}-\dot{r}_{\text{icb}}$ and is eroded by phase change, such that we can write that $\text{D}h/\text{D}t=u_{r}-\dot{r}_{\text{icb}}+V_{r}$, where $V_{r}$ is the velocity of phase change in the radial direction. $V_{r}$ at first order is $-h/\tau_{\phi}$, where $\tau_{\phi}$ is a typical timescale for the phase change. Considering a dynamical equilibrium for the topography, we obtain $u_{r}-\dot{r}_{\text{icb}}=h/\tau_{\phi}$ and the continuity of normal stress is written as

[TABLE]

9.2 Non-dimensionalisation and final set of equations

The governing equations are made dimensionless using characteristic scales for time, length, velocity, pressure and $\chi$ (potential temperature or composition) as, respectively, the diffusion time $r_{\text{icb}}^{2}/(6\kappa_{\chi})$, its radius $r_{\text{icb}}$, $\kappa_{\chi}/r_{\text{icb}}$, $\eta\kappa_{\chi}/r^{2}_{\text{icb}}$ and $S_{\chi}r_{\text{icb}}^{2}/(6\kappa)$. Using the same symbols for dimensionless quantities, the non-dimensional set of equations is

[TABLE]

The dimensionless number $Ra$ is a Rayleigh number expressed as

[TABLE]

The momentum equation (171) can also be written for the poloidal decomposition in spherical harmonics as

[TABLE]

The last step is to express the boundary conditions in term of the poloidal decomposition and in non-dimensional form. Noting that $u_{r}=L^{2}P/r$ and that the horizontal integration of the momentum equation taken at $r=1$ gives $-p^{\prime}+\partial(r\nabla^{2}P)/\partial r=\text{cste}$, the stress free condition takes the form

[TABLE]

and the normal stress balance is

[TABLE]

where $\mathcal{P}$ is a dimensionless number comparing the timescale of viscous relaxation of the boundary $\eta/\Delta\rho g_{\text{icb}}r_{\text{icb}}$ and the time scale of phase change $\tau_{\phi}$, de facto quantifying the permeability of the inner core boundary. It is defined as

[TABLE]

10 Unstable or stable stratification in the inner core?

The core crystallises from the center outward because the solidification temperature of the core mixture increases with depth faster than the (adiabatic) core geotherm (Jacobs, 1953). One consequence of this solidification mode is that the inner core is cooled from above, a configuration which is potentially prone to thermal convection. Thermal convection further requires the inner core temperature profile to be superadiabatic, which depends on a competition between extraction of the inner core internal heat by diffusion and advection, and cooling at the ICB. Equation (158) shows that superadiabaticity in the inner core (i.e. a potential temperature increasing with depth) requires $S_{T}$ to be positive. Fast cooling and a low inner core thermal diffusivity ($S_{T}>0$) promotes superadiabaticity; slow cooling and high thermal diffusivity ($S_{T}<0$) results in a stable thermal stratification.

The expression of $S_{T}$ (equation 159) can be rewritten by writing the time derivative of the temperature at the ICB as a function of the rate of inner core growth (Deguen and Cardin, 2011). This gives

[TABLE]

where $\text{d}T_{s}/\text{d}T_{\text{ad}}$ is the ratio of the Clapeyron slope to the adiabat. It is then straightforward that the term $S_{T}$ is positive only if

[TABLE]

Unfortunately, the uncertainties on the thermal conductivity and inner core growth rate are such that the sign of $S_{T}$ is not known with much certainty. The high value of the thermal conductivity currently favoured (de Koker et al., 2012; Pozzo et al., 2012; Gomi et al., 2013, 2016) results in a stable thermal stratification.

The inner core may also have developed a compositional stratification. The concentration in light elements of newly crystallised solid, $c^{s}_{\text{icb}}$, is linked through the partition coefficient $D$ to the concentration in the liquid from which it crystallises, $c^{l}_{\text{icb}}$, as

[TABLE]

while its derivative with respect to inner core size is

[TABLE]

A stable compositional stratification would develop if $c^{s}_{\text{icb}}$ increases with increasing inner core size (more light elements in the upper part of the inner core); conversely, an unstable stratification would develop if $c^{s}_{\text{icb}}$ decreases with increasing inner core size. The first term on the right-hand-side term of equation (182) is very likely positive, due to the gradual enrichment of the outer core in light elements expelled during crystallisation. The second term depends on how $D$ varies with pressure and temperature along the $(P,T)$ path defined by the evolution of the position of the inner core boundary. Ab initio calculations (Gubbins et al., 2013) suggest that it is negative, and of the same order of magnitude as the first term on the r.h.s. of equation (182). The relative importance of the two terms depend on the exact composition of the inner core (Gubbins et al., 2013; Labrosse, 2014), which is not very well constrained. Again, it is difficult to be definitive: given our current knowledge of the composition of the core and of the partitioning behaviour of its light elements, stable and unstable compositional stratifications seem equally plausible.

Natural (thermal or compositional) convection in the inner core has been studied in details (e.g. Weber and Machetel, 1992; Wenk et al., 2000; Alboussière et al., 2010; Monnereau et al., 2010; Deguen and Cardin, 2011; Cottaar and Buffett, 2012; Deguen et al., 2013; Mizzon and Monnereau, 2013; Deguen et al., 2018). In the limit $\mathcal{P}\rightarrow 0$, the convection instability takes the form of a translation of the inner core, with melting on one hemisphere and solidification on the other (Alboussière et al., 2010; Monnereau et al., 2010; Deguen et al., 2013; Mizzon and Monnereau, 2013; Deguen et al., 2018). We will focus here on the case of neutral or stable stratification, and consider the flow forced by the Lorentz force associated with the magnetic field diffused in the inner core from the outer core.

11 Deformation induced by the Lorentz force

As discussed in section 6, the flow in the outer core sustains a magnetic field extending upward to the surface of the Earth but also inward inside the inner core. The magnetic Reynolds number of the inner core is likely small: assuming for example a velocity $10^{-10}$ m.s*-1* gives a magnetic Reynolds number on the order of $10^{-5}$. This shows that the magnetic field is only diffused inside the inner core, with no net advection or generation of the field. A diffused magnetic field in the inner core will add two terms in the set of equations: the Lorentz force in the momentum equation, and Joule heating in the energy equation. We are interested here in the flow driven by the Lorentz force in the inner core (Lasbleis et al., 2015).

Geodynamo simulations often exhibits a strong toroidal magnetic field close to the inner core boundary. As we are interested in the largest effect on the inner core dynamics, we consider here only low-order toroidal components of the magnetic field at the ICB, which have the largest penetration length scale.

We thus add in the momentum equation the Lorentz force due to a purely toroidal and axisymmetric magnetic field of degree two at the ICB $\mathbf{B}|_{\text{icb}}=B_{0}\sin\theta\cos\theta\mathbf{e}_{\phi}$ (Buffett and Bloxham, 2000). Imposing this field at the ICB and solving for its diffusion inside the inner core assuming it does not vary with time ($\boldsymbol{\nabla}^{2}\mathbf{B}=0$), we obtain $\mathbf{B}=B_{0}r^{2}/r_{\text{icb}}^{2}\cos\theta\sin\theta\mathbf{e}_{\phi}$. This field is associated to an electric current density $\mathbf{j}=\frac{1}{\mu_{0}}\boldsymbol{\nabla}\times\mathbf{B}$, with $\mu_{0}$ the magnetic permeability. The associated Lorentz force is $\mathbf{f}_{L}=\mathbf{j}\times\mathbf{B}$ which non-potential part (magnetic tension) can be written as

[TABLE]

with $f_{r}$ and $f_{\theta}$ two functions of $\theta$ expressed as

[TABLE]

Injecting this force into the dimensionless Stokes equation, we obtain

[TABLE]

where

[TABLE]

is similar to a Hartmann number, quantifying the ratio of the Lorentz force to the viscous force, and $\mathbf{\tilde{f}}_{L}$ is defined as in (183) without the prefactor $\frac{B_{0}^{2}}{\mu_{0}r_{\text{icb}}}$.

Equation (186) is solved using a poloidal decomposition and horizontal spherical harmonics decomposition. The term corresponding to the Lorentz force gives

[TABLE]

where $Y_{2}^{0}=\sqrt{5}(3\cos^{2}\theta-1)/2$. The momentum equation can thus be written as an equation for the spherical harmonics components $P_{l}^{m}$ and $t_{l}^{m}$of respectively the poloidal component of the velocity and the temperature as

[TABLE]

where $\delta$ is the Kronecker symbol.

11.1 Neutral stratification

We first consider the neutral stratification end-member where $Ra=0$. In that case, the only force driving flows in the system is the Lorentz force, and we do not need to solve for the temperature or composition fields. The flow results from a balance between the Lorentz and viscous forces. Since the characteristic length scale of velocity variations must be the size of the inner core (1 in dimensionless form), we have

[TABLE]

which implies that the magnitude of the velocity field should be proportional to $M$.

We can now solve analytically the flow field for a neutral stratification. With $Ra=0$, equation (189) reduces for $l=2$ and $m=0$ to

[TABLE]

which we solve with the boundary conditions at $r=1$ described in section 9.1:

[TABLE]

Equation (191) is a fourth order non-homogeneous differential equation, which solution can be obtained by solving the homogeneous equation and noticing that

[TABLE]

is one solution of the complete equation. Searching for a polynomial solution, we find that $r^{\alpha}$ is solution of the homogeneous equation $D_{2}^{2}P_{2}^{0}=0$ if $\alpha$ is a zero of the polynomial expression $[\alpha(\alpha+1)-6][(\alpha-2)(\alpha-1)-6]$. We then obtain the general solution of equation (191) as

[TABLE]

$A$ and $B$ must be equal to 0 for the velocity to remain finite at $r=0$. $C$ and $D$ are obtained from the boundary conditions at $r=1$, and we finally obtain

[TABLE]

From the expression of $P_{2}^{0}$, we can now obtain the expressions for the velocity field from

[TABLE]

Defining the r.m.s. velocity as

[TABLE]

we obtain

[TABLE]

This expression for the RMS velocity gives insight on the effect of $\mathcal{P}$ on the global dynamics. As predicted at the beginning of the subsection, the velocity is indeed a linear function of $M$, with the boundary conditions modifying the prefactor. For $\mathcal{P}\rightarrow 0$, $V_{\text{RMS}}\sim 0.066M$, and for $\mathcal{P}\rightarrow\infty$, $V_{\text{RMS}}\sim 0.008M$. Permeable boundary conditions ($\mathcal{P}\rightarrow 0$) give velocities about one order of magnitude higher than impermeable boundary conditions ($\mathcal{P}\rightarrow\infty$).

11.2 Stable stratification

If the inner core is stably stratified ($Ra<0$), the buoyancy forces resulting from the deformation of constant density surfaces would tend to oppose further deformation, and inhibit vertical motions. We can thus anticipate that the flow obtained in the limit of $Ra=0$ can be significantly altered by a strong stratification.

To estimate to what extent the neutral stratification solution can be altered by the presence of a stable density stratification, we consider the vorticity equation (obtained by taking the curl of equation (186)), which is

[TABLE]

where $\boldsymbol{\zeta}=\boldsymbol{\nabla}\times\mathbf{u}$ is the vorticity. Because the form of the magnetic field considered here forces a degree 2 flow, and because in non-dimensional form $\chi$ varies between 0 and 1, the $\theta$ derivative of $\chi$ is $\lesssim 1$. The magnitude of the baroclinic vorticity production $-Ra\frac{\partial\chi}{\partial\theta}$ is therefore $\lesssim-Ra$. Since the vorticity production associated with the Lorentz force is on the order of $M$, a stable stratification can affect significantly the flow only if $-Ra\gtrsim M$. As one of the main unknown for inner core dynamics is the viscosity, it is interesting to note that the ratio $M/Ra$,

[TABLE]

does not depend on the viscosity, so that the boundary between a strongly stratified regime and a neutral stratification does not depends on the viscosity.

Figure 16 shows the temperature and vorticity fields obtained by solving numerically equations (172) and (189) at $M=10^{4}$ and $Ra=-10^{4}$, $-6\times 10^{4}$, $-10^{5}$, and $-10^{6}$. The stratification has a negligible effect on the flow at $Ra=-10^{4}$, but already has a significant effect at $Ra=-6\times 10^{4}$, which is consistent with the criterion we just derived ($-Ra\gtrsim M$ for a strong effect of the stratification). At $Ra=-10^{5}$ and $-10^{6}$, the flow induced by the Lorentz force is essentially confined to a thin layer below the inner core boundary, in which the flow direction is essentially horizontal. Lasbleis et al. (2015) found that the thickness of this layer is $\propto(-Ra)^{-1/6}$ and the strain rate in the layer is $\propto M(-Ra)^{-1/3}$. For application to Earth’s inner core, we take $\Delta\rho\sim-1$ kg.m*-3* and $B_{0}\sim 4$ mT (Gillet et al., 2010), and find $M/(-Ra)\sim 10^{6}$, which would imply that stratification strongly affects the flow forced by the Lorentz force.

Chapter 7 Core formation

The final section of this chapter concerns the formation of Earth’s core. The Earth formed about 4.56 Gy ago through the accretion of solar nebula material, a process which is estimated to have taken a few tens of millions of years. Accretion in the solar system went through different dynamical phases which involved increasingly energetic and catastrophic impacts and collisions (Lunine et al., 2011; Walsh et al., 2011). The last phase of accretion, in which most of the Earth mass was accreted, involved extremely energetic collisions between already differentiated planetary embryos ($\sim 100-1000$ km size), which resulted in widespread melting and the formation of magma oceans.

The basic ingredients of a terrestrial planet – an iron-rich metal and silicates – are immiscible, and can separate under the action of gravity to form an iron-rich core surrounded by a rocky mantle. The timescale of phase separation is much smaller than the accretion timescale if at least one of the two phases is liquid, and it is thus believed that core formation has been concomitant with Earth’s accretion. Importantly, the metal added to the Earth by each impact had no reason to be in thermodynamic equilibrium with the silicate mantle of the growing Earth, which implies that heat and chemical elements would have been transferred from one phase to the other when the two phases were in contact. Chemical elements and heat released during accretion have thus been partitioned between the core and mantle, in a way which depends on the exact physical mechanism by which metal and silicate have separated.

This has important implications for the state of the planet at the end of its accretion, and its subsequent evolution. The partitioning between core and mantle of the heat released during accretion has set the initial temperature contrast between the mantle and core, a key parameter for the early dynamics of the planet, with implications for the possibility of forming a basal magma ocean (Labrosse et al., 2007), the existence of an early dynamo (Williams and Nimmo, 2004; Monteux et al., 2011), and the subsequent thermal and magnetic evolution of the planet. The partitioning of chemical elements between the core and mantle has been used to constrain the timing of differentiation and the physical conditions under which it occurred (e.g. Yin et al., 2002; Halliday, 2004; Wood et al., 2006; Rudge et al., 2010). It has also important implications for a number of geodynamical issues. One example is the identity and abundance of light elements and radioactive elements in the core (Corgne et al., 2007; Badro et al., 2015), which depend on the conditions (pressure $P$, temperature $T$, oxygen fugacity $f{O_{2}}$) at which metal and silicate have interacted for the last time.

The large impacts which dominated the last stages of Earth formation injected enormous amounts of kinetic energy into the magma oceans, creating highly turbulent environments in which it has been conjectured that the cores of the bodies impacting the Earth would fragment down to centimetre scale, at which metal and silicates can efficiently exchange chemical elements and heat (Stevenson, 1990; Karato and Murthy, 1997; Rubie et al., 2003). The metal is envisioned to disperse in the magma ocean and equilibrate with the silicates, before raining out. It would then collect at the base of the magma ocean, and finally migrate toward the core as large diapirs (Stevenson, 1990; Karato and Murthy, 1997; Monteux et al., 2009; Samuel et al., 2010) or by the propagation of iron dykes (Stevenson, 2003), at which point further chemical equilibration is unlikely to be significant.

Most geochemical models of core formation are based on this so-called iron rain scenario, but the fluid dynamics involved is actually poorly understood, even at a qualitative level. Efficient chemical exchange requires a high metal-silicates interfacial area-to-volume ratio, which requires fragmentation or stretching of the metal down to $\sim$ cm size. Differentiation of terrestrial bodies started early, and it is now recognized that most of the mass of the Earth was accreted from already differentiated planetary bodies, with cores of their own. Whether or not these large volumes of iron ($\sim 100$ to 1000 km) would indeed fragment down to cm scale at which chemical equilibration can occurs therefore remains an open question, and a matter of much speculation (Dahl and Stevenson, 2010; Deguen et al., 2011; Samuel, 2012; Deguen et al., 2014; Wacheul et al., 2014).

We have yet no well-tested and self-consistent theory of fragmentation and thermal or compositional equilibration in the context of metal-silicate separation in a magma ocean. For this reason we have here made the choice of focusing on the basic concepts and mechanisms which we believe are important to know and understand to tackle this problem, rather than trying to give a definitive answer to this question.

12 Problem set-up and non-dimensional numbers

To keep things reasonably simple, we will ignore here the dynamics of the impact itself and consider the fate of an initially spherical mass of molten iron falling into silicates, which can be either solid or molten. This molten iron mass may be either the core of a planetary body impacting the Earth, or a fragment of the core if it has been significantly dispersed by the impact. As an additional simplifying hypothesis, we will even assume that the metal mass has no initial velocity.

The volume of molten iron is assumed to be close to spherical, and has a radius $R_{0}$. We denote by $\rho_{m}$ and $\rho_{s}=\rho_{m}-\Delta\rho$ the densities of metal and silicates, and by $\eta_{m}$ and $\eta_{s}$ their viscosities. The metal and silicates phases are immiscible and we denote by $\gamma$ the interfacial tension of the metal-silicates interface. We denote by $g$ the acceleration of gravity.

The evolution of the metal mass depends on 7 dimensional parameters ($R_{0}$, $\rho_{m}$, $\rho_{s}$, $\eta_{m}$, $\eta_{s}$, $\gamma$) involving 3 fundamental units (length, weight, and time). According to Vashy-Buckingham’s theorem, the number of independent non-dimensional numbers to be used to describe the problem is equal to $7-3=4$. One possible set is the following:

[TABLE]

The Bond number is a measure of the relative importance of buoyancy and interfacial tension. The Grashof number is basically a Reynolds number obtained by taking Stokes’ velocity as a velocity scale. Additionnal useful numbers include Reynolds and Weber numbers, which can be defined as

[TABLE]

where $U$ is a velocity scale to be defined. The Reynolds and Weber numbers compare inertia to viscous forces and interfacial tension, respectively. The density ratio is close to 2.

If in addition we consider heat or mass transfer between metal and silicates, then two additional dimensional parameters enter the problem: the diffusivities (thermal or compositional) $\kappa_{m}$ and $\kappa_{s}$ in the metal and silicates. Since there is no additional fundamental units, Vashy-Buckingham’s theorem implies that two additional non-dimensional numbers must be used. One possible choice is to use the ratio of the diffusivities and a Péclet number:

[TABLE]

In what follows, we will assume that $\kappa_{s}=\kappa_{m}$ for the sake of simplicity.

13 Preliminary considerations

13.1 Terminal velocity

Since we have chosen to focus on the case of metal mass falling with no initial velocity, a relevant velocity scale is its terminal velocity, reached when the buoyancy force ($\sim\Delta\rho gR_{0}^{3}$) is balanced by the drag on the surface of the metal mass. Two different scalings can be obtained depending on whether the drag is dominated by viscous stresses ($\sim\eta_{s}U/R_{0}$) or dynamic pressure ($\sim\rho_{s}U^{2}$). The ratio of the dynamic pressure and viscous stresses contributions to the total drag is

[TABLE]

The drag on the metal mass is obtained by multiplying the dominant stress by the surface area of the metal mass: it is on the order of $\eta_{s}UR_{0}$ if the drag is dominated by the contribution of viscous stresses (low $Re$), and on the order of $\sim\rho_{s}U^{2}R_{0}^{2}$ if it is dominated by the contribution of dynamic pressure (high $Re$). The force balance on the metal mass can thus be written as

[TABLE]

The terminal velocity obtained from this balance is

[TABLE]

The first scaling gives Stokes’ settling velocity, and corresponds to the low $Re$ limit; the second scaling is the so-called newtonian scaling, and corresponds to the high $Re$ limit. The two velocities are on the same order of magnitude when

[TABLE]

which defines the boundary between the two scalings (here $\nu_{s}=\eta_{s}/\rho_{s}$). The terminal velocity is thus

[TABLE]

With these scalings, the Reynolds and Weber numbers based on the terminal velocity are given by

[TABLE]

13.2 Maximal stable size of a falling drop

Interfacial tension (unit J.m*-2*) can be interpreted as an energy per unit area. Deforming an interface in a way which results in an increase of the interfacial area costs energy, and interfacial tension effects will tend to minimise the surface area of the interface. If no other force acts on a drop, interfacial tension would keep it spherical, hence minimising its surface area.

Interfacial tension can also be seen as a force per unit of length (it can be verified that J.m${}^{-2}=$N.m*-1*): if a piece of an interface is divided into two parts, the force imparted by one part of the surface on the other is parallel to the interface and has a magnitude given by the product of the interfacial tension with the length of the curve separating the two parts of the surface. If integrated over a curved surface, one can also show that interfacial tension induces a pressure jump across the interface equal to

[TABLE]

where $R_{1}$ and $R_{2}$ are the principal radii of curvature. This pressure jump is called Laplace’s pressure. Across a spherical interface (the surface of a drop or bubble of radius $R_{0}$), Laplace’s pressure is equal to $2\gamma/R_{0}$.

A falling drop can be deformed by the stresses imparted by the surrounding fluid onto the drop, or in other words by the fluid drag. If the total drag on the drop is $F_{\mathrm{drag}}$, the mean stress on the surface of the drop is $\sim F_{\mathrm{drag}}/R_{0}^{2}$. One can expect significant deformation of the drop if the hydrodynamic stress variations due to the drag exceed Laplace’s pressure:

[TABLE]

If the drop reached its terminal velocity, then the drag must be equal to the total buoyancy of the drop, $F_{\mathrm{drag}}\sim\Delta\rho gR_{0}^{3}$. Combining this with equation (220), we find that strong deformation of the drop will happen if its radius is larger than a critical radius $R_{c}$ given by

[TABLE]

This length is also known as a capillary length: it is the length scale over which buoyancy effects dominate over surface tension effects. Interfacial tension will keep the drop close to spherical if its radius is smaller than $R_{c}$. Equation (221) is equivalent to writing that deformation is significant if the Bond number of the drop is large compared to 1. The interfacial tension between metal and silicates is on the order of 1 J.m*-2* and $\Delta\rho\sim 4000$ kg.m*-3*. With $g\sim 10$ m.s*-2*, we thus have $R_{c}\sim 5$ mm.

Strong drop deformation may happen before reaching terminal velocity, and in this case the above scaling will not be the most relevant. If drag is dominated by viscous effects ($F_{\mathrm{drag}}\sim\eta_{s}UR_{0}$), then we find that deformation of the drop may happen if its velocity is larger than a critical velocity

[TABLE]

This criterion is of limited use since if the drag is dominated by viscous effect (which means that $Re\ll 1$), then the drop velocity will very quickly reach a terminal velocity equal to Stokes’ settling velocity. Using Stokes’ velocity for $U$ in equation (222) gives equation (220).

If instead drag is dominated by the contribution of dynamic pressure ($F_{\mathrm{drag}}\sim\rho_{s}U^{2}R_{0}^{2}$), then we find that strong deformation requires the drop Weber number is large:

[TABLE]

This criterion reduces to equation (221) when the drop reaches its terminal velocity ($U\sim\sqrt{({\Delta\rho}/{\rho_{s}})gR_{0}}$ since in this limit $Re\gg 1$).

14 The low Reynolds limit: diapirism

A first relevant limit of the problem described in section 12 corresponds to molten metal diapirs travelling through a solid, or partially molten, silicate layer. The radius of these diapirs may be similar to the size of the core of the impactors, say somewhere between 1 km and 1000 km. The acceleration of gravity is smaller or equal to its current value in Earth’s mantle, $g\sim 10$ m.s*-2*. The viscosity $\eta_{s}$ of the silicates is a strong function of temperature, and can also be significantly decreased if the silicate layer is partially molten. A reasonable range is $10^{15}$ to $10^{21}$ Pa.s. The viscosity of molten iron is $\sim 10^{-2}$ Pa.s.

With this parameter values, $Bo\gtrsim 10^{10}$, $Gr\lesssim 10^{-4}$, $\rho_{m}/\rho_{s}\sim 2$, $\eta_{m}/\eta_{s}\lesssim 10^{-18}$. The Grashof number being small, we are well into the low Reynolds regime: viscous forces dominate over inertia in the silicates. The limit of low $Re$ and high $Bo$ has been studied numerically in the context of core formation (Samuel and Tackley, 2008; Monteux et al., 2009; Samuel et al., 2010), and experimentally in other contexts (e.g. Ribe, 1983; Bercovici and Kelly, 1997). In this limit surface tension is unimportant, but the volume of metal is kept roughly spherical because $Re\ll 1$. Since viscous forces are so important in the silicates, the flow around the metal mass is limited to spatial scales on the order of $R_{0}$, which limits the deformation of the metal mass (in other words, small scale perturbations of the metal-silicate interface shape are damped viscously).

The falling velocity is thus simply given by Stokes’ velocity (equation (215)). The law of heat or mass transfer between the diapir and its surrounding is also well known (e.g. Clift et al., 1978; Ribe, 2007; Ulvrová et al., 2011): the heat flux is given by

[TABLE]

where $k_{s}$ is the thermal conductivity in the silicates, $\bar{T}$ the mean temperature in the diapir, $T_{s}(z)$ the temperature of the surrounding silicate layer at depth $z$, and $a$ is a constant on the order of 1. A similar expression can be written for chemical elements transfer.

The heat balance of the diapir writes

[TABLE]

where $c_{p,m}$ is the heat capacity of the metal phase. Transforming the time derivative into a derivative with respect to the distance $z$ travelled by the diapir (using $d(...)/dt=Ud(...)/dz$) and re-arranging gives

[TABLE]

where the characteristic length $\ell$ is given by

[TABLE]

$\ell$ is the characteristic distance over which the temperature of the diapir responds to changes of the surrounding temperature, the thermal equilibration distance. The general solution of equation (226) is

[TABLE]

where $T_{0}$ is the initial temperature of the diapir. In practice the amount of heat transfer is small because $\ell$ is typically larger than the mantle thickness. The compositional diffusivity in the solid silicates being perhaps four orders of magnitude smaller that the thermal diffusivity, exchange of chemical elements would be even smaller (the equilibration distance being $\propto\,Pe^{1/2}\propto\,\kappa^{-1/2}$). Diapirs migrating through a solid part of the mantle would therefore exchange a negligible amount of heat and chemical elements with the surrounding mantle.

15 The high Reynolds limit: metal-silicates separation in a magma ocean

Let us now consider the case of a volume of molten metal falling into molten rocks - a magma ocean. The parameter values are similar to what we have considered when discussing the case of diapirism, except that the viscosity of the silicates is much smaller, on the order of $10^{-1}$ Pa.s. With these parameter values, $Bo\gtrsim 10^{10}$, $Gr\gtrsim 10^{22}$, $\rho_{m}/\rho_{s}\sim 2$, $\eta_{m}/\eta_{s}\lesssim 10^{-18}$. Since $Gr\gg 1$, the relevant velocity scale is the newtonian scaling given by equation (216). This gives $Re\sim Gr^{1/2}\gtrsim 10^{11}$ and $We\sim Bo\gtrsim 10^{10}$. The very large values of $Re$, $Bo$, and $We$ imply that neither viscous forces nor interfacial tension can keep the metal volume spherical: a molten mass of metal falling into a magma ocean should suffer significant deformation, possibly resulting in its fragmentation into drops.

15.1 Observations from laboratory experiments

The very large values of the Bond and Weber numbers imply that interfacial tension should be unimportant for the large scale dynamics. This suggests to first look at the case of infinite Bond and Weber numbers, which corresponds to the limit of miscible fluids.

Figure 17a shows snapshots from an experiment in which a negatively buoyant volume of an aqueous solution of sodium iodide is released into fresh water. The volume of the (dyed) negatively buoyant fluid is seen to increase as it falls, which indicates that it entrains and incorporates ambient fluid, resulting in its gradual dilution. Measurements show that the mean radius of the dyed mixture increases linearly with the distance from the point of release. This is what is known as a turbulent thermal in the fluid mechanics and atmospheric science communities (e.g. Batchelor, 1954; Scorer, 1957; Woodward, 1959). The name thermal is inherited from the usage of glider pilots, for whom a thermal is an isolated mass of warm air rising through the lower atmosphere. Though the buoyancy in atmospheric thermals is due to temperature differences, the nature of the source of buoyancy (thermal or compositional) happens to be of secondary importance and introduces no qualitative difference. The term thermal has since been used to denote an isolated buoyant mass of a fluid rising or falling (depending on the sign of the buoyancy), irrespectively of the nature of the source of buoyancy. Here we will also often use the term buoyant cloud instead of thermal.

Figure 17b shows snapshots from a similar experiment in which a negatively buoyant volume of sodium iodide is now released into a low-viscosity silicon oil. The NaI solution and the silicon oil are immiscible, so we are one important step closer to the core-mantle differentiation configuration. The experimental fluids and configuration have been chosen so as to maximize the values of the Bond and Reynolds numbers, which are $Bo=3.4\times 10^{4}$ and $\mathrm{Re}=5.5\times 10^{4}$. The density ratio is $\rho_{m}/\rho_{s}=1.95$ (close to metal-silicate), and the viscosity ratio is $\eta_{m}/\eta_{s}=1.2$. The large-scale evolution of the negatively buoyant volume is strikingly similar to what has been observed in the miscible experiment: the volume of the negatively buoyant fluid increases linearly with distance, which indicates that it entrains and incorporates silicon oil. PIV measurements on a similar experiment (figure 18) show that the velocity field has a vortex ring structure, with most of the entrainment of silicon oil probably occurring from the rear of the cloud.

In contrast, the small scale structure in the immiscible experiment is qualitatively different from what we can observe in the miscible experiment. In the miscible experiment, the negatively buoyant solutions mixes with the entrained water, diffusion of the NaI salt allowing homogenisation at the molecular scale. In the immiscible experiment, the NaI solution of course does not mix with the entrained silicon oil since the two liquids are immiscible. A close inspection of the last snapshot of the immiscible experiment reveals that the dense phase has been fragmented into droplets.

15.2 Large-scale dynamics: turbulent entrainment model

We consider the evolution of a mass of negatively buoyant fluid falling into another one (figure 19). The mass of negatively buoyant fluid has an initially spherical shape and an initial radius $R_{0}$, and has a density $\rho_{a}+\Delta\rho$, where $\rho_{a}$ is the density of the surrounding fluid. We denote by $u_{z}$ the vertical velocity of the center of mass of the negatively-buoyant cloud, and by $R(z)$ its spatial extension. Following Batchelor (1954), we assume that far from the source $u_{z}$ and $R$ only depend on either distance $z$ or time $t$, and on its total amount of buoyancy defined as

[TABLE]

where $V_{0}$ is the initial volume of the released buoyant mass. We thus assume that surface tension has no effect of the evolution of $u_{z}$ and $R$. Dimensional analysis then shows that

[TABLE]

which predicts that the spatial extension of the cloud increases linearly with $z$: the cloud must therefore entrain ambient fluid. The prediction that the mean velocity decreases as $z^{-1}$, or, equivalently, as $R^{-1}$, is consistent with the fact that the total buoyancy of the cloud is conserved, but not its volume. The buoyancy is “diluted” by the incorporation of neutrally buoyant ambient fluid to the cloud.

A more physical (and more general) way to obtain the evolution of a turbulent thermal has been given by Morton et al. (1956), who based their analysis on the assumption that the rate of entrainment of ambient fluid within the buoyant cloud is simply proportional to the mean vertical velocity $u_{z}$ and to the surface area of the cloud (figure 19). This is the basic assumption of the turbulent entrainment models used to describe the dynamics of turbulent clouds, plumes, and jets.

This assumption implies that the time derivative of the cloud volume is given by

[TABLE]

where $\alpha$ is the entrainment coefficient. Noting that $d(...)/dt=u_{z}d(...)/dz$, where $z$ is the vertical position of the center of mass of the cloud, integration of equation (229) gives

[TABLE]

which is consistent with the prediction of dimensional analysis.

Conservation of momentum allows to obtain a predictive law for the vertical velocity $u_{z}$ of the center of mass of the cloud. Ignoring fluid drag on the cloud and a possible loss of buoyancy in the wake of the cloud, conservation of momentum simply states that

[TABLE]

where $\bar{\rho}$ is the mean density of the cloud, given by

[TABLE]

Conservation of mass implies that

[TABLE]

Using this relation in equation (231), using the transformation $d(...)/dt=u_{z}d(...)/dz=\alpha u_{z}d(...)/dR$, and re-arranging gives

[TABLE]

a linear first order ordinary differential equation with varying coefficients. The solution, written for $u_{z}$ as a function of $z$, is

[TABLE]

where

[TABLE]

The function $\mathcal{F}$ tends toward 1 at large $\alpha z/R_{0}$. The full solution is thus consistent with the dimensional analysis prediction when far from the source. Close to the source, the velocity is given (at first order in $\alpha z/R_{0}$), by

[TABLE]

Though the turbulent thermal model just described has been developed to model miscible flows, experiments suggest that it can also be applied to immiscible fluids in situations where the Bond and Weber numbers are large (Deguen et al., 2014; Landeau et al., 2014; Wacheul et al., 2014; Wacheul and Le Bars, 2017). A qualitative comparison of figures 17a and 17b suggests that it is indeed the case. The similarity between the miscible and immiscible experiments also holds on a quantitative level: from a series of experiments similar to that presented on figure 17, we have been able to show that the evolution of both $R$ and $u_{z}$ are very well described by the turbulent entrainment model with $\alpha=0.25\pm 0.05$, similar to miscible flows (Deguen et al., 2014; Landeau et al., 2014). This demonstrates that there is indeed no effect of surface tension on the large scale part of the flow.

15.3 Fragmentation

15.3.1 Qualitative observations from experiments

In the experiment shown on figure 17b and in similar experiments, most of the fragmentation of the dense liquid occurs during a relatively short time span. In figure 17b, the dense phase is essentially continuous until the third snapshot, and almost entirely fragmented into drops at the fourth snapshot. The analysis of images obtained with a high-speed camera (1kHz) shows that drops formation results from two mechanisms:

the fragmentation of stretched cylindrical ligaments of aqueous solution through the Rayleigh-Plateau capillary instability, as shown on figure 20a (the mechanism of the Rayleigh-Plateau instability is explained below), 2. 2.

the fragmentation of thin liquid films, as shown in figure 20b. In this regime, thin films of aqueous solutions are stretched by the flow before eventually being punctured. The film then quickly retracts, the liquid forming the film gathering into ligaments which then fragment into drops due to the Rayleigh-Plateau instability.

These two modes of fragmentation are classically observed in fragmentation problems in a variety of contexts. In fact, liquid fragmentation necessitates a capillary instability, irrespectively of the nature of the flow (Villermaux, 2007). What varies from one problem to another is the sequence of mechanisms resulting in the formation of ligaments which can fragment as a result of the Rayleigh-Plateau capillary instability. In experiments such as shown on figure 17b, the observed sequence is the following: (i) the interface is destabilised and deformed by the combined effect of shear and Rayleigh-Taylor instabilities; (ii) three-dimensional structures generated by the destabilisation of the interface are stretched and stirred by the mean flow and velocity field fluctuations; (iii) stirring produces ligaments and films, which will then break up and produce a population of drops.

15.3.2 Rayleigh-Plateau capillary instability

In many situations, the deformation of an interface results in an increase of its surface area, and hence of its interfacial energy. In this case deformation is not energetically favoured, and mechanical work therefore has to be provided to deform the interface. This is for example the case of initially planar interface: any perturbation of the interface results in an increase of its surface area and energy.

In contrast, the deformation of a cylindrical interface can, under certain conditions, result in a decrease of its surface area, and hence of its interfacial energy. Take a cylinder of one liquid into another, of length $L$ and radius $R_{0}$. Its surface area is $2\pi R_{0}L$ and its interfacial energy is $2\pi R_{0}L\gamma$. It is easy to see that the cylindrical shape is not very favourable from an energetic point of view: if the liquid of the cylinder is re-arranged to form a sphere of the same volume ($\pi R_{0}^{2}L$), the sphere will have a radius equal to $\left[(3/4)R_{0}^{2}L\right]^{1/3}$, and a surface area equal to $4(3/4)^{2/3}\pi R_{0}^{4/3}L^{2/3}$, which is smaller than the cylinder surface area if the length of the cylinder is larger than $(9/2)R_{0}$. This shows that the fragmentation into drops of liquid cylinder is energetically favoured if the ratio of its length and radius is larger than $9/2$.

Is fragmentation dynamically possible? To see if it is, let us consider again a liquid cylinder of radius $R_{0}$, and assume now that its surface is perturbed from its initial shape as

[TABLE]

where $x$ is the coordinate along the axis of symmetry of the cylinder, and $\epsilon$ and $\lambda$ the amplitude and wavelength of the perturbation (figure 21). Note that conservation of mass implies that $\bar{R}<R_{0}$: comparing the volume of a section of length $\lambda$ of the unperturbed and perturbed states indeed shows that

[TABLE]

Calculating the surface area of the perturbed cylinder shows that the perturbation induces a decrease of the surface area (and hence of energy) if $\lambda>2\pi R_{0}$, which suggests that the cylindrical shape may be unstable to perturbations with wavelengths larger than $2\pi R_{0}$.

To see how the instability works, let us consider the two principal curvatures of the interface (figure 21). One is the curvature associated with the radius of the cylinder, $1/R(x)$, and the other is the curvature associated with the longitudinal perturbation, $1/R_{\lambda}(x)$, which is equal to the divergence in the $x-$direction of the normal $\mathbf{n}$ of the interface. The Laplace pressure jump across the interface is equal to

[TABLE]

The contribution of $1/R_{\lambda}$ is positive where $R>R_{0}$, and negative where $R<R_{0}$. It thus produces a pressure gradient from regions of large $R$ to regions of small $R$, which can drive a flow from large to small $R$ that would decrease the amplitude of the radius perturbation. It thus has a stabilizing effect. In contrast, the pressure jump associated with $1/R$ is larger in regions of small $R$. It thus produces a pressure gradient from small to large $R$, which may drive a flow that would increase the amplitude of the perturbation. The amplitude of the radius perturbation can therefore grow if the pressure gradient associated with the curvature $R(x)$ is larger in magnitude that the pressure gradient associated with the curvature $R_{\lambda}(x)$.

At first order in $\epsilon$, one find

[TABLE]

and

[TABLE]

which gives a pressure jump across the interface given by

[TABLE]

Its gradient along $x$ is given by

[TABLE]

This shows that the pressure gradient inside the liquid cylinder is from small to large $R$ if $\frac{4\pi^{2}}{\lambda^{2}}<\frac{1}{R_{0}^{2}}$, and from large to small $R$ instead. In other words, the initial perturbation will grow if

[TABLE]

which will eventually lead to the fragmentation of the cylinder into drops. This dynamical criterion is slightly more restrictive that the energy criterion, since $2\pi\simeq 6.26>9/2=4.5$.

Things get more complicated in situations where liquid ligaments are deformed and stretched by the ambient flow. Fragmentation can be significantly protracted by stretching effects (Taylor, 1934; Tomotika, 1936; Mikami et al., 1975; Eggers and Villermaux, 2008), which can be understood as follows. Let us consider a stretched ligament, the surface of which is modulated by a longitudinal perturbation of wavelength $\lambda$. The stretching will affect the disturbance, which will see its wavelength increase in proportion to the amount of stretching. If the perturbation wavelength initially corresponds to the optimal wavelength for the growth of the capillary instability, increasing the wavelength will decrease the rate of growth of the disturbance.

15.4 Chemical and heat transfer at the drop scale

If the metal phase ends up being fragmented into drops of size equal or smaller than the maximal stable size $R_{c}\sim\sqrt{\gamma/(\Delta\rho g)}$ (equation (221)), then thermal and chemical equilibration of the metal phase with the surrounding silicates is not an issue (e.g. Stevenson, 1990; Karato and Murthy, 1997; Rubie et al., 2003; Ulvrová et al., 2011), as shown below.

One can show (Lherm and Deguen, 2018) that the timescale of chemical equilibration of a falling metal drop with its surrounding is given by

[TABLE]

where $\kappa_{s}$ and $\kappa_{m}$ are the compositional diffusivities in the silicates and metal, $D_{\text{m/s}}$ the metal/silicates partitioning coefficient of a given chemical element, and $Pe_{s}=UR/\kappa_{s}$. The distance $\ell_{\text{eq}}$ fallen by the drop during a time $\tau_{\text{eq}}$ is given by

[TABLE]

In the case of siderophile elements ($D_{\text{m/s}}\gg 1$), this is approximated by

[TABLE]

For a metal drop falling into a magma ocean, $R_{c}$ is about 5 mm. The corresponding Grashof number is

[TABLE]

Magma ocean viscosity is estimated to be in the range $10^{-3}$-$10^{-2}$ Pa.s, which implies that the drop is in the newtonian regime and has a terminal velocity given by $\sim\sqrt{(\Delta\rho/\rho_{s}gR_{c}}$. With a compositional diffusivity $\kappa_{s}\sim 10^{-9}$ m.s*-2*, this gives a Péclet number around $10^{6}$ and we thus have $(R/6)\mathrm{Pe_{s}}^{1/2}\simeq 1$ m. Siderophile elements like Nickel or Cobalt have partitioning coefficients $D_{\text{m/s}}$ around $10^{3}$ at low pressure, and as low as $\sim 10$ when approaching the pressure of the core-mantle boundary. This would give $\ell_{\text{eq}}\sim 1$ km at low pressure, and $\ell_{\text{eq}}\sim 10$ m at high pressure. This is in both cases much smaller than typical magma ocean depth, which implies that drops of metal a few mm in size will readily equilibrate with the surrounding molten silicates.

However, the above conclusion rests on the assumption that the metal phase fragments into drops of a few mm in radius. We have no well tested fragmentation model that can be used in the context of core formation, so whether the metal phase would fragment or not is still an open question. It is also possible that efficient equilibration does not require fragmentation of the metal phase. The metal phase would necessarily be intensely stirred and stretched before fragmentation, and this may allow for efficient chemical transfer between the metal and silicates phase (Lherm and Deguen, 2018).

15.5 Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No 716429).