Subcritical thermal convection of liquid metals in a rapidly rotating sphere

TL;DR

This study investigates nonlinear thermal convection in rapidly rotating liquid metal spheres, revealing subcritical bifurcation, hysteresis, and the dominance of advection over wave-like convection near planetary core conditions.

Contribution

It demonstrates that in rapidly rotating planetary cores, convection can persist below linear onset through subcritical bifurcation and hysteresis, highlighting the role of Reynolds stress.

Findings

Convection onsets as thermal Rossby waves near critical Rayleigh number.

At higher Ra, convection transitions to advection-dominated states.

Strong branch convection persists below linear onset, indicating subcriticality.

Abstract

Planetary cores consist of liquid metals (low Prandtl number $Pr$) that convect as the core cools. Here we study nonlinear convection in a rotating (low Ekman number $Ek$) planetary core using a fully 3D direct numerical simulation. Near the critical thermal forcing (Rayleigh number $Ra$), convection onsets as thermal Rossby waves, but as the $Ra$ increases, this state is superceded by one dominated by advection. At moderate rotation, these states (here called the weak branch and strong branch, respectively) are smoothly connected. As the planetary core rotates faster, the smooth transition is replaced by hysteresis cycles and subcriticality until the weak branch disappears entirely and the strong branch onsets in a turbulent state at $Ek < 10^{-6}$. Here the strong branch persists even as the thermal forcing drops well below the linear onset of convection ($Ra=0.7Ra_{crit}$ in this study). We highlight the importance of the Reynolds stress, which is required for convection to subsist below the linear onset. In addition, the P\'eclet number is consistently above 10 in the strong branch. We further note the presence of a strong zonal flow that is nonetheless unimportant to the convective state. Our study suggests that, in the asymptotic regime of rapid rotation relevant for planetary interiors, thermal convection of liquid metals in a sphere onsets through a subcritical bifurcation.

Figures

Tables

Table 1. Table 1: Critical parameters at the onset of convection

$Ek$	$Pr$	$R{a}_{crit}$	$m_c$	${\omega }_c/\Omega$
${10}^{-5}$	$0.1$	$8.440\times {10}^6$	11	$-0.04024$
$3\times {10}^{-6}$	$0.03$	$2.336\times {10}^7$	12	$-0.04275$
${10}^{-6}$	$0.01$	$5.475\times {10}^7$	11	$-0.03895$
$3\times {10}^{-7}$	$0.003$	$1.255\times {10}^8$	12	$-0.04287$
${10}^{-7}$	$0.01$	$1.01574\times {10}^9$	27	$-0.01930$
${10}^{-7}$	$0.01$	$1.01567\times {10}^9$	28	$-0.01941$

Equations

$${\partial }_t\mathbf{u}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}+\frac{2}{Ek}\hat{\mathbf{z}}\times \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$${\partial }_t\Theta +\mathbf{u}\cdot \nabla \Theta -\frac{2}{Pr}r{u}_r$$

Full text

Self-sustaining magnetic fields of terrestrial planets are generated in their liquid metal cores. Left on their own they would decay from ohmic dissipation, but the slow cooling of the planets drives the convective flows thought to maintain the fields Jones08 . Numerical models of planetary dynamos Glatzmaier95 ; Miyagoshi10 are widely used to study these strongly nonlinear systems. They solve the Navier-Stokes equation coupled to a temperature equation and to the induction equation that governs the magnetic field. For simplicity, most direct numerical simulations (DNS) of the dynamos have set the Prandtl number $Pr=1$. However, liquid metals have $Pr\lesssim 0.1$ and the nature of their convection differs from $Pr=1$ chandrasekhar1961 ; Guervilly.JFM.2016 ; calkins2016 .

