How WebGCM works
Part 4: The primitive equations

The equations at the heart of the model, taken apart one term at a time.

Part 3 built the grid and sliced the air above it into layers. The rules that move the air from cell to cell and layer to layer are the primitive equations, the heart of every weather and climate model, and this part takes them apart, one rule at a time. Each rule is an equation, and each equation is a sum of a few terms, each with a job of its own. For every rule we'll write the equation with its terms in colors, play with each term on its own until it makes sense, then put the terms back together in one figure, and finally look at what the model actually computes in their place. At the end the four rules come together again inside a storm. Drag any of the 3D figures to turn them.

What the model stores

Before the rules, the state they act on. The model keeps very little. In each layer it stores one potential temperature θ at the center of every cell and one wind speed across every edge, the across-edge wind of Part 3. Under each column it stores one surface pressure, π (on this page π is always a pressure, never 3.14159). That is all. Temperature, the pressure at every level, the height of every layer, the vertical motion, the wind along each edge: the model works out every one of them again from these few numbers at every step.

The patch here has six layers; the climate model has 36. At N = 128 that is 5.9 million values of θ, 17.7 million winds and 164 thousand surface pressures, about 24 million numbers for the dry atmosphere alone, with water vapor, cloud and the ocean on top. Peel the layers off from the top and watch θ grow with height and fall toward the pole. That is the stable atmosphere of Part 1, colder toward the pole as the real one is, stored one number at a time.

That list of numbers is the heart of the difference between the equations and the model. The equations are written for smooth fields: π, θ and the wind have a value at every point of the sphere and at every instant, and a derivative like ∂/∂t is a limit, the change over an interval divided by the interval's length, as the interval shrinks to nothing. A computer can hold none of that. It holds a finite list of numbers, so every derivative in space has to become a difference between neighboring numbers and every integral a sum. That step is called discretization, and there is more than one way to take it; the choices are where a model's character comes from. For each rule, once its terms make sense, we'll look at the model's choice.

Air in, air out

The first rule is the continuity equation: air is never created or destroyed, only moved. It is written for one thin layer of air, between two levels of σ, the pressure divided by the pressure at the ground, the sigma coordinate of Part 3:

\class{term-a}{\frac{\partial \pi}{\partial t}} + \class{term-b}{\nabla\cdot(\pi\mathbf{u})} + \class{term-c}{\frac{\partial (\pi\dot{\sigma})}{\partial \sigma}} = 0

Here π is the surface pressure, u the horizontal wind, and σ̇, said “sigma dot,” how fast air moves across the levels of σ, positive downward, since σ grows toward the ground. Every ∇ is taken along the layer, at constant σ. There are three terms: the change in the layer's air, the air it sheds sideways, and the air that leaves through its top and bottom. Each can be positive or negative, and since the three add up to zero, their signs trade off: where the air leaving sideways or through the top and bottom is positive, the layer is losing air, and its change is negative. The model works them out in a different order from the one they are written in: first the walls, since the winds through them are what it stores, then their total down the column, which is the change in surface pressure, and only then the flow between the layers. The sections below take them in that order, with a figure each.

Through the walls\nabla\cdot(\pi\mathbf{u})

The model starts with the divergence, the spread of Part 3: what flows out through the walls of a cell, minus what flows in. The figure is one cell of the climate model, seen from above, 60 km across, with the air crossing each of its six walls drawn in blue. Drag any wall's arrow to change the wind through it.

A uniform wind of 10 m/s carries 80 percent of the cell's air in through one side every hour and the same 80 percent out through the other, and changes nothing. Converging at only 1 m/s through every wall moves less than a third as much air, yet the layer gains 24 percent of its air an hour. In “a front” the wind slows from 10 to 9 m/s across the cell: 80 percent comes in and 72 percent goes out, so the layer gains 8 percent an hour. And in “spinning” the air races around the cell along its walls and none of it crosses them at all. Only the difference between in and out matters, and it is usually a small difference of big flows. Each meter per second through one wall moves 4 percent of this cell's air an hour.

The sign is Part 3's, positive where the layer drains and negative where it fills, and the readout shows it: “diverging” gives a positive divergence of 6.7 × 10⁻⁵ per second, and the layer loses 24 percent of its air an hour; “converging” gives −6.7 × 10⁻⁵, and the layer gains 24 percent.

The amount of air\frac{\partial \pi}{\partial t}

Add up what every layer gains or loses through its walls and you have what the whole column gains or loses, and that shows up at the ground, as surface pressure. Why there? Because the pressure at the ground is the weight of everything above it. A layer between two levels of σ holds π Δσ/g of air over each square meter, where g is the strength of gravity, so every layer is a fixed share of the column, and the column's air is measured by π itself.

At 1000 hPa the whole column holds 10.2 metric tons of air over each square meter. Press “Pile air in” and π climbs to 1010 hPa: 102 kg more over each square meter, shared out among the layers in their fixed proportions, the top layer still holding 20 percent of it. Every hectopascal of surface pressure is 10.2 kg of air over each square meter, so ∂π/∂t of 1 hPa an hour piles about 10 kg onto each square meter every hour. Drag the column from a deep low of 960 hPa to a strong high of 1040 hPa: the whole range of the weather's surface pressure is only about 4 percent of the column's air either way.

The sign of this term is the plainest of the three. Positive, the column has taken in more air through its walls than it let out, and its surface pressure is rising, as under a building high; negative, it has let out more than it took in, and the pressure is falling, as under a deepening low.

Between the layers\frac{\partial (\pi\dot{\sigma})}{\partial \sigma}

The last term, and the last one the model works out, is the difference between the air crossing a layer's top and the air crossing its bottom. That flow is πσ̇: multiplying σ̇ by π turns a speed across the levels into an amount of air. The figure is a column of six layers seen from the side, with the flow across each level drawn in green. Drag the arrows.

