Formulation of the (magneto)-hydrodynamics problem

The general equations describing thermal convection and dynamo action of a rotating compressible fluid are the starting point from which the Boussinesq or the anelastic approximations are developed. In MagIC, we consider a spherical shell rotating about the vertical $z$ axis with a constant angular velocity $\mathrm{\Omega }$. Equations are solve in the corotating system.

The conservation of momentum is formulated by the Navier-Stokes equation

(1)$$\rho ({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}\overrightarrow{u})=-\overrightarrow{\mathrm{\nabla }}p+{\displaystyle \frac{1}{{\mu }_0}}(\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B})\times \overrightarrow{B}+\rho \overrightarrow{g}-2\rho \overrightarrow{\mathrm{\Omega }}\times \overrightarrow{u}+\overrightarrow{\mathrm{\nabla }}\cdot \mathsf{S},$$

where $\overrightarrow{u}$ is the velocity field, $\overrightarrow{B}$ the magnetic field, and $p$ a modified pressure that includes centrifugal forces. $\mathsf{S}$ corresponds to the rate-of-strain tensor given by:

$$\begin{array}{r}\begin{array}{rl}{S}_{ij}& =2\nu \rho [e_{ij}-{\displaystyle \frac{1}{3}}{\delta }_{ij}\overrightarrow{\mathrm{\nabla }}\cdot \overrightarrow{u}],\\ e_{ij}& ={\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}).\end{array}\end{array}$$

Convection is driven by buoyancy forces acting on variations in density $\rho$. These variations have a dynamical part formulated by the continuity equation describing the conservation of mass:

(2)$${\displaystyle \frac{\partial \rho }{\partial t}}=-\overrightarrow{\mathrm{\nabla }}\cdot \rho \overrightarrow{u}.$$

In addition an equation of state is required to formulate the thermodynamic density changes. For example the relation

(3)$${\displaystyle \frac{1}{\rho }}\partial \rho =-\alpha \partial T+\beta \partial p-\delta \partial \xi$$

describes density variations caused by variations in temperature $T$, pressure $p$, and composition $\xi$. The latter contribution needs to be considered for describing the effects of light elements released from a growing solid iron core in a so-called double diffusive approach.

To close the system we also have to formulate the dynamic changes of entropy, pressure, and composition. The evolution equation for pressure can be derived from the Navier-Stokes equation, as will be further discussed below. For entropy variations we use the so-called energy or heat equation

(4)$$\rho T({\displaystyle \frac{\partial s}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}s)=\overrightarrow{\mathrm{\nabla }}\cdot (k\overrightarrow{\mathrm{\nabla }}T)+{\mathrm{\Phi }}_{\nu }+{\displaystyle \frac{\lambda }{{\mu }_0}}{(\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B})}^2+ϵ,$$

where ${\mathrm{\Phi }}_{\nu }$ corresponds to the viscous heating expressed by

$${\mathrm{\Phi }}_{\nu }=2\rho [e_{ij}{e}_{ji}-{\displaystyle \frac{1}{3}}{(\overrightarrow{\mathrm{\nabla }}\cdot \overrightarrow{u})}^2]$$

Note that we use here the summation convention over the indices $i$ and $j$. The second last term on the right hand side is the Ohmic heating due to electric currents. The last term is a volumetric sink or source term that can describe various effects, for example radiogenic heating, the mixing-in of the light elements or, when radially dependent, potential variations in the adiabatic gradient (see below). For chemical composition, we finally use

(5)$$\rho ({\displaystyle \frac{\partial \xi }{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}\xi )=\overrightarrow{\mathrm{\nabla }}\cdot (k_{\xi }\overrightarrow{\mathrm{\nabla }}\xi )+{ϵ}_{\xi },$$

The induction equation is obtained from the Maxwell equations (ignoring displacement current) and Ohm’s law (neglecting Hall effect):

(6)$${\displaystyle \frac{\partial \overrightarrow{B}}{\partial t}}=\overrightarrow{\mathrm{\nabla }}\times (\overrightarrow{u}\times \overrightarrow{B}-\lambda \overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B}).$$

When the magnetic diffusivity $\lambda$ is homogeneous this simplifies to the commonly used form

(7)$${\displaystyle \frac{\partial \overrightarrow{B}}{\partial t}}=\overrightarrow{\mathrm{\nabla }}\times (\overrightarrow{u}\times \overrightarrow{B})+\lambda \overrightarrow{\mathrm{\Delta }}\overrightarrow{B}.$$

