PLUME

What Is Plume?

Plume simulates Rayleigh-Bénard convection — the spontaneous flow that develops when a horizontal fluid layer is heated from below and cooled from above. The setup is the simplest possible: flat hot plate at the bottom (T = 1, amber), flat cold plate at the top (T = 0, teal), fluid in between. At first, the fluid conducts heat — a linear temperature gradient and no movement. Then, if the temperature difference is large enough relative to the fluid’s viscous and thermal resistance, a phase transition occurs: the motionless conducting state becomes unstable, and any infinitesimal disturbance grows into organized convection rolls.

In Plume, the transition happens when you drag the Ra slider past the critical value Ra_c = 1707.8. Below it: any perturbation decays, the Nusselt number stays at 1, and the temperature field holds its smooth linear gradient. Above it: within seconds, amber plumes rise and teal sheets sink, locking into counter-rotating rolls whose width is approximately twice the layer height. The Nusselt number climbs above 1, measuring how much more efficiently convection moves heat than conduction alone.

The Governing Equations

Plume uses the stream-function / vorticity formulation of the 2D incompressible Boussinesq equations. Instead of tracking velocity components directly, the velocity field is encoded in the stream function ψ (whose contours are streamlines) and the vorticity ω (the curl of velocity). The stream function automatically satisfies incompressibility (∇·u = 0) by construction, so the mass-conservation equation disappears. Only three scalar fields — T, ω, ψ — are needed to describe the entire flow.

// Non-dimensional Boussinesq equations  (H=1, κ=1, ΔT=1, Pr=1)
// Vorticity  ω = ∂v/∂x − ∂u/∂y,  stream function  ∇²ψ = −ω
// Velocity   u = ∂ψ/∂y,  v = −∂ψ/∂x

∂ω/∂t = Pr·∇²ω − u·∇ω + Ra·Pr·∂T/∂x   // buoyancy source
∂T/∂t = ∇²T   − u·∇T                    // heat equation
∇²ψ   = −ω                               // Poisson (solved by SOR)

// Boundary conditions:
T(y=0) = 1  (hot bottom),  T(y=1) = 0  (cold top)
ψ = 0 at walls  (no penetration)
ω_wall = −2ψ_interior / h²  (Thom no-slip condition)
Periodic in x

// Nusselt number  (heat transport ratio):
Nu = −⟨∂T/∂y⟩|_{y=0}  =  1 + convective contribution
Nu = 1 → pure conduction,  Nu > 1 → convection enhances heat flux

The buoyancy term Ra · Pr · ∂T/∂x is the engine of Rayleigh-Bénard instability. Wherever the temperature has a horizontal gradient — hotter on the right — the fluid there is lighter (lower density) than its surroundings, and the curl of the resulting buoyancy force generates counterclockwise vorticity. This vorticity drives the fluid upward at that location. The rising warm fluid enhances the horizontal temperature gradient, which drives more vorticity, which drives more upward flow: a positive feedback that amplifies any seed perturbation once Ra exceeds the critical threshold.

The Critical Rayleigh Number

The Rayleigh number is the dimensionless ratio of buoyancy force to the product of viscous and thermal dissipation:

Ra = g α ΔT H³ / (ν κ)

  g  = gravitational acceleration
  α  = thermal expansion coefficient
  ΔT = temperature difference (bottom − top)
  H  = layer height
  ν  = kinematic viscosity
  κ  = thermal diffusivity

Ra_c = 1707.8  (no-slip horizontal plates, onset of steady rolls)

// At Ra_c: most unstable horizontal wavenumber  k_c = π/√2 ≈ 2.22
// Critical wavelength: λ_c = 2π/k_c ≈ 2.83H  ≈  2H  (approximate)
// In Plume: domain width 4H → fits ≈ 2 roll pairs at onset

The value Ra_c = 1707.8 was computed analytically by Lord Rayleigh in 1916, solving the linear stability problem for the motionless conducting state. It is one of the most precisely known numbers in fluid mechanics. Above Ra_c, the spectrum of modes with k near k_c all grow exponentially, with growth rate σ ∝ (Ra − Ra_c). The initial pattern of the simulation — before nonlinear saturation — reflects whichever of these modes first emerges from the noise.

The Prandtl number Pr = ν/κ sets the relative importance of viscous versus thermal diffusion. Plume uses Pr = 1 (common for gases; water is Pr ≈ 7; Earth’s mantle has Pr → ∞). At Pr = 1, momentum and heat diffuse at the same rate, and the stability threshold is determined purely by the competition between buoyancy and the combined dissipation.

The Nusselt Number and Heat Transport

The Nusselt number Nu is the central observable of Rayleigh-Bénard convection. It measures how much faster convection moves heat than pure conduction: Nu = 1 means the fluid is conducting only; Nu = 5 means convection moves heat five times faster than conduction alone. In laboratory experiments, the first measurable quantity after the onset of convection is the jump in Nu at Ra_c — from exactly 1 to a value just above 1, then growing as Nu ∝ (Ra − Ra_c)^(1/2) near onset, and eventually following Nu ∝ Ra^(2/7) in the turbulent regime (Kraichnan–Spiegel scaling).