The onset of convection in a full sphere is relevant to the early history of the Earth’s Olson.Science.2013 or Moon’s Laneuville2014 core, and has thus received a great deal of attention (e.g. Jones.JFM.2000 ). A significant thermal forcing is required to overcome the stabilizing rotational constraint and drive convective instabilitiesBusse70 ; Jones.JFM.2000 . At and near this threshold, convection onsets in the form of columnar vortices aligned with the axis of rotation, drifting in the azimuthal direction Busse70 ; Dormy04 . The nonlinear regime has been extensively studied for $Pr=1$ gastine2016 . The $Pr\ll 1$ regime is more difficult to tackle Zhang92 , but its nonlinear state was recently described using a quasigeostrophic model, which relies on a two dimensional description of the axial vorticity. In their simplified model, the authors of Guervilly.JFM.2016 found first clues of subcritical convection—that is convection below the linear onset of instability—anticipated by weakly nonlinear theoretical predictions Soward77 ; Plaut08 .

In the present work on the rotating convection problem, we use three-dimensional direct numerical simulations to describe the nature of the weak and strong convective branches, especially when the strong branch becomes subcritical. We discuss the insights gained from the differences and similarities between the two branches.

I Formulation of the Model

We study thermal Boussinesq convection driven by internal heating in a sphere rotating at constant angular velocity $\Omega\mathbf{\hat{z}}$. The acceleration due to gravity is radial and increases linearly, as expected for a constant density medium, $\mathbf{g}=g_{0}r\mathbf{\hat{r}}$. The radius of the sphere is $r_{o}$. The fluid has kinematic viscosity $\nu$, thermal diffusivity $\kappa$, density $\rho$, heat capacity at constant pressure $C_{p}$, and thermal expansion coefficient $\alpha$, all of which are constant. We consider an homogeneous internal volumetric heating S.

The dynamic system is nondimensionalized by scaling lengths with $r_{o}$, times with $r_{o}^{2}/\nu$, and temperature with $\nu Sr_{o}^{2}/(6\rho C_{p}\kappa^{2})$Guervilly.JFM.2016 . The system is governed by the incompressible Navier-Stokes equation and an advection-diffusion equation of the temperature perturbation,

[TABLE]

where $\mathbf{u}$ is the velocity field, $p$ is the modified pressure, which includes the centrifugal potential, and $\Theta$ is the temperature perturbation relative to $T_{s}$. The dimensionless numbers are the Ekman number $Ek=\nu/r_{o}^{2}\Omega$, the Prandtl number $Pr=\nu/\kappa$, and the Rayleigh number $\left(Ra={\alpha g_{0}Sr_{o}^{6}}/{6\rho C_{p}\nu\kappa^{2}}\right)$. At $r=r_{o}$, the boundary condition for the velocity is no-slip and impenetrable ($\mathbf{u}=\mathbf{0}$) and the temperature is fixed ($\Theta=0$).

This dynamic system is modeled in a 3 step process. First, for each set of $(Ek,Pr)$, the linear onset of convection $Ra_{crit}$ and the associated eigenmode are computed precisely by using the singe code Vidal.GJI.2015 to find the value of $Ra$ for which the least damped eigenmode of the linearized equations has almost zero growth rate (see supplementary table 1). The linear eigenmodes are then used as initial values of fully 3d simulations, which are then run until they reach a statistical equilibrium. Hysteresis cycles and subcriticality are explored by changing the thermal forcing or rescaling the final flow state and evolving the system to its new equilibrium.

The fully three-dimensional simulations are run with our spherical code xshells schaeffer2013 ; marti2014 . Most runs employ a form of hyper-viscosity affecting only the 20% highest spherical harmonics kaplan2016 to speed up computations, but we have checked that this does not alter the solutions. Similarly, careful spatial and temporal convergence checks have been made, and we were surprised to find that the amplitude of the zonal winds (axisymmetric azimuthal flow) is very sensitive to the bulk radial grid spacing, although no sharp gradients were seen. We suspect that the balance between the small viscous stress and the small Reynolds stress requires a high resolution to be computed accurately with second order finite differences. As a result, simulations were run using up to 576 cores for 1152 radial levels and spherical harmonics up to degree 199.

