Vertical Coordinate Systems: The Foundation of Earth System Models

Modeling the Earth's atmosphere and oceans requires translating the continuous laws of physics onto a discrete computational grid. A foundational choice in this process is the selection of a vertical coordinate system, a decision that dictates how we slice the fluid world into layers for simulation. This choice is far from a simple technicality; it presents a core dilemma for modelers. How can one accurately represent both the complex, rugged topography of the Earth's surface and the subtle physical forces that drive motion, all while avoiding the pitfalls of numerical errors? Different coordinate systems offer different solutions, each with its own set of powerful advantages and weaknesses.

This article explores the critical trade-offs inherent in this decision. The following chapters will guide you through this complex landscape. The "Principles and Mechanisms" chapter dissects the fundamental concepts behind geopotential, pressure, terrain-following, hybrid, and isentropic coordinate systems, revealing the physical reasoning and numerical challenges associated with each. Following this, the "Applications and Interdisciplinary Connections" chapter illustrates the real-world consequences of these choices, showing how the abstract geometry of a grid can profoundly impact simulations of everything from deep ocean currents to global climate phenomena.

To simulate the grand, swirling dance of the atmosphere and oceans, we must first write its rules. These rules are the laws of physics, expressed as mathematical equations. But how do we solve these equations for a fluid sloshing around on a spinning, bumpy sphere? The first, most fundamental choice a modeler must make is how to draw the map, how to lay down the grid upon which the drama of weather and climate will unfold. This is the choice of a ​​vertical coordinate system​​. It may sound like a drab technicality, but as we shall see, it is a choice fraught with deep physical meaning, clever tricks, and devilish pitfalls. The story of vertical coordinates is a journey into the heart of how we translate the elegant laws of nature into the practical art of prediction.

Imagine the most straightforward way to build a model of the atmosphere. You might create a three-dimensional grid, much like the floors of a skyscraper, where each level is perfectly flat and parallel to the one below it. This is the ​​geopotential height coordinate​​, or ​​$z$-coordinate​​, system. Here, the vertical coordinate is simply the geometric height above sea level.

The beauty of this "grid of perfect flatness" is its mathematical simplicity. The force that drives all winds—the ​​pressure gradient force​​, which arises from differences in air pressure—is calculated by simply measuring how pressure changes along these perfectly horizontal grid surfaces. The equations look clean and elegant.

But the Earth, inconveniently, is not flat. Mountains, some of the most dramatic drivers of weather, come crashing through this pristine grid. In a $z$-coordinate model, a mountain range is not a smooth, sloping surface but a crude, jagged staircase of blocked-out grid cells. This "stair-step" representation makes it notoriously difficult to accurately model crucial processes that happen near the ground, like friction in the boundary layer, which feels the terrain's true shape. It's like trying to describe the curve of a sculpture using only Lego bricks.

Instead of imposing our simple grid on the complex Earth, perhaps we could be cleverer. What if we chose a coordinate system that listens to the physics of the atmosphere itself? For large-scale weather systems—the cyclones and anticyclones that span continents—the atmosphere is in a state of exquisite equilibrium known as ​​hydrostatic balance​​. This is a profound statement: the upward-pushing pressure gradient force is almost perfectly cancelled by the relentless downward pull of gravity.

A scale analysis reveals just how superb this approximation is. For a typical weather system, the vertical acceleration of air is millions of times smaller than the acceleration due to gravity. The atmosphere is not shooting up and down; it is in a grand balancing act. This balance leads to a magical consequence. If we integrate the hydrostatic equation, $\frac{\partial p}{\partial z} = -\rho g$, we find that the total mass of air in a column between any two pressure surfaces is directly proportional to the pressure difference between them.

Suddenly, pressure isn't just a force; it's a proxy for ​​mass​​. This insight allows for a brilliant change of coordinates. What if we abandon height and use ​​pressure ($p$)​​ as our vertical coordinate? When we rewrite the governing equations in this new language, a wonderful simplification occurs. The mass continuity equation, which describes how density changes as air moves, sheds its explicit dependence on density, $\rho$. This is a huge victory for numerical modelers. By aligning our mathematics with a fundamental physical balance, we've made the problem simpler to solve. We've chosen the right language to describe the phenomenon.