In the default, 1 hPa an hour of air rises through the middle of the column: the top two layers gain, the lowest two lose, and the two in the middle pass the air through without changing at all. Whatever you set, the whole column's total stays exactly zero. Each arrow moves air from one layer into its neighbor, and no air can cross the top or the ground: σ̇ = 0 at σ = 0 and at σ = 1. Rules like these, which say what happens at the edges of the region an equation covers, are called boundary conditions, and these two are why this term rearranges the column's air but can never change its weight. Only the walls can do that. The same 1 hPa an hour is a different speed at different heights: about 10 m an hour near 830 hPa, 15 m at 500 hPa and 38 m near 170 hPa, because the air up there is thinner.

The sign works the same way as the divergence's, in the vertical. Since σ grows downward, πσ̇ is positive where air sinks, and the term is the flow out through the layer's bottom minus the flow in through its top. Positive, more air leaves the layer through its top and bottom than arrives, and the layer drains; negative, the layer fills. The figure shows each layer's gain, which is the term with its sign flipped: in the default, the term is positive in the lowest two layers, which the rising air leaves, and negative in the top two, where it arrives.

All three together

In the real equation the three terms are tied together, because each layer must keep its fixed share of the column. In this figure each slider sets how much air flows out of, or into, one layer through its walls, and the model works out the rest: the flows between layers in green and, in orange, the change in surface pressure, which lights up the ground and raises or lowers an arrow beside the column whenever the column gains or loses air.

On the page the three terms look like equals, any one of them fixed by the other two. In the model they are not. The winds and the surface pressure are what it stores, so the divergence is the only term it can work out directly; the change in surface pressure and the flow between the layers are unknowns it has to solve for, and the equation alone does not pin them down. In this column of six layers there are six copies of the equation, one per layer, and eight unknowns: the change in surface pressure and the flow across each of the seven levels. The two boundary conditions, no flow through the top and none through the ground, remove two of them, and the equations then fix the remaining six exactly. The model works them out in the order of the sections above, each term built from the one before:

  1. The divergence of every layer, from the winds on its walls. These are what the sliders set.
  2. Their total down the column. Add the six equations together and the flows between layers cancel in pairs, each one leaving a layer as it enters the next, while the two at the ends are zero by the boundary conditions. What survives says that the change in surface pressure is minus the column's average divergence, which is also what the weight of the column demands: whatever piles in anywhere adds to it.
  3. The flow between the layers. With the first two terms known, each layer's copy of the equation has one unknown left in it. Each layer has to keep its fixed share of the column, a layer between σ = 0.8 and 0.9 always holding a tenth of its air, so whatever its own divergence differs from the column's average has to leave or arrive through its top and bottom. The flow across each level is the running total of those differences.

That running total can be added up from either end. Start from zero at the top and add downward, or from zero at the ground and add upward, and the flows come out the same, because the differences from the column's average add up to exactly zero over the whole column. What cannot be skipped is the second step. Leave out the surface pressure and add from the top, and the sum would reach the ground with the column's total still left over, air leaking through the ground; add from the ground, and it would leak out of the top. The change in surface pressure is what lets both boundary conditions hold at once. Nothing here steps the vertical motion forward in time, either: like every hydrostatic model, WebGCM never stores it, and works it out from the divergence afresh each time. The continuity equation itself is done once the surface pressure is known; the flows between the layers are worked out for the other rules. Both the wind and the heat equations have a term that carries their quantity up and down the column, and πσ̇ is what does the carrying.

The rising column, with air flowing in near the ground and out near the top, leaves the surface pressure steady and lifts air up through the middle at about 0.6 hPa an hour: the warmed columns of Part 1, and the inside of every storm. Air piling in at every level raises the surface pressure by a hectopascal an hour and moves no air between layers at all. The whole column gets heavier, every layer keeps its share, and the layers simply thicken together. “Air draining out” is the mirror image: the same hectopascal an hour lost, the arrow pointing down, and still no air crossing between the layers.

Can one layer take in air from every side at once, through its walls and across both its top and its bottom? Its share of the column is fixed, but its mass is that share times the column's, so it can, whenever the whole column is gaining. “Squeezed from every side” has air converging in every layer, six times harder in the other five than in the third, so the third falls short of the column's average and its neighbors push their surplus into it: the surface pressure rises by 2.58 hPa an hour, 0.14 hPa an hour sinks into the third layer from above and 0.21 rises into it from below, on top of the 0.08 its own walls bring in. Its gain is 0.43 hPa an hour, a sixth of the column's, exactly its share. What the fixed share forbids is a layer gaining while the column does not: with the surface pressure steady, air can arrive through a layer's top or through its bottom, never both. “Emptied from every side” is the mirror image, the third layer shedding air through its walls and to both of its neighbors while the whole column drains.

Press “Rising and sinking, out of step” to watch the column swing slowly from rising to sinking and back, each layer a little behind the one above it, so that the change runs down the column like a wave and the lowest layer is a quarter of a cycle behind the top. Whenever the air spreading out aloft and the air converging below are out of step, the difference shows at the ground: the surface pressure falls while the outflow aloft runs ahead, by up to 0.7 hPa an hour, and rises while the inflow below catches up. A surface low deepens in just this way, when the air spreading out aloft takes away more than the air converging below brings in.

In symbols, since the boundary conditions make σ̇ = 0 at the top, σ = 0, and at the ground, σ = 1, integrating the equation down the whole column removes the third term and leaves the change in surface pressure, minus the column's average divergence, D̄. Put that back into the equation and each layer's vertical term is minus its departure from the average, and integrating that from either end gives the same flow across each level:

\begin{aligned} \overline{D} &= \int_0^1 \nabla\cdot(\pi\mathbf{u})\,d\sigma, \qquad \class{term-a}{\frac{\partial \pi}{\partial t}} = -\overline{D} \\ \class{term-c}{\frac{\partial (\pi\dot{\sigma})}{\partial \sigma}} &= -\left(\class{term-b}{\nabla\cdot(\pi\mathbf{u})} - \overline{D}\right) \\ \pi\dot{\sigma} &= -\int_0^{\sigma} \left(\nabla\cdot(\pi\mathbf{u}) - \overline{D}\right)d\sigma' = \int_{\sigma}^{1} \left(\nabla\cdot(\pi\mathbf{u}) - \overline{D}\right)d\sigma' \end{aligned}