The simulations output several useful diagnostic values over the course of the run. We present the velocity using the dimensionless Péclet number ($Pe=Ur_{o}/\kappa$) which is the ratio between the rate of advection to the rate of diffusion of temperature. $U$ is defined by the square root of the volume averaged kinetic energy; the zonal Péclet number $\left(Pe_{zon}\right)$, uses only the azimuthal component of the flow averaged in longitude ($m=0$); the convective Péclet number $\left(Pe_{conv}\right)$ uses the $m\neq 0$ velocity. We also use the Nusselt number $\left(Nu\equiv\frac{1}{Pr\Theta_{c}}\right)$ which measures the convective heat transfer from the core to the surface. ($\Theta_{c}$ is the temperature at the origin).

II Weak and Strong Branches

Our studies seek to verify that the subcritical behavior found in the quasigeostrophic approximation Guervilly.JFM.2016 is also seen in the fully three dimensional system, for $Ek\in\left[10^{-7},10^{-5}\right]$ and $Pr\in\left[0.003,0.1\right]$.

Our results are summarized in Fig. 1, which shows the Péclet number $\left(Pe\right)$ vs the Rayleigh number $\left(Ra\right)$ for the simulations we’ve run. Two branches are clearly visible. a strong branch where the $Pe$ number is larger than 10 and a weak branch with lower $Pe$ numbers. The advantage of representing the kinetic energies in terms of $Pe$ is that it collapses all of the simulations onto a single scale with a clear divide between the strong and weak branches. A more typical Reynolds number representation sees values of kinetic energy that are in either the weak branch or the strong branch depending on $Pr$; an example of this is in the supplementary figure S1. The $Pe$ is also important because one of the markers of the strong branch is a significant cooling of the sphere’s core; $Pe>10$ indicates that there is enough convective power to draw thermal energy away from the center. At lower $Ek$, the strong branch persists below the linear critical Rayleigh number.

The weak branch arises from a supercritical bifurcation at $Ra_{crit}$. For Rayleigh numbers near this value, convection onsets as a thermal Rossby wave Busse70 . The thermal anomaly $\Theta$ in the equatorial plane for a typical case is visible in Fig. 2a Zhang92 . The azimuthal velocity $v_{\phi}$ in a single meridional slice is shown in Fig. 2b. The $z$ invariance implied by the quasigeostrophic approximation is well displayed. Near onset, the flow driven by the thermal gradient scales roughly with $\left(Ra-Ra_{crit}\right)^{1/2}$ Gillet06 , as shown in Fig. 1.

At large Rayleigh numbers, the weak branch becomes unstable Guervilly.JFM.2016 and the flow evolves towards a state where the velocity is an order of magnitude greater. A snapshot of a typical strong branch flow is shown in Fig. 3. The difference with Fig. 2 is stark. There is strong cooling near the origin of the sphere, and a noticeable prograde zonal flow near the axis of rotation. There is some departure from $z$ invariance, but the flow is still mostly columnar.

At $Ek=10^{-5}$ the strong branch is smoothly connected with the weak branch. Below this $Ek$ the strong branch onsets with a discontinuity, and is characterized by $Pe\gtrapprox 10$. This state onsets as a subcritical bifurcation at $Ek\leq 3\times 10^{-6}$, allowing a small hysteresis cycle. The only green open-faced symbol in Fig. 1 represents a simulation that persists in the strong branch for $1.5\tau_{\kappa}$ ($\tau_{\kappa}=r_{o}^{2}/\kappa$ is the thermal diffusion time) before suddenly decaying to the weak branch. At $Ek=10^{-6}$ the hysteresis persists below the onset of the weak branch instability.

At lower Ekman numbers ($Ek<10^{-6}$), the weak branch is absent and the thermal convection operates on the strong branch always at $Pe>10$. Indeed, we have looked very close above the onset for $Ek=10^{-7}$ and $Pr=0.01$: even with an initial perturbation of very low amplitude ($Re\ll 1$) the kinetic energy quickly increased to reach the strong branch. This demonstrates that the saturation mechanism leading to the weak branch at moderate Ekman numbers is lost at higher rotation rates.

