How WebGCM works
Part 3: The grid

How to cover a sphere with cells that a computer can march across, and how to slice the air above them.

Part 1 ended with a row of columns of air standing side by side, and Part 2 set the planet they stand on turning. The real model needs columns all over the globe, and a rule for which ones are neighbors. That is a grid, and on a sphere a grid is surprisingly hard to get right. This part is about the one WebGCM uses and the handful of ideas that make it work. Every globe can be turned by dragging it.

Squares of latitude and longitude

The obvious grid is the one on every classroom globe: lines of latitude and longitude, with a cell wherever they cross. It has two problems, and the figure shows both. Color the cells by their area and the ones along the equator are many times bigger than the ones beside the poles, because the meridians converge. A cell that is 1100 km wide at the equator is 100 km wide at 85° north and shrinks to a sliver at the pole itself. Slide the cell size and watch the ratio.

The second problem is the time step. A model marches forward in steps, and a step cannot be much longer than the time the fastest wave needs to cross the smallest cell, or the arithmetic feeds on its own errors and blows up. That rule is the Courant–Friedrichs–Lewy condition, and Part 4 meets it properly. The slivers at the poles force the whole globe to take steps ten times shorter than its typical cell would need, so most of the work is wasted. Older models dealt with this by smoothing the poles or bolting on special treatments. The modern answer is a different grid.

An equal-area grid

Start with an icosahedron, the twenty-sided die, inflated onto the sphere. Its twenty triangles are identical, so whatever you do to one you can do to all of them. Divide each edge into N parts, fill the triangle with N² smaller triangles, and carry the result onto the sphere with a map projection that keeps areas equal: the one John Snyder published in 1992 for exactly this purpose, which is the basis of the ISEA grids that geographers use for global datasets. The model does not use the triangles themselves but the cell around each of their corners, which is a hexagon, except at the twelve corners of the original icosahedron, where it is a pentagon. There are always exactly twelve pentagons, whatever N is, and 10N² + 2 cells in all.

The cells are drawn on the same color scale as the latitude–longitude grid, and they barely register on it: the largest is about a fifth bigger than the smallest, with the pentagons smallest of all, against a factor of twenty or more before. There is no pole to speak of, and the time step is set by a typical cell rather than a freak one. This kind of mesh is called a geodesic grid. The model runs at N = 128 on a desktop, which is 163,842 cells with centers about 60 km apart, and at N = 64 on a phone, with cells twice as wide.

Nudging the cells into shape

The raw grid has one flaw that matters later. The method the model uses for the wind, which a later section explains, wants each cell edge to be cut at its midpoint by the line joining the two cell centers on either side of it. On the raw grid that line can miss the midpoint by 16 percent of the edge's length or more, depending on the resolution, and the worst misses are always next to the pentagons. So the model nudges the grid before it starts, in passes. One pass moves every cell's center to the center of mass of its cell and then redraws the cells around the new centers, as the set of points closer to each center than to any other. That changes the centers of mass a little, so the pass is applied again, and again, thirty times in all. This is Lloyd's algorithm, and what it converges toward is a centroidal Voronoi tessellation, a grid in which every center sits exactly at its cell's center of mass. The figure looks straight down on the north pole, one of the twelve pentagons, where the misses are largest. Apply the passes yourself and watch the numbers fall, and slide N to see how they depend on the resolution.

At N = 10 the first pass does most of the work, taking the worst miss from 16 percent to 12. The gains shrink with every pass after that: eight passes bring the worst edge to 9 percent and the typical one to 1.4, thirty bring them to 8.7 and 1.3, and a hundred more would barely move them. Thirty is where the model stops, because in a standard test of steady flow around the sphere the height error keeps falling until about then and only creeps lower afterward. The two numbers behave differently as the grid gets finer. The worst edge always touches a pentagon, and on the raw grid it gets worse with resolution, reaching a 21 percent miss at the model's N = 128, while the relaxation holds it near 9 or 10 percent at every N. The typical edge improves on its own, because the distortion is confined to a few rings of cells around each of the twelve pentagons, and those rings are a shrinking share of the grid: at N = 128 the typical miss is 0.2 percent. Small numbers, but the wind method is exact only when they are zero, and in that same test the error shrinks sevenfold between the raw grid and the relaxed one.

Where the numbers live