The two integrals in the last line are equal because the integral of the departures over the whole column, from 0 to 1, is zero.

What the model computes

Now the discrete version. It needs a little notation, which the other rules use too:

The model replaces the divergence with a sum over the cell's walls: the air crossing each wall, times the wall's length, divided by the cell's area, exactly what the one-cell figure in “Through the walls” adds up. That is the divergence theorem applied to a single cell: the cell's average divergence is exactly the net flow through its walls. The approximation is in the flow through each wall, which the model takes from the one wind it stores there and the average of the two cells' π. The integrals over σ become sums over layers:

\begin{aligned} F_{e,k} &= \pi_e\,u_{e,k}, \qquad \pi_e = \tfrac{1}{2}(\pi_i + \pi_j) \\ \class{term-b}{D_{i,k}} &= \frac{1}{A_i}\sum_{e\in i} n_{e,i}\,l_e\,F_{e,k} \\ \class{term-a}{\frac{d\pi_i}{dt}} &= -\sum_{k} D_{i,k}\,\Delta\sigma_k \\ \class{term-c}{(\pi\dot{\sigma})_{i,k+1/2}} &= -\sum_{m\le k} D_{i,m}\,\Delta\sigma_m - \sigma_{k+1/2}\,\frac{d\pi_i}{dt} \end{aligned}

F_{e,k} is the flux of air across the edge, with π on the edge the average of its two cells, and D_{i,k} is the divergence that the sliders in “All three together” set. The equations run in the same order as the three steps: divergences, then their column total, then the flows between layers, which the model adds up from the top and, at the ground, sets to exactly zero to clear the last digit of rounding. Two things have changed on the way. The derivatives in space are gone: the equations are now a finite list of numbers, each with its own rate of change. And ∂/∂t has become d/dt, an ordinary derivative, because each number now depends on time alone. Time itself is still continuous; the section on stepping takes care of it. Every flux leaves one cell and enters its neighbor with the opposite sign, so the sum of A_i\,d\pi_i/dt over the whole globe is zero to the last digit the computer keeps, which is why the mass of the atmosphere never drifts.

Warm air takes up room

The second rule is hydrostatic balance, from Part 1: the pressure drops with height just fast enough to hold up the weight of the air, ∂p/∂z = −ρg, where ρ is the air's density. The model never uses height directly, though. It writes the balance for the geopotential Φ, the height of a level times g, in terms of the Exner function Π = (p/p₀)κ, which turns potential temperature into temperature, T = θΠ.

Here p₀ is a fixed reference pressure, about 1000 hPa and the same everywhere, κ = R/cp, and cp and R are the specific heat capacity and gas constant of air, per kilogram. With the ideal gas law p = ρRT the balance takes a remarkably compact form:

d\Phi = g\,dz = -\frac{dp}{\rho} = -R\,T\,d\ln p = -\frac{R\,T}{\kappa}\,d\ln\Pi = -c_p\,\theta\,d\Pi
\class{term-a}{\frac{\partial \Phi}{\partial \Pi}} = \class{term-b}{-c_p\,\theta}

It has only two parts: how fast the height climbs as Π falls, and the potential temperature that sets it.

Every layer a straight line\class{term-a}{\frac{\partial \Phi}{\partial \Pi}} = \class{term-b}{-c_p\,\theta}

Plot the height of the air against Π instead of against pressure, and the equation says something you can see. The model's θ is one number per layer, so in this plot every layer is a straight line, and its steepness is its θ. Drag a layer to warm or cool it.

The bottom layer, at 290 K, climbs 2.96 km for every 0.1 that Π falls; the top one, at 379 K above the tropopause, climbs 3.87 km. Roughly, the steepness in kilometers per 0.1 of Π is θ divided by 100 K. Warm the layer from 500 to 300 hPa by 20 K and it grows from 3,590 m to 3,818 m thick, and everything above it rises by the same 228 m, while nothing below it moves. Warm the lowest layer by the same 20 K instead and the surfaces above rise only 93 m, because that layer spans less than half as much Π: across it the pressure falls by only 15 percent, against 40 percent from 500 to 300 hPa, so it starts only 1,347 m thick, and the same warming stretches every layer by about the same share, 6 to 7 percent.

Warming a patch

Now the same rule over a whole region. Warm a patch of air and every surface of equal pressure above it lifts.

Six degrees of warming through the whole column lifts the 850 hPa surface by 29 m and the 200 hPa surface by 283 m: each layer adds its own extra thickness to everything above it, so the bulge grows with height. The ground does not move, because no air has gone anywhere yet; this is the instant before the wind responds. Warm only the air below 700 hPa and every surface above the warm layer rises by the same 63 m. The air up there is no warmer, but it has been lifted. The slopes of those raised surfaces are what push the air outward aloft, the first step of the circulation in Part 1, and they are the next rule's business.

What the model computes

As in the figure, each layer is a straight line in Π, and the integral from the ground becomes a sum of each layer's contribution, cpθ times the layer's span in Π:

\begin{aligned} \Pi_{i,k+1/2} &= \left(\frac{\pi_i\,\sigma_{k+1/2}}{p_0}\right)^{\kappa}, \qquad \overline{\Pi}_{i,k} = \frac{\Pi_{i,k+1/2}\,\sigma_{k+1/2} - \Pi_{i,k-1/2}\,\sigma_{k-1/2}}{(1 + \kappa)\,\Delta\sigma_k} \\ \Phi_{i,k} &= \Phi_{s,i} + \sum_{m>k} \class{term-b}{c_p\,\theta_{i,m}}\left(\Pi_{i,m+1/2} - \Pi_{i,m-1/2}\right) + \class{term-b}{c_p\,\theta_{i,k}}\left(\Pi_{i,k+1/2} - \overline{\Pi}_{i,k}\right) \end{aligned}