The physical properties determining above equations are rotation rate $\mathrm{\Omega }$, the kinematic viscosity $\nu$, the magnetic permeability ${\mu }_0$, gravity $\overrightarrow{g}$, thermal conductivity $k$, Fick’s conductiviy $k_{\xi }$, magnetic diffusivity $\lambda$. The latter connects to the electrical conductivity $\sigma$ via $\lambda =1/({\mu }_0\sigma )$. The thermodynamics properties appearing in (3) are the thermal expansivity at constant pressure (and composition)

(8)$$\alpha =-{\displaystyle \frac{1}{\rho }}{\left({\displaystyle \frac{\partial \rho }{\partial T}}\right)}_{p,\xi },$$

the compressibility at constant temperature

$$\beta ={\displaystyle \frac{1}{\rho }}{\left({\displaystyle \frac{\partial \rho }{\partial p}}\right)}_{T,\xi }$$

and an equivalent parameter $\delta$ for the dependence of density on composition:

(9)$$\delta =-{\displaystyle \frac{1}{\rho }}{\left({\displaystyle \frac{\partial \rho }{\partial \xi }}\right)}_{p,T},$$

Sketch of the spherical shell model and its system of coordinate.

The reference state

The convective flow and the related processes including magnetic field generation constitute only small disturbances around a background or reference state. In the following we denote the background state with a tilde and the disturbance we are interested in with a prime. Formally we will solve equations in first order of a smallness parameters $ϵ$ which quantified the ratio of convective disturbances to background state:

(10)$$ϵ\sim {\displaystyle \frac{T^{\prime }}{\tilde{T}}}\sim {\displaystyle \frac{p^{\prime }}{\tilde{p}}}\sim {\displaystyle \frac{{\rho }^{\prime }}{\tilde{\rho }}}\sim ...\ll 1.$$

The background state is hydrostatic, i.e. obeys the simple force balance

(11)$$\mathrm{\nabla }\tilde{p}=\tilde{\rho }\tilde{\overrightarrow{g}}.$$

Convective motions are supposed to be strong enough to provide homogeneous entropy (and composition). The reference state is thus adiabatic and its gradients can be expressed in terms of the pressure gradient (11):

(12)$${\displaystyle \frac{\mathrm{\nabla }\tilde{T}}{\tilde{T}}}={\displaystyle \frac{1}{\tilde{T}}}{\left({\displaystyle \frac{\partial T}{\partial p}}\right)}_s\mathrm{\nabla }p={\displaystyle \frac{\alpha }{c_p}}\tilde{\overrightarrow{g}},$$

(13)$${\displaystyle \frac{\mathrm{\nabla }\tilde{\rho }}{\tilde{\rho }}}={\displaystyle \frac{1}{\tilde{\rho }}}{\left({\displaystyle \frac{\partial \rho }{\partial p}}\right)}_s\mathrm{\nabla }p=\beta \tilde{\rho }\tilde{\overrightarrow{g}}.$$

The reference state obviously dependence only on radius. Dimensionless numbers quantifying the temperature and density gradients are called dissipation number $Di$ and compressibility parameter $Co$ respectively:

$$Di={\displaystyle \frac{\alpha d}{c_p}}\tilde{g},$$

and

$$Co=d\beta \tilde{\rho }\tilde{g}.$$

Here $d$ is a typical length scale, for example the shell thickness of the problem. The dissipation number is something like an inverse temperature scale hight while the compressibility parameters is an inverse density scale hight. The ratio of both numbers also helps to quantify the relative impact of temperature and pressure on density variations:

(14)$${\displaystyle \frac{\alpha \mathrm{\nabla }T}{\beta \mathrm{\nabla }\rho }}\approx \alpha \tilde{T}{\displaystyle \frac{Di}{Co}}.$$

As an example we demonstrate how to derive the first order continuity equation here. Using $\rho =\tilde{\rho }+{\rho }^{\prime }$ in (2) leads to

$${\displaystyle \frac{\partial \tilde{\rho }}{\partial t}}+{\displaystyle \frac{\partial {\rho }^{\prime }}{\partial t}}=-\overrightarrow{\mathrm{\nabla }}\cdot \left(\tilde{\rho }\overrightarrow{u}\right)-\overrightarrow{\mathrm{\nabla }}\cdot \left({\rho }^{\prime }\overrightarrow{u}\right).$$

