Subcritical thermal convection of liquid metals in a rapidly rotating sphere

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.

💡 Research Summary

This paper investigates thermal convection of liquid metals in a rapidly rotating spherical shell, a configuration relevant to the interiors of terrestrial planets. The fluid is modeled as a Boussinesq liquid with low Prandtl number (Pr ≪ 1) and low Ekman number (Ek ≪ 1), heated internally. The governing equations are the incompressible Navier‑Stokes equations with Coriolis force, a temperature equation with internal heating, and appropriate non‑dimensional parameters: Ek = ν/(r₀²Ω), Pr = ν/κ, and Rayleigh number Ra = αg₀Sr₀⁶/(6ρCₚνκ²). No‑slip, impermeable velocity boundary conditions and fixed temperature at the outer sphere are imposed.

The authors first compute the linear onset of convection for each (Ek, Pr) pair using a dedicated eigenvalue solver, obtaining the critical Rayleigh number Ra_crit and the associated eigenmode (thermal Rossby wave). These linear modes are then used as initial conditions for fully three‑dimensional direct numerical simulations performed with the XSHELLS code. Hyper‑viscosity is applied only to the highest spherical harmonic degrees to accelerate the runs, and convergence tests confirm that the results are insensitive to this treatment. Resolutions up to ℓ = 199 and 1152 radial points are employed, requiring up to 576 CPU cores.

Two distinct nonlinear solution branches emerge from the simulations:

1. Weak branch – Appears just above Ra_crit. The flow consists of columnar thermal Rossby waves, with a Péclet number Pe that scales as (Ra − Ra_crit)¹ᐟ². The kinetic energy is modest, and the Nusselt number is only slightly above unity. This branch is the continuation of the linear instability.

2. Strong branch – At slightly higher Ra the weak branch loses stability and the system jumps to a state with much larger velocities (Pe > 10) and vigorous heat transport (Nu ≫ 1). The flow remains largely columnar but exhibits a pronounced prograde zonal jet near the rotation axis. The strong branch persists even when Ra is reduced below Ra_crit, i.e., it is subcritical.

The transition between the branches depends critically on Ek. For moderate rotation (Ek ≈ 10⁻⁵) the two branches are smoothly connected; the strong branch can be reached continuously as Ra increases. When Ek ≤ 3 × 10⁻⁶, the transition becomes discontinuous, producing a hysteresis loop: the system can remain on the strong branch for Ra < Ra_crit, or revert to the weak branch (or to a quiescent state) if perturbed sufficiently. The hysteresis width shrinks as Ek is further reduced, but the strong branch remains subcritical down to Ra ≈ 0.69 Ra_crit in the most extreme cases (Ek = 10⁻⁷, Pr = 0.01).

A key finding is that the Reynolds stress generated by the non‑axisymmetric convective motions is essential for sustaining subcritical convection. Numerical experiments that remove the nonlinear advection term (u·∇)u cause the strong branch to disappear, whereas eliminating the temperature advection term (u·∇)Θ actually enhances convection. This demonstrates that the self‑sustaining mechanism is kinetic rather than thermal. Moreover, the authors identify a threshold Pe ≈ 10: as long as the convective Péclet number stays above this value, the strong branch can survive below Ra_crit; once Pe falls below 10, the branch collapses.

Zonal flows, while prominent in the strong branch, are shown to be a by‑product rather than a driver of the subcritical transition. When the axisymmetric (m = 0) velocity components are artificially set to zero at each time step, the convective vigor and heat transport remain essentially unchanged. The scaling analysis reveals that the zonal Péclet number Pe_zon scales with the convective Péclet number Pe_conv as Pe_zon ∝ Pe_conv³ᐟ², a relationship that holds across all Ek, Pr, and both branches.

The paper concludes that in the asymptotic regime of rapid rotation and low Prandtl number, thermal convection in a sphere can onset via a subcritical bifurcation, sustained by Reynolds stresses rather than by the mean zonal flow. This behavior contrasts with earlier predictions for moderate Ekman numbers, where mean flows can weaken the stabilizing Coriolis effect. The results have important implications for the early thermal evolution of planetary cores and for the onset of dynamo action, suggesting that vigorous convection may be present even when linear theory predicts stability.