\Phi_{s,i} is the ground's own geopotential, g times its height, and \overline{\Pi}_{i,k} is the exact average of Π through layer k, the level where the model places the layer's geopotential and its temperature, T_{i,k} = \theta_{i,k}\overline{\Pi}_{i,k}. The sum adds the full thickness of every layer below and the part of the layer itself that lies below the level \overline{\Pi}_{i,k}: about half of it in most layers, but only a fifth in the top one, which reaches up to σ = 0.

What pushes the wind

The third rule is the momentum equation, the wind's own budget. It has the most going on, four terms in all:

\frac{\partial \mathbf{u}}{\partial t} = \class{term-a}{-(\zeta + f)\,\hat{\mathbf{z}}\times\mathbf{u}} \class{term-b}{- \nabla(\Phi + K)} \class{term-c}{- R\,T\,\nabla\ln\pi} \class{term-d}{- \dot{\sigma}\,\frac{\partial \mathbf{u}}{\partial \sigma}}

They are a turning, by the spin of the planet, f, and of the flow, ζ, the slope of the geopotential plus the kinetic energy per unit mass, K = ½|u|², a second pressure term that comes from the layers following the ground, and the wind carried up and down between layers. Here \hat{\mathbf{z}} points straight up, so \hat{\mathbf{z}}\times\mathbf{u} is the wind turned a quarter turn to the left. We'll take them in the order that makes each easiest to see: first the push of the pressure, then the turning, then the wind carrying itself, and last the vertical.

Pushed along a slope-\nabla\Phi

The pressure gradient force, the push toward low pressure of Part 1, is simplest on a surface of equal pressure. There it is, as Part 1's second view showed, the downhill slope of the surface's height times g, −∇Φ, with ∇ the slope of Part 3: air on a tilted surface of equal pressure is pushed downhill along it, toward where the surface is lower. The figure is a cross-section 2,000 km wide between a warm column and a cold one.

Warm air takes up room, so the surfaces of equal pressure stand higher over the warm side, and the higher up, the more they tilt. With 10 °C between the two sides the 300 hPa surface stands 352 m lower on the cold side, a slope of only 1 in 5,700, yet g times that slope speeds the air up by 6.2 m/s every hour. Step down through the surfaces and the push weakens, 3.6 m/s an hour at 500 hPa and 0.8 at 850, until at the ground, where both columns share the same pressure, there is none: warm air aloft spreads outward, as it did in Part 1. Let the air go and in six hours the 300 hPa air would be moving at 37 m/s, if nothing else acted. Something does: at 45° the Coriolis force would already have turned it 64° to its right. That is the turning term, which comes after the mountain.

Over a mountain\class{term-b}{-\nabla\Phi}\;\class{term-c}{-RT\,\nabla\ln\pi}

The model's layers are not surfaces of equal pressure, though. They are surfaces of constant σ, and over a mountain they bend up with the ground. Moving the gradient off a surface of constant height, ∇z, and onto the sloping layer, ∇σ, picks up the change in height along it, and since the pressure on the surface is p = σπ, the change in pressure along it, σ∇π:

-\frac{1}{\rho}\,\nabla_{z}\,p = \class{term-b}{-\nabla_{\sigma}\Phi} \class{term-c}{- R\,T\,\nabla_{\sigma}\ln\pi}

Over flat ground both are ordinary pushes of a few meters per second an hour, and they share the work between them; near the ground the second is often the larger, as the figure of all four terms shows below. Over a mountain both are huge, and they nearly cancel.

This figure evaluates the model's own versions of the two terms on the climate model's 36 layers, with columns 100 km apart, for air at rest over a mountain, with the temperature falling with height as in the standard atmosphere. The right answer is zero: the dotted surfaces of equal pressure are flat, so nothing should push. Each term on its own, low down on the flank of a 3 km mountain, would speed the air up by about 350 m/s in an hour. What is left over is 0.02 m/s an hour, enough to make a breeze of 0.4 m/s in a day where none should be.

The continuous terms cancel exactly; their discrete versions are each a difference standing in for a derivative, and their errors do not cancel. Move up through the layers and the leftover shrinks, almost to nothing near 350 hPa, then jumps at the tropopause, where the temperature profile bends sharply, to as much as 25 m/s a day. This is the price Part 3 warned of, the classic weakness of layers that follow the ground. It is why most climate models switch to layers of constant pressure higher up, and it is one reason WebGCM smooths its mountains a little before it uses them.

Turned by the spin-(\zeta + f)\,\hat{\mathbf{z}}\times\mathbf{u}

The first term in the equation turns the wind. Its strength is the absolute vorticity, ζ + f: f is the Coriolis parameter of Part 2, 2Ω sin φ, twice the rate at which the ground turns beneath the air; ζ is the spin of the flow itself, its vorticity, the curl of Part 3. The push is always at right angles to the wind, to its right when ζ + f is positive. Suppose a parcel felt this term and nothing else.

Let it go and it traces the inertial circle of Part 2: at 45° north a 10 m/s wind goes once around a circle 97 km in radius in 16.9 hours, at a speed that never changes, since a push at right angles to the motion can only turn it. This term does no work. Change the speed and the circle grows or shrinks with it, but the time around stays 16.9 hours, set by ζ + f alone. Move the parcel toward the pole and the turn tightens to 12 hours, as in Part 2; at the equator f is zero and the parcel goes straight. Now give the surrounding flow a spin: ζ simply adds to f, so 4.3 × 10⁻⁵ per second of counterclockwise spin at 45° turns the parcel as fast as the pole does, and spin of the opposite sign can cancel f altogether. The ζ half never acts alone in the real equation; it is half of the wind carrying itself, which is next.

The wind carries itself\class{term-a}{-\zeta\,\hat{\mathbf{z}}\times\mathbf{u}}\;\class{term-b}{-\nabla K}