The zero order term vanishes since the background density is considered static (or actually changing very slowly on very long time scales). The second term in the right hand side is obviously of second order. The ratio of the remaining two terms can be estimated to also be of first order in $ϵ$, meaning that the time derivative of $\rho$ is actually also of second order:

$${\displaystyle \frac{\left[\partial \rho /\partial t\right]}{[\overrightarrow{\mathrm{\nabla }}\cdot \rho \overrightarrow{u}]}}\approx {\displaystyle \frac{{\rho }^{\prime }}{\tilde{\rho }}}\approx ϵ.$$

Square brackets denote order of magnitude estimates here. We have used the fact that the reference state is static and assume time scale of changes are comparable (or slower) ${\rho }^{\prime }$ than the time scales represented by $u$ and that length scales associated to the gradient operator are not too small. We can then neglect local variations in ${\rho }^{\prime }$ which means that sound waves are filtered out. This first order continuity equation thus simply reads:

(15)$$\overrightarrow{\mathrm{\nabla }}\cdot \left(\tilde{\rho }\overrightarrow{u}\right)=0.$$

This defines the so-called anelastic approximation where sound waves are filtered out by neglecting the local time derivative of density. This approximation is justified when typical velocities are sufficiently smaller than the speed of sound.

Boussinesq approximation

For Earth the dissipation number and the compressibility parameter are around $0.2$ when temperature and density jump over the whole liquid core are considered. This motivates the so called Boussinesq approximation where $Di$ and $Co$ are assumed to vanish. The continuity equation (2) then simplifies further:

$${\displaystyle \frac{1}{\tilde{\rho }}}\overrightarrow{\mathrm{\nabla }}\cdot \tilde{\rho }\overrightarrow{u}={\displaystyle \frac{\overrightarrow{u}}{\tilde{\rho }}}\cdot \mathrm{\nabla }\tilde{\rho }+\mathrm{\nabla }\cdot \overrightarrow{u}\approx \mathrm{\nabla }\cdot \overrightarrow{u}=0.$$

When using typical number for Earth, (14) becomes $0.05$ so that pressure effects on density may be neglected. The first order Navier-Stokes equation (after to zero order hydrostatic reference solution has been subtracted) then reads:

(16)$$\tilde{\rho }({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}\overrightarrow{u})=-\overrightarrow{\mathrm{\nabla }}{p}^{\prime }-2\rho \overrightarrow{\mathrm{\Omega }}\times \overrightarrow{u}+\alpha {\tilde{g}}_o{T}^{\prime }{\displaystyle \frac{\overrightarrow{r}}{r_o}}+\delta {\tilde{g}}_o{\xi }^{\prime }{\displaystyle \frac{\overrightarrow{r}}{r_o}}+{\displaystyle \frac{1}{{\mu }_0}}(\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B})\times \overrightarrow{B}+\tilde{\rho }\nu \mathrm{\Delta }\overrightarrow{u}.$$

Here $u$ and $B$ are understood as first order disturbances and $p^{\prime }$ is the first order non-hydrostatic pressure and $T^{\prime }$ the super-adiabatic temperature and $\xi$ the super-adiabatic chemical composition. Above we have adopted a simplification of the buoyancy term. In the Boussinesq limit with vanishing $Co$ and a small density difference between a solid inner and a liquid outer core a linear gravity dependence provides a reasonable approximation:

$$\tilde{\overrightarrow{g}}={\tilde{g}}_o{\displaystyle \frac{\overrightarrow{r}}{r_o}},$$

where we have chosen the gravity ${\tilde{g}}_o$ at the outer boundary radius $r_o$ as reference.

The first order energy equation becomes

(17)$$\tilde{\rho }({\displaystyle \frac{\partial T^{\prime }}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}{T}^{\prime })=\kappa \mathrm{\Delta }{T}^{\prime }+ϵ,$$

where we have assumed a homogeneous $k$ and neglected viscous and Ohmic heating which can be shown to scale with $Di$ as we discuss below. Furthermore, we have used the simple relation

$$\partial s\approx {\displaystyle \frac{\tilde{\rho }{c}_p}{\tilde{T}}}\partial T,$$

defined the thermal diffusivity

$$\kappa ={\displaystyle \frac{k}{\tilde{\rho }{c}_p}},$$