Using pressure coordinates helps, but it doesn't fully solve our mountain problem, as pressure surfaces themselves can intersect the ground. So, scientists came up with another ingenious idea: a ​​terrain-following coordinate​​, often called a ​​sigma ($\sigma$) coordinate​​. The idea is to create a "stretchy" grid that is squashed to fit the Earth's surface. A common definition is $\sigma = \frac{p}{p_s}$, where $p$ is the pressure and $p_s$ is the pressure at the ground.

By this definition, the surface of the Earth, no matter how rugged, is always the single coordinate surface $\sigma=1$. The top of the atmosphere might be $\sigma=0$. The entire grid of the model elegantly drapes over the mountains and valleys. This has a fantastic benefit for boundary conditions. The physical rule that air cannot flow through the ground becomes the trivially simple mathematical statement that the vertical velocity in sigma-coordinates, $\dot{\sigma}$, is zero at the bottom boundary. This smooth representation is far better for modeling the boundary layer than the ugly stair-steps of a $z$-coordinate model.

But this cleverness comes at a terrible price. The villain of our story now enters with a vengeance: the infamous ​​Pressure Gradient Force (PGF) error​​. Recall that the PGF is what drives the wind, and it's physically defined by pressure differences on a constant height surface. But in our new $\sigma$-system, the coordinate surfaces are no longer horizontal; they slope steeply over mountains. To calculate the true horizontal PGF, a model must compute the pressure gradient along its sloped grid line and then subtract a large correction term involving the slope of the grid surface.

Imagine trying to determine if a ping-pong ball will roll on a table that's sitting on the deck of a violently pitching ship. You measure the slope of the table relative to the deck, but to find out which way the ball will actually roll, you must subtract the enormous, rapidly changing slope of the entire ship. This is what the computer must do. It has to calculate the PGF as the small difference between two very large, opposing numbers. Digital computers have finite precision, and tiny rounding errors in the two large terms can lead to a massive error in their small difference. This ​​spurious force​​ can be so large that it generates unrealistic winds out of thin air, especially high in the atmosphere over steep mountain ranges. The PGF error is a direct, and dangerous, consequence of the misalignment between our chosen coordinate surfaces and the true horizontal.

So, we have a dilemma. Height coordinates have a simple PGF but a terrible lower boundary. Sigma coordinates have a simple lower boundary but a terrible PGF. What can we do? We can compromise. This is the engineer's solution, and it is the foundation of most modern weather and climate models: the ​​hybrid coordinate ($\eta$)​​.

The idea is to have the best of both worlds. We want our grid to follow the terrain near the surface, but we want it to flatten out and become parallel to pressure surfaces at high altitudes, where the PGF error is most severe. A hybrid coordinate is designed to do just this. It is defined by a formula like $p(\eta) = A(\eta) + B(\eta) p_s$, where the functions $A$ and $B$ are carefully chosen. Near the surface, the coordinate behaves like $\sigma$, slavishly following the terrain. As we move up through the atmosphere, it smoothly transitions, and high in the stratosphere, it becomes a pure pressure coordinate, with coordinate surfaces that are nearly flat.

This elegant blend mitigates the PGF error aloft while retaining the beautifully simple boundary representation at the ground. It is a pragmatic and powerful solution to a very difficult problem.

So far, our choices have been guided by geometry (mountains) and force balances (hydrostatics). But there is another, perhaps more profound, way to think. What if we let the thermodynamics of the fluid define our grid?

Let's consider a property called ​​potential temperature ($\theta$)​​. You can think of it as the temperature an air parcel would have if you moved it from its current pressure to a standard reference pressure, without letting it exchange any heat with its surroundings. For the large-scale, nearly heat-free (adiabatic) motions that dominate the atmosphere, potential temperature is conserved. A parcel of air, once set in motion, is forever trapped on the surface of its initial potential temperature. These surfaces are called ​​isentropic surfaces​​.

This is a deep physical constraint. It suggests a radical choice of coordinate: use $\theta$ itself as the vertical coordinate! What happens if we do this? The vertical velocity in this system is simply the rate of change of a parcel's potential temperature, $\frac{D\theta}{Dt}$. For adiabatic motion, this is zero! In an isentropic coordinate system, the complex three-dimensional motion of the atmosphere becomes a series of stacked, purely two-dimensional flows. Air parcels glide along these $\theta$-surfaces, never crossing them.