The wind carries its own motion along with it, so the wind at a point changes as faster or slower air arrives from upwind. That is advection, −(u·∇)u: the wind, times how fast the wind changes along its own direction. The model never computes it in that form. An identity of vector calculus splits it into the spin of the flow times the wind turned sideways, the ζ half of the turning term, and the slope of the flow's kinetic energy per unit mass, the K in the second term:

-(\mathbf{u}\cdot\nabla)\,\mathbf{u} = \class{term-a}{-\zeta\,\hat{\mathbf{z}}\times\mathbf{u}} \class{term-b}{-\nabla K}, \qquad \zeta = \hat{\mathbf{z}}\cdot\nabla\times\mathbf{u}

In a flow along straight lines, with the same wind all the way across it, there is no spin, and the kinetic-energy part is the whole of the advection. The figure is a straight channel 1,000 km long, with a wind you can shape by dragging its points.

Where the wind speeds up from 10 to 30 m/s along the channel, K climbs downstream and the push down its slope points back upstream: at the middle, where the wind is 20 m/s, it slows the wind by 1.4 m/s every hour. Run it and watch why. The ring follows one bit of air, which keeps its 20 m/s all the way down the channel; it is the air arriving at the middle from upstream that is slower, so the wind at that fixed point falls, to 18.7 m/s after an hour and 14 m/s after six. Make it slow down instead and the wind at the middle rises; make it steady and nothing at any point changes, though the air is moving just as fast.

When the flow curves, or the wind changes across it, the spin part joins in. Stir the flow and look at both pieces.

In the vortex, the whole push points outward: a spinning patch of air flings itself apart unless a pressure gradient holds it in, which is the gradient wind of Part 2. Its two parts disagree. Near the center, where the air turns like a wheel, the vortex part points outward and is stronger than the whole, and the kinetic-energy part points inward, taking the difference back. Out past the fastest winds the vortex's spin reverses, and the two parts swap roles. In the straight jet the whole push is exactly zero, since air arriving from upwind is moving just as fast, yet the two parts are each large and cancel each other. On this coarse demonstration grid the two parts add up to the whole only approximately, because the figure computes each with simple differences; on paper they agree exactly.

The split is worth the trouble for two reasons. The vortex part has the same shape as the Coriolis force, so the model adds the two together and handles them as one term, the absolute vorticity ζ + f times the sideways wind. That is the term TRiSK was built for in Part 3, the one that needs the reconstructed along-edge wind, and its weights make it neither create nor destroy energy. And the kinetic-energy part is a slope, just like the pressure gradient force, so the model adds it to Φ and takes one slope of the sum, on the edges, at almost no cost.

Up and down-\dot{\sigma}\,\frac{\partial \mathbf{u}}{\partial \sigma}

The last term is the same idea in the vertical: air rising or sinking carries its wind with it into the layers above or below. In height it reads −w ∂u/∂z, the vertical speed times how fast the wind changes with height. The figure is a column with a west wind that strengthens with height, as it does under the jet stream.

With the wind growing 3 m/s for every kilometer and the air rising 5 cm/s, the term slows the wind by 0.54 m/s an hour at every height up to 10 km, where the shear stops: the rising air brings slower wind up from below. Run it for six hours and the wind at 5 km drops from 20 to 16.8 m/s. Sinking air does the opposite, bringing the fast wind of the upper air down. Reverse the shear and rising air speeds the wind aloft instead; make the wind the same at every height and nothing happens, however fast the air rises. The term needs both: vertical motion and a wind that changes with height. It only moves wind between levels and never creates any.

All four together

Now all four terms at once, computed by the model's own code from a standard test set out by Jablonowski and Williamson. It starts from a balanced jet in each hemisphere, over air that is warm at the equator and cold at the poles, and a small bump in the northern jet grows into a chain of storms called a baroclinic wave, which Part 5 runs in full. This is its ninth day, as the lows deepen. Pick a layer and a term, and tap anywhere to read all four.

The turning and the two pressure terms nearly cancel. Near 850 hPa, at the ring, the turning is 11.7 m/s an hour and points south-southeast, a quarter turn to the right of the wind, which blows from the southwest; the slope and the pressure term, 4.9 and 6.0, point north and northwest, toward lower pressure, and together almost exactly against it. What is left, all four together, is only 1.1 m/s an hour, drawn five times longer than the rest so that it shows. At the typical cell of the map the pressure terms and the turning cancel to within 9 percent, and higher up the balance tightens: to 3 percent at 500 hPa and 1.4 percent at 250 hPa. This near-cancellation is the geostrophic balance of Part 2, the wind blowing along the isobars, and here it emerges from the model's own terms rather than being assumed. The small leftover is the weather changing. The vertical term is smaller still, at most 0.3 m/s an hour. Switch between “b alone” and “c alone” to see the pressure push split between the two pieces the σ coordinates make of it.

What the model computes