and adjusted the definition of $ϵ$. Finally the first order equation for chemical composition becomes

(18)$$\tilde{\rho }({\displaystyle \frac{\partial {\xi }^{\prime }}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}{\xi }^{\prime })={\kappa }_{\xi }\mathrm{\Delta }{\xi }^{\prime }+{ϵ}_{\xi },$$

where we have assumed a homogeneous $k_{\xi }$ and adjusted the definition of ${ϵ}_{\xi }$.

MagIC solves a dimensionless form of the differential equations. Time is scaled in units of the viscous diffusion time $d^2/\nu$, length in units of the shell thickness $d$, temperature in units of the temperature drop $\mathrm{\Delta }T=T_o-T_i$ over the shell, composition in units of the composition drop $\mathrm{\Delta }\xi ={\xi }_o-{\xi }_i$ over the shell and magnetic field in units $(\mu \lambda \tilde{\rho }\mathrm{\Omega }{)}^{1/2}$. Technically the transition to the dimensionless form is achieved by the substitution

$$r\to rd,t\to (d^2/\nu )t,T\to \mathrm{\Delta }TT,\xi \to \mathrm{\Delta }\xi \xi ,B\to {\left(\mu \lambda \tilde{\rho }\mathrm{\Omega }\right)}^{1/2}B$$

where $r$ stands for any length. The next step then is to collect the physical properties as a few possible characteristic dimensionless numbers. Note that many different scalings and combinations of dimensionless numbers are possible. For the Navier-Stokes equation in the Boussinesq limit MagIC uses the form:

(19)$$({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}\overrightarrow{u})=-\overrightarrow{\mathrm{\nabla }}{p}^{\prime }-{\displaystyle \frac{2}{E}}\overrightarrow{e_z}\times \overrightarrow{u}+{\displaystyle \frac{Ra}{Pr}}{T}^{\prime }{\displaystyle \frac{\overrightarrow{r}}{r_o}}+{\displaystyle \frac{R{a}_{\xi }}{Sc}}{\xi }^{\prime }{\displaystyle \frac{\overrightarrow{r}}{r_o}}+{\displaystyle \frac{1}{EPm}}(\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B})\times \overrightarrow{B}+\mathrm{\Delta }\overrightarrow{u},$$

where ${\overrightarrow{e}}_z$ is the unit vector in the direction of the rotation axis and the meaning of the pressure disturbance $p^{\prime }$ has been adjusted to the new dimensionless equation form.

Anelastic approximation

The anelastic approximation adopts the simplified continuity (15). The background state can be specified in different ways, for example by providing profiles based on internal models and/or ab initio simulations. We will assume a polytropic ideal gas in the following.

Analytical solution in the limit of an ideal gas

In the limit of an ideal gas which follows $\tilde{p}=\tilde{\rho }\tilde{T}$ and has $\alpha =1/\tilde{T}$, one directly gets:

$$\begin{array}{r}\begin{array}{rl}{\displaystyle \frac{d\tilde{T}}{dr}}& =-Di\tilde{g}(r),\\ \tilde{\rho }& ={\tilde{T}}^{1/(\gamma -1)},\end{array}\end{array}$$

where $\gamma =c_p/c_v$. Note that we have moved to a dimensionless formulations here, where all quantities have been normalized with their outer boundary values and $Di$ refers to the respective outer boundary value. If we in addition make the assumption of a centrally-condensed mass in the center of the spherical shell of radius $r\in [r_i,r_o]$, i.e. $g\propto ,1/r^2$, this leads to

$$\begin{array}{r}\begin{array}{rl}\tilde{T}(r)& =Di\frac{r_o^2}{r}+(1-Di{r}_o),\\ \tilde{\rho }(r)& ={\tilde{T}}^m,\\ Di& ={\displaystyle \frac{r_i}{r_o}}(\exp {\displaystyle \frac{N_{\rho }}{m}}-1),\end{array}\end{array}$$

where $N_{\rho }=\ln ({\tilde{\rho }}_i/{\tilde{\rho }}_o)$ is the number of density scale heights of the reference state and $m=1/(\gamma -1)$ is the polytropic index.

Warning

• The relationship between $N_{\rho }$ and the dissipation number $Di$ directly depends on the gravity profile. The formula above is only valid when $g\propto 1/r^2$.

• In this formulation, when you change the polytropic index $m$, you also change the nature of the fluid you’re modelling since you accordingly modify $\gamma =c_p/c_v$.

