Earth’s Inner Core dynamics induced by the Lorentz force

Earth’s Inner Core dynamics induced by the Lorentz force

[
N

M. Lasbleis, R. Deguen, P. Cardin, S. Labrosse]M. Lasbleis, R. Deguen, P. Cardin, S. Labrosse
LGL-TPE, Université Lyon 1-ENS de Lyon-CNRS, Lyon, France, ISTerre, UGA-CNRS, Grenoble, France

{summary}

Seismic studies indicate that the Earth’s inner core has a complex structure and exhibits a strong elastic anisotropy with a cylindrical symmetry. Among the various models which have been proposed to explain this anisotropy, one class of models considers the effect of the Lorentz force associated with the magnetic field diffused within the inner core. In this paper we extend previous studies and use analytical calculations and numerical simulations to predict the geometry and strength of the flow induced by the poloidal component of the Lorentz force in a neutrally or stably stratified growing inner core, exploring also the effect of different types of boundary conditions at the inner core boundary (ICB). Unlike previous studies, we show that the boundary condition that is most likely to produce a significant deformation and seismic anisotropy is impermeable, with negligible radial flow through the boundary. Exact analytical solutions are found in the case of a negligible effect of buoyancy forces in the inner core (neutral stratification), while numerical simulations are used to investigate the case of stable stratification. In this situation, the flow induced by the Lorentz force is found to be localized in a shear layer below the ICB, which thickness depends on the strength of the stratification, but not on the magnetic field strength. We obtain scaling laws for the thickness of this layer, as well as for the flow velocity and strain rate in this shear layer as a function of the control parameters, which include the magnitude of the magnetic field, the strength of the density stratification, the viscosity of the inner core, and the growth rate of the inner core. We find that the resulting strain rate is probably too small to produce significant texturing unless the inner core viscosity is smaller than about Pa.s.

umerical solutions; Seismic anisotropy; Composition of the core

1 Introduction

The existence of structures within the inner core was first discovered by Poupinet et al. (1983), who discussed the possibility of lateral heterogeneity from the observation of P-waves travel time anomalies. These were then attributed to the existence of seismic anisotropy (Morelli et al., 1986; Woodhouse et al., 1986), with P-waves travelling faster in the north-south direction than in the equatorial plane. Since then, more complexities have been discovered in the inner core: a slight tilt in the fast axis of the anisotropy, radial variations of the anisotropy with a nearly isotropic upper layer, hemispherical variations of the thickness of the upper isotropic layer, an innermost inner core with different properties in anisotropy or attenuation, and anisotropic attenuation (See Souriau et al., 2003; Tkalčić & Kennett, 2008; Deguen, 2012; Deuss, 2014, for reviews, and references therein).

The seismic anisotropy can be explained either by liquid inclusions elongated in some specific direction (shape preferred orientation, SPO) (Singh et al., 2000) or by the alignment of the iron crystals forming the inner core (lattice preferred orientation, LPO). In the case of LPO, the orientation is acquired either during crystallization (e.g. Karato, 1993; Bergman, 1997; Brito et al., 2002) or by texturing during deformation of the inner core. Several mechanisms have been proposed to provide the deformation needed for texturing: solid state convection (Jeanloz & Wenk, 1988; Weber & Machetel, 1992; Buffett, 2009; Deguen & Cardin, 2011; Cottaar & Buffett, 2012; Deguen et al., 2013), or deformation induced by external forcing, due to viscous adjustment following preferential growth at the equator (Yoshida et al., 1996, 1999; Deguen & Cardin, 2009), or Lorentz force (Karato, 1999; Buffett & Bloxham, 2000; Buffett & Wenk, 2001).

Thermal convection in the inner core is possible if its cooling rate, related to its growth rate, or radiogenic heating rate is large enough to maintain a temperature gradient steeper than the isentropic gradient. In other words, the heat loss of the inner core must be larger than what would be conducted down the isentrope. However, the thermal conductivity of the core has been recently reevaluated to values larger than 90 W.m.K at the core mantle boundary and in excess of 150 W.m.K in the inner core (de Koker et al., 2012; Pozzo et al., 2012; Gomi et al., 2013; Pozzo et al., 2014), and this makes thermal convection in the inner core unlikely (Yukutake, 1998; Deguen & Cardin, 2011; Deguen et al., 2013; Labrosse, 2014). Inner core translation, that has been proposed to explain the hemispherical dichotomy of the inner core (Monnereau et al., 2010), results from a convection instability (Alboussière et al., 2010; Deguen et al., 2013; Mizzon & Monnereau, 2013) and is therefore also difficult to sustain.