The payoff is spectacular. Imagine tracking a plume of pollution. In a pressure-coordinate model, as the plume drifts along a sloping isentropic surface, it must be numerically moved from one pressure-based grid level to another. This process of interpolation inevitably smears the plume out, an artifact we call ​​spurious numerical diffusion​​. In an isentropic-coordinate model, the plume stays on its original grid level. There is no vertical transport to simulate, and therefore no spurious vertical mixing. This makes isentropic coordinates extraordinarily powerful for applications like tracking chemical tracers or understanding the transport of ozone.

Of course, no choice is without its own challenges. Isentropic surfaces can become steeply sloped or even fold over. In the ocean, where density surfaces can intersect the surface in a process called "outcropping," a pure isopycnal (constant-density) model would fail. This again necessitates hybrid schemes that blend isopycnal coordinates in the deep ocean with height-based coordinates near the turbulent surface. Even in these sophisticated models, tiny misalignments between the coordinate grid and the true density surfaces can introduce spurious mixing, and modelers must carefully ensure this numerical error does not overwhelm the small but crucial physical mixing that occurs in the real ocean.

The choice of a vertical coordinate is a beautiful microcosm of the entire field of geophysical modeling. It is a constant dialogue between the immutable laws of physics, the complex geometry of the Earth, and the clever, pragmatic art of computation. There is no single "best" answer, only a series of trade-offs, each choice illuminating a different facet of the fluid's intricate behavior.

We have spent some time understanding the various ways one might choose to slice up the world—be it the ocean or the atmosphere—into a stack of layers for a computer model. You might be tempted to think this is a mere technicality, a dry bookkeeping exercise for the programmer. But nothing could be further from the truth. The choice of a vertical coordinate system is not just a choice of grid; it is the choice of a lens through which our model views the world. And just as different lenses can reveal, distort, or hide features of a landscape, so too can different coordinate systems illuminate, obscure, or even create phantom physics in our simulations. This choice has profound, far-reaching, and often surprising consequences, rippling out from the tiniest boundary layer to the grandest climate patterns. Let us now explore this fascinating landscape where the abstract geometry of coordinates meets the tangible reality of the fluid Earth.

Imagine dense, salty water from the Mediterranean Sea spilling over the submarine sill at Gibraltar and plunging into the Atlantic. Or picture a cold, dense current cascading down the continental slope. How do we capture this dramatic, plunging flow in a model? Here we face a classic dilemma, a kind of "original sin" of ocean modeling.

If we use simple, horizontal $z$-levels, our smoothly sloping seabed becomes a crude staircase. As the dense water tries to flow down this staircase, our advection schemes inevitably cause it to mix with the lighter water in the steps below. This is not real, physical mixing; it is a numerical artifact, a "spurious diapycnal mixing" that can fatally erode the very dense water plume we are trying to simulate. It is as if our simulation is haunted by a ghost that constantly stirs the water, weakening the stratification wherever the grid has steps.

To avoid these clumsy staircases, we might be tempted to use a coordinate system whose layers bend and follow the terrain, the so-called $\sigma$-coordinates. This seems much more elegant; the bottom boundary is now perfectly smooth. But we have traded one devil for another. In a stratified ocean at rest, the surfaces of constant pressure (isobars) are perfectly flat. Our new, sloping coordinate surfaces are not. To calculate the horizontal pressure force that drives the flow—which is all about the slope of the isobars—we now have to compute the difference between two very large, nearly equal numbers. In the continuous, perfect world of mathematics, they cancel out perfectly to give the correct, zero force for a fluid at rest. But in the finite, approximate world of a computer, tiny truncation errors in this subtraction leave behind a residual, a phantom force. This "pressure gradient error" can be large enough to spin up spurious currents from nothing, especially over steep slopes. The model, in a sense, sees a hill where there is none.

This trade-off is fundamental. Do we accept a blurry, overly mixed world with $z$-levels, or a world of phantom forces with $\sigma$-coordinates? The answer depends on the problem, and much of the art of numerical oceanography lies in navigating this compromise, for example by designing "hybrid" coordinates that try to get the best of both worlds.