The phase trajectories of several simulations at $Pr=0.03,Ek=3\times 10^{-6}$ are shown in Fig. 4. A simulation at $Ra=1.03Ra_{crit}$, initialized from a weak branch state is shown in Fig. 4a. The limit cycle is simple and stable; as the convective power of the flow (indicated by $Pe$) increases, the core cools (indicated by an increase in $Nu$), weakening the convection, letting the core heat up again. Fig. 4b shows a simulation at $Ra=1.1Ra_{crit}$, initialized from a weak branch state. The limit cycle here is more complex. The transition to the strong branch happens when other growing modes besides the critical one reach similar amplitude. Simulations at the same two Rayleigh numbers, but initialized in the strong branch, are shown in Fig. 4c & d. Here we see that the phase trajectories are stochastic rather than approaching a limit cycle. At $Ra=1.03Ra_{crit}$, the system persists in the strong branch for 1.5$\tau_{\kappa}$ before suddenly jumping to the weak branch. A second observation from Fig. 1 is that the fluctuation levels in the strong branch actually decrease as the thermal forcing is increased. These fluctuations are categorically different from the fluctuations of the weak branch. Here, the fluctuations are stochastic in the phase spaces of the various diagnostic parameters.

III Subcriticality

As Fig. 1 shows, subcritical convection exist in our system, as convection on the strong branch can occur below the linear onset of convection, at $Ra<Ra_{crit}$. Interestingly at $E/Pr=10^{-4}$ the subcriticality is small and fragile: a small perturbation can kick the system back to rest. This is indicated by the open symbols below $Ra<Ra_{crit}$ (Fig. 1), where the flow transitions to rest because of its intrinsic fluctuations, sometimes after several thermal diffusion times.

A much more robust subcriticality is observed at $E/Pr=10^{-5}$ and $Pr=0.01$. There, the strong branch can be violently perturbed without loosing the convection. As an example, we can divide flow and temperature anomalies by 10 at $Ra=0.94Ra_{crit}$ and the convection will quickly recover, while dividing instead by 100 sends the system back to rest. By gradually lowering $Ra$, we have observed robust convection at decreasing $Ra$ down to $Ra=0.69Ra_{crit}$ and lost the convection at $Ra=0.64Ra_{crit}$, when the Péclet number dropped below $Pe\simeq 10$. Further stress tests were applied at $Ra=0.94Ra_{crit}$. We artificially removed all $m=0$ (axisymmetric) components, including the zonal flow and zonal temperature anomaly: the convection stayed firm. We also kept only the most unstable mode at onset and its harmonics: the convection endured. When removing the $(\mathbf{u}.\nabla)\Theta$ term, the kinetic energy increases significantly. Instead, when removing the $(\mathbf{u}.\nabla)\mathbf{u}$ term, the subcritical convection dies out.

In the 2D planar geometry studied by Chandrasekhar chandrasekhar1961 , Veronis predicted a possible subcriticality in a window of low rotation rates veronis1959 , recently confirmed numerically beaume2013 . The subcritical behavior is associated with the presence of large mean flows which reduce locally the effective rate of rotation and consequently, the rotational constraint on the flow and the critical Rayleigh number toward its non rotating value. This mechanism implies a Rossby number $Ro$ close to unity. In contrast, it is very small here ($Ro=Re.Ek<10^{-3}$) because of the high rotation rate ($Ek=10^{-7}$) and despite the large Reynolds numbers ($Re\sim 10^{4}$) observed near the convection onset. It is thus not surprising that the zonal flow is not important for subcritical convection. Furthermore, outside the boundary layers and in the subcritical regime, the root-mean-square local vorticity fluctuations never exceed $5\%$ of the background vorticity $2\Omega$. Not only did Veronis find no subcritical motion in this rapidly rotating regime veronis1968 , but the mechanism of lowering the effective (local) rotation rate cannot explain the large amount of subcriticality we found (convection down to $Ra=0.69Ra_{crit}$).