Compositional convection is possible if the partition coefficient of light elements at the inner core boundary (ICB) decreases with time (Deguen & Cardin, 2011; Gubbins et al., 2013) or if some sort of compositional stratification develops in the outer core (Alboussière et al., 2010; Buffett, 2000; Gubbins & Davies, 2013; Deguen et al., 2013) so that the concentration of the liquid that crystallizes decreases with time. However, the combination of both thermal and compositional buoyancy does not favor convection in the inner core (Labrosse, 2014).

The strong thermal stability of the inner core resulting from its high thermal conductivity (Labrosse, 2014) is a barrier to any vertical motion and other forcing mechanisms need to work against it. This situation has already been considered in the case of deformation induced by preferential growth in the equatorial belt (Deguen & Cardin, 2009), and has been shown to produce a layered structure. Deguen et al. (2011) and Lincot et al. (2014) evaluated the predictions of anisotropy from this model and found that although it can induce significant deformation, it is difficult to explain the strength and geometry of the anisotropy observed in the inner core.

In this paper, we consider another major external forcing that was proposed, Maxwell stress. This was first proposed by Karato (1999) who considered the action of the Lorentz force assuming the inner core to be neutrally buoyant throughout. This situation is rather unlikely and, as discussed above, we expect the inner core to be stably stratified. Buffett & Bloxham (2000) have shown that in this case the flow is confined in a thin layer at the top of the inner core, similar to the case discussed above for a flow driven by preferential growth at the equator. However, the growth of the inner core gradually buries the deformed iron and this scenario may still produce a texture in the whole inner core. All these previous studies considered a fixed inner core size and infinitely fast phase change at the ICB. The moving boundary brings an additional advection term in the heat balance which can influence the dynamics. In the context of inner core convection Alboussière et al. (2010) and Deguen et al. (2013) have proposed a boundary condition at the ICB that allows for a continuous variation from perfectly permeable boundary conditions, that was considered in previous studies, to impermeable boundary conditions when the timescale for phase change is large compared to that for viscous adjustment of the topography.

In this paper, we investigate the dynamics of a growing inner core subject to electromagnetic forcing, and include the effects of a stable stratification, of the growth of the inner core, and different types of boundary conditions. We propose a systematic study of the dynamics induced by a poloidal Lorentz force in the inner core and develop scaling laws to estimate the strain rate of the flow.

In Section 2 , we develop a set of equations taking into account the Lorentz force and a buoyancy force from either thermal or compositional origin. Analytical and numerical results are presented in Section 3, scaling laws for the maximum velocity and strain rate are developed in Sections 4 and 5 and compared to numerical solutions. In Section 6, we use our results to predict the instantaneous strain rates and cumulative strain in the Earth’s inner core due to the Lorentz force.

2 Governing equations

2.1 Effect of an imposed external magnetic field

The magnetic field produced by dynamo action in the liquid outer core extends up to the surface of the Earth, but also to the center-most part of the core. Considering for example a flow velocity of the order of the growth rate of the inner core gives a magnetic Reynolds number (comparing advection and diffusion of the magnetic field) of the inner core of about . This shows that the magnetic field in the inner core is only maintained by diffusion from its boundary.

Two dynamical effects need to be taken into account: the Lorentz force and Joule heating. The Lorentz force acts directly on the momentum conservation, while Joule heating is part of the energy budget and modifies the temperature distribution, inducing a flow through buoyancy forces.

In this paper, we will discuss the effect of the Lorentz force in the case of a purely toroidal axisymmetric magnetic field with a simple mathematical form. The effect of Joule heating in the case of a non growing inner core was studied by Takehiro (2010) and will not be investigated further here.

The poloidal magnetic field intensity at the core mantle boundary (CMB) can be inferred from surface observations of the field at the Earth’s surface, but both poloidal and toroidal components are poorly known deeper in the core. The root mean square (RMS) strength of the field at the ICB has been estimated using numerical simulations to be around a few milliteslas (e.g. Glatzmaier & Roberts, 1996; Christensen & Aubert, 2006). It can be also constrained by physical observations: for example, Koot & Dumberry (2013) give an upper bound of 9-16mT for the RMS field at the ICB by looking at the dissipation in the electromagnetic coupling, while Gillet et al. (2010) suggest 2-3mT from the observation of fast toroidal oscillations in the core. Buffett (2010) obtains similar values from measurements of tidal dissipation. Numerical simulations also predict a strong azimuthal component at the vicinity of the inner core, possibly one order of magnitude higher than the vertical component (Glatzmaier & Roberts, 1996), though this depends on the magnitude of inner core differential rotation.