A cell has to hold three things: how much air it has, which Part 4 measures by the pressure at the ground, how warm it is, and how fast it moves. Where exactly on the cell they live turns out to matter as much as anything in this part. The simplest choice puts all of them at the center. Akio Arakawa, who classified these arrangements in the 1970s, called that the A-grid. His C-grid keeps pressure at the center but stores the wind on the edges, and only the part of it that crosses each edge: one number per edge, not a vector.

The difference shows up in the two operations a model does most. To move air, it needs the pressure difference across each edge, which is what pushes the wind that lives there. To change a cell's pressure, it adds up the air flowing in through the edges. On the C-grid both use neighbors that touch. On the A-grid the push on a cell has to be estimated from the pressures of its neighbors, and the cell's own pressure drops out, because its contributions to the six edges cancel. So a pattern that alternates from one cell to the next can be invisible: every cell's neighbors agree with each other, the estimated push is zero, and the model cannot feel the pattern, let alone remove it. The figure starts with such a pattern. Switch grids, then paint patterns of your own with the brush and see which ones each grid can feel.

Nobody has to paint these patterns. Every step a model takes rounds its numbers, and the tiny noise that leaves behind contains a little of every pattern there is, including the ones the grid cannot feel. Those are the dangerous ones. A pattern the grid can feel gets pushed on, spread out and damped away like any other disturbance; a pattern it cannot feel sits unopposed while the model's other arithmetic, the terms that carry air and heat around, keeps feeding it. It grows a little every step, which is to say exponentially, and after enough steps it is no longer tiny. The first version of WebGCM was an A-grid model, and its runs died the same way every time: alternating bands of pressure along the equator that the grid could not feel, growing out of nothing until they outgrew everything else. Rewriting the core on the C-grid is what fixed it. The checkerboard preset is the purest such pattern on a hexagonal grid, with every cell surrounded by six of the opposite sign; try it on both grids.

Slope, spread and spin

The two operations of the last section have names, and there is a third that the grid answers just as easily. The pressure difference across each edge measures how the pressure slopes; the air added up through a cell's walls measures how the wind spreads out of the cell; and the wind added up around a corner where three cells meet measures how the air there spins. Mathematics calls them the gradient, the divergence and the curl, written ∇, ∇· and ∇×. The gradient is the one in Part 1's pressure gradient force, and every sideways change in Part 4's equations is one of the three. The figure is a flat patch of cells the size of the model's, with the pressure at their centers and the wind across their edges, and it runs the model's own code for all three. Choose what to show, and tap a cell, or a triangle in the spin view, to read it.

Start with the slope. Subtract the pressure in the cell on one side of an edge from the pressure in the cell on the other, divide by the 60 km between their centers, and you have how steeply the pressure changes across that edge. That is the gradient as the model keeps it, one number per edge, and it is exactly what the wind stored on that edge needs: Part 1's push, the change in pressure over a distance, divided by the density, with the distance measured across the edge. In the steady slope the pressure falls 1 hPa every 100 km toward the east. The two walls of each cell that face east and west see all of it, 0.6 hPa between centers 60 km apart, and the four slanting walls see half, the part of the slope that crosses them. Put a cell's six edges back together, each along its own direction, and the whole slope returns as an arrow on the cell: 1.0 hPa every 100 km, 1.00 × 10⁻³ Pa/m, and a push of 3.0 m/s an hour toward the east, the numbers of Part 1. The arrow ∇p points up the slope, toward higher pressure, and the push, −∇p/ρ, points down it. The model never pushes the wind with that cell arrow. It keeps the six edge numbers, because summing them back to the cell is the A-grid's estimate, the one that cannot feel the checkerboard. In the low every push points in, strongest about 100 km out, 2.02 × 10⁻³ Pa/m and 6.1 m/s an hour, where the pressure falls off fastest, and at the center, at 996 hPa, the pressure is level and nothing pushes. This is the ∇ of Part 4, and the model takes every slope in its equations this way, one edge at a time.