Anelastic MHD equations

In the most general formulation, all physical properties defining the background state may vary with depth. Specific reference values must then be chosen to provide a unique dimensionless formulations and we typically chose outer boundary values here. The exception is the magnetic diffusivity where we adopt the inner boundary value instead. The motivation is twofold: (i) it allows an easier control of the possible continuous conductivity value in the inner core; (ii) it is a more natural choice when modelling gas giants planets which exhibit a strong electrical conductivity decay in the outer layer.

The time scale is then the viscous diffusion time $d^2/{\nu }_o$ where ${\nu }_o$ is the kinematic viscosity at the outer boundary. Magnetic field is expressed in units of $({\rho }_o{\mu }_0{\lambda }_i\mathrm{\Omega }{)}^{1/2}$, where ${\rho }_o$ is the density at the outer boundary and ${\lambda }_i$ is the magnetic diffusivity at the inner boundary.

This leads to the following sets of dimensionless equations:

(20)$$({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}\overrightarrow{u})=-\overrightarrow{\mathrm{\nabla }}\left({\displaystyle \frac{p^{\prime }}{\tilde{\rho }}}\right)-{\displaystyle \frac{2}{E}}\overrightarrow{e_z}\times \overrightarrow{u}+{\displaystyle \frac{Ra}{Pr}}\tilde{g}{s}^{\prime }\overrightarrow{e_r}+{\displaystyle \frac{R{a}_{\xi }}{Sc}}\tilde{g}{\xi }^{\prime }\overrightarrow{e_r}+{\displaystyle \frac{1}{PmE\tilde{\rho }}}(\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B})\times \overrightarrow{B}+{\displaystyle \frac{1}{\tilde{\rho }}}\overrightarrow{\mathrm{\nabla }}\cdot \mathsf{S},$$

(21)$$\overrightarrow{\mathrm{\nabla }}\cdot \tilde{\rho }\overrightarrow{u}=0,$$

(22)$$\overrightarrow{\mathrm{\nabla }}\cdot \overrightarrow{B}=0,$$

(23)$$\tilde{\rho }({\displaystyle \frac{\partial {\xi }^{\prime }}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}{\xi }^{\prime })={\displaystyle \frac{1}{Sc}}\overrightarrow{\mathrm{\nabla }}\cdot ({\kappa }_{\xi }(r)\tilde{\rho }\overrightarrow{\mathrm{\nabla }}{\xi }^{\prime })$$

(24)$${\displaystyle \frac{\partial \overrightarrow{B}}{\partial t}}=\overrightarrow{\mathrm{\nabla }}\times (\overrightarrow{u}\times \overrightarrow{B})-{\displaystyle \frac{1}{Pm}}\overrightarrow{\mathrm{\nabla }}\times (\lambda (r)\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B}).$$

Here $\tilde{g}$ and $\tilde{\rho }$ are the normalized radial gravity and density profiles that reach one at the outer boundary.

Entropy equation and turbulent diffusion

The entropy equation usually requires an additional assumption in most of the existing anelastic approximations. Indeed, if one simply expands Eq. (4) with the classical temperature diffusion an operator of the form:

$$ϵ\overrightarrow{\mathrm{\nabla }}\cdot \left(K\overrightarrow{\mathrm{\nabla }}{T}^{\prime }\right)+\overrightarrow{\mathrm{\nabla }}\cdot \left(K\overrightarrow{\mathrm{\nabla }}\tilde{T}\right),$$

will remain the right-hand side of the equation. At first glance, there seems to be a $1/ϵ$ factor between the first term and the second one, which would suggest to keep only the second term in this expansion. However, for astrophysical objects which exhibit strong convective driving (and hence large Rayleigh numbers), the diffusion of the adiabatic background is actually very small and may be comparable or even smaller in magnitude than the $ϵ$ terms representing the usual convective perturbations. For the Earth core for instance, if one assumes that the typical temperature fluctuations are of the order of 1 mK and the temperature contrast between the inner and outer core is of the order of 1000 K, then $ϵ\sim {10}^{-6}$. The ratio of the two terms can thus be estimated as

(25)$$ϵ{\displaystyle \frac{T^{\prime }/{\delta }^2}{T/d^2}},$$

where $d$ is the thickness of the inner core and $\delta$ is the typical thermal boundary layer thickness. This ratio is exactly one when $\delta =1\text{ m}$, a plausible value for the Earth inner core.

