Chaos That Begins at the Bottom of a Pot — Rayleigh–Bénard Convection and the Lorenz Equations

In 1900 the French physicist Henri Bénard heated a thin layer of oil from below and saw something strange. A liquid that had been sitting perfectly still suddenly cracked itself into a honeycomb of hexagons. No one stirred it, yet it began to move in a regular pattern. Sixteen years later Lord Rayleigh explained the threshold with a single number, and forty-seven years after that Edward Lorenz found chaos hiding in the same equations.

This post follows what happens when you heat a pot from below, in three steps. There's the Rayleigh number (the ratio of buoyancy to viscous and diffusive braking), the critical Rayleigh number of 1708 where convection switches on, and the Lorenz chaos that pops out once you shrink that convection down to three equations. We integrate the Lorenz system in Python to watch all three faces of convection directly.

Why a Still Liquid Starts to Move#

Picture a liquid layer that is hot at the bottom and cold at the top. The warm fluid below expands and becomes lighter. Light stuff below, heavy stuff above — an upside-down density arrangement. Buoyancy wants to flip it over.

Yet the fluid does not immediately surge upward. Two things hold buoyancy back. Viscosity resists the flow. Thermal diffusion bleeds the heat of a rising warm blob into its surroundings, killing the buoyancy that lifted it. Both are brakes on convection.

So whether convection switches on comes down to a tug of war. If buoyancy wins, the fluid overturns and a circulation (a convection roll) forms. If the brakes win, the fluid stays put and heat moves by conduction alone.

The Rayleigh Number — Overturning Force over Restraining Force#

That tug of war is compressed into a single dimensionless number: the Rayleigh number.

Here $g$ is gravitational acceleration, $\beta$ is the thermal expansion coefficient (how much volume grows as temperature rises), $\Delta T$ is the top-to-bottom temperature difference, $L$ is the layer thickness, $\nu$ is the kinematic viscosity, and $\alpha$ is the thermal diffusivity. The numerator is the overturning buoyancy; the denominator is the restraining viscosity and diffusion.

The key is $L^3$. Thicken the layer and the Rayleigh number grows as the cube. A thin oil film rarely convects, but the deep ocean or the Earth's mantle convects at a tiny temperature difference. That's why a slightly deeper pot switches convection on far more easily.

Critical Rayleigh Number 1708 — the Threshold Where Convection Turns On#

Rayleigh's linear stability analysis gives a remarkably clean answer. With idealized free top and bottom boundaries, convection begins at $\mathrm{Ra}_c = 27\pi^4/4 \approx 657.5$. With the rigid walls common in real experiments, $\mathrm{Ra}_c \approx 1708$.

Below $\mathrm{Ra}_c$, any small wobble is eaten by viscosity and diffusion and dies out. The fluid returns to rest. The moment you cross $\mathrm{Ra} > \mathrm{Ra}_c$, a wobble of a particular wavelength starts to grow on its own. That is the seed of a convection roll.

Try it yourself in the simulation below. Raise the Rayleigh number and you can watch the threshold where convection turns on.

Keep Ra below 1708 and the particles stay nearly frozen while the temperature forms smooth horizontal layers (conduction). Cross the threshold and counter-rotating rolls appear, with hot columns rising and cold columns sinking. Push Ra higher and the rolls spin faster.

The Shape of a Roll and the Nusselt Number#

Just past the threshold, convection appears as rolls of a particular wavelength. The critical wavenumber is $k_c \approx 3.117/L$ — one roll is about as wide as the layer is deep. The rolls are roughly square. That shape is no accident: the mode that grows most easily (switches on at the lowest Ra) has exactly that wavelength.

How well convection carries heat is measured by the Nusselt number (total heat transfer over pure conductive heat transfer).

Here $q_\text{total}$ is the actual heat flux and $q_\text{conduction}$ is the flux with no convection. With conduction alone $\mathrm{Nu}=1$. As convection strengthens, experiments give roughly $\mathrm{Nu} \sim \mathrm{Ra}^{1/3}$. Convection doesn't just move a little more heat — it shovels it at a scale conduction can never reach.

The Convection Lorenz Reduced to Three Equations#

In 1963 the meteorologist Edward Lorenz simplified this convection to the extreme. Keep only the coarsest few modes of the temperature and flow fields (a Galerkin truncation) and the partial differential equations collapse to just three ordinary differential equations.

Here $x$ is the circulation strength of the roll, $y$ is the temperature contrast between rising and sinking columns, and $z$ is the distortion of the vertical temperature profile. $\sigma$ is the Prandtl number (momentum diffusivity over thermal diffusivity), $\beta$ is a geometric constant, and $\rho$ is exactly $\mathrm{Ra}/\mathrm{Ra}_c$ — the very Rayleigh ratio we have been discussing.

Raising $\rho$ is the same as heating the pot harder. For $\rho<1$ every trajectory heads to the origin (conduction). For $1<\rho<24.74$ it converges to an off-origin fixed point (a steady convection roll). Cross $\rho>24.74$ and the trajectory never settles onto any fixed point; it wanders forever between two wings. That is chaos.

Python — Integrating the Three Faces of Convection#

We integrate the Lorenz system with fourth-order Runge–Kutta to see how convection changes with $\rho$.

import numpy as np
 
SIGMA, BETA = 10.0, 8.0 / 3.0
 
def lorenz_deriv(state, rho):
    x, y, z = state
    return np.array([
        SIGMA * (y - x),        # change in roll circulation strength
        x * (rho - z) - y,      # change in rising/sinking temperature contrast
        x * y - BETA * z,       # change in vertical temperature distortion
    ])
 
def integrate_rk4(rho, dt=0.01, n=20000, s0=(0.6, 0.4, 0.2)):
    s = np.array(s0, dtype=float)
    traj = np.empty((n, 3))
    for i in range(n):
        k1 = lorenz_deriv(s, rho)
        k2 = lorenz_deriv(s + 0.5 * dt * k1, rho)
        k3 = lorenz_deriv(s + 0.5 * dt * k2, rho)
        k4 = lorenz_deriv(s + dt * k3, rho)
        s = s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        traj[i] = s
    return traj
 
def convection_regime(rho):
    tail = integrate_rk4(rho)[-4000:]          # tail after dropping the transient
    spread = tail[:, 0].max() - tail[:, 0].min()
    if abs(tail[-1, 0]) < 1e-3 and rho < 1.0:
        return "conduction"                     # converges to the origin
    if spread < 0.1:
        return "steady roll"                    # converges to a fixed point
    return "chaos"                              # oscillates between two wings
 
for rho in (0.8, 15.0, 28.0):
    print(f"rho={rho:5.1f} -> {convection_regime(rho)}")

The output is:

rho=  0.8 -> conduction
rho= 15.0 -> steady roll
rho= 28.0 -> chaos

The same equations swing between rest, steady convection, and chaos with a single heating knob. Try it yourself in the simulation below.

Keep $\rho$ below 24.74 and the trajectory is pulled into one of the yellow fixed points (a steady convection roll). Cross the threshold and the famous butterfly wings appear. Hit reset with a slightly different starting point and two trajectories that overlap at first soon split completely — the butterfly effect.

What to Remember#

• The Rayleigh number is the ratio of overturning buoyancy to restraining viscosity and diffusion. It scales as $L^3$, so deeper layers convect more easily.
• Cross 1708 and a still liquid organizes itself into convection rolls. Convection is a threshold phenomenon that switches on and off.
• Lorenz's three equations are a miniature of this convection, and turning up the heating $\rho$ collapses steady convection into chaos. The limits of weather prediction began at the bottom of this very pot.