The spread comes from the walls. Take the wind across each of a cell's six walls, counted positive outward, multiply each by the wall's length, 34.6 km, add them all the way around, and divide by the cell's area, 3,118 km². That is the divergence, ∇·u: how fast air leaves the cell, as a share of the cell's air per second. It is positive where the air spreads out and the cell drains, and negative where it converges and the cell fills. In the converging wind the figure picks the cell 120 km east of the middle. Counterclockwise from its east wall, the wind out through its six walls is −5, −3, +1, +3, +1 and −3 m/s: 5 m/s in through its outer wall, 3 m/s out through its inner one, and around all six walls 6 m/s more in than out. Its divergence is −6.67 × 10⁻⁵ per second, and the cell gains 24 percent of its air an hour. The middle cell takes in the same 6 m/s more simply, 1 m/s through every wall, and every other cell gains the same 6 m/s net, though the wind grows stronger outward. The other winds move air around and pile it up nowhere, the spinning disk included: their divergence is zero, in the lone vortex to within two parts in a thousand of the converging wind's, a small grid error of the kind the spin shows below. Part 4 meets the middle cell again, one wall at a time: the divergence of the air the wind carries is the heart of its first rule.

The spin needs a different shape. Around each corner where three cells meet, the lines joining their centers make a triangle, and each side of it crosses one edge at right angles, so the wind the grid stores on that edge is the wind along that side. Add the three, each times its side's length, 60 km, going counterclockwise around the triangle, and divide by the triangle's area, 1,559 km². That is the curl, ∇×u, how fast the air there is spinning about the vertical. Like the gradient it is an arrow, and for a level wind it points straight up or straight down, so one signed number says it all: positive counterclockwise seen from above. The number comes out at twice the rate at which a tiny paddle wheel set in the air would turn. For a wind it is called the vorticity, the ζ of Part 4. The spinning disk turns once every 23 hours 56 minutes, like the ground at the North Pole, and every triangle reads 1.46 × 10⁻⁴ per second, twice the disk's own rate of turning, 7.29 × 10⁻⁵ radians per second. The triangle the figure starts on, with one corner at the center, shows how. Two of its sides run straight out from the center, so the disk's wind blows square across them and none of it runs along them; the third carries 3.79 m/s, and 3.79 m/s times 60 km, divided by 1,559 km², is the whole of it. That factor of two is the 2 in Part 2's f = 2Ω sin φ. The ground turns at Ω sin φ about the local vertical, so air at rest on it, seen from space, already has a curl of f, which is why Part 4 adds f to ζ. The breaking jet of Part 5 is colored by the curl of its flow, computed this way at every corner of the grid and averaged around each cell.

Spin does not need anything to go around. The shear blows straight from the west, 12 m/s through the middle of the patch and 5 m/s faster every 100 km to the north, and every triangle reads −5.00 × 10⁻⁵ per second, clockwise: the stronger wind on the north side of a paddle wheel would turn it. The lone vortex is the opposite case. The twelve triangles nearest its center turn like a solid disk and read 2.00 × 10⁻⁴ per second; the ring that straddles the core's edge at 100 km reads 1.42 × 10⁻⁴, part disk and part not. Beyond that the air still circles the center, but its speed falls off as one over the distance, just fast enough that its shear cancels its turning, and its curl is zero. The triangles out there read up to 7 percent of the core's value, of either sign, the error a grid of 60 km makes on a wind that changes this much across a single triangle. None of the three needed more than the grid stores: on the C-grid each comes from neighbors that touch. The Coriolis force needs one number more, and it is the price of the C-grid, the subject of the next section.

The wind a cell does not store

One number per edge, the wind across it, is enough to move air around. It is not enough for the Coriolis force of Part 2, which pushes sideways, at right angles to the wind. To apply it at an edge you need the wind along the edge: along the wall itself, not along the line between centers, which the last section used. That number lives nowhere. The scheme the model uses, called TRiSK after the four people who worked it out (Thuburn, Ringler, Skamarock and Klemp, in 2009 and 2010), reconstructs it as a weighted sum of the across-edge winds on the ten edges of the two cells that share the edge. The weights come from the geometry alone, and they are chosen so that the reconstruction conserves energy exactly and treats the spin of the flow consistently. That is the step which needed every cell edge to be cut at its midpoint, and the reason the grid was nudged.

On the relaxed grid the reconstruction is accurate to better than one percent of the wind speed for a smooth flow, whatever its direction, and the typical error shrinks as the grid gets finer; the worst edges are the ones next to pentagons. The model does this for every edge at every step. It might seem simpler to store the along-edge wind too and skip the sum, but there is nothing new to store: the across-edge winds of a cell's neighbors already fix it, and a separately stored copy would drift away from them and take the scheme's conservation of energy with it. The weights are chosen for something else: the Coriolis force turns the wind without ever speeding it up or slowing it down, as the inertial circles of Part 2 showed, so it hands the wind no energy. The weights make the grid's version keep that property exactly, so that the reconstruction can never create or destroy energy however it errs at a given edge, and over a long run that matters more than a fraction of a percent at a few edges.

