Onset of convection (2-D)

A layer heated from below convects only if buoyancy wins over viscous resistance and thermal diffusion. Below a critical Rayleigh number every temperature perturbation decays and heat crosses the shell by conduction alone; above it the perturbation grows exponentially until convection reaches a finite amplitude. This example finds that threshold in the spherical annulus by measuring the growth rate of one perturbation at six Rayleigh numbers.

Ra = 4000, temperature

Temperature (300 to 1300 K) at Ra = 4000, Nr = 64 (64 × 512 cells). The wavenumber-4 perturbation grows on the conductive profile, saturates after about 3.2 Gyr (0.1 diffusion times) and settles into four steady cells.

Set-up

shell	$R_{\mathrm{in}}=1220$ km, $R_{\mathrm{out}}=2220$ km ($f=0.5495$, $D=1000$ km), spherical annulus
boundaries	free slip; 1300 K at the core-mantle boundary, 300 K at the surface
fluid	Boussinesq, constant viscosity set by eta_for_Ra
initial state	the conductive profile of the spherical annulus plus a 1 K perturbation of azimuthal wavenumber $m=4$ near the core-mantle boundary
transport	temperature on the grid with the TVD scheme, so no marker noise seeds other modes
run	0.3 thermal diffusion times, $0.3D^2/\kappa \approx 9.5$ Gyr

The conductive state of the spherical annulus is not linear in $r$. With the measure $r^{2}dr$, steady conduction $\frac{1}{r^2}\frac{d}{dr}(r^2\frac{dT}{dr})=0$ gives

$$T_{\mathrm{c}}(r)=T_{\mathrm{s}}+(T_{\mathrm{b}}-T_{\mathrm{s}})\frac{1/r-1/R_{\mathrm{out}}}{1/R_{\mathrm{in}}-1/R_{\mathrm{out}}},$$

which is what temp.initial_profile = "conductive" sets. Starting from the exact conductive state of the discrete problem means a run below the threshold has nothing to do but let the perturbation decay.

Measuring the growth rate

Near the conductive state the equations are linear in the perturbation, so each mode grows or decays as $e^{st}$. After the initial perturbation has projected onto the slowest mode of wavenumber 4, the rms velocity follows the same law,

$$V_{\mathrm{rms}}(t)\propto e^{st},\ln V_{\mathrm{rms}}=st+\mathrm{const},$$

and $s$ is the slope of $\ln V_{\mathrm{rms}}$ against time. The script docs/scripts/fig_onset.jl fits it by least squares from each run's diagnostics.csv. It uses the last two thirds of a run, or, for a run that grows by more than a factor 100 and saturates, the linear phase only (from a tenth of the run to the time $V_{\mathrm{rms}}$ first reaches 5 % of its maximum). The critical Rayleigh number is where $s$ changes sign, interpolated linearly in $\ln \mathrm{Ra}$ between the two runs that bracket it.

A: rms velocity of the six runs, logarithmic scale. Straight lines are exponential growth or decay; Ra = 2000 and 4000 bend over where they reach finite amplitude. B: fitted growth rate in units of the diffusion time $D^2/\kappa$; the red line marks where it vanishes.

Ra	600	800	1000	1400	2000	4000
$s$ (per $D^2/\kappa$)	−2.73	+1.52	+5.79	+14.3	+27.2	+70.0

The growth rate changes sign between Ra = 600 and 800, giving ${\mathrm{Ra}}_c\approx 722$ for wavenumber 4.

How it compares

For a plane layer with free-slip, isothermal walls the linear-stability threshold is ${\mathrm{Ra}}_c=27{\pi }^4/4\approx 657.5$, reached at the wavelength $2\sqrt{2}D\approx 2.83D$ (Chandrasekhar 1961; Turcotte & Schubert 2014). Wavenumber 4 in this shell has a wavelength of $2\pi \cdot 1720\mathrm{km}/4\approx 2.7D$ at mid-depth, close to that most unstable plane-layer wavelength. The annulus threshold is about 10 % higher than the plane-layer value. In the curved shell the 1/r conductive profile puts most of the temperature drop in the lower half, where the boundary is shorter, so the plane-layer number is a guide and not a target. The measurement is the threshold of one fixed wavenumber; other wavenumbers have their own, and the lowest of them is the threshold of the shell.

Running it

for ra in 600 800 1000 1400 2000 4000; do
  MILEVA_OUTDIR=output/onset_Ra$ra julia --project=. -t 2 \
      examples/annulus_onset_2d.jl $ra 4 64
done
julia --project=. docs/scripts/fig_onset.jl output

Arguments are the Rayleigh number, the wavenumber $m$ and Nr. Each run took about 10 minutes on 2 threads.

References

Chandrasekhar, S. (1961), Hydrodynamic and Hydromagnetic Stability, Oxford, Clarendon Press. Turcotte, D. L. & Schubert, G. (2014), Geodynamics, 3rd edn., Cambridge University Press, section 6.19, doi:10.1017/CBO9780511843877.