The challenge deepens when we look closer at the physics we need to resolve. The bottom boundary layer, where the water feels the friction of the seabed, is a region of critical importance. How well does our grid "see" this layer? We can quantitatively assess whether our coordinate layers are well-aligned with the boundary and if we have enough layers packed inside the thin boundary region to resolve its structure. We quickly find that for a sloping bottom, $\sigma$-coordinates naturally concentrate layers near the boundary and align with it, while $z$-coordinates do a poor job on both counts. Again, the choice is not trivial; it determines whether we can even represent the physics we care about.

The consequences of this boundary interaction are not just local. Consider the tides. The great, basin-spanning barotropic tide, a depth-uniform sloshing of the entire water column, contains immense energy. When this flow encounters a bump on the seafloor—a seamount or a mid-ocean ridge—it must rise and fall. This vertical motion, forced by the boundary condition $w \approx \mathbf{U} \cdot \nabla h$, pushes and pulls on the stable stratification of the ocean's interior. This disturbance then radiates away as a new kind of wave: an internal tide, a depth-dependent, or baroclinic, motion. A vast amount of energy is transferred from the barotropic tide to these internal waves, which carry this energy thousands of kilometers before breaking and mixing the ocean. This process is a key driver of the global ocean circulation. A $z$-coordinate model captures this energy conversion naturally, because the interaction between the flow $\mathbf{U}$ and the bottom topography $\nabla h$ is explicit. The staircase grid might be clumsy, but it faithfully transmits the vertical push that generates these crucial waves.

You might think that these numerical artifacts, however annoying, are small-scale issues. Surely they average out and do not affect the big picture, like the climate? This is a dangerous assumption. Let us travel to the equatorial Pacific, the home of the El Niño–Southern Oscillation (ENSO), a climate pattern whose influence is felt across the globe.

The physics of ENSO is a delicate dance between the atmosphere and the ocean, critically dependent on the temperature of the sea surface and the structure of the thermocline—the sharp transition layer between the warm surface waters and the cold abyss. In a $z$-level model, the spurious mixing we discussed earlier is constantly at work, especially where the thermocline is tilted. Along the equator, this mixing artificially weakens and deepens the thermocline. A weaker thermocline means that when the ocean upwells, it brings less-cold water to the surface, damping the temperature anomalies that drive ENSO. Furthermore, the speed of the equatorial Kelvin wave, a key messenger that carries signals across the Pacific, is set by the strength of the stratification. A spuriously weakened stratification slows down the wave, altering the timing and period of the entire ENSO cycle. The end result? A seemingly small, technical choice in how we draw our vertical grid can lead to a model that simulates ENSO with the wrong amplitude, the wrong frequency, and the wrong timing, crippling our ability to predict this vital climate phenomenon.

The influence of coordinates extends to how we represent physics that the grid cannot see at all. Processes like vertical convection—the turbulent overturning of an unstable water column—happen at scales far too small for a global model to resolve. We must parameterize it, writing a recipe for the model to follow whenever it finds denser water on top of lighter water. But the recipe itself depends on the coordinate system! In a $z$-level model, the recipe is simple: find unstable layers and mix them vertically. In a $\sigma$-coordinate model, one must be careful not to mix along the sloping coordinate surfaces, which would create enormous spurious mixing; instead, one must map the water column back to true vertical ($z$) space, perform the mixing, and then map it back. In an isopycnal model, where layers are already defined by density, instability appears as a bizarre "folding" or re-ordering of layers in physical space, and fixing it requires a careful exchange of mass between layers. Each coordinate system demands its own bespoke algorithm to represent the very same physical process.

So far, our story has been one of troublesome trade-offs. But sometimes, a clever choice of coordinates can make complexity melt away, revealing a beautiful, underlying simplicity. Nowhere is this clearer than in the atmosphere, with the use of pressure as a vertical coordinate.

The atmosphere, for large-scale motions, is in a state of near-perfect hydrostatic balance: the pressure at any point is simply the weight of the air above it. This leads to a wonderfully simple relationship: a small change in pressure, $dp$, is directly proportional to the mass of a thin slice of air, $dm$. Specifically, $dm = -dp/g$, where $g$ is the acceleration of gravity.

