Evolution of the inner core of the earth: consequences for geodynamo

TL;DR

This paper models Earth's core evolution over billions of years to understand changes in geodynamo parameters and their impact on the Earth's magnetic field generation.

Contribution

It introduces a comprehensive evolutionary model of Earth's core that links thermal and compositional convection to geomagnetic field regimes over geological time.

Findings

Evolution of energy sources for convection over Earth's history.

Changes in spatial scales of convective zones.

Variation in dimensionless parameters affecting geodynamo regimes.

Abstract

Using models of the Earth's core evolution and the length of the day observations the change of the dimensionless geodynamo parameters is considered. The evolutionary model includes cooling of the liquid adiabatic core, growing solid core, and the region in the outer part of the core with a subadiabatic temperature gradient. The model covers time period 4.5Ga in the past till 1.5Gy to the future, and produce evolution of the energy sources of the thermal and compositional convection, spatial scales of the convective zone. These quantities are used for Ekman, Rayleigh and Rossby numbers estimates. So far these numbers determine regime of the geomagnetic field generation, we discuss evolution of the geomagnetic field over Earth's evolution.

Tables

Table 1. Table 1:

Parameter	Denotation	Value
Gravitational constant	$G$	${\mathrm{6.6873\hspace{0.17em}10}}^{-11}{m}^3/\left(\mathrm{kg}{s}^2\right)$
Thermal diffusivity	$k$	${\mathrm{7\hspace{0.17em}10}}^{-6}{m}^2/s$
Kinematic viscosity	$\nu$	${10}^{-6}{m}^2/s$
Light element diffusivity	$\lambda$	${10}^{-9}{m}^2/s$
Coefficient of chemical expansion	$\beta$	$1$
Gruüneisen parameter	$\gamma$	$1.5$
Core radius	$r_b$	$3480\mathrm{km}$
Entropy of crystallization	$\delta S$	$118J/\left(\mathrm{kg}K\right)$
Density at zero pressure	${\rho }_{\circ }$	$7500\mathrm{kg}/m^3$
Density jump at ICB	$\delta \rho$	$500\mathrm{kg}/m^3$
Temperature of solidi¢cation at the cente	$T_s^{\circ }$	$5270K$
Initial temperature in the center of the core	$T_{\circ }$	$6000K$
Incompressibility of the core	$K_{\circ }$	${\mathrm{4.76\hspace{0.17em}10}}^{11}\mathrm{Pa}$
Specific heat	$C_p$	$860J/\left(\mathrm{kg}K\right)$

Equations

$$\nabla P=-\rho g,g(r)=\frac{4\pi G}{r^2}\underset{0}{\overset{r}{\int }}\rho (u)u^{2}du,$$

$$P=K_{\circ }\frac{\rho }{{\rho }_{\circ }}\ln \frac{\rho }{{\rho }_{\circ }},$$