In numerical simulations however, the over-estimated diffusivities restrict the computational capabilities to much lower Rayleigh numbers. As a consequence, the actual boundary layers in a global DNS will be much thicker and the ratio (25) will be much smaller than unity. The second terms will thus effectively acts as a radial-dependent heat source or sink that will drive or hinder convection. This is one of the physical motivation to rather introduce a turbulent diffusivity that will be approximated by

$$\kappa \tilde{\rho }\tilde{T}\overrightarrow{\mathrm{\nabla }}s,$$

where $\kappa$ is the turbulent diffusivity. Entropy diffusion is assumed to dominate over temperature diffusion in turbulent flows.

The choice of the entropy scale to non-dimensionalize Eq. (4) also depends on the nature of the boundary conditions: it can be simply the entropy contrast over the layer $\mathrm{\Delta }s$ when the entropy is held constant at both boundaries, or $d(ds/dr)$ when flux-based boundary conditions are employed. We will restrict to the first option in the following, but keep in mind that depending on your setup, the entropy reference scale (and thus the Rayleigh number definition) might change.

(26)$$\tilde{\rho }\tilde{T}({\displaystyle \frac{\partial s^{\prime }}{\partial t}}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla }}{s}^{\prime })={\displaystyle \frac{1}{Pr}}\overrightarrow{\mathrm{\nabla }}\cdot (\kappa (r)\tilde{\rho }\tilde{T}\overrightarrow{\mathrm{\nabla }}{s}^{\prime })+{\displaystyle \frac{PrDi}{Ra}}{\mathrm{\Phi }}_{\nu }+{\displaystyle \frac{PrDi}{P{m}^{2}ERa}}\lambda (r){(\overrightarrow{\mathrm{\nabla }}\times \overrightarrow{B})}^2,$$

A comparison with (20) reveals meaning of the different non-dimensional numbers that scale viscous and Ohmic heating. The fraction $Pr/Ra$ simply expresses the ratio of entropy and flow in the Navier-Stokes equation, while the additional factor $1/EPm$ reflects the scale difference of magnetic field and flow. Then remaining dissipation number $Di$ then expresses the relative importance of viscous and Ohmic heating compared to buoyancy and Lorentz force in the Navier-Stokes equation. For small $Di$ both heating terms can be neglected compared to entropy changes due to advection, an limit that is used in the Boussinesq approximation.

Dimensionless control parameters

The equations (20)-(26) are governed by four dimensionless numbers: the Ekman number

(27)$$E=\frac{\nu }{\mathrm{\Omega }{d}^2},$$

the thermal Rayleigh number

(28)$$Ra=\frac{{\alpha }_o{g}_o{T}_o{d}^3\mathrm{\Delta }s}{c_p{\kappa }_o{\nu }_o},$$

the compositional Rayleigh number

(29)$$R{a}_{\xi }=\frac{{\delta }_o{g}_o{d}^3\mathrm{\Delta }\xi }{{\kappa }_{\xi }{\nu }_o},$$

the Prandtl number

(30)$$Pr=\frac{{\nu }_o}{{\kappa }_o},$$

the Schmidt number

(31)$$Sc=\frac{{\nu }_o}{{\kappa }_{\xi }},$$

and the magnetic Prandtl number

(32)$$Pm=\frac{{\nu }_o}{{\lambda }_i}.$$

In addition to these four numbers, the reference state is controlled by the geometry of the spherical shell given by its radius ratio

(33)$$\eta =\frac{r_i}{r_o},$$

and the background density and temperature profiles, either controlled by $Di$ or by $N_{\rho }$ and $m$.

In the Boussinesq approximation all physical properties are assumed to be homogeneous and we can drop the sub-indices $o$ and $i$ except for gravity. Moreover, the Rayleigh number can be expressed in terms of the temperature jump across the shell:

(34)$$Ra=\frac{\alpha g_o{d}^3\mathrm{\Delta }T}{\kappa \nu }.$$

See also

In MagIC, those control parameters can be adjusted in the &phys_param section of the input namelist.

Variants of the non-dimensional equations and control parameters result from different choices for the fundamental scales. For the length scale often $r_o$ is chosen instead of $d$. Other natural scales for time are the magnetic or the thermal diffusion time, or the rotation period. There are also different options for scaling the magnetic field strength. The prefactor of two, which is retained in the Coriolis term in (20), is often incorporated into the definition of the Ekman number.