Buffett & Wenk (2001) have considered the effect of the azimuthal component of the Lorentz force resulting from the combination of the and components of the magnetic field. We will focus here on the effect of the azimuthal component of the magnetic field, for which the associated Lorentz force is poloidal and axisymmetric. The flow calculated by Buffett & Wenk (2001) is decoupled from the flow induced by the azimuthal component of the magnetic field, and thus the total axisymmetric flow can be obtained by simply summing the two flows.

One of the most intriguing feature of the Earth’s magnetic field is the existence of reversals. However, since the Lorentz force depends quadratically on the magnetic field, its direction is not modified by a reversal of the field. For simplicity, we will consider that the magnetic field is constant in time.

The magnetic field inside the inner core is calculated by diffusing the field from the ICB. The magnetic Reynolds number for the inner core being very small, is not advected by the flow. Because the seismic observation of anisotropy is of large scale, and also because low-order toroidal component penetrates deeper inside the inner core, only the lowest order of the azimuthal component of the magnetic field is taken into account, following the work of Karato (1999) and Buffett & Bloxham (2000).

We consider a purely toroidal axisymmetric field of degree two in the vicinity of the ICB, of the form (Buffett & Bloxham, 2000). Solving , the field inside the inner core is

(1)

in spherical coordinates, which is associated to an electric current density , where is the radius of the inner core and is the magnetic permeability.

The Lorentz force is a volume force given by . The Lorentz force can be decomposed as the sum of the gradient of a magnetic pressure and a non-potential part as , which is a unique Helmholtz decomposition for . With the magnetic field as defined in Eq. (1), we find that and are given by

(2)

and

(3)
Figure 1: Meridional cross sections showing the intensity of the magnetic field (a), the Lorentz force (b) and the non-potential part of the Lorentz force as defined in equation (3) (c).

The potential part of the Lorentz force can only promote a new equilibrium state but no persisting flow. We are thus only interested in the non potential part of the Lorentz force, shown in Fig. 1. Eq. (3) provides a characteristic scale for the force as .

Karato (1999) investigated the effect of the Maxwell stress by applying a given normal stress on the inner core boundary. This is different from our study, where, as in Buffett & Bloxham (2000), we consider a volumetric forcing, as shown on Fig. 1, and not a forcing on the surface of the inner core.

2.2 Conservation equations

2.2.1 Conservation of mass, momentum and energy

We consider an incompressible fluid in a spherical domain, with a newtonian rheology of uniform viscosity , neglecting inertia. Volume forces considered here are the buoyancy forces, with density variations due to temperature or compositional variations, and the Lorentz force as defined above.

The equations of continuity and conservation of momentum are written as

(4)
(5)

where is the velocity, the dynamic pressure that also includes the magnetic pressure, the density difference compared to the reference density profile, and the acceleration of gravity with the acceleration of gravity at .

The density depends on both the temperature and the light element concentration . We define a potential temperature as , with the isentropic temperature profile anchored to the liquidus at the ICB, and introduce a potential composition , where is the composition of the solid at the ICB. We will consider separately the effects of composition and temperature, but both can induce a density stratification, which is quantified through a variation of density which is either or , where is the reference density, and and the coefficients of thermal and compositional expansion, respectively. Because the potential temperature and composition are solutions of mathematically similar equations, we will use a quantity which is either the potential temperature or composition . In this paper, quantities that apply for both cases will have no subscript, whereas we will use for quantities referring to the thermal stratification, and for compositional stratification.

The momentum conservation equation (5) is thus written as

(6)

The equations for the evolution of the potential temperature (energy conservation) and of light element concentration (solute conservation) have a common form, which will be written as

(7)

where is the diffusivity of either heat () or composition () and a source term given by

(8)

and

(9)

As discussed in Deguen & Cardin (2011), the inner core is stably stratified when the source term is negative, and no convective instability can develop. In this paper, we will focus on this case, with either or negative.

2.2.2 Growth of the inner core

To take into account the growth of the inner core, we use a front fixing approach to solve the moving boundary problem (Crank, 1984) by scaling lengths with the inner core radius at time . We define a new coordinate system with . This modifies slightly the spatial derivatives by bringing a factor to radial derivatives, but also adds a radial advection term in the equations where the time derivative is present. In the new coordinate system, we obtain

(10)

where is the instantaneous growth rate of the inner core. Eq. (7) becomes

(11)

where and are now spatial derivative operators in the new coordinate system , with and the colatitude and longitude.

2.3 Dimensionless equations and parameters

2.3.1 Definition of the dimensionless quantities

