A slice through the Earth with a live thermal-convection solver in the mantle. Tune the Rayleigh number from just below the onset of convection to vigorous plume-and-slab flow, inject a hot plume or a cold sinking slab, and check the solver against the analytic onset and a published benchmark.
Runs four small steady convection cells (Ra = 10³ to 3·10⁴), measures the Nusselt number and fits Nu ∝ Raβ.
The science
Equations. The mantle is a very viscous fluid, so inertia is negligible (infinite Prandtl number). In units of the layer depth H, the thermal diffusion time H²/κ and the temperature difference ΔT: ∂T/∂t + u·∇T = ∇²T + Hh and Stokes flow −∇p + ∇²u + Ra T ŷ = 0, ∇·u = 0. With a stream function ψ (u = ∂ψ/∂y, w = −∂ψ/∂x) this becomes ∇²ω = −Ra ∂T/∂x, ∇²ψ = −ω, solved by a Fourier transform in x and a tridiagonal solve in y. The Rayleigh number Ra = ρ g α ΔT H³ / (κ η) measures buoyancy against viscous and thermal resistance.
Onset. For free-slip boundaries the linear growth rate of a roll of wavenumber k is σ = −(π²+k²) + Ra k²/(π²+k²)², which vanishes first at Rac = 27π⁴/4 ≈ 657.5, kc = π/√2. The solver reproduces this growth rate and the critical value (see the checks).
Heat transport. The Nusselt number Nu is the surface heat flux divided by the purely conductive value; in steady state it is the same at the top and the bottom. Vigorous convection concentrates the temperature drop into thin thermal boundary layers, and boundary-layer theory gives Nu ∝ Ra1/3; the sweep measures the exponent on a coarse grid.
Benchmark. Blankenbach et al. (1989) published steady-state results for an isoviscous 1 × 1 box with free-slip walls at Ra = 10⁴ (Nu = 4.884, vrms = 42.865). A periodic domain of width 2 containing two mirror-image rolls is equivalent, so the same numbers apply. These published values are quoted from memory and should be spot-checked by a human reader; the solver lands within about 1–2 % of them, with the remaining difference due to the finite grid.
Scale to the real Earth. With H = 2900 km and κ = 10-6 m²/s one unit of time is 266 Gyr and one unit of velocity is 1.09×10-3 cm/yr. The real mantle has Ra of roughly 107 (the value depends on the assumed viscosity), beyond what this real-time solver can resolve; flow speeds scale roughly as Ra2/3, so the cm/yr values shown here are those of a much less vigorous mantle. The Gyr readout is the diffusion time scale, not a prediction for Earth.
Simplifications. Cartesian 2D box wrapped around the ring for display (geometry in a real annulus differs), constant viscosity (no plates, no stagnant lid, no viscosity-temperature dependence), Boussinesq approximation (no compressibility), free-slip boundaries, fixed temperature at top and bottom, no phase changes. The tracers are advected with the model velocity but are only a display aid. The core, crust and continents are drawn for orientation and are not part of the model. Models built with Blender, scene by Three.js.