In Plume, Nu is computed as the average vertical temperature gradient at the hot bottom wall, normalized by the conductive gradient ΔT/H = 1. Drag Ra slowly past 1708 and watch the Nusselt number climb from 1 — this is the bifurcation made visible in real time.

Regimes: Rolls → Oscillatory → Turbulent

As Ra increases beyond Ra_c, Rayleigh-Bénard convection passes through a cascade of increasingly complex regimes. Between Ra_c and roughly 10 × Ra_c, the flow settles into steady convection rolls — time-independent, geometrically ordered counter-rotating cells. These are the Bénard cells that Henri Bénard first photographed in 1900, heating a thin layer of spermaceti wax.

At higher Ra (≳ 20 × Ra_c), the rolls become time-dependent: they oscillate, merge, split, and drift. At still higher Ra, the flow becomes spatiotemporal chaos — fully turbulent convection where thermal plumes detach from the boundary layers and traverse the layer unpredictably. In Plume, drag Ra to 30000–60000 to watch this transition. The Nusselt number fluctuates rapidly, the temperature field loses its spatial order, and the amber and teal plumes become intermittent and chaotic.

Where This Physics Appears

Rayleigh-Bénard convection is one of the most pervasive instabilities in nature. Earth’s mantle convects on geological timescales: Ra ≈ 10⁷. The slow creep of solid rock driven by radiogenic heat from the core drives plate tectonics — the mantle rolls are what move continents. The plates are the cool top boundary layer detaching and sinking (subduction) while hot upwellings (mantle plumes, hotspots like Hawaii) rise from below.

The Sun’s convection zone extends from 0.7 R_☉ to the surface. Ra ≈ 10²⁰. Granules visible in solar images are the tops of convection cells — rising hot plasma (bright granule center), radiating, cooling, sinking at the dark intergranular lanes. The 5-minute oscillations of the solar surface are partly driven by the convective overturn. The same physics governs other stars: the violent convective envelopes of red giants, the magnetoconvection in white dwarfs.

Ocean thermohaline circulation (the "great ocean conveyor belt") is driven by density differences from both temperature and salinity (double-diffusive convection, a generalization of Rayleigh-Bénard). The deep water formation in the North Atlantic — cold, salty water sinking and driving a global circulation pattern — is the same instability operating at planetary scale with a 1000-year timescale.

Earth’s atmosphere convects daily: the boundary layer heating from the sun-warmed surface drives cumulus clouds (rising thermals — Plume’s amber columns), with clear downdrafts between them. The cumulus convective parameterization in climate models is essentially a Rayleigh-Bénard model at coarse resolution. The tradewind inversion is the stable conductive layer above — the top boundary condition of the Bénard problem imposed by large-scale atmospheric dynamics.

The Simulation Method

Plume solves the 2D stream-function/vorticity equations on a 120 × 30 staggered grid (aspect ratio 4:1, non-dimensional domain Lx = 4, Ly = 1). Horizontal boundaries are periodic (allowing the rolls to form without preferential start positions). Vertical boundaries are no-slip walls with Thom’s boundary condition for vorticity: ω_wall = −2ψ_interior/h² — derived from applying the no-slip condition u = 0 at the wall to the stream function equation via a ghost-cell argument.

// Time stepping: explicit forward Euler, dt = 1.9×10⁻⁴
// Stability: dt < h²/4 = 2.78×10⁻⁴  (von Neumann for Pr=1 diffusion)
// Steps per frame: 50  →  Δt ≈ 0.0095 per animation frame
//
// Poisson solve: SOR (successive over-relaxation, ω = 1.65)
//   22 inner iterations per time step, warm-started from previous step
//   Fast convergence: psi changes by O(dt) each step → small residual
//
// Color map:
//   T = 0 (cold): rgb(16, 168, 196)  — teal
//   T = 0.5:      rgb(8, 10, 22)     — near-black
//   T = 1 (hot):  rgb(240, 140, 20)  — amber
//   Vorticity view: amber = CCW (ω>0), teal = CW (ω<0)
//   Speed view: black → amber by |u|

The Poisson equation ∇²ψ = −ω must be solved at every time step to update the stream function from the new vorticity field. Plume uses warm-started SOR: since ψ only changes by O(dt) between steps, the previous ψ is an excellent initial guess, and 22 SOR iterations reduce the residual to negligible levels. The over-relaxation factor ω_SOR = 1.65 accelerates convergence significantly over pure Gauss-Seidel (ω = 1). The entire simulation — 50 time steps of vorticity transport plus 50 Poisson solves — runs comfortably at 60 fps in the browser using Float32Arrays and minimal JavaScript overhead.