$$\rho (r)=\{\begin{array}{ll}\rho (r),& \text{if }r\le c\\ \rho (r)-\delta \rho ,& \text{if }r>c.\end{array}$$

$$T_{ad}(r)=T_c(c)exp\left(-\underset{c}{\overset{r}{\int }}\frac{\alpha (u)g(u)}{C_p}du\right),$$

$$\begin{array}[]{l}\displaystyle\alpha(r)={\gamma C_{p}\rho_{\circ}\over K_{\circ}\Big{(}1+\displaystyle\ln{\rho\over\rho_{\circ}}\Big{)}},\end{array}$$

$$4\pi r_1^2{q}_1=-4\pi \underset{0}{\overset{r_1}{\int }}\rho C_p\frac{\partial T_{ad}}{\partial t}{r}^{2}dr=-\frac{\partial T_{\circ }S}{\partial t},S=4\pi \underset{0}{\overset{r_1}{\int }}\rho C_{p}exp\left(-\underset{0}{\overset{r}{\int }}\frac{\alpha g}{C_p}\right)r^{2}dr,$$

$$\begin{array}[]{l}\displaystyle T_{s}(r)=T_{s}^{\circ}\left(\displaystyle{\rho(r)\over\displaystyle\rho(c)}\right)^{2(\gamma-{1\over 3})},\end{array}$$

$$\begin{array}[]{l}\displaystyle r_{b}^{2}q_{b}-c^{2}q_{c}=\dot{c}\Big{(}c^{2}(P_{L}+P_{G})+P_{C}\Big{)},\end{array}$$

$$\begin{array}[]{l}\displaystyle P_{L}(c)=\rho(c)\delta S\,T_{s}(c),\end{array}$$

$$\begin{array}[]{l}\displaystyle E_{G}={2\pi\over 5}GM_{\circ}\delta\rho{c^{3}\over c_{b}}\left(1-\left({c\over r_{b}}\right)^{2}\right),\end{array}$$

$$\begin{array}[]{l}\displaystyle\dot{E}_{G}=P_{G}\dot{c},\qquad P_{G}={12\pi\over 5}{GM_{\circ}\delta\rho\over r_{b}}c\left(1-{2c^{2}\over r_{b}^{2}}\right).\end{array}$$

$$\begin{array}[]{l}\displaystyle P_{C}=-\int\limits_{c}^{r_{1}}\displaystyle\rho C_{p}{\partial T_{ad}\over\partial c}r^{2}\,dr,\qquad dc\equiv dr,\end{array}$$

$$\begin{array}[]{l}\displaystyle{\partial T\over\partial t}=k\Delta T,\end{array}$$

$$T_{\mathrm{d}}(t)=\{\begin{array}{ll}21.435+0.974(t-4.5),& 0\le t\le 3.86,\\ 24+4.98(t-4.5),& 3.86\le t\le 6,\end{array}$$

$${\mathrm{Ra}}_{\mathrm{C}}=\frac{\beta g_0{\dot{\chi }}_{\circ }{D}_{\mathrm{C}}^5}{\lambda {\nu }^2},$$

$${\dot{\chi }}_{\circ }=\frac{3{\chi }_{\circ }}{1-{\chi }_{\circ }}\frac{c^2\dot{c}}{r_b^3-c^3}.$$

$${\chi }_{\circ }(t)=Cexp\left(3{\left(\frac{c(t)}{r_b}\right)}^3\right),$$

Full text

Introduction

Dynamo theory is the most promising candidate to explain existence of the geomagnetic field observed at the surface of the Earth [1]. If the direct three-dimensional modeling of the magnetic field generation over the times comparable to the age of the Earth ($4.5$Gy) is still beyond the modern computer capacities, because of the wide range of the temporal-space scales of the MHD turbulence in the core, then study of influence of the evolution of the core on the geodynamo parameters is already quite realistic task [2].

The straightforward application of the evolutionary scenarios to the geodynamo models is complicated by the different representation of parameters. Thus, the evolutionary models produce evolution of the inner core and subadiabatic region boundaries, heat fluxes. In other words the models give the behavior of the physical quantities. On the other hand, the geodynamo models have deal with a set of the dimensionless parameters, which are very useful for analysis of the physical state in the core. As a result translation from one language to another is needed.

Below we repeat the well-known and accepted by the geodynamo community modeling of the thermal evolution model of the core, originated to the pioneer papers [3, 4, 5, 6, 7] and references therein, which include cooling of the liquid core, evolution of the geometry of the liquid core, concerned with the growth of the inner core and subadiabatic region, change of the energy sources. So far, the Earth’s angular velocity, which contributes to the dimensionless numbers, evolves on the geological times due to the tidal forces, information on variations of the length of the day from observations was used as well [8]. As a result we get evolution of the Ekman, Rayleigh, Rossby numbers, the critical Rayleigh number for the thermal and compositional convection starting from the origin of the liquid core about 4.5Ga, and extrapolate it to a close, on a geological scales, future 1.5Gy. We also discuss the known paradox [9], concerned with the fact that compositional convection, started with the origin of the inner core, should change the geomagnetic field generation in the core.

1 Thermodynamic model of the core

Following [4, 5, 7] we consider scenario of the Earth’s evolution, where soon after the end of the accretion process, the Earth’s core of radius $r_{b}$ was fully convective. Then, it cooled due to the thermal flux density $q_{b}$ at the core-mantle boundary (CMB) $r=r_{b}$, and as a result, depending on the amplitude of $q_{b}$, two regions could appear: the solid inner core ($0\leq r\leq c$, region I) and subadiabatic layer in the outer part of the core ($r_{1}\leq r\leq r_{b}$, region III). The rest convective part of the core $c\leq r\leq r_{1}$ here and after is denoted as region II.

Radial distributions of density $\rho(r)$, pressure $P(r)$ and gravity $g(r)$ satisfy to the hydrostatic balance equations:

[TABLE]

with $G$ the gravitational constant. To close system of equations for $(P,\,\rho,\,g)$ the logarithmic equation of state is used:

[TABLE]

where $K_{\circ}$, $\rho_{\circ}$ are incompressibility and density at zero pressure, respectively. The optional in the model jump of the density, observed at the surface of the inner core, and which effect on the evolution of the core is quite small, is introduced as follows:

[TABLE]

Eqs(1–3) with given $c$ can be solved numerically. Then, with known $(P,\,\rho,\,g)$, adiabatic temperature profile can be derived:

[TABLE]

where $T_{c}(c)$ is the temperature at $r=c$, thermal expansion coefficient

[TABLE]

with $C_{p}$ specific heat, and $\gamma$ for Grüneisen parameter.

If the inner core is still absent, $c=0$, then $T_{c}(c)=T_{\circ}$, where the temperature in the center of the Earth $T_{\circ}$ can be found from the heat balance equation:

[TABLE]

with $q_{1}$ for heat flux density at $r_{1}$. The growth of the inner core starts, when temperature of the liquid core is equal to the temperature of solidification:

[TABLE]

where $T_{s}^{\circ}$ is the temperature of solidification in the center of the Earth. Solidification process starts in the core’s center, i.e. $T_{c}=T_{\circ}=T_{s}^{\circ}$, $r=c=0$. Then, for $c>0$, $T_{s}$ defines adiabatic temperature at the boundary $c$ in (4): $T_{c}(c)=T_{s}(c)$.

Position of the inner core boundary $c$ can be derived from the heat flux equation:

[TABLE]

where on the left side is the total heat flux in the region II, and on the right one are the cooling sources, and the dot over $c$ stands for the time derivative.

The latent heat source is defined as

[TABLE]

with $\delta S$ entropy of crystallization.

Estimate of the release of the gravitational energy due to the growth of the inner core has the form [11]:

[TABLE]

with mass of the core $\displaystyle M_{\circ}={4\over 3}\pi\int\limits_{0}^{r_{b}}\rho(r)r^{2}\,dr$ constant in the model. Then it leads to

[TABLE]

The main term, concerned with adiabatic cooling, has the form:

[TABLE]

with $q_{c}$ heat flux density through the boundary $c$. Eq(8) was resolved with respect to $\dot{c}$ and then integrated in time. This defines evolution of the inner core boundary $c$ in time.

From condition of continuity of the temperature at the boundary $c$, follows that $T_{s}(c)$ is the boundary condition for the thermal-diffusion equation in the region $0<r<c$, I, with a moving boundary $c(t)$ [10]:

[TABLE]

where $k$ is the thermal diffusivity. The second boundary condition in the center $r=0$ is $T^{\prime}=0$, where ′ is a derivative on $r$. The joined system $(\ref{c1}-\ref{c16})$ defines evolution of the fields in the regions I and II.

If the adiabatic heat flux density $\displaystyle q_{ad}(r)=-\kappa T^{\prime}_{ad}(r)$, with thermal conductivity $\kappa=k\rho C_{P}$, becomes larger than the heat flux density at the outer boundary $r_{b}$: $\displaystyle q_{ad}(r)<\left({r_{b}\over r}\right)^{2}q_{b}$, the subadiabatic stably stratified thermal region III develops at the outer part o the core, where the heat flux density is smaller. The temperature profile in the region III can be derived from Eq(9) with the moving boundary $r_{1}(t)$, and two boundary conditions: $T(r_{1})=T_{ad}(r_{1})$ at the inner boundary, and given heat flux density $q_{b}(t)$ at the outer boundary $r_{b}$. In the general case Eqs(1–13) in regions I–III are solved numerically, using iterative methods with under-relaxation method to provide numerical stability.

2 Evolution of the boundaries $c$ and $r_{1}$

We consider two regimes similar to that ones in [5], proposed by the authors, to get the realistic size of the inner core at the present time. These regimes differ in decrease of the heat flux density $q_{\beta}$ in the linear dependency $q_{b}(t)=q_{b}(0)-q_{\beta}t$ for the heat flux intensity at CMB. Let at $t=0$ $q_{b}(0)=0.075$ W/m2, then in the case $A$ $q_{\beta}=10^{-19}$ W/(s m2), and $q_{\beta}=3\,10^{-19}$ W/(s m2) in the case $B$. These values correspond to decrease of the heat flux density of 4.2% and 12.6% during the first 1Gy, correspondingly. Eqs(1–13), with parameters summarized in Table 1, were solved numerically for the time interval 6Gy. This time interval corresponds to 4.5Ga to the past, starting from the origin of the liquid core, till 1.5Gy to the future.

Evidently that cases $A$ and $B$ differ in evolution of the boundaries $c$ and $r_{1}$, see Fig.1.

In the case $A$ region III is absent during the whole time interval, and the inner core appears earlier ($t=3$Gy) than for $B$ ($t=3.6$Gy). In the case $B$ firstly appears region III ($t=3.2$ Gy), and after that the inner core (I). The velocity $\dot{r}_{1}$ in case $B$ is approximately constant in time and greater than $\dot{c}$. The age of the inner core in cases $A$ and $B$ is $1.5$Ga and $0.9$Ga, correspondingly. In the case $A$ and $B$ at the present time the inner core radius $c$ is equal to $1230$km and $810$km, correspondingly. Note, that seismological estimate of $c$ is $1220$km. In the case $B$, where the heat flux at CMB decreases faster, the inner core radius is smaller.

3 Evolution of the length of the day

So far convection and magnetic field generation in the liquid core are closely related to the angular rotation velocity $\Omega$ of the Earth, evolution of $\Omega$ in time should be also considered. Due to the tidal forces Earth’s length of the day (LOD) $\displaystyle T_{\rm d}={2\pi\over\Omega}$ increases in time. The estimate of this increase over the last three centuries is $(2.0\pm 0.2)$ms cy*-1*. This process can be described analytically, using astronomical predictions on the evolution of the Earth’s orbit, and the results are in agreement with observations. However astronomical methods can not be extrapolated to the geological times, because they do not include evolution of the Earth itself. This is the reason why the new palaeontological and palaeosedimental observations should be used, see review [8]. The relevant information comes from fossils (bivalves, brachiopods and corals) of the Phanerozoic (time interval from 542Ma to the present time), from stromatolites mainly of the Proterozoic (from 2500Ma to the Phanerozoic), and from palaeodeposits of the Proterozoic. The main idea of these approaches is that one can resolve as the daily marks of the growth of these organisms, as well as the annual, and even seasonal ones. Then finding the number of days per year is the trivial procedure. To the moment information on LOD is available for the last 2.5Gy. To my knowledge such a long time series of LOD were never used in the geodynamo theory. Here we extrapolate data from [8] to the last 6Gy in the following form:

[TABLE]

where LOD $T_{\rm d}$ is measured in hours, and $t$ in Gy. The present time corresponds to $t=4.5$. There is the break at $t=3.86$, when $T_{\rm d}$ starts to decrease faster in factor 5. However the certain estimates of accuracy of the data hardly can be done, authors of [8] demonstrate that the break is statistically significant. One of the possible explanation of such a break, considered by authors, is the change of the plate tectonics of the Earth.

4 Dimensionless numbers

The Ekman number $\displaystyle\mathop{\rm E}\nolimits={\nu\over 2\Omega D^{2}}$, with $D$ typical scale of convection, characterizes ratio of the viscous and Coriolis forces. Definition of the scale $D$ depends on the model of convection. For the thermal convection, denoted by index $\rm T$, the stably stratified layer III is excluded, and $D=D_{\rm T}=r_{1}-c$. For the compositional convection (index $\rm C$), the full convective liquid core is considered $D=D_{\rm C}=r_{b}-c$. If region III is absent, then $D_{\rm T}=D_{\rm C}$.

Evolution of $\mathop{\rm E}\nolimits_{\rm T}$ is presented in Fig.2a. There is the slight increase of $\mathop{\rm E}\nolimits_{\rm T}$ before the inner core origin, caused by change of $T_{\rm d}$ (14). After regions I and III have developed, evolution of $\mathop{\rm E}\nolimits_{\rm T}$ is controlled by the scale of convection $D_{\rm T}(t)$. For compositional convection for the case $A$ $\mathop{\rm E}\nolimits_{\rm T}=\mathop{\rm E}\nolimits_{\rm C}$ because $D_{\rm T}=D_{\rm C}$. In the regime $B$ the growth of $\rm E_{\rm C}$ is less, because $D_{\rm C}$ does not include $r_{1}(t)$. As a result, curves do not intersect.

The measure of the heat sources intensity is the Rayleigh number: $\displaystyle\mathop{\rm Ra}\nolimits_{\rm T}=\displaystyle{\alpha gq_{r_{1}}D_{\rm T}^{4}\over\kappa k\nu}$, with \displaystyle q_{r_{1}}(t)=\Big{(}{r_{b}\over r_{1}}\Big{)}^{2}q_{b}(t). It is instructive to consider $\displaystyle\mathop{\rm Ra}\nolimits_{\rm T}$ in units of its critical value $\displaystyle\rm Ra^{cr}$: $\rm\widehat{\mathop{\rm Ra}\nolimits}=\mathop{\rm Ra}\nolimits/\mathop{\rm Ra}\nolimits^{cr}$. For $\mathop{\rm E}\nolimits_{\rm T}\ll 1$ $\displaystyle\rm Ra^{cr}\sim{\mathop{\rm Pr}\nolimits^{4/3}\over(\mathop{\rm Pr}\nolimits+1)^{1/3}}E_{T}^{-4/3}$, and in the case of the thermal convection, where the Prandtl number $\displaystyle\mathop{\rm Pr}\nolimits={\nu\over\kappa}\approx 0.1$ is small, $\displaystyle\rm Ra^{cr}\sim{\mathop{\rm Pr}\nolimits^{4/3}}E_{T}^{-4/3}$. Before the inner core appears, $\rm\widehat{Ra_{T}}$, due to evolution of $q_{b}(t)$ and $\Omega(t)$, evolves slightly different in cases $A$ and $B$, Fig. 2b. However, after the start of the inner core solidification, when $\dot{c}$ is large, there is sharp decrease of $\rm\widehat{Ra_{T}}$ in the both cases. Duration of such decrease is limited by the time moment, when the length of the day changes in the model ($t=3.86$Gy). After that $\rm\widehat{Ra_{T}}$ continues to grow in the case $A$, and decrease in $B$. Note, that if the origin of the inner core would coincide with the break in dependence (14), then the jump in $\rm\widehat{Ra_{T}}$ will disappear. This fact can be an implicit justification that the break in (14) corresponds to the inner core emergency.

The Rossby number $\mathop{\rm Ro}\nolimits$ is the ratio of $T_{\rm d}$ to the typical convective time. The larger is $\mathop{\rm Ro}\nolimits$ the more negligible is the Coriolis force in the system. The value of $\mathop{\rm Ro}\nolimits$ can be derived from the following scaling-law $\displaystyle\mathop{\rm Ro}\nolimits_{\rm T}=\left({\mathop{\rm Ra}\nolimits_{\rm T}\,\mathop{\rm E}\nolimits_{\rm T}^{2}}\right)^{0.41}$ [12]. The behavior of $\mathop{\rm Ro}\nolimits_{\rm T}$ is similar in the both cases $A$ and $B$, see Fig. 2c: the slow change till the break in $T_{\rm d}$, and after that increase for 50% during 2Gy. For the present time one has to divide 50% on factor 4. Summarizing, we conclude that increase of $\mathop{\rm Ro}\nolimits_{\rm T}$, based on the thermal convection model, after the origin of the inner core is quite small.

For the compositional convection the following estimate of the Rayleigh number is used [13]:

[TABLE]

with $\lambda$ for the light element diffusivity, $\displaystyle\beta=-{1\over\rho}{\partial\rho\over\partial\chi}$ coefficient of the chemical expansion, and time derivative of the light element concentration $\chi_{\circ}$ is related to $\dot{c}$ as [13, 2]:

[TABLE]

In the limit of $\chi_{\circ}\ll 1$, and $c\ll r_{b}$, Eq(16) has exact solution:

[TABLE]

where $C$ is the constant, defined by the initial condition that in the present time $\chi_{\circ}\sim 0.1$ [2]. Using calculated dependence $c(t)$ one estimates that $\chi_{\circ}$ in cases $A$ and $B$ is in the range of $[0.89\div 0.11]$, $[0.97\div 0.1]$, correspondingly.

Estimate of $\rm\mathop{\rm Ra}\nolimits_{\rm C}$ gives value three orders of magnitude larger than for the thermal convection [13], however this difference for the normalized values $\widehat{\mathop{\rm Ra}\nolimits}$ is only factor 20, see Fig. 2b and 3b. It is because for the compositional convection $\displaystyle\mathop{\rm Pr}\nolimits={\nu\over\lambda}\sim 10^{3}$, and $\mathop{\rm Ra}\nolimits^{\rm cr}\sim\mathop{\rm Pr}\nolimits\mathop{\rm E}\nolimits^{-4/3}_{\rm C}$. The time dependence of $\rm\widehat{\mathop{\rm Ra}\nolimits_{\rm C}}$ is similar to that ones for the thermal convection regime in Fig. 2b. Again, if the origin of the inner core coincides with the break in $T_{\rm d}$, the jump in curve $\widehat{\mathop{\rm Ra}\nolimits_{\rm C}}$ disappears.

By analogy with the thermal convection regime, estimate of the Rossby number for the compositional convection $\displaystyle\rm\mathop{\rm Ro}\nolimits_{C}=\left({\mathop{\rm Ra}\nolimits_{\rm C}\,\mathop{\rm E}\nolimits_{C}^{2}}\right)^{0.41}$ is plotted in Fig. 3b. Firstly note, that this estimate leads to the quite large $\displaystyle\rm\mathop{\rm Ro}\nolimits_{C}\gg 1$, which hardly corresponds to the geostrophic state in the liquid core. Possibly, the additional normalization for $\rm\mathop{\rm Ra}\nolimits_{C}$ is needed. Let focus our attention at the time behavior of the curves, which differs also in the cases $A$ and $B$: for $A$ contribution of \dot{\mathop{\rm E}}\nolimits>0 is more essential than of \dot{\mathop{\rm Ra}}\nolimits_{C}<0, and for $B$ situation is quite opposite. It follows from 3D simulations [12] that Rossby number controls the frequency of the geomagnetic reversals: the smaller is $\rm Ro$ the faster geomagnetic dipole changes its polarity. From this point of view the case $A$ looks more attractive, because of the general increase of the number of reversals after the Cretaceous superchrone.

5 Conclusions

The evolutionary models predict existence of the young inner core $\sim 1$Ga. Switch on of the compositional convection, concerned with the inner core solidification, leads to increase of convection intensity in the core. In units of its critical value the Rayleigh number increases in factor 20 that corresponds to increase of the intensity of the geomagnetic field times $(20)^{1\over 3}\approx 2.7$ [12]. Moreover, influence of the inner core on the magnetic field generation can be even smaller. The reason is that compositional convection generates magnetic field located deeper in the liquid core, and as a result, intensity of the magnetic field at the surface of the planet is weaker than in the case of the thermal convection. In some sense our result confirms recent results [14] that switch on of the compositional convection does not change generation of the magnetic field essentially, as it was supposed earlier [9]. However the absolute values of the Rayleigh numbers are still under consideration, and their comparison for the compositional and thermal convection models is a tricky procedure, the time evolution of these parameters is not so questionable: the considered above scenarios predict increase of the geomagnetic reversals frequency with the inner core growth in agreement with observations.

It is worth noting that our analysis includes two independent physical data: concerned with the thermal evolution of the core, and evolution of the length of the day, based on observations. Of course, these data are interrelated, and the more complicated physical model is required for description. The hint is that identification of the break in the curve of length of the day with the origin of the inner core simplifies the general evolution scenario for the considered model. In this case the young inner core is more likely.