Slicing the column

Everything in this part so far lies on the surface of the sphere, but the air above each cell reaches up tens of kilometers, most of it in the lowest ten, and it has to be sliced too. Part 1 stacked its columns into layers, twelve to a column, without saying where one layer ends and the next begins. The obvious answer is height, a level every kilometer like the floors of a building, and it has two problems. The ground is not flat: a level 1 km up runs into every mountain range taller than that, so the cells over them lose their lowest levels, and a cell 60 km across with valleys and ridges in it, has levels that are partly underground. And a slice between two heights holds whatever the weather has put there, while the pressure at any level is the weight of the air above it, the lesson of Part 1's column of air, so a slice between two pressures always holds the same amount of air.

Slicing by pressure fixes the second problem but not the first, because the ground itself stands at different pressures in different places. In the standard atmosphere, taken with 1000 hPa at sea level as in Part 1, the ground in Denver, 1.6 km up, stands at about 820 hPa. The Tibetan Plateau, 4.5 to 5 km up, stands at 570 to 530 hPa. The summit of Everest, 8.85 km up, stands at 310 hPa, with more than two thirds of the air below it; the real summit, in air warmer than standard, measures about 337 hPa. The 900 hPa surface, 880 m above the sea, is underground in Denver, and the 700 hPa surface, at 2.9 km, is underground across Tibet. The figure draws both kinds of level over a mountain whose height you can set, and dashes every stretch the mountain cuts.

At 3.2 km the mountain cuts three height levels, and switched to pressure levels it cuts the three from 900 to 700 hPa; the summit itself stands at 675 hPa. The way out is to measure each column against its own ground. Divide the pressure at a point by the pressure at the ground beneath it and call the ratio σ. That is the sigma coordinate. At the ground σ is 1 everywhere, at sea level and on a summit alike, and at the top it is 0, where the pressure is zero, so no air is left out above the model. A layer between two values of σ holds the same share of its column wherever it is: between 0.9 and 1, a tenth, which is 100 hPa of air over the sea and 57 hPa over a plateau 4.5 km up. Switch the figure to σ levels and they drape over the mountain, and no level is ever cut, however high the mountain grows. Against height they ride up over the mountain almost together: where the summit rises 3.2 km, the σ = 0.1 surface rises 2.5 km, because a surface's displacement in pressure, σ times the drop in surface pressure, shrinks toward the top about as fast as the air thins. Put pressure up the page instead and the same surfaces fan out from the ground to a flat top at zero pressure, each displaced in proportion to its σ. The columns of Part 1 were sliced this way all along: each of their twelve layers held a twelfth of the air in its column.

The model's own layers, the figure's last setting, are slices of σ as well, but not evenly spaced. There are 36 of them, thinnest in meters near the ground, in the boundary layer, where the air trades heat, water and drag with the surface and changes fastest with height. The lowest is 40 m thick over the sea and holds less than half a percent of its column, and nine layers fit in the lowest kilometer, as the figure's inset shows. Above that they spread out, and the top layer holds everything above 0.3 hPa, more than 50 km up, which is three hundredths of a percent of the air and reaches all the way to σ = 0.

Layers that follow the ground come at a price. Along a layer that climbs a mountain, moving sideways also means moving up, so the pressure push along the layer comes in two large pieces that nearly cancel: one from the slope of the layer and one from the fall in pressure along it. Their small difference is where σ models make their best-known error, and Part 4 shows it over a mountain. In return the bookkeeping is easy. Because every layer's share is fixed, one number per column, the pressure at the ground, says how much air is in every layer above it, and keeping count of the air comes down to following that number. That is the subject of Part 4's first rule.

That is the whole kit: a grid of nearly equal cells, pressure at their centers and wind across their edges, the slope, spread and spin read off them, a way to recover the wind along an edge, and columns sliced into layers that follow the ground. Part 4 writes down the rules that move the air through them, and Part 5 puts the pieces on the sphere, on the turning planet of Part 2, and lets them run.