Why is this so magical? Many of the most complex physical processes in the atmosphere, like the absorption and emission of radiation or the condensation of water in clouds, are fundamentally framed in terms of mass. For example, a heating rate is most naturally expressed in joules per second per kilogram of air. In a $z$-coordinate model, to find the total heating in a layer, we would need to integrate this rate multiplied by the density $\rho$, which varies with height in a complicated way. But in a pressure-coordinate model, this integral over mass becomes a simple integral over pressure! A mass-weighted average is identical to a pressure-weighted average. The complex, spatially varying density field vanishes from the calculation, replaced by the simple pressure interval $\Delta p$ of the layer. This makes coupling the dynamical core of the model to the complex physics modules incredibly simple and elegant. Furthermore, the total mass of the atmospheric column is directly proportional to the surface pressure. This ensures that the model conserves mass, energy, and water perfectly, a task that is much clumsier in other coordinate systems. The pressure coordinate system is a triumph of finding the "natural" variable for the problem, a choice that makes the physics clearer and the numerics simpler.

There is no free lunch, however. The choice of coordinates also has a very practical price tag: computational cost. Explicit numerical models are bound by a stability constraint known as the Courant-Friedrichs-Lewy (CFL) condition, which states that information cannot be allowed to travel more than one grid cell per time step. For vertical motion, this means the time step $\Delta t$ must be less than $\Delta s / |c_s|$, where $\Delta s$ is the vertical grid spacing and $c_s$ is the vertical velocity. If we choose a $z$-coordinate system with very fine vertical resolution (small $\Delta z$) to resolve a boundary layer, we may be forced to take incredibly tiny time steps, making the simulation prohibitively expensive. An isopycnal model, on the other hand, might allow for a much larger time step for vertical processes, because the "velocity" across density surfaces is typically very small. The choice of coordinates thus directly impacts the speed and feasibility of our simulations.

To get around these limitations, modelers have developed sophisticated tricks. One such trick is "mode splitting." The ocean's motion can be separated into a fast, depth-averaged (barotropic) part, and a slower, depth-dependent (baroclinic) part. It turns out that when one derives the equations for the fast barotropic mode, which involves the sea surface height, the complexities of the 3D vertical coordinate system get "integrated out." The problem collapses to a much simpler 2D equation that is the same whether you started with $z$-levels or $\sigma$-coordinates. Modelers can then solve this 2D problem with a short time step, and the full 3D problem for the slow internal motion with a much longer, more efficient time step. The vertical coordinate choice still affects the coefficients in the 2D problem, but its fundamental structure is unchanged. It is a clever way of acknowledging that different physical phenomena "see" the grid in different ways.

Finally, our journey takes us to the frontier where the idealized world of the model meets the messy, chaotic real world of observations. Our weather and ocean forecasts are not just the product of a free-running model; they are the result of constantly nudging the model back toward reality using a process called data assimilation. What happens when our observations live in a different coordinate system from our model?

Imagine a robotic ocean glider that is programmed to measure temperature not at fixed depths, but along a surface of constant density (an isopycnal). We now want to assimilate this measurement into our $z$-level model. The task is no longer a simple interpolation. To find the model's equivalent of the observation, we first have to find the depth $z$ where the model's density field matches the observed density. This depth depends on the model's temperature and salinity fields, making the "observation operator" $\mathcal{H}$ that maps the model state to the observation a complex, nonlinear function.

Worse, the real isopycnal surface is not a smooth sheet; it is corrugated with unresolved small-scale waves and eddies. The glider measures temperature on this wavy, true surface. The model, however, only knows about its own smooth, large-scale representation of that surface. The difference in temperature between the true, wavy location and the model's smooth location is not an error in the instrument, but an error in representation. This "representativeness error" is a real headache. A vertical displacement $\eta$ of the isopycnal creates an apparent temperature error of $\eta \times (\partial T / \partial z)$. If we know something about the statistics of the unresolved waves (the variance of $\eta$), we can estimate the variance of this error and tell the assimilation system to give this observation a little less weight. This shows that the choice of coordinates has profound implications not just for how we simulate the world, but for how we fuse our simulations with measurements to create the best possible picture of our planet. The conversation between model and data is spoken in the language of coordinates, and understanding that language is essential to understanding our world.