See also

Those references timescales and length scales can be adjusted by several input parameters in the &control section of the input namelist.

Usual diagnostic quantities

Characteristic properties of the solution are usually expressed in terms of non-dimensional diagnostic parameters. In the context of the geodynamo for instance, the two most important ones are the magnetic Reynolds number $Rm$ and the Elsasser number $\mathrm{\Lambda }$. Usually the rms-values of the velocity $u_{rms}$ and of the magnetic field $B_{rms}$ inside the spherical shell are taken as characteristic values. The magnetic Reynolds number

$$Rm=\frac{u_{rms}d}{{\lambda }_i}$$

can be considered as a measure for the flow velocity and describes the ratio of advection of the magnetic field to magnetic diffusion. Other characteristic non-dimensional numbers related to the flow velocity are the (hydrodynamic) Reynolds number

$$Re=\frac{u_{rms}d}{{\nu }_o},$$

which measures the ratio of inertial forces to viscous forces, and the Rossby number

$$Ro=\frac{u_{rms}}{\mathrm{\Omega }d},$$

a measure for the ratio of inertial to Coriolis forces.

$$\mathrm{\Lambda }=\frac{B_{rms}^2}{{\mu }_0{\lambda }_i{\rho }_o\mathrm{\Omega }}$$

measures the ratio of Lorentz to Coriolis forces and is equivalent to the square of the non-dimensional magnetic field strength in the scaling chosen here.

See also

The time-evolution of these diagnostic quantities are stored in the par.TAG file produced during the run of MagIC.

Boundary conditions and treatment of inner core

Mechanical conditions

In its simplest form, when modelling the geodynamo, the fluid shell is treated as a container with rigid, impenetrable, and co-rotating walls. This implies that within the rotating frame of reference all velocity components vanish at $r_o$ and $r_i$. In case of modelling the free surface of a gas giant planets or a star, it is preferable to rather replace the condition of zero horizontal velocity by one of vanishing viscous shear stresses (the so-called free-slip condition).

Furthermore, even in case of modelling the liquid iron core of a terrestrial planet, there is no a priori reason why the inner core should necessarily co-rotate with the mantle. Some models for instance allow for differential rotation of the inner core and mantle with respect to the reference frame. The change of rotation rate is determined from the net torque. Viscous, electromagnetic, and torques due to gravitational coupling between density heterogeneities in the mantle and in the inner core contribute.

See also

The mechanical boundary conditions can be adjusted with the parameters ktopv and kbotv in the &phys_param section of the input namelist.

Magnetic boundary conditions and inner core conductivity

When assuming that the fluid shell is surrounded by electrically insulating regions (inner core and external part), the magnetic field inside the fluid shell matches continuously to a potential field in both the exterior and the interior regions. Alternative magnetic boundary conditions (like cancellation of the horizontal component of the field ) are also possible.

Depending on the physical problem you want to model, treating the inner core as an insulator is not realistic either, and it might instead be more appropriate to assume that it has the same electrical conductivity as the fluid shell. In this case, an equation equivalent to (24) must be solved for the inner core, where the velocity field simply describes the solid body rotation of the inner core with respect to the reference frame. At the inner core boundary a continuity condition for the magnetic field and the horizontal component of the electrical field apply.

See also

The magnetic boundary conditions can be adjusted with the parameters ktopb and kbotb in the &phys_param section of the input namelist.

Thermal boundary conditions and distribution of buoyancy sources

In many dynamo models, convection is simply driven by an imposed fixed super-adiabatic entropy contrast between the inner and outer boundaries. This approximation is however not necessarily the best choice, since for instance, in the present Earth, convection is thought to be driven by a combination of thermal and compositional buoyancy. Sources of heat are the release of latent heat of inner core solidification and the secular cooling of the outer and inner core, which can effectively be treated like a heat source. The heat loss from the core is controlled by the convecting mantle, which effectively imposes a condition of fixed heat flux at the core-mantle boundary on the dynamo. The heat flux is in that case spatially and temporally variable.

See also

The thermal boundary conditions can be adjusted with the parameters ktops and kbots in the &phys_param section of the input namelist.

Chemical composition boundary conditions

They are treated in a very similar manner as the thermal boundary conditions

See also

The boundary conditions for composition can be adjusted with the parameters ktopxi and kbotxi in the &phys_param section of the input namelist.