The model keeps only the wind across each edge, so it needs each term's component along \mathbf{n}_e. These formulas work within one layer, so they leave off its index k: u_e is u_{e,k}, F_{e'} is F_{e',k}, and K_i, \zeta_v and Q_v are the layer's own. The slope and the pressure term are each a difference across the edge, divided by the distance between the two cell centers: a derivative along the edge's normal, estimated from the two nearest values the grid has. The kinetic energy of each cell is built from the winds on its walls:

K_i = \frac{1}{A_i}\sum_{e\in i} \frac{l_e\,d_e}{4}\,u_e^2, \qquad \class{term-b}{-\frac{(\Phi + K)_j - (\Phi + K)_i}{d_e}} \class{term-c}{- R\,\frac{\theta_i\overline{\Pi}_i + \theta_j\overline{\Pi}_j}{2}\,\frac{\ln\pi_j - \ln\pi_i}{d_e}}

The vorticity lives at the vertices, as the circulation around the triangle joining three cell centers divided by the triangle's area A_v, with t_{e,v} = +1 where the edge's normal \mathbf{n}_e runs counterclockwise around the triangle, seen from above, and −1 where it runs clockwise. Added to the Coriolis parameter at the vertex, f_v, and divided by \pi_v, the average of the three surrounding cells' π weighted by how much of the triangle lies in each, it becomes the absolute vorticity per unit of air, the potential vorticity Q_v. The component across the edge of the turned wind −\hat{\mathbf{z}}\times\mathbf{u} is just the wind along the edge, so the turning is TRiSK's sum over the other edges of the edge's two cells, ten of them between hexagons, with the weights w_{e,e'} of Part 3 and the mass fluxes F_{e'} of the continuity rule:

\zeta_v = \frac{1}{A_v}\sum_{e\in v} t_{e,v}\,d_e\,u_e, \qquad Q_v = \frac{\zeta_v + f_v}{\pi_v}, \qquad \class{term-a}{\frac{1}{d_e}\sum_{e'} w_{e,e'}\,l_{e'}\,F_{e'}\,\frac{\hat{Q}_e + \hat{Q}_{e'}}{2}}

\hat{Q}_e is the potential vorticity averaged over the edge's two vertices. Averaging it inside the sum, pair by pair, is what makes the term do no work, as the circling parcel in “Turned by the spin” showed it must: the weights are antisymmetric, w_{e,e'} = -w_{e',e}, so in the budget of the total kinetic energy every pair of edges cancels. The vertical term is a centered difference between the layers above and below, with πσ̇ on the edge, (\pi\dot{\sigma})_{e,k\pm 1/2}, the average of its two cells like \pi_e, arranged so that air moving between layers carries its kinetic energy without creating or destroying any:

\class{term-d}{-\frac{(\pi\dot{\sigma})_{e,k+1/2}\,(u_{e,k+1} - u_{e,k}) + (\pi\dot{\sigma})_{e,k-1/2}\,(u_{e,k} - u_{e,k-1})}{2\,\pi_e\,\Delta\sigma_k}}

Heat goes where the air goes

The fourth rule is the thermodynamic equation, and in this model it is short: θ is carried along with the air, and changes only where heat is added. It has the same budget as the air itself, with θ riding along in every flux, and a source H for the heating, measured as a rate of change of θ:

\class{term-a}{\frac{\partial (\pi\theta)}{\partial t}} + \class{term-b}{\nabla\cdot(\pi\theta\mathbf{u})} + \class{term-c}{\frac{\partial (\pi\dot{\sigma}\,\theta)}{\partial \sigma}} = \class{term-e}{\pi H}

The terms mirror the continuity equation's: the change in the layer's θ, weighted by its air, the θ it sheds sideways and the θ carried through its top and bottom, balanced by the heating. Subtracting θ times the continuity equation and dividing by π leaves an equation for θ itself, which says that following the air, θ changes only by heating:

\frac{\partial \theta}{\partial t} + \class{term-b}{\mathbf{u}\cdot\nabla\theta} + \class{term-c}{\dot{\sigma}\,\frac{\partial \theta}{\partial \sigma}} = \class{term-e}{H}

Carried past\mathbf{u}\cdot\nabla\theta

At a fixed point θ changes as air of a different θ is carried past it, at the rate −u·∇θ. The map is 1,500 km across, with a front where θ climbs 10 K across 300 km, and a wind you set with the dial.

With 15 m/s from the southwest, the probe, a fixed spot on the map, warms by 1.3 K an hour, as warmer air arrives from upwind. Turn the wind to blow straight across the front, from the south, and it warms at 1.8 K an hour; blow along the front, from the west, and nothing changes, though the wind is just as strong. Press “Play”: the probe climbs from 286 to 290 K in about three hours, while the parcel that started there, moving with the wind, keeps its 286 K the whole time. This term is only the wind carrying θ; following the air, nothing changes.

Rising, sinking and heating\class{term-c}{\dot{\sigma}\,\frac{\partial \theta}{\partial \sigma}}\;\class{term-e}{H}

The same in the vertical. In a stable atmosphere θ grows with height, so rising air brings lower θ up to each level and cools it, and sinking air warms it. Heating, H, can push the other way.

Air rising 2 cm/s at 6 km, through θ that climbs 3.3 K per kilometer, cools the air there by 5.7 K a day: the vertical term. Run it for six hours and the lines of equal θ ride up with the air. Now heat it by 5.7 K a day at 6 km and the two balance at every height: the air keeps rising through the lines of equal θ, warming as it goes, while the lines themselves stay put. That is the balance of the tropics, where water condensing in the rising air heats it. Sinking air with cooling makes the mirror image, the balance of the deserts and the subtropical oceans, where the air sinks and cools by radiating to space.

Riding the isentropes

Together the terms say that with no heating every parcel keeps its θ, which means it stays on its own surface of equal potential temperature, its isentropic surface, or isentrope, wherever it goes. In an atmosphere warmer at the equator than at the poles, as the Earth's is near the ground in every season, those surfaces climb toward the pole.

The parcel starts 700 m up at 34° north, on the 300 K surface, and drifts poleward with the jet carrying it east. Because it keeps its θ, it has to climb as it goes: by 63° north, nearly five days later, the 300 K surface is near 490 hPa, about 5.6 km up, and the parcel rides it all the way. Nothing pushed it up. It rose because the surface it lives on rises. Heat it and it climbs across the surfaces instead, faster than the slope alone would carry it, the way air does when its water vapor condenses. Air moving poleward and up, or equatorward and down, along these slanting surfaces is how heat travels poleward in the middle latitudes, and it is what the storm figure after the next section shows.

What the model computes

The model computes the tendency of θ in the flux form, with the same fluxes of air as the continuity equation, so every evaluation keeps the global total of πθ to the computer's last digit. θ on an edge is the average of its two cells, and θ at an interface between layers, \theta_{i,k\pm 1/2}, is interpolated linearly in Π between the layers above and below. In the runs of this part H is zero; in the forced runs of Part 5 it is a nudging toward a target temperature.

\frac{d\theta_{i,k}}{dt} = \class{term-b}{-\frac{1}{\pi_i}\left[\frac{1}{A_i}\sum_{e\in i} n_{e,i}\,l_e\,F_{e,k}\,\theta_{e,k} - \theta_{i,k}\,D_{i,k}\right]} \class{term-c}{- \frac{(\pi\dot{\sigma})_{i,k+1/2}\,(\theta_{i,k+1/2} - \theta_{i,k}) - (\pi\dot{\sigma})_{i,k-1/2}\,(\theta_{i,k-1/2} - \theta_{i,k})}{\pi_i\,\Delta\sigma_k}}

Stepping forward in time

Every equation has now become a list: the rate of change of each of the model's numbers, computed from the current values of its neighbors. Collect all of them, every \pi_i, \theta_{i,k} and u_{e,k}, into one long vector y, and all the right-hand sides into one function f, and the whole atmosphere is a single system of ordinary differential equations:

\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y})

Space has been discretized; time has not. The model still has to get from now to a moment Δt later, and the obvious way is to replace the derivative with a difference, assuming the rate of change stays what it is now for the whole step. That is the forward Euler method:

\mathbf{y}^{n+1} = \mathbf{y}^{n} + \Delta t\,\mathbf{f}(\mathbf{y}^{n})

Forward Euler has a lot going for it. It needs one evaluation of f per step, as cheap as a step can be. It is explicit: the new state comes straight from the old one, with no equations to solve. And it is easy to reason about. Each step's error is of order Δt², so over a fixed stretch of time, which takes 1/Δt steps, the error is of order Δt: halving the step halves the error. That is called first-order accuracy, and it is the method's first weakness, since being ten times more accurate costs ten times as many steps.

The second weakness is fatal. Most of what the atmosphere does is oscillate: gravity waves, the inertial circles of Part 2, the Rossby waves of Part 5, the wind carrying a pattern past a point. The simplest oscillation is dy/dt = iωy, with y a complex number and i here the square root of −1 rather than a cell, whose exact solution goes around a circle at a steady speed, ω radians per unit of time. One forward Euler step multiplies y by 1 + iωΔt, and for every step size

\left|1 + i\,\omega\Delta t\right| = \sqrt{1 + \omega^2\Delta t^2} > 1

Every oscillation grows, every step, however small the step. A smaller step only slows the growth. The figure runs the inertial circle of Part 2, the wind turning at the rate f with nothing else acting on it, so that ω is the Coriolis parameter f, once with forward Euler and once with Runge–Kutta 4, the method the model actually uses, which comes next.

Forward Euler steps along the tangent to the circle, always a little outward, and spirals away. At the figure's starting step, 1/25 of a turn, it speeds the wind up by a factor of 2.1 every turn; at 1/13 of a turn, f Δt = 0.5, by 4.1. At N = 128 the model's step turns the inertial oscillation at 45° by only 0.017 radians, and even so forward Euler would make it grow by 8 percent a day. For the fastest waves the grid can carry it would be hopeless: a factor of 2.8 every step.

The model instead uses the classical fourth-order Runge–Kutta method, RK4. It samples the rate of change four times across the step, once at the start, twice at the middle and once at the end, each sample taken at a trial state built from the one before, and combines them:

\begin{aligned} \mathbf{s}_1 &= \mathbf{f}(\mathbf{y}^{n}), & \mathbf{s}_2 &= \mathbf{f}\left(\mathbf{y}^{n} + \tfrac{\Delta t}{2}\,\mathbf{s}_1\right), \\ \mathbf{s}_3 &= \mathbf{f}\left(\mathbf{y}^{n} + \tfrac{\Delta t}{2}\,\mathbf{s}_2\right), & \mathbf{s}_4 &= \mathbf{f}\left(\mathbf{y}^{n} + \Delta t\,\mathbf{s}_3\right) \end{aligned}
\mathbf{y}^{n+1} = \mathbf{y}^{n} + \frac{\Delta t}{6}\left(\mathbf{s}_1 + 2\,\mathbf{s}_2 + 2\,\mathbf{s}_3 + \mathbf{s}_4\right)

It costs four evaluations of f per step instead of one, and it buys two things. It is fourth-order accurate: halving the step cuts the error sixteenfold. And it handles oscillations. For dy/dt = iωy one RK4 step multiplies y by a factor G, the |G| of the figure's chart, whose size is

\left|G(i\,\omega\Delta t)\right|^2 = 1 - \frac{(\omega\Delta t)^6}{72} + \frac{(\omega\Delta t)^8}{576}

which stays at or below 1 all the way up to ωΔt = 2√2 ≈ 2.83. Well below that limit RK4 damps a wave only very slightly, by about (ωΔt)⁶/144 per step, so that a wave of about twelve steps a cycle keeps 99.99 percent of its amplitude each step; close to the limit the damping is strong. Beyond the limit it grows. Drag the time step past the marked limit in the figure to watch it go.

That limit sets the model's time step. The fastest waves the model's equations allow are Lamb waves: hydrostatic balance removes the sound waves that travel up and down, but keeps pressure waves that run along the ground at about the speed of sound, 310 to 340 m/s, faster in warmer air. The shortest wave the grid can hold, a pattern that flips from cell to cell, oscillates at about 2.7c/d̄, where c is the wave speed and d̄ the mean distance between cell centers, 60 km at N = 128. For a wave at 340 m/s, the speed of sound at 15 °C, the model's step of 169 s turns that pattern 2.6 radians per step, 92 percent of the way to RK4's limit, which is the mark on the figure. This is the form the Courant–Friedrichs–Lewy condition of Part 3 takes for this grid and RK4: the fastest wave may cross only about one cell per step, here cΔt/d̄ ≈ 0.95. It is why the step is 1350 s × 16/N, halving every time the grid is refined, so that each doubling of the resolution costs eight times as much computing: four times the cells and twice the steps. Push past the limit and the model blows up. At N = 16, where the climate model would step every 1350 s, the model's dynamical core, the four rules alone with none of the climate model's physics, runs the baroclinic wave stably with steps up to 1600 s and blows up within three days at 1625 s. That is a little beyond the 340 m/s estimate, because the fastest wave in that test's colder air runs at about 315 m/s.

At the model's step the grid-scale wave keeps only 58 percent of its amplitude each step, and its phase comes out wrong. That is no loss. A wave two cells long cannot be represented faithfully on any grid, and damping it keeps it from feeding back on the weather.

RK4 steps the four rules, the dynamics. Everything else is applied after it, one process at a time: the ocean, which takes a longer step of its own every few steps, heating by radiation and the surface, a smoothing of the shortest waves, a sponge near the model's top that absorbs waves before they reflect off it, mixing in the boundary layer near the ground, condensation and convection, the drag of the ground, of mountains and of breaking gravity waves on the wind, and finally the kinetic energy that the smoothing, mixing and drag removed, returned to the air as heat:

\mathbf{y}^{*} = \mathrm{RK4}_{\Delta t}(\mathbf{y}^{n}), \qquad \mathbf{y}^{n+1} = \mathrm{P}_{\Delta t}(\mathbf{y}^{*})

Here RK4Δt is one RK4 step of the four rules and PΔt applies the other processes in turn. This is called operator splitting. Several of those pieces use forward Euler after all, and that is fine, because they damp, or change slowly, rather than oscillate. For a decay, dy/dt = −λy, forward Euler multiplies y by 1 − λΔt each step and is stable as long as λΔt stays below 2. The smoothing is a hyperdiffusion, a diffusion with ∇⁴ in place of ∇², which damps each wave at a rate proportional to the fourth power of its wavenumber, so it removes the shortest waves while barely touching the weather:

\frac{\partial \mathbf{u}}{\partial t} = -\nu_4\,\nabla^4\mathbf{u}, \qquad \frac{\partial \theta}{\partial t} = -\frac{\nu_4}{\Pi}\,\nabla^4 T

with ν₄ chosen so that the wave two grid spacings long fades by a factor of e in three hours, and at N = 128 its fastest-damped pattern has λΔt of only about 0.12. Splitting makes the step as a whole only first-order accurate in time, but the physics changes slowly compared with the waves, and it is the waves that need RK4's care. The runs behind the next two figures use that bare dynamical core at N = 16, with 27 layers rather than the climate model's 36, a 450 s step and the smoothing included in every RK4 stage instead of after it.

All four in a storm

Here are all four rules working together, inside the baroclinic wave that the momentum budget above was taken from. On the sixth day, as the lows begin to deepen, 192 parcels are released at two heights across the storm's path, and the model's winds carry them for five days. Heights are exaggerated so the climbs show.

Ahead of each low, east of it, the parcels move poleward and rise, roughly along the isentropes; behind it they move toward the equator and sink. That pairing of warm air rising and cold air sinking flattens the sloping isentropes and lowers the air's center of mass, which releases energy, and the next figure follows it. On this coarse grid the motions are gentle: after five days 171 of the 192 parcels are within 150 hPa of where they started, and the highest climb is about 290 hPa. All but two of the parcels that rose more than 50 hPa also moved poleward, and all but one of those that sank as far moved toward the equator. In a real storm the rising stream ahead of the low, called the warm conveyor belt, can carry air from near the ground to the upper troposphere in a day or two. Color the paths by potential temperature to see how well they keep it. The parcels released near 700 hPa do well: half of them end within 0.7 K of the θ they started with. Those released near 900 hPa, down where the fronts form, drift by about 4 K, the rising ones losing θ and the sinking ones gaining it. With cells 450 km across the model mixes warm and cold air together at the fronts, and the winds the parcels follow are blurred as well.

Where the energy goes

The storm's winds come from somewhere. The chart follows the energy of the whole atmosphere over the twelve days of the test, as the change in each kind since the start: the kinetic energy of the wind, and the potential and internal energy of the air, its height in gravity and its heat, which in a hydrostatic atmosphere travel together. Hover over it to read the values.

For the first week almost nothing changes, apart from a small zigzag: the starting state is not perfectly balanced on this coarse grid, and as it settles a few kJ per square meter slosh back and forth between the wind and the air. The two curves mirror each other, so the sloshing cancels in the total, which drifts down only slowly. Then, as the storm grows, kinetic energy rises by 90 kJ per square meter by day twelve, about 5 percent more than the jets started with, while the potential and internal energy fall by 155 kJ. That is the conversion that drives every storm: the sinking of cold air under warm turns stored energy into motion.

The total falls too, by 65 kJ per square meter, and that part is the model's doing rather than the atmosphere's. As the storm grows, its fronts sharpen and its swirls tighten until they reach the scale of the grid, where the model's smoothing, the ν₄ term of the stepping section, removes them. Real turbulence does the same job but turns the motion it destroys into heat, so the total is kept. The climate model does that too, handing the energy its smoothing removes back to the air as heat; this bare run of the core does not, so its total falls. The mass of the atmosphere is a different matter: it changes by about one part in a million billion, the limit of the computer's arithmetic, because every flux of air leaves one cell and enters the next. The total amount of θ changes by a few parts in a million over the twelve days. The fluxes keep it exactly at every evaluation, but the smoothing, which in this run is part of every evaluation, works on temperature rather than on the amount of θ, and each RK4 step updates θ and π separately, so neither keeps their product exactly.

The whole machine

Four rules: air moves from cell to cell and layer to layer without being lost; warm columns stand taller; the wind is pushed by the slope of the pressure, turned by the spin of the planet and the flow, and carries itself; and the heat goes where the air goes. In the computer each becomes a sum over a cell's walls, a difference across an edge or a sum down a column, and the model evaluates them for every cell and edge of every layer four times a step, 512 steps a simulated day at N = 128.

What they leave out is everything that makes weather wet and the climate the one we live in: sunlight and radiation that warm the ground and cool the sky, water vapor and the heat it carries, clouds, an ocean, ice, soil and the real shape of the land. Part 5 first lets the four rules run on their own: the Rossby waves of a single layer, then the storms, then a whole dry climate with one simple rule standing in for sunlight. The parts after it add the rest, one piece at a time.