The set of equations (4), (6), (7) is now made dimensionless, using , the age of the inner core , , and as characteristic scales for length, time, velocities, pressure, and density variations. The density scale is the difference of density across the inner core due to either thermal or compositional stratification. The quantity is scaled by . The characteristic velocity scale is defined using the diffusion time scale rather than the inner core growth rate, to make it usable in both the growing and non-growing inner core cases. The quantity is made dimensionless using . Using the same symbols for the dimensionless quantities (including using now for the dimensionless radius defined in the last subsection), we obtain

(12)
(13)
(14)

with four dimensionless parameters defined as

(15)
(16)
(17)
(18)

The last term in Eq. (14) comes from the time evolution of the scale for , .

and characterize the growth of the inner core. The Péclet number compares the apparent advection from the moving boundary to diffusion. A large Péclet number thus corresponds to a fast inner core growth compared to the diffusion rate. In the case of a non-growing inner core, , and the relevant timescale is no longer but the diffusion time scale, which gives . This approach allows us to treat both non-growing and growing cases with the same set of dimensionless parameters.

is an effective Hartmann number, which quantifies the ratio of the Lorentz force to the viscous force, using thermal diffusivity in the velocity scale. This effective Hartmann number is related to the Hartmann number often used in magnetohydrodynamics (Roberts, 2007), , through , where is the magnetic diffusivity.

defined in equation (15) is the Rayleigh number that characterizes the stratification, and is negative since is negative for a stable stratification. The density stratification depends on the importance of diffusion, and on the time-dependence of the inner core radius. Expressions for and will be given in section 2.3.2.

To solve numerically the momentum equation (13), the velocity field is decomposed into poloidal and toroidal components. The complete treatment of this equation and the expression of the Lorentz force in term of poloidal and toroidal decomposition are described in appendix A.

2.3.2 Simplified growth of the inner core

A realistic model for the inner core thermal evolution can be obtained by calculating the time evolution of the source term and the radius from the core energy balance (Labrosse, 2003, 2015), as done by Deguen & Cardin (2011). The result is sensitive to the physical properties of the core. To focus on the effect of the Lorentz forces, we choose a simpler growth scenario and assume that the inner core radius increases as the square root of time (Buffett et al., 1992). Using with the present radius of the inner core, the growth rate is thus .

This leads to the following expressions for the control parameters:

(19)
(20)
(21)
(22)

where the subscript corresponds to values for the present inner core, and is dimensionless.

The Péclet number is constant and equal to , and the parameter is proportional to . We are left with only three independent dimensionless parameters: the Rayleigh number characterizes the density stratification, the effective Hartman number the strength of the magnetic field, and the Péclet number the the relative importance of advection from the growth of the inner core and diffusion.

The value and time dependence of depends on whether a stratification of thermal or compositional origin is considered:

  • In the thermal case, the source term for thermal stratification defined in Eq. (8) can also be written as

    (23)

    where is the ratio of the Clapeyron slope to the adiabat gradient, , the Gruneisen parameter, and the isentropic bulk modulus (Deguen & Cardin, 2011). With , the product is constant, and so is .

    Solving the energy conservation equation for the potential temperature () assuming , , and constant gives

    (24)

    in dimensional form (see appendix B for the derivation). If , then the potential temperature difference across the inner core is , which corresponds to a balance between effective heating () and diffusion. In contrast, tends toward if diffusion is ineffective and . From Eq. (24), we obtain

    (25)

    and

    (26)

    With and , this gives and

    (27)
  • We estimate the density stratification due to composition from the equation of solute conservation, assuming that the outer core is well-mixed and that the partition coefficient is constant. The compositional Péclet number is large ( with a diffusivity m.s) and solute diffusion in the inner core is therefore neglected.

    The composition of the solid crystallized at time at the ICB is estimated as

    (28)

    (see Appendix C), from which the density difference across the inner core is

    (29)
    (30)

    We take advantage of the smallness of to approximate as

    (31)

    which gives

    (32)

    With and , this gives and

    (33)
