Implication of the lopsided growth for the viscosity of Earth’s inner core
Two main seismic features characterize the Earth’s inner core: a North-South polar anisotropy and an East-West asymmetry of P-wave velocity and attenuation. Anisotropy is expected if shear deformation is induced by convective motions. Translation has recently been put forward as an important mode of convection of the inner core. Combined with a simple diffusive grain growth model, this mechanism is able to explain the observed seismic asymmetry, but not the bulk anisotropy. The source of anisotropy has therefore to be sought in the shear motions caused by higher modes of convection. Using a hybrid finite-difference spherical harmonics Navier-Stokes solver, we investigate the interplay between translation and convection in a 3D spherical model with permeable boundary conditions at the inner core boundary. Three parameters act independently: viscosity, internal heating and convection velocity in the outer core. Our numerical simulations show the dominance of pure translation for viscosities of the inner core higher than Pas. Translation is almost completely hampered by convective motions for viscosities lower than Pas. Between these values, translation and convection develop, but convective downwellings are restricted to the coldest hemisphere where crystallization occurs. On the opposite side, shear is almost absent, thereby allowing grain growth. We propose that the coexistence of translation and convection observed in our numerical models leads to a seismic asymmetry but localizes deformation only in one hemisphere.
keywords:Earth’s inner core, numerical modeling, thermal convection, translation
The image of the inner core growing slowly at the center of the Earth by gradual cooling and solidification of the surrounding liquid outer core is being replaced by the more vigorous image of a “deep foundry” (Buffett2011), where melting and crystallization rates exceed by many times the net growth rate (Monnereau2010; Alboussiere2010; Gubbins2011)
During seventy years, the analysis of compressive waves (P waves) and free oscillations excited after large earthquakes have been depicting a more and more complex structure of the inner core. It appears anisotropic, with a fast axis parallel to Earth’s spin axis. It is also asymmetric: within the outermost 100 km, velocity and attenuation of P waves increase from the hemisphere facing America (West) to the one facing Asia (East) (Tanaka1997) ; anisotropy also seems stronger in the Western hemisphere than in the Eastern one (Deuss2010).
In the 80’s and 90’s, anisotropy was thought to be the prominent feature. This is commonly attributed to the preferred orientation of iron crystals, possibly acquired during solidification but most probably resulting from creep flow. Thermal convection, developing a flow characterized by a spherical harmonic degree (a clementine shape), orientated along the spin axis of the Earth, was one of the first candidates to account for anisotropy (Jeanloz1988). In this model and the following ones (Weber1992; Buffett2009), the nature of the inner core boundary (i.e. a phase change) was not considered, a classical impermeable boundary condition being preferred. With permeable boundary conditions — a phase change does not prevent material transfer —, the expression of convection at the first harmonic degree is a constant velocity field across the inner core, that obviously does not produce any deformation, but implies melting on one side and crystallization on the opposite side. This peculiar situation, called translation of the inner core, has recently been put forward to explain the hemispherical asymmetry of velocity and attenuation (Monnereau2010).
If the inner core grows in a superadiabatic regime, which is the condition for the onset of convection, an unstable thermal stratification develops, so that any infinitesimal thermal heterogeneity of harmonic degree (i.e. one side colder than the other one) will be amplified. Such a heterogeneity will induce a displacement of the inner core to maintain its center of mass — shifted toward the denser/colder hemisphere — at the center of the Earth. The inner core acquires a positive topography on the hotter and lighter side and a negative one on the opposite side. This topography is thermodynamically unstable: the side emerging above the phase change melts, bringing hotter material up to the surface, while the sinking side allows crystallization. The phase change acts to remove the topography, which is continuously rebuilt by isostatic equilibrium. This feedback results in a permanent drift from the crystallizing side to the melting side, the drift velocity being controlled by the ability of the outer core to restore the adiabatic condition at the surface of the inner core.
Translation is consistent with multiple scattering models of wave propagation. If they do not experience deformation, iron crystals grow as they transit from one hemisphere to the other. Larger crystals constituting a faster and more attenuating medium, a translation velocity of some cm/yr (about ten times the growth rate) is enough to account for the superficial asymmetry observed for P-wave velocity and attenuation, with grains of a few hundred meters on the crystallizing side (West) growing up to a few kilometers before melting on the East side, and a drift direction located in the equatorial plane.
Translation was also proposed to be responsible for the formation of a dense layer at the bottom of the outer core, since the high rate of melting and crystallization would release a liquid depleted in light elements at the surface of the inner core (Alboussiere2010). This would explain the anomalously low gradient of P wave velocity in the lowermost 200 km of the outer core (Poupinet1983).
Clearly, translation cannot account for anisotropy. However, convective modes developing at higher harmonic degrees () will necessarily induce shear deformation that could be a source of anisotropy. The development of these modes depends on the Rayleigh number, that controls the vigor of the convection, and thus mainly on the viscosity of the solid portion of the core. At high Rayleigh number (low viscosity), these modes can be dominant and dissipate the degree of the thermal heterogeneities: the source of the translation. Thus a viscosity threshold may be expected below which translation would not take place. This may constrain the viscosity of iron at the conditions of the inner core, based on seismological observation. In this paper, we present dynamics model of inner core taking into account the phase change boundary and study the interaction of translation and convection.
Inner core dynamics obey similar governing equations to those used for mantle convection. The specificity lies in the boundary conditions required to treat the phase change at its surface.
Inertial forces can be neglected because of the high viscosity of solid iron at the inner core temperature and pressure conditions, which is at least Pas (Yoshida1996). Conservation of momentum just expresses the balance between buoyancy forces and viscous dissipation. It is time independent.
is the deviatoric stress tensor, the pressure, the density and g the gravity. The permeable surface condition is introduced. It describes the balance between the radial stress and buoyancy forces induced by topography:
where is the deviatoric radial stress, the inner core radius, the density of solid iron, the density of liquid iron, the gravity at the surface of the inner core and the topography. In addition, the surface is considered as tangential stress free:
At the first spherical harmonic degree , the momentum equation written for the full sphere with the above conditions comes down to a simple isostasy equilibrium where thermal heterogeneities are balanced by a topography ( in Eq. 1) (Monnereau2010). This singularity requires a particular treatment. The three components of the first harmonic degree of the topography are directly related to the position of the center of mass anomaly caused by the presence of the inner core in the outer core:
with the thermal expansion coefficient and the radial profile of the temperature heterogeneity at degree and order (see Appendix A for more details).
Surface heat exchange with the outer core
The topography is thermodynamically unstable and eroded at a rate depending on the vigor of convection within the outer core. In the mantle, the position of mineral phase transitions is mainly related to the ambient temperature, latent heat exchanges having almost no effect. For instance, a transition with a positive Clapeyron slope like olivine to spinel occurs deeper in ascending (hot) currents than in dipping slabs. The reverse situation happens for the ICB. The turbulent flow in the liquid maintains the temperature above the surface close to the adiabat so that the topography only depends on the temperature variations induced by the latent heat effects and not on temperature anomalies within the inner core. Topography is thus positive where material exits the inner core because of the cooling induced by the latent heat consumption. The local thermodynamical equilibrium is achieved when the rate at which the inner core consumes or releases latent heat equals the rate at which the outer core brings or takes the energy to maintain the adiabatic temperature. This can be written as (Alboussiere2010):
is the radial velocity across the phase change, the melting temperature of iron, the adiabatic temperature, the amplitude of the outer core convective flow at the surface of the inner core, the specific heat and L the latent heat of the phase change. We studied the kinetics of this equilibrium in a simple 1D model (see Appendix B) in order to check that it can be considered as instantaneous. The temperature adjustment is reached in a few thousands years which is at least one order of magnitude smaller than the characteristic time of the convective fluctuations at the highest Rayleigh number we considered.
The inner core dynamics is a moving boundary problem that can be treated in all its complexity (Deguen2011a). For the sake of simplicity, we may neglect the variation of the radius with time, and focus on the present time dynamics. This assumption is plainly justified since the translation velocity required to account for seismic properties of the inner core, but also for the formation of a dense layer above ICB, should exceed the inner core growth rate (by one or two orders of magnitude) (Monnereau2010; Alboussiere2010). The energy equation is thus written in terms of the temperature relative to the adiabat anchored at the ICB (), in which the decrease of ICB temperature with time plays the role of an internal heating:
where v is the velocity vector, the conductivity, the internal heating rate and the Laplacian operator. The surface is assumed isothermal, . This comes to neglect the small temperature perturbation due to the latent heat (K) compared to the one involved in the dynamics (1K). Indeed, the latent heat effects play a much more important role in the momentum equation than in the energy equation.
, the internal heating, results from a competition between the heat lost by conduction along the adiabat and the secular cooling of the core:
|Thermal expansion coefficient||K|
|Heat capacity||800||J/K kg|
|Adiabatic temperature anchored at ICB||5500||K|
|thermal conductivity||34||W/K m|
|melting curve slope||K/Pa|
|solid iron density||11800||kg/m|
|liquid iron density||11200||kg/m|
|inner core radius||1220||km|
|gravity at the ICB||4||m/s|
|Outer core velocity||m/s|
The superadiabatic regime, required for the onset of translation or convection, is defined by . Uncertainties on both the conductivity and the cooling rate contribute to the uncertainty on . The former was evaluated around 60 W/mK (Stacey2001) ; an estimate that has been recently revisited and reduced by a factor of 2: 36 W/mK (Stacey2007). is proportional to the inner core growth rate :
being itself proportional to the heat flux at the core-mantle boundary (CMB) (Labrosse2003), whith walues for the latter estimated to be 6 TW and 14 TW (Turcotte2002). It results that can be as large as W/m.
The vigor of convection is controlled by the Rayleigh Number:
written here for pure internal heating in a spherical system. At first sight, the variation of radius appears of primary importance for the evolution of the dynamics of the inner core along its history, notably because of the successive development of convective modes at higher and higher degree as the inner core grows (Deguen2011a). However, as mentioned before, we focus on the present day dynamics. In this case, the least constrained parameters are the internal heating, and the viscosity, , both subject to uncertainties of several orders of magnitude.
A minimum viscosity of 10 Pas was deduced from attenuation of seismic waves in the inner core and a maximum of 10 Pas was inferred from conjectures on iron rheology close to its melting point for pressures at the center of the Earth (Yoshida1996). Since translation does not involve deformation, it is not sensitive to , but only to . Their effects have to be studied independently.
The third and last parameter is the amplitude of the convective velocity within the outer core, , that acts only on the translation velocity, the effect of the topography remaining very small on the convective circulation. Other parameters in the expression of the factor F (Eq. 5) remain much more constrained than . The average convection velocity in the outer core was assessed to be 10 m/s with 10 m/s maxima, from the secular variations of the geomagnetic field (Amit2006). It is possible that a dense layer forms at the base of the outer core (Alboussiere2010; Buffett2011), convection velocity could be lower at the surface of the inner core. We will also test a velocity ten times lower than the average.
The momentum equation is solved using a classical spherical harmonics expansion together with a radial finite difference solver for each harmonic. The conservation of energy is solved with a second order finite volume method. More details on the code may be found in Monnereau2002b. For all calculations we have employed 128 spherical harmonics and 128 points in the radial direction.
Table 1 summarizes the range of investigation for the parameters , and . Figure 1 presents the velocity fields of experiments with an intermediate value for the outer core velocity, m/s. A rapid overview shows that translation dominates at high viscosity and high internal heating (Fig. 1c, e-i), and is absent from experiments at low viscosity and low internal heating (Fig. 1a, b & d). Translation is characterized by a surface radial velocity, i.e. the radial velocity through the ICB, positive on one hemisphere and negative on the opposite, inducing melting and crystallization respectively.
At high viscosity, Pas, translation seems to be the sole active mode of heat transfer. Nevertheless, the velocity field does not coincide with a pure solid motion, except at very low internal heating rate (Fig. 1g). It also contains a weak degree component, attesting the development of convective modes at degree . This feature is not visible on the surface radial velocity maps, but is revealed by the non-collinearity of the velocity vectors depicted on the cross section of Fig. 1h. At lower viscosity, , translation does not develop for the lowest internal heating rate Fig. 1d. A degree dominates in this case, where a sheet-like downwelling draws a great circle. When translation still takes place (Fig. 1e & f), it is associated with one or several planar dipping currents. These short wavelengths structures only form within the crystallizing hemisphere. At low viscosity, Pas, these 2D currents give way to cold spots. Translation disappears at low and intermediate internal heating rates (Fig. 1a & b). It persists at the highest heating rate, but, contrary to the situations at Pas, high degree convective structures are present on both crystallizing and melting hemispheres even if an asymmetry in the number of structures is still present.
In summary, a low viscosity enhances the vigor of the convection and favors the development of short wavelength structures, a result expected since it corresponds to an increase in Rayleigh number.
Internal heating effect
Internal heating also enhances convection, as evidenced by the increase of the number of structures. However, when degree dominates the circulation, the major effect is an amplification of the translation, whose velocity seems only dependent on the internal heating rate. For instance, the comparison of situations depicted in Fig. 1e & h and Fig. 1f & i. shows that translation reaches the same velocity at Pas and Pas, irrespective of the convective pattern. This highlights the fact that the Rayleigh number is not an independent parameter of the coupled system. Besides, this system expresses differently at equal Rayleigh numbers as in Fig. 2a where point like cold currents inhibit the translation and in Fig. 2f where the reverse prevails. Internal heating has a promoting effect on the translation.
Translation vs convection
Translation induces an asymmetry of the thermal boundary layer with a thickening on the crystallizing side and a thinning on the melting side. When this thinning is enough to reduce the local Rayleigh number of the thermal boundary layer below the critical value, no convective instability can begin. This leads to situations where convective structures are restricted to the crystallizing side as in Fig. 1e & f. This also explains the large difference between Fig. 1d & i, both corresponding to experiments with a Rayleigh number more than 200 times the critical value for the onset of convective structure at a spherical harmonic degree (2607 with impermeable surface conditions (Chandrasekhar1961)). At high Rayleigh number, when the conditions of this thinning are not fulfilled, higher harmonic degree structures, such as coldspots, develop. These uniformly distributed downwellings suppress the thermal anomaly at harmonic degree that translation requires to develop. Thus a real competition is playing between convection and translation.
These are two distinct modes of heat transfer, which can be verified by examining global quantities. For instance, horizontally averaged temperature profiles exhibit very different shapes depending on whether translation operates or not. The classical temperature profile obtained for convection in a full sphere is marked by a top thermal boundary layer and a maximum located just below. When the translation develops, this shape evolves towards an elliptic integral profile where the maximum is now reached at the center (see in Fig. 2 changes from profiles b to e, then f). This is the analytical solution of the temperature equation when conduction and time dependence are neglected. An additional major effect is the strong cooling produced by the translation as illustrated by the comparison of profiles of experiments run at the same Rayleigh number (a and f, or d and i). The permeable boundary that allows advective heat transfer across the surface is responsible for this cooling. This result was already shown in a former study, where such permeable conditions were used to mimic the cooling due to intense volcanism on Io (Monnereau2002a). In this case, the exponent of the power scaling relation between averaged temperature and Rayleigh number is with permeable conditions, instead of the classical value found with impermeable conditions (McKenzie1974).
The averaged temperature decreases proportionally with only when translation dominates the circulations (see Fig. 3). Elsewhere the power low exponent is . Actually, the non dimensional averaged temperature is just proportional to , since viscosity has no effect on translation. This relation can be easily understood by considering that, when translation dominates, the energy equation comes down to a balance between internal heating and advection, so that the product of velocity and temperature is proportional to . This also expresses the fact that all the energy produced internally is evacuated by advection. As the translation velocity is proportional to the temperature (see Eqs  & ), the temperature is proportional to . Finally, the non dimensional temperature is obtained by dividing the temperature by the temperature scale . When the geometry does not allow the translation mode, as in the case of convection between spherical shells, or for all convective modes at degree , the velocity is proportional to , so that the non-dimensional temperature varies in proportion to .
The origin of the difference between these modes lies in the fact that the permeable boundary conditions act differently at degree and degrees . In the case of translation, the radial velocity across the ICB is simply proportional to the topography through Eq. . For convection (), the topography is also proportional to the vertical stress Eq.  and thus to . This last relation expresses the reverse of Eq. : in absence of phase change, i.e. when , the topography is maximal when the velocity across the ICB is zero. The antagonism of both relations reduces the amplitude of the radial velocity at degree . Further, since the stress is proportional to the viscosity, the ICB becomes impermeable at low viscosity, as observed in Fig. 4.
Most of the circulations shown in Fig. 1 are time dependent except those run with a viscosity Pas. In a relatively limited range of viscosities, we observe periodic regimes. For instance at Pas the energy is exchanged between degree 2 to degree 1 in an oscillating regime of period of 1.5 Gyr. Such a period comparable with the age of the inner core age is irrelevant because inner core growth is not modeled. At Pas translation evolves pseudo-steadily while higher degrees develop chaotically. Circulation becomes completely chaotic below Pas. The time scale of variation are, in these cases, smaller than 100 Myr and thus appears rapid at the scale of the inner core age.
Effect of convection velocity in the outer core
Convection velocity drives material transfer through the phase change, owing to its capacity to extract or bring latent heat. It consequently drives surface radial velocity and mainly affects the translation. For instance at low viscosity ( Pas), we observed that increasing has almost no impact on convection, except on the surface radial velocities, which are proportional to . But neither the shape or the number of structures are affected. On the other hand, its enhancing effect on the translation is clear. This one is absent at m/s, just appears at m/s for an internal heating rate W/m, with convective structure equally distributed on both hemisphere (Fig. 1c). At m/s, translation is present whatever the internal heating rate, but dominates for an internal heating rate W/m, as revealed by the asymmetric repartition of convective structures,
Fig. 5 summarizes experiments performed in this study. Globally, it highlights the fact that i) viscosity weakens the convection, ii) the liquid velocity at the inner core surface promotes translation and iii) internal heating promotes both convection and translation, but with a larger effect on the latter. We can see on this figure that most of the situations found at given values of , and are equivalent to the ones found at , , and
Implication for the viscosity of the inner core
The seismic asymmetry of the inner core and particularly the correlation of velocity with attenuation of compressive wave traveling in the uppermost part of the inner core has been interpreted as the signature of iron crystal growth during their transit from a freezing side (facing Peru) to a melting side (facing Indonesia). This model requires that translation be dynamically possible. A first condition is that the inner core be in a superadiabatic state, corresponding here to (see Eq. ). A second condition concerns the viscosity of solid iron at inner core conditions, . If the viscosity is too low, high harmonic degree convection takes place, preventing the development of the translation. Whatever the internal heating considered here (i.e. the core cooling rate), almost pure translation prevails, above a viscosity threshold that only depends on the vigor of the outer core convection. In Fig. 5, this threshold appears where the constant product of the average liquid iron velocity above the ICB and the viscosity exceeds Pam. Since the most accepted value for lies around m/s, this threshold corresponds to Pas.
Just below, at Pas, harmonic degree establishes or alternates with the harmonic degree (translation) in a periodic or quasi periodic regime. Because of their time scale exceeding the age of the inner core, these situations appear irrelevant and would necessitate a modeling that takes into account inner core growth. In any case, a circulation dominated by a degree would not satisfy the observed seismic asymmetry.
Translation is recovered at viscosity one order of magnitude lower, Pas, but embedded with short wavelength convective structures that develop on the freezing side only (see cross section of Fig. 1e). The resulting deformation precludes crystal growth, but may imprint a texture that could be responsible for the seismic anisotropy of the inner core. Seismic anisotropy may be subject to a hemispherical variation, with the western hemisphere displaying much stronger anisotropy than the eastern hemisphere (Tanaka1997; Deuss2010). On Fig. 1e, the hotter hemisphere is not sheared by convective currents, material being simply advected toward the melting side. Annealing during this transit would erase the texture and allow crystal growth (Bergman2010). This asymmetric convection is still present at lower viscosity, Pas, but only for values of ten times larger than the classically expected at the base of the outer core. Such a vigorous circulation appears in opposition with the possible presence of a stable high density layer above the ICB, produced by the fusion of inner core material depleted in light elements, for which the translation of the inner core has also been invoked (Alboussiere2010). For the expected value, a degree remains for Pas, but short wavelength convective structures stir both hemispheres preventing any asymmetry of crystal size and orientation and the formation of a dense layer at the surface of the inner core.
In this study, the external environment was considered as homogeneous. Global circulation in the liquid core may induce variations of the thermodynamic conditions at the surface of the inner core, which are able to significantly disturb its growth rate (Aubert2008; Gubbins2011). External forcing (Yoshida1996; Aubert2008; Deguen2011b) may act in conjunction with thermal convection and may stabilize degree 2 and degree 1 together.
If translation generates the seismic asymmetry of the inner core and if its anisotropy results from convection, our numerical simulations show that viscosity is comprised between and Pas for reasonable values of conductivity and age of the inner core. For viscosities lower than Pas the ”sluggish” inner core is impermeable, meaning that exchanges with the outer core are small, it is deformed and has no hemispherical pattern. For viscosities higher than Pas the rigid inner core is permeable, non-deformed but hemispherical. For intermediate viscosities ( Pas) particular conditions implying a boundary layer close to its Rayleigh number, allows convection to develop over a bulk translation. This configuration possibly makes thermal convection a process able to explain seismic dichotomy and large-scale variation of anisotropy.
Appendix A Force balance of the inner core
Convection within the inner core affects the density field. The harmonic degree of these variations will shift the center of mass of the inner core. To visualize this shift, let us consider the inner core as a solid sphere of average density , entirely covered with a liquid ocean (the outer core) of density (see Fig. 6). The density of the solid, marked by a variation at a spherical harmonic degree due to a temperature variation, is:
where are the components of degree of the temperature field. More simply, we may write:
where is the angle made by the radius originating at the center of the figure of the solid and the direction of the axis defined by and the center of mass of the system , being the radial profile of the temperature heterogeneity at degree . The ocean surface is an equipotential of the gravity field and thus is centered on the center of mass . The position of the center of mass is the one that reduces to zero the first torque of the mass anomaly that represents the inner core within the outer core, the contribution of the outer core being zero:
Introducing the radial vector leads to:
The distance also corresponds to the amplitude of the harmonic degree of the topography of the solid sphere referenced at ; this allows to relate to through the coordinate of relatively to .
Appendix B Phase change kinetics 1D model
A delicate point in solving this problem concerns the computation of the inner core topography. At the harmonic degree , the topography is given by the offset position of center of mass, whereas it results from the convection velocity, at degrees . The solution would consist in writing an integral boundary condition via a topography differential equation:
where is the surface horizontal velocity, the radial velocity across the ICB and the solidification/melting velocity. If the kinetics of the phase change is fast, the above equation comes down to a simple balance between and , so that the topography can be calculated in the spectral domain. In order to verify that adjustment of melting front is almost instantaneous compared to inner core convection time scale, we developed a 1D model of the time evolution of the melting/solidification front.
We solve the energy equation in which the latent heat, , has been introduced as in Christensen1985 :
is a function describing the solid fraction:
is the pressure offset due the phase change, normalized by its pressure thickness . is the Clapeyron slope of the pause change, p the pressure, a reference pressure. S(z), the source/sink term, corresponds to the capacity of outer core to provide or extract the phase change latent heat:
and introducing the temperature dependence of latent heat:
allows to write the specific heat and the thermal expansion coefficient modified by latent heat (Christensen1985):
Finally conservation of heat (16) becomes:
We use a Backward Euler scheme to solve differential equation (23). Boundary conditions are differently defined when material enters or leave the box. Solid Iron (the inner core) is supposed to occupy the bottom half of the domain. In the case of solidification (), the top temperature is ascribed to the outer core temperature, which is . At the bottom, we impose the temperature gradient to be equal to the adiabatic gradient. In the case of melting (), the bottom temperature is set to the adiabatic temperature plus an anomaly resulting from internal heating: . At the top, we impose the temperature gradient to be equal the adiabatic gradient.
We performed five simulations varying the vertical velocities — supposed to be the radial velocity across the Inner core Bounadry (ICB)—, the convective velocity within the liquid outer core just above the ICB and the temperature within the inner core. Table 2 summarizes the characteristics of these experiments.
|Exp.||(cm/yr)||(m/s)||Initial T(z)||(h) (m)|
Fig. 7 show the time evolution of temperature profiles in case of melting by decompression (). We see that the consumption of latent heat cools the inner core by a few hundredths of kelvin, shifting up the phase change by a hundred meters. This state is reached in less that years, which is shorter than time scale involved in convective process. A symmetric situation is observed in the case of solidification when . The time necessary to reach a steady state depends on the difference between the temperature of the inner core material and the adiabat. This time is four time larger for 1K of difference. On the other hand, we find, as expected, that the topography reached is non dependent on this parameter and that it is just proportional to and to
Equilibrium specific time, on the order of to years, remains short enough when compared to convection time scales to consider thermal equilibrium of the topography as instantaneous.