Our numerical experiments highlight the key role of the Reynolds stress to sustain convection below the linear onset. However, it is the Péclet number that is more conveniently used to characterize the strong branch and subcriticality. If $Pe$ is much larger than $10$ at the linear onset ($Pe(Ra_{crit})\simeq 50$), it is possible to have convection at $Ra$ well below $Ra_{crit}$, as long as $Pe\gtrsim 10$. Below that threshold, the strong branch of convection cannot survive.

Our study also shows the trend of lower and lower subcritical $Ra/Ra_{crit}$ as $Ek/Pr$ is reduced while keeping $Pr\sim 0.01$. This suggests an even larger effect at planetary core conditions ($Ek/Pr<10^{-10}$), currently out of reach for numerical models.

IV Conclusions

Subcriticality, like hysteresis, implies the presence of active nonlinearities. The zonal flow in this system is produced by nonlinear interaction of the convective velocity (Reynolds stress). In Fig. 5 we plot the $Pe_{zon}$ as a function of $Pe_{conv}$ over the full set of simulations carried out. The data points seem to align to an inertial scaling law with a power of 3/2 Gillet06 . This scaling law persists not only between the different $Ek$ and $Pr$ numbers, but between the strong and weak branches as well. A set of subcritical runs were rerun with their zonal components artificially canceled (set to zero at every time step) with no significant change in the convective flow or core temperatures. These two factors combined show that the zonal flow is only ever a byproduct of the convective flows and not a driver of any of the dynamics. This is fundamentally different from the subcriticality predicted for moderate Ekman numbers ($Ek\sim 1$), for which the mean flow weakens the stabilizing effect of global rotation veronis1959 ; beaume2013 .

V Acknowledgements

We thank three anonymous reviewers for their constructive comments. The SHTns, SINGE and XSHELLS codes are freely available at https://bitbucket.org/nschaeff/. We acknowledge GENCI for awarding us access to resource Occigen (CINES) and Turing (IDRIS) under grant x2015047382 and x2016047382. Part of the computations were also performed on the Froggy platform of CIMENT (https://ciment.ujf-grenoble.fr), supported by the Rhône-Alpes region (CPER07_13 CIRA), OSUG@2020 LabEx (ANR10 LABX56) and Equip@Meso (ANR10 EQPX-29-01). This work was partially supported by the French Agence Nationale de la Recherche under grants ANR-13-BS06-0010 (TuDy) and ANR-14-CE33-0012 (MagLune).

Appendix A Supplementary Materials

A.1 Linear results

The critical parameters of the linear onset of convection shown in table 1 have been determined using the freely available singe code Vidal.GJI.2015 , which has been benchmarked with the results of Jones et al Jones.JFM.2000 .

For $Ek=10^{-7}$, $Pr=0.01$, the onset of modes $m=27$ and $m=28$ are very close to each-other, and we list both for completeness.

A.2 Complementary nonlinear results

A more typical measure of the convective velocity is the Reynolds number ($Re=Ur_{o}/\nu$) which is in our case related to the Péclet number by $Re=Pe/Pr$. Fig. S1a shows the convective velocity versus thermal forcings $Ra$. Near the linear onset of convection, the strong branch flows seem to align with a single trendline for each ratio $Ek/Pr$ that seems to scale with $Ra^{3/2}$. As the $Ra$ increases, the flows depart from this ’scaling law’, most likely as a result of Ekman pumping. Fig. S1b shows the convective heat transfer as a function of convective velocity. Large forcings (strong branch) seems to get a lower exponent than the 2 exponent expected by a weakly non linear analysis at the onset of convection Gillet06 . These trendlines should be seen as suggestions to guide the eyes, rather than hard and fast scaling laws that the data must conform to.

All the data behind Figs. 5 & S1 and a jupyter notebook reproducing the figures is available at https://dx.doi.org/10.6084/m9.figshare.4540846 for anybody to plot data as needed.

Finally, to complement the fields shown in figure 2 and 3, we represent thermal anomaly in the meridional plane and vertical vorticity in the equatorial plane in figure S2 and S3.