Parameter Symbol Typical value Typical range Magnetic field  T  T Thermal diffusivity  m.s  m.s Chemical diffusivity  m.s  m.s Viscosity Pa.s Pa.s Age of IC Gyrs Gyrs Density stratification (thermal case)  kg.m  kg.m Density stratification (compositional case)  kg.m  kg.m Phase change timescale  yrs - yrs Inner core radius  km Acceleration of gravity (  m.s Density of the solid phase  kg.m Density difference at the ICB  kg.m Thermal expansivity  K Permeability  H.m

from Gubbins et al. (2013)

from PREM Dziewoński & Anderson (1981)

obtained using  W.m.K,  J.K.kg (Pozzo et al., 2012; Gomi et al., 2013)

assuming  K.Gyrs (Deguen & Cardin, 2011)

from Deguen & Cardin (2011)

Table 1: Typical values for the parameters used in the text, and typical range of values when useful.
Dimensionless parameter Symbol Thermal Compositional
Rayleigh number
effective Hartmann number
Péclet number
Phase change number

See the definitions of the dimensionless parameters in the text.

Table 2: Typical values of the dimensionless parameters discussed in the text for the present inner core, using typical values from table 1.

2.4 Boundary conditions

The Earth’s inner core boundary is defined by the coexistence of solid and liquid iron, at the temperature of the liquidus for the given pressure and composition. By construction, the potential temperature and concentration are both 0 at the ICB : . The mechanical boundary conditions are tangential stress-free conditions and continuity of the normal stress at the inner core boundary.

When allowing for phase change at the ICB, the condition of continuity of the normal stress gives

(34)

in dimensionless form. The parameter was introduced by Deguen et al. (2013) to characterize the phase change, and is the ratio of the phase change timescale to the viscous relaxation timescale ,

(35)

where is the density difference between liquid and solid iron at the inner core boundary. has been estimated to be  years (Alboussière et al., 2010; Deguen et al., 2013). The limit corresponds to perfectly permeable boundary conditions where the phase change occurs instantaneously, and corresponds to perfectly impermeable boundary conditions with no phase change allowed at the boundary.

With and constant, the parameter is expressed using the current value as

(36)

2.5 Numerical modelling

The code is an extension of the one used in Deguen et al. (2013), adding the effect of the magnetic forcing as in Eq. (6). The system of equations derived in appendix A in term of poloidal/toroidal decomposition is solved in axisymmetric geometry, using a spherical harmonic expansion for the horizontal dependence and a finite difference scheme in the radial direction. The nonlinear part of the advection term in the temperature (or composition) equation is evaluated in the physical space at each time step. A semi-implicit Crank-Nicholson scheme is implemented for the time evolution of the linear terms and an Adams-Bashforth procedure is used for the nonlinear advection term in the heat equation.

The boundary conditions are the same as in Deguen et al. (2013), but for most of the runs we use , which correspond to impermeable boundary conditions as discussed in section 2.4.

When keeping the inner core radius constant, the code is run until a steady state is reached. Otherwise, the code is run from to .

3 Flow description

3.1 Neutral stratification

In this subsection, we investigate the effect of the boundary conditions on the geometry and strength of the flow by solving analytically the set of equations in the case of neutral stratification. The analytical solution for neutral stratification has also been used to benchmark the code for .

In the case of neutral stratification, with , the equations for the temperature or composition perturbation (14) and momentum conservation (37) are no longer coupled. The diffusivity is no longer relevant and the problem does not depend on the Péclet number. Eq. (37) is solved in appendix D using the boundary conditions presented in the previous section. The flow velocity is found to be proportional to the effective Hartman number times a sigmoid function of . The solution is shown on Fig. 2, with dimensionless maximum horizontal velocity and root mean square velocity as functions of the phase change number , as well as streamlines corresponding to the two extreme cases, (fully permeable boundary conditions) and (impermeable boundary conditions).

In the limit , corresponding to permeable boundary conditions, the streamlines of the flow cross the ICB, which indicates significant melting and freezing at the ICB. In contrast, the streamlines in the limit are closed lines which do not cross the ICB, which indicates negligible melting or freezing at the ICB. The velocity is proportional to the effective Hartmann number whereas the dependence is more complex. The velocities reaches two asymptotic values for low and large values, separated by a sharp kink. The discontinuity in the derivative of the maximum horizontal velocity slightly above corresponds to a change of the spatial position of the maximum, when the streamlines become closed and the maximal horizontal velocity is obtained at the top of the cell and no longer at its bottom. The change of behavior of the boundary from permeable to impermeable induces a significant decrease of the strength of the flow, since the velocity magnitude in the regime is one order of magnitude smaller than when permeable boundary conditions () are assumed.

Fig. 2.b shows the maximum of the velocity, now given in m.s, as a function of the viscosity, using typical values of the parameters given in Table 1 and five different values for the phase change timescale , from zero to infinite. In term of dimensionless parameters, a high viscosity corresponds to small Hartmann number and phase change number . Changing the timescale translates the position of the transition between the two regimes, the viscosity value corresponding to the transition being proportional to , but does not change the general trend of the curve, which is a linear decrease of the velocity magnitude in log-log space, except for the kink between the two regimes. The linear decrease is due to the viscosity dependence of the Hartmann number . For typical values of the phase change timescale between 100 years and years, the kink between the two regimes occurs at a viscosity in the range Pa.s.

In what follows, we will focus on the conditions which are the most favorable to deformation due to the poloidal component of the Lorentz force, and therefore focus on the case of low viscosity and large . The ICB would act as a permeable boundary only if (see Fig. 2), corresponding to Pa.s. Under these conditions, the typical flow velocity would be m.s, i.e. two orders of magnitude or more smaller than the inner core growth rate, and would be unlikely to result in significant texturing. For this reason, we will let aside the high viscosity/low regimes to focus on low viscosity/high cases, for which the ICB acts as an impermeable boundary. This gives boundary conditions very different from previous studies, where perfectly permeable boundary conditions were assumed (Karato, 1999; Buffett & Bloxham, 2000). In particular, this implies that the flow velocity estimated by (Karato, 1999) was overestimated by one order of magnitude.

Figure 2: Analytic solution for . (a) Evolution of the dimensionless velocity as a function of the phase change number , with streamlines for (left) and (right). The RMS velocity and the maximum of the horizontal velocity are plotted. (b) Evolution of the RMS velocity as a function of , with velocity in m.s. Except for the viscosity and the phase change time scale , the parameters used for definition of and are given in Table 1. The kink in the curves corresponds to the change in regime between large (low viscosity) and low (large viscosity), and the corresponding viscosity value is a function of the phase change timescale .

According to Eq. (36), the parameter varies linearly with time, which means that must have been small early in inner core’s history. However, this is true for a very short time, when the inner core radius was very small, of the order , and this episode is unlikely to have observable consequences in the present structure of the inner core.

3.2 Zero growth rate

Figure 3: Snapshots of meridional cross section of the temperature and the vorticity fields for and a constant inner core radius, for four different values of the Rayleigh number (from top right, going clockwise: , , , ). When the stratification is large enough (), the flow is confined at the top of the inner core and the temperature field has a spherical symmetry. When the stratification is weak (), the flow is similar to the one in Fig. 2 for and the temperature is almost uniform. The vorticity is scaled by and the temperature by . For , reduces to .

We first investigate the effect of the Lorentz force without taking into account the secular growth of the inner core (). Fig. 3 shows the vorticity and temperature fields obtained for different values of the Rayleigh number, at a given effective Hartmann number , for a thermally stratified inner core.

When the Rayleigh number is small, the vorticity field is organized in two symmetric tores wrapped around the N-S axis.The stratification is too weak to alter the flow induced by the Lorentz force and the temperature field is advected and mixed by the flow. The velocity field is equal to the one calculated analytically for (see subsection 3.1 and appendix D).

However, when the Rayleigh number is larger, the flow is altered by the stratification and is confined in an uppermost layer, as found by Buffett & Bloxham (2000). The velocity is smaller than in the case of neutral stratification. The temperature field is strongly stratified and the perturbations due to radial advection are small. The flow obtained here is similar to the flow induced by differential inner core growth with a stable stratification (Deguen et al., 2011), with a notable difference: we impose a large implying a near zero radial flow across the ICB, whereas Deguen et al. (2011) impose a given as the driving force. The confinement of the flow in a thin layer is likely to concentrate the deformation and thus we may expect higher strain rates for a highly stratified inner core, but a different spatial distribution of the deformation.

Figure 4: Maximum velocity (normalized by the Hartmann number ) in the upper vorticity layer, for a zero growth rate. The velocity is scaled by the diffusion velocity . The maximum size of the dots corresponds to the value for computed analytically. For some values of , the vorticity field is plotted in the meridional cross section. The red line with a slope of shows the limit between the two regimes.

To explore the parameter space in terms of Rayleigh and effective Hartmann numbers, we computed runs with Rayleigh numbers from to and effective Hartman number from to . The maximum velocity (normalized by ) is plotted in Fig. 4 as a proxy to determine the regime. The largest velocity coincides with the flow velocity obtained for neutral stratification. The vorticity field corresponding to some of the points in the regime diagram are also shown in Fig. 4.

The systematic exploration of the parameter space reveals two different dynamical regimes, which domains of existence in a (,) space are shown in Fig. 4. In the upper left part of the diagram (large effective Hartmann number, low Rayleigh number), the flow is very similar (qualitatively and quantitatively) to the analytical solution for a neutral stratification, and deformation extends deep in the inner core. This regime is characterized by a negligible effect of the buoyancy forces, and will therefore be referred to as the weakly stratified regime. In the lower right part, the flow is confined in a shallow layer which thickness depends on the Rayleigh number only (not on ) and in which the velocity is smaller than for neutral stratification. This regime will be referred to as the strongly stratified regime.

3.3 Growing inner core

To investigate the effect of inner core growth, we compute several runs with a given set of parameters , with the time between and . Unlike in subsections 3.1 and 3.2, the dimensionless numbers evolve with time, as described by Eqs (19) to (22).

Fig. 5 shows the evolution of the vorticity field in six simulations, for a thermal stratification, with the same Rayleigh and effective Hartman numbers, , , but different values of the Péclet number, which corresponds to increasing diffusivity from left to right. For each run, snapshots of the vorticity field corresponding to four time steps are shown, from top-right and going clockwise.

Fig. 5 shows that the thickness of the upper layer increases with increasing Péclet number. The transition between the two regimes of strong and weak stratification is shifted toward larger Rayleigh numbers when the Péclet number is increased. At low or moderate Péclet numbers ( in the cases presented here), the magnitude of vorticity is almost constant time, implying that the deformation rate in the uppermost layer is also constant.

Figure 5: Snapshots of the vorticity field for simulations with dimensionless parameters , , and , 1, 10, , and (from left to right), with . Each panel corresponds to one simulation, with four time steps represented : , , and dimensionless time, from top-right and going clockwise. See Fig. 8 for strain rates of corresponding runs.

4 Scaling laws

In this section, we determine scaling laws for the thickness of the shallow shear layer and the maximum velocity in the layer in the strongly stratified regime from the set of equations developed in section 2. We will first discuss the transition between the strongly stratified and weakly stratified regimes discussed in Fig. 4. We will then focus on the strongly stratified regime and estimate the deformation in the uppermost layer. Thermal and compositional stratification are discussed separately. The flow in the weakly stratified regime is given by the analytical model discussed in section 3.1 and appendix D for neutral stratification.

4.1 Balance between magnetic forcing and stratification

We start here by discussing the transition between the strongly-stratified and weakly-stratified regimes. We base our analysis on the vorticity equation obtained by taking the curl of the momentum conservation equation (Eq. (13)),

(37)

where is the vorticity. The quantity (denoting either potential temperature or composition) is split into two parts, , where is the reference radial profile corresponding to . The vorticity equation then writes

(38)

In the vorticity equation, the three terms must balance if the effect of stratification is important. Starting from a state with no perturbations, , the flow velocity is initially set by a balance between the Lorentz force and the viscosity force. Isosurfaces of are deformed by the resulting flow, and the buoyancy force increases, eventually balancing the magnetic force if the stratification is strong enough. In this case, further radial motion is prevented and the flow tends to be localized in a layer below the ICB, as found in our numerical simulations. Denoting by the thickness of the shear layer and and the horizontal and vertical velocity, respectively, the vorticity is . We therefore have, from Eq. (38),

(39)

The perturbation thus scales as

(40)

if the stratification is strong enough for the induced buoyancy forces to balance the Lorentz force.

The effect of the stratification is negligible if the buoyancy forces, which are , cannot balance the Lorentz force, which is . Since is necessarily smaller than , which by construction is equal to 1, the effect of the stratification will therefore be negligible if . This is consistent with the boundary between the two regimes found from our numerical calculations, as shown in Fig. 4, as well as with the results of Buffett & Bloxham (2000) who found that the Lorentz force can displace isodensity surfaces by . This estimate is valid for both a growing or non-growing inner core.

4.2 Scaling laws in the strongly stratified regime

The reference profile is solution of

(41)

Subtracting Eqs (41) to (14), and assuming that , we obtain

(42)

Three of these terms depend on the growth rate: , , and . With our assumption of , we have and thus . Thus, the largest term among the growth rate-dependent terms is .

Comparing the effect of the diffusion term, which is , with the inner core growth term, which is , we find that the effect of the inner core growth is negligible if

(43)

This suggest the existence of two different regimes depending on whether is small or large. We develop below scaling laws for these two cases.

4.2.1 Small limit

Neglecting the growth terms, we have

(44)

The conservation of mass implies that , and with , we obtain

(45)

We now assume that the advection of the perturbation is small compared to the vertical advection of the reference state, which requires that . Balancing the advection and diffusion terms, we obtain

(46)

Combining this expression with the relation obtained from the vorticity equation (Eq. (39)), we have

(47)
(48)

4.2.2 Large limit

In this limit, the diffusion time is larger than the age of the inner core, which allows us to neglect the diffusion term. Keeping only the largest growth rate-dependent term, and using Eq. (39), we have

(49)

Assuming again that the advection of the perturbation is small compared to the vertical advection of the reference state, the main balance is between the second and third terms, which gives

(50)

Combined with from Eq. (39), we find that the thickness and maximum velocity of the upper layer are

(51)
(52)

Two conditions have to be fulfilled for these scaling laws to be valid. First, we must have , which is using Eq. (51). Also, we have assumed that the upper layer is thin () and that the horizontal advection is small compared to the vertical one.

4.2.3 Time-dependence

Our derivation does not make any assumption on the form of the inner core growth , and the scaling law validity should not be restricted to the case assumed in the numerical simulations.

These scaling laws are valid at all time during the growth of the inner core, provided that and and that the time dependence of the control parameters is properly taken into account. For example, under the assumption of and with , and given by Eqs. (20), (27) and (22), we obtain

(53)

for the thermal case (small ), and

(54)

for the compositional case (large ).

4.3 Comparison with numerical results

Fig. 6 shows the thickness of the uppermost vorticity layer and the maximum horizontal velocity obtained in numerical simulations with a constant inner core radius, which corresponds to the small Péclet number limit. When , the flow has the geometry and amplitude predicted by our analytical model for . When and the thickness of the upper layer is smaller than , the data points align on straight lines in log-log scale, with slopes close to the predictions of the scaling laws developed in subsection 4.2.1 (Eqs. (47) and (48)).

Figure 6: Results from simulations with a constant inner core radius. Evolution of the thickness of the uppermost vorticity layer (a) and maximum horizontal velocity (b) with the absolute value of the Rayleigh number . Colors correspond to different values of the effective Hartmann number , and dashed lines to the corresponding line. The velocity has been scaled with the effective Hartman number and the extreme value for low corresponds to the analytical model with no stratification (black horizontal dotted line). The solid black lines are the best fit for , and . The orange lines show the slopes predicted in subsection 4.2.1.
Figure 7: Thickness of the uppermost vorticity layer (a) and maximum velocity (b) as functions of . 50 runs are plotted, with from to and from to 5000, with 10 time steps for each runs, and filtered by . Solid green lines are the best fit for , which is and . The red lines are expected scaling laws developed in subsection 4.2.2.

Fig. 7 shows the vorticity layer thickness and maximum velocity as functions of , in log-log scale, for runs in the large Péclet limit. The thickness of the upper layer and the maximum velocity align on slopes close to the 1/5 and 3/5 slopes predicted in subsection 4.2.2. Fig. 7 has been constructed from runs with , but we have checked that, as long as the condition is verified, the geometry of the flow does not depend on and that the velocity is proportional to .

5 Strain rate

Fig. 8 shows the von Mises equivalent strain rate for the runs corresponding to Fig. 5, highlighting regions of high deformation. The von Mises equivalent strain rate is the second invariant of the strain rate tensor, measuring the power dissipated by deformation (Tome et al., 1984; Wenk et al., 2000; Tackley, 2000). Comparing Figs 5 and 8, we see that the deformation and vorticity fields have a similar geometry when the flow is organized in several layers, whereas the location of the regions of high deformation and high vorticity differ when the effect of stratification is small and the flow is organized in one cell only. In this case, which is similar to that studied by Karato (1999), the maximum deformation is at the edges of the cells, whereas in the case of large stratification, the strain is confined in the uppermost layer. In the strongly stratified regime, the deformation can be predicted from the scaling laws discussed in section 4.2, as in dimensionless form, or in dimensional form.

In the small number case, relevant for the Earth’s inner core with a thermal stratification, using the scaling laws (47) and (48) for the velocity and shear layer thickness gives

(55)

In the large number case, relevant for the Earth’s inner core with a compositional stratification, using the scaling laws (51) and (52) for the velocity and shear layer thickness gives

(56)

Notice that is proportional to and in the thermal and compositional cases, and therefore increases with decreasing viscosity in spite of the fact that the velocity in the shear layer decreases with decreasing . This is because the thickness of the shear layer decreases with viscosity faster than the velocity.

These scaling laws give upper bounds for the actual strain rate in the inner core, which evolves with time because of the time dependence of the parameters. The quantity is the time needed to deform the shear layer to a cumulated strain .

To estimate the cumulated deformation in the inner core, we assume that the strain rate found above is applied only on the uppermost shear layer of thickness . The simplest is to assume that both and are evolving slowly with time, on a timescale long compared to the time over which a layer of thickness is crystallized. Assuming that the strain rate is given by within the shear layer of thickness and is negligible elsewhere, the cumulated deformation is then

(57)

with the dimensionless maximum horizontal velocity given by Eqs. (48) or (52) depending on the value of the Péclet number. With dimensional, the cumulated strain can also be written as

(58)

The deformation magnitude below the upper shear layer is given by the ratio between the horizontal velocity induced by the Lorentz force and the growth rate of the inner core.

A more elaborated method of estimating is discussed in appendix E. The results are close to what Eq. (57) predicts for , but are more accurate for smaller . The validity of both estimates is restricted to conditions under which the strong stratification scaling laws applies, which requires that . This is only verified if the inner core radius is larger than in the thermal case, and in the compositional case.

6 Application to the inner core

To determine in which regime is the inner core, the first step is to estimate the ratio