Intro to weather modeling

19 May, 2026

Edit history:

2026-05-20: Published first draft

Modeling weather is an extremely challenging problem. If you’re good at it, you can make a lot of money. Polymarket has weather markets. Citadel’s commodity trading desk has its own proprietary weather models. It is rumored that Ken Griffin would call his internal team for the weather report rather than check a weather app. So how do people model weather?

Spherical approximations

As we know, the Earth is NOT flat. It is, in fact, an ellipsoid, which we approximate as a sphere. In spherical coordinates, we can define the mapping as follows:

$$f\begin{pmatrix} \phi\\ \theta \end{pmatrix} = R\begin{pmatrix} \sin \phi \cos \theta \\ \sin \phi \sin \theta \\ \cos \phi \end{pmatrix}$$

where $R$ is the radius of the Earth, approximately 6371 km. As we know, we compute the Jacobian of $f$, whose columns are valid basis vectors for the tangent space. In other words:

$$f_\theta= R\begin{pmatrix} -\sin \phi \sin \theta \\ \sin \phi \cos \theta \\ 0 \end{pmatrix}, \quad f_\phi = R\begin{pmatrix} \cos \phi \cos \theta \\ \cos \phi \sin \theta \\ -\sin \phi \end{pmatrix}$$

One thing to note is that not only do these vectors serve as a basis for the tangent space, they are also orthogonal to each other. Specifically, $f_\theta \cdot f_\phi = 0$. We also have that:

$$||f_\phi||_2 = R, \quad ||f_\theta||_2 = R \sin \phi$$

Atmospheric fields

Because we have made this approximation, it is natural to define scalar atmospheric fields as functions on the sphere. For example, an atmospheric field $A(\phi, \theta)$ is a function that maps a point on the sphere to a scalar value. The coordinate-based total derivative is given by:

$$dA = A_\phi d\phi + A_\theta d\theta$$

This, however, is not ideal because we are relating changes in the two-dimensional system of $(\phi, \theta)$ to changes in the atmospheric field. Instead, it would be more intuitive to relate small physical displacements on the sphere to changes in the atmospheric field. For this to make sense, we need:

$$dA = \nabla_{\text{surface}} A \cdot d\vec{x}$$

Recall that a small displacement vector on the sphere can be written as:

$$d\vec{x} = R e_\phi d\phi + R \sin \phi e_\theta d\theta$$

And because $e_\phi$ and $e_\theta$ are unit vectors that span the tangent space, we also have that:

$$\nabla_{\text{surface}} A = g_\phi e_\phi + g_\theta e_\theta$$

For consistency, we must have:

$$\nabla_{\text{surface}} A \cdot (R e_\phi d\phi + R \sin \phi e_\theta d\theta) = A_\phi d\phi + A_\theta d\theta$$

This implies that:

$$g_\phi = \frac{A_\phi}{R}, \quad g_\theta = \frac{A_\theta}{R \sin \phi}$$

In other words, the surface gradient in physical dimensions is:

$$\nabla_{\text{surface}} A = \frac{A_\phi}{R} e_\phi + \frac{A_\theta}{R \sin \phi} e_\theta$$

The Laplacian

Just as with Fourier modes, we want to define a set of basis functions on the sphere. We can take them to be the eigenfunctions of the spherical Laplacian. The Laplacian is also known as the divergence of the gradient. Let’s compute the divergence of the gradient field above. To do so, consider a rectangle with corners $[\theta, \theta + \Delta \theta] \times [\phi, \phi + \Delta \phi]$. This is mapped to a parallelogram of area $R^2 \sin \phi \, \Delta \theta \, \Delta \phi$. The outward normal on the $\theta + \Delta \theta$ face is $e_\theta$, and the outward normal on the $\phi + \Delta \phi$ face is $e_\phi$.

Now suppose we have a vector field on the sphere, $V(\phi, \theta) = V_\phi e_\phi + V_\theta e_\theta$. To compute the divergence, consider the face at $\theta + \Delta \theta$. The flux through it is:

$$V(\phi, \theta + \Delta \theta) \cdot e_\theta \cdot R \Delta \phi$$

On the face at $\theta$, the flux is:

$$-V(\phi, \theta) \cdot e_\theta \cdot R \Delta \phi$$

The net flux through the $\theta$-faces is therefore:

$$(V_\theta(\phi, \theta + \Delta \theta) - V_\theta(\phi, \theta)) \cdot R \Delta \phi = \frac{\partial V_\theta}{\partial \theta} R \Delta \phi \Delta \theta$$

Similarly, for the $\phi$-faces, the flux through the top face is:

$$V(\phi + \Delta \phi, \theta) \cdot e_\phi \cdot R \sin (\phi + \Delta \phi) \Delta \theta$$

On the face at $\phi$, the flux is:

$$-V(\phi, \theta) \cdot e_\phi \cdot R \sin \phi \Delta \theta$$

The net flux through the $\phi$-faces is therefore:

$$(V_\phi(\phi + \Delta \phi, \theta)\sin(\phi + \Delta \phi) - V_\phi(\phi, \theta)\sin(\phi)) \cdot R \Delta \theta$$

In the limit, this simplifies to:

$$R\frac{\partial \sin\phi V_\phi}{\partial \phi} \Delta \phi \Delta \theta$$

The total flux is then:

$$R \frac{\partial V_\theta}{\partial \theta} \Delta \phi \Delta \theta + R\frac{\partial \sin\phi V_\phi}{\partial \phi} \Delta \phi \Delta \theta$$

Thus the divergence is:

$$\nabla \cdot V = \frac{1}{R \sin \phi} \left( \frac{\partial V_\theta}{\partial \theta} + \frac{\partial \sin\phi V_\phi}{\partial \phi} \right)$$

Applying this to the surface gradient, we have $V_\theta = \frac{A_\theta}{R\sin \phi}$ and $V_\phi = \frac{A_\phi}{R}$. The divergence of the gradient of the atmospheric field is therefore:

$$\nabla \cdot \nabla A = \frac{1}{R \sin \phi} \left( \frac{\partial}{\partial \theta}\left(\frac{A_\theta}{R\sin \phi}\right) + \frac{\partial}{\partial \phi}\left(\frac{A_\phi}{R}\right) \right)$$

In other words:

$$\nabla^2_{\text{surface}} A = \frac{1}{R^2 \sin^2 \phi} \frac{\partial^2 A}{\partial \theta^2} + \frac{1}{R^2 \sin \phi} \frac{\partial}{\partial \phi}\left(\sin \phi \frac{\partial A}{\partial \phi}\right)$$

Spherical harmonics

Since we have defined what the spherical Laplacian is, we can define a set of basis functions on the sphere given by the eigenfunctions of the Laplacian.

$$\nabla_{\text{surface}}^2 f = c f $$

Substituting the definition, we get:

$$\frac{1}{R^2 \sin^2 \phi} \frac{\partial^2 f}{\partial \theta^2} + \frac{1}{R^2 \sin \phi} \frac{\partial}{\partial \phi}\left(\sin \phi \frac{\partial f}{\partial \phi}\right) = c f$$

We can define further:

$$\begin{align*} \frac{1}{\sin^2 \phi}\frac{\partial^2 f}{\partial \theta^2} + \frac{1}{\sin \phi} \left( \cos \phi \frac{\partial f}{\partial \phi} + \sin \phi \frac{\partial^2 f}{\partial \phi^2}\right) &= c R^2 f \\ \frac{1}{\sin^2 \phi}\frac{\partial^2 f}{\partial \theta^2} + \frac{\cos \phi}{\sin \phi} \frac{\partial f}{\partial \phi} + \frac{\partial^2 f}{\partial \phi^2} &= c R^2 f \\ \end{align*}$$

Let’s now provide an ansatz for $f(\phi, \theta) = g(\phi) h(\theta)$. Substituting this into the equation, we get the following where:

$$\begin{align*} \frac{\partial^2 f}{\partial \theta^2} = g(\theta)h''(\theta) \\ \frac{\partial f}{\partial \phi} = g'(\phi) h(\theta) \\ \frac{\partial^2 f}{\partial \phi^2} = g''(\phi) h(\theta) \\ \end{align*}$$

Substituting these into the above, we get:

$$\begin{align*} \frac{1}{\sin^2 \phi} g(\theta)h''(\theta) + \frac{\cos \phi}{\sin \phi} g'(\phi) h(\theta) + g''(\phi) h(\theta) &= c R^2 g(\phi) h(\theta) \\ \end{align*}$$

Dividing by $g(\phi) h(\theta)$ and rearranging, we get:

$$\begin{align*} \frac{1}{\sin^2 \phi} \frac{h''(\theta)}{h(\theta)} + \frac{\cos \phi}{\sin \phi} \frac{g'(\phi)}{g(\phi)} + \frac{g''(\phi)}{g(\phi)} &= c R^2 \\ cR^2 - \frac{\cos \phi}{\sin \phi} \frac{g'(\phi)}{g(\phi)} - \frac{g''(\phi)}{g(\phi)} &= \frac{1}{\sin^2 \phi} \frac{h''(\theta)}{h(\theta)} \\ \sin^2 \phi \left( cR^2 - \frac{\cos \phi}{\sin \phi} \frac{g'(\phi)}{g(\phi)} - \frac{g''(\phi)}{g(\phi)} \right) &= \frac{h''(\theta)}{h(\theta)} \end{align*}$$

The left side only depends on $\phi$ and the right side only depends on $\theta$. This means that both sides are equal to a constant, which by convention we call $-m^2$. We get the following two equations:

$$\begin{align*} \frac{h''(\theta)}{h(\theta)} &= -m^2 \\ \sin^2\phi \left( cR^2 - \frac{\cos \phi}{\sin \phi} \frac{g'(\phi)}{g(\phi)} - \frac{g''(\phi)}{g(\phi)} \right) &= -m^2 \\ \end{align*}$$

The first equation gives us a simple differential equation:

$$h'' = -m^2 h \implies h(\theta) = e^{i m \theta}$$

The second equation is a more complex differential equation, which after rearranging we get:

$$\frac{g''(\phi)}{g(\phi)} + \frac{\cos \phi}{\sin \phi} \frac{g'(\phi)}{g(\phi)} = \left( \frac{m^2}{\sin^2 \phi} + cR^2 \right)$$

$$g'' + \frac{\cos \phi}{\sin \phi} g' = \left( \frac{m^2}{\sin^2 \phi} + cR^2 \right) g$$

Now let $x = \cos \phi$. We will simplify this equation down, where $\phi = \arccos x$. Define $G(x) = g(\arccos x)$. For the first derivative, we have:

$$G'(x) = -g'(\arccos x) \frac{1}{\sqrt{1-x^2}} \implies g'(\phi) = -G'(x) \sqrt{1-x^2}$$

Now for the second derivative, we have:

$$\begin{align*} g''(\phi) &= \frac{d}{d\phi}\left(-G'(x) \sin \phi \right) \\ &= G''(x) \sin^2 \phi - G'(x) \cos \phi \\ &= G''(x) (1 - x^2) - G'(x) x \\ \end{align*}$$

Now plugging this back into the original equation, we get:

$$(1 - x^2)G''(x) - 2xG'(x) = \left( \frac{m^2}{(1 - x^2)} + cR^2 \right) G(x)$$

Since we are operating in spherical harmonics, $c = -\frac{l(l+1)}{R^2}$. So we have:

$$(1 - x^2)G''(x) - 2xG'(x) = \left( \frac{m^2}{1 - x^2} - l(l+1) \right) G(x)$$

The solution to this associated Legendre differential equation are given by the associated Legendre polynomials. Thus, the spherical harmonics are given by:

$$\nabla^2_{\text{surface}} Y_l^m(\phi, \theta) = -\frac{l(l+1)}{R^2} Y_l^m(\phi, \theta)$$

The physics of weather

Coriolis effect

Modeling the Coriolis effect is a standard mechanics problem. We must consider the motion of a particle from an inertial and rotating frame. From the inertial reference frame, the particle’s position is given by (using the rotating basis as seen by the inertial frame):

$$r_{\text{inertial}}(t) = x(t) e_x(t) + y(t) e_y(t) + z(t) e_z(t)$$

From the rotating reference frame, the particle’s position is given by:

$$r_{\text{rotating}}(t) = x(t) e_x + y(t) e_y + z(t) e_z$$

In other words, for some rotation matrix $R(t)$, we have:

$$r_{\text{inertial}}(t) = R(t) r_{\text{rotating}}(t)$$

If we differentiate with time w.r.t to the inertial frame, we get:

$$v_{\text{inertial}}(t) = R'(t) r_{\text{rotating}}(t) + R(t) v_{\text{rotating}}(t)$$

But we know that:

$$R'(t)r_{\text{rotating}}(t) = x(t) e_x'(t) + y(t) e_y'(t) + z(t) e_z'(t) = \Omega \times r_{\text{inertial}}(t)$$

Thus, the above becomes:

$$v_{\text{inertial}}(t) = R(t) v_{\text{rotating}}(t) + \Omega \times r_{\text{inertial}}(t)$$

Now, we compute the acceleration in the inertial frame:

$$\begin{align*} a_{\text{inertial}}(t) &= \frac{d^2}{dt^2} r_{\text{inertial}}(t) \\ &= \frac{d}{dt} \left( R(t) v_{\text{rotating}}(t) + \Omega \times r_{\text{inertial}}(t) \right) \\ &= R'(t) v_{\text{rotating}}(t) + R(t) a_{\text{rotating}}(t) + \Omega \times v_{\text{inertial}}(t) \\ &= R(t) a_{\text{rotating}}(t) + \Omega \times v_{\text{inertial}}(t) + \Omega \times (R(t) v_{\text{rotating}}(t)) \\ &= R(t) a_{\text{rotating}}(t) + 2 \Omega \times (R(t) v_{\text{rotating}}(t)) + \Omega \times (\Omega \times r_{\text{inertial}}(t)) \\ \end{align*}$$

We can see from the above that:

$$a_{\text{inertial}} = R(t) a_{\text{rotating}} + \underbrace{2\Omega \times R(t) v_{\text{rotating}}}_{\text{Coriolis}} + \underbrace{\Omega \times (\Omega \times r_{\text{inertial}})}_{\text{centrifugal}}$$

In the inertial frame Newton’s second law reads $F = m a_{\text{inertial}}$. To get the equation of motion in the rotating frame, multiply through by $m$ and move the inertial pseudo-forces to the right-hand side:

$$m a_{\text{rotating}} = R^{-1}(t) F - 2m\,\Omega_{\text{rotating}} \times v_{\text{rotating}} - m\,\Omega_{\text{rotating}} \times (\Omega_{\text{rotating}} \times r_{\text{rotating}})$$

where we used $R^{-1}(t) = R^T(t)$ since rotation matrices are orthogonal, and $\Omega_{\text{rotating}}$ is expressed in the rotating frame.

Centrifugal force and effective gravity

The centrifugal term $-\Omega \times (\Omega \times r)$ can be simplified using the BAC-CAB identity $A \times (B \times C) = B(A \cdot C) - C(A \cdot B)$:

$$-\Omega \times (\Omega \times r) = r\,|\Omega|^2 - \Omega(\Omega \cdot r) = \Omega^2 r_\perp$$

where $r_\perp = r - \frac{\Omega(\Omega \cdot r)}{|\Omega|^2}$ is the component of $r$ perpendicular to the rotation axis. The centrifugal acceleration points radially outward from the Earth’s spin axis with magnitude $\Omega^2 d$, where $d = R\sin\phi$ is the distance to the axis.

For the Earth, $\Omega \approx 7.3 \times 10^{-5}\ \text{rad/s}$, so the maximum centrifugal acceleration (at the equator) is:

$$a_\text{cf} = \Omega^2 R \approx 0.034\ \text{m/s}^2 \approx 0.3\%\, g$$

Because it is static (independent of the parcel’s velocity) and small, the centrifugal acceleration is simply added to the Newtonian gravitational acceleration $g_0$ to define an effective gravity:

$$\mathbf{g}_\text{eff} = -g_0\,\hat{r} + \Omega^2 r_\perp\,\hat{r}_\perp$$

This is exactly what a plumb bob measures: not the pull toward the Earth’s center, but the vector sum of gravity and centrifugal acceleration. The Earth’s oblate shape is itself a consequence of this, and the surface bulges until it is everywhere perpendicular to $\mathbf{g}_\text{eff}$. From this point on we treat $g = |\mathbf{g}_\text{eff}|$ as a known constant and drop the centrifugal term from the equations of motion.

Hydrostatic balance

For most weather modeling, we assume a particularly useful property of the atmosphere that may or may not be true in reality. We assume that the atmosphere is in hydrostatic balance. This does not mean that vertical velocity is zero. It means that the vertical pressure force nearly balances gravity, so we do not explicitly model vertical acceleration using a vertical momentum equation.

To see where the equation comes from, consider a thin slab of air with horizontal area $A$ and height $dz$. The pressure at the bottom of the slab is $p(z)$ and the pressure at the top is $p(z+dz)$. The bottom face is pushed upward by the air below it, while the top face is pushed downward by the air above it. Gravity also pulls the slab downward.

Force balance gives:

$$p(z)A - p(z+dz)A - \rho g A\,dz = 0$$

Rearranging and dividing by $A\,dz$, then taking $dz\to 0$:

$$\frac{\partial p}{\partial z} = -\rho g$$

This is the hydrostatic balance identity. It says that pressure decreases upward at exactly the rate needed for the pressure difference across a thin layer to support the weight of that layer.

Equivalently, integrating upward from height $z$ to the top of the atmosphere gives

$$p(z)=\int_z^\infty \rho(z')g\,dz'$$

Thus pressure at a given height is the weight per unit area of the air above that height. This is a vertical statement, but it is still extremely important for horizontal motion: if neighboring columns have different temperature, moisture, surface pressure, or topography, then their vertical pressure profiles differ, which creates horizontal pressure gradients.

Equation of state and virtual temperature

Hydrostatic balance by itself relates pressure and density:

$$\frac{\partial p}{\partial z} = -\rho g$$

But in weather modeling we also need to know how density relates to temperature and moisture. This is where the ideal gas law enters. For dry air, the local equation of state is

$$p = \rho R_d T$$

where $R_d$ is the gas constant for dry air. Equivalently,

$$\rho = \frac{p}{R_dT}$$

or, if we define the specific volume

$$\alpha = \frac{1}{\rho}$$

then

$$\alpha = \frac{R_dT}{p}$$

For moist air, the model often uses virtual temperature $T_v$. The point of virtual temperature is to rewrite the moist-air equation of state in the same form as the dry-air equation:

$$p = \rho R_d T_v$$

Thus

$$\alpha = \frac{R_dT_v}{p}$$

A rough approximation is

$$T_v \approx T(1 + 0.61q_v - q_l - q_i)$$

where $q_v$ is water vapor, $q_l$ is cloud liquid, and $q_i$ is cloud ice. Water vapor tends to make air lighter, while condensed water and ice add mass without contributing to gas pressure in the same way.

This equation-of-state step is what connects temperature and moisture back to the dynamics. The chain is:

$$T, q_v, q_l, q_i \Rightarrow T_v \Rightarrow \rho,\alpha \Rightarrow \Phi \Rightarrow -\nabla_p\Phi$$

So temperature and moisture matter for winds because they determine density and geopotential, which determine the pressure-gradient acceleration.

Hydrostatic geopotential diagnosis

It is often convenient to use geopotential rather than height. Define

$$\Phi = gz$$

Near the surface, this is gravitational potential energy per unit mass. Its units are $\mathrm{m^2/s^2}$, so a horizontal gradient of geopotential has units of acceleration.

Starting from hydrostatic balance,

$$\frac{\partial p}{\partial z} = -\rho g$$

we can rewrite it in pressure coordinates. Since

$$\Phi = gz$$

we have

$$d\Phi = g\,dz$$

Hydrostatic balance gives

$$dp = -\rho g\,dz$$

Substituting $d\Phi=g\,dz$,

$$dp = -\rho\,d\Phi$$

$$d\Phi = -\frac{1}{\rho}dp$$

Using $\alpha=1/\rho$, we get

$$\frac{\partial \Phi}{\partial p} = -\alpha$$

Using the equation of state,

$$\alpha = \frac{R_dT_v}{p}$$

$$\frac{\partial \Phi}{\partial p} = -\frac{R_dT_v}{p}$$

This is important because it tells us how to diagnose geopotential from temperature and moisture. Given $T_v$ in a vertical column, we can integrate vertically to get $\Phi$.

In sigma coordinates, where

$$\sigma = \frac{p}{p_s}$$

we have

$$p = \sigma p_s$$

Inside a fixed column, $p_s$ is fixed while integrating over $\sigma$, so

$$d\log p = d\log \sigma$$

Thus the hydrostatic integral has the schematic form

$$\Phi(\sigma) = \Phi_s + \int_{\sigma}^{1} R_dT_v(\sigma')\,d\log\sigma'$$

Here $\Phi_s$ is the surface geopotential, which represents topography:

$$\Phi_s = gh_s$$

Thus mountains enter the dynamics because they change the lower boundary geopotential, and this affects the geopotential field above the surface.

Material derivative

There are two natural perspectives on a fluid. In the Eulerian picture we sit at a fixed point in space and watch fluid flow past, so $\partial_t q$ measures how a field $q$ changes at that point. In the Lagrangian picture we ride with a fluid parcel, and $Dq/Dt$ measures how $q$ changes for that parcel as it moves.

Consider a scalar field $q(t, \mathbf{x})$ and a parcel whose trajectory is $\mathbf{x}(t)$. Differentiating along the trajectory by the chain rule:

$$\frac{Dq}{Dt} = \partial_t q + \frac{\partial q}{\partial x}\dot{x} + \frac{\partial q}{\partial y}\dot{y} + \frac{\partial q}{\partial z}\dot{z} = \partial_t q + \mathbf{v}\cdot\nabla q$$

The term $\mathbf{v}\cdot\nabla q$ is the advection term. It accounts for the fact that even if $q$ is steady in time, a parcel moving into a region of higher or lower $q$ will experience a change. The material derivative is not a new physical law. It is just the chain rule applied along the path of a moving parcel.

On the sphere, with velocity

$$\mathbf{v}=u e_\theta + v_m e_\phi + w\hat r$$

we need to be careful because $u$ and $v_m$ are physical speeds in meters per second, while $\dot\theta$ and $\dot\phi$ are angular speeds in radians per second.

A small physical displacement is

$$d\mathbf{x} = R\sin\phi\,d\theta\,e_\theta + R\,d\phi\,e_\phi + dz\,\hat r$$

Dividing by $dt$ gives

$$\frac{d\mathbf{x}}{dt} = R\sin\phi\,\dot\theta\,e_\theta + R\dot\phi\,e_\phi + \dot z\,\hat r$$

But this is the physical velocity:

$$\frac{d\mathbf{x}}{dt} = \mathbf{v} =u e_\theta + v_m e_\phi + w\hat r$$

Therefore

$$\dot\theta = \frac{u}{R\sin\phi}$$$$\dot\phi = \frac{v_m}{R}$$$$\dot z = w$$

Substituting these into the chain rule gives

$$\frac{Dq}{Dt} = \partial_t q + \frac{u}{R\sin\phi}\frac{\partial q}{\partial \theta} + \frac{v_m}{R}\frac{\partial q}{\partial \phi} + w\frac{\partial q}{\partial z}$$

This is the material derivative in spherical coordinates under the shallow-atmosphere approximation $r\approx R$.

Continuity of mass

Mass cannot be created or destroyed. For a fixed control volume $\mathcal{V}$ with boundary $\partial\mathcal{V}$, the rate of change of enclosed mass equals the negative of the net outward mass flux:

$$\frac{d}{dt}\int_{\mathcal{V}} \rho\,dV = -\oint_{\partial\mathcal{V}} \rho \mathbf{v}\cdot \hat n\,dA$$

The object $\rho\mathbf{v}$ is the mass flux vector. It has units of mass per area per time. The divergence $\nabla\cdot(\rho\mathbf{v})$ therefore has units of mass per volume per time, so it measures the net outward mass flux per unit volume.

Applying the divergence theorem,

$$\oint_{\partial\mathcal{V}} \rho\mathbf{v}\cdot \hat n\,dA = \int_{\mathcal{V}} \nabla\cdot(\rho\mathbf{v})\,dV$$

Thus

$$\int_{\mathcal{V}}\left(\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})\right)dV=0$$

Since this must hold for every fixed control volume,

$$\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0$$

Expanding with the product rule,

$$\nabla\cdot(\rho\mathbf{v}) = \mathbf{v}\cdot\nabla\rho + \rho\nabla\cdot\mathbf{v}$$

$$\frac{\partial \rho}{\partial t} + \mathbf{v}\cdot\nabla\rho + \rho\nabla\cdot\mathbf{v}=0$$

The first two terms are the material derivative of density, so

$$\frac{D\rho}{Dt}+\rho\nabla\cdot\mathbf{v}=0$$

This says that a moving parcel’s density changes if the parcel expands or compresses. If $\nabla\cdot\mathbf{v}>0$, the parcel expands and density decreases. If $\nabla\cdot\mathbf{v}<0$, the parcel compresses and density increases.

Sigma coordinates

For weather models, we often do not use physical height $z$ as the vertical coordinate. Instead, we use a pressure-like vertical coordinate. The simplest version is the sigma coordinate:

$$\sigma = \frac{p}{p_s}$$

where $p_s=p_s(\theta,\phi,t)$ is surface pressure. Thus

$$p=\sigma p_s$$

At the surface, $\sigma=1$. At the top of the atmosphere, $\sigma$ is close to $0$. The point of sigma coordinates is that vertical levels follow the mass of the atmospheric column.

If we change $\sigma$ inside a fixed column, then $p_s$ is fixed during that vertical differential, so

$$dp = p_s\,d\sigma$$

Combining this with hydrostatic balance,

$$dp=-\rho g\,dz$$

we get

$$p_s\,d\sigma = -\rho g\,dz$$

$$d\sigma = -\frac{\rho g}{p_s}\,dz$$

This gives the vertical Jacobian between physical height and sigma coordinate.

The more important point is that sigma coordinates change the material derivative. In sigma coordinates, vertical motion is described by

$$\dot\sigma = \frac{D\sigma}{Dt}$$

and the material derivative becomes

$$\frac{D}{Dt} = \partial_t + \mathbf{U}_h\cdot\nabla_h + \dot\sigma\partial_\sigma$$

where $\mathbf{U}_h$ is the horizontal wind. Thus vertical transport still exists, but it appears through $\dot\sigma\partial_\sigma$ rather than through $w\partial_z$.

Mass continuity in sigma coordinates

In a hydrostatic atmosphere, surface pressure is closely related to total column mass. The column mass per unit horizontal area is approximately

$$\frac{p_s}{g}$$

Thus mass continuity does two things in sigma coordinates. It evolves the column mass, equivalently $p_s$, and it diagnoses the vertical coordinate velocity $\dot\sigma$.

A useful schematic sigma-coordinate continuity equation is

$$\frac{\partial \log p_s}{\partial t} + \mathbf{U}_h\cdot\nabla_h\log p_s + \nabla_h\cdot\mathbf{U}_h + \frac{\partial \dot\sigma}{\partial \sigma} =0$$

The terms have the following meaning. The first term is the local change of column mass. The second term says that air can advect surface-pressure variations horizontally. The third term is horizontal convergence or divergence of the wind. The last term is the vertical divergence of mass flux in sigma coordinates.

If the horizontal flow converges into a column, mass continuity requires either an increase in column mass or vertical motion through sigma levels. Thus, in a hydrostatic model, vertical motion is not solved from vertical acceleration. Instead, it is diagnosed from mass continuity.

This is why we should not say that the model has no vertical motion. It has no independently prognosed vertical momentum equation, but it still has vertical transport through $\dot\sigma$.

Variables of interest

In the atmosphere, there are a select number of quantities that we are concerned with modeling. The hydrostatic primitive-equation models prognose variables such as vorticity, divergence, temperature, water vapor, cloud liquid, cloud ice, and surface pressure. In this writeup, we explicitly do not model vertical acceleration. We use hydrostatic balance to remove the vertical momentum equation and use mass continuity to infer vertical coordinate motion.

Horizontal momentum

Newton’s second law for a unit-mass air parcel in the rotating frame is:

$$\frac{D\mathbf{v}}{Dt} = \underbrace{-2\Omega\times\mathbf{v}}_{\text{Coriolis}} - \underbrace{\frac{1}{\rho}\nabla_z p}_{\text{pressure gradient}} + \underbrace{\mathcal{F}}_{\text{friction/physics}}$$

where $D/Dt$ is the material derivative following the parcel. The centrifugal term has already been absorbed into effective gravity as discussed above.

For the horizontal momentum equation, we only keep the horizontal part of the pressure-gradient acceleration:

$$\mathbf{a}_{\text{pressure},h} = -\frac{1}{\rho}\nabla_z p$$

This term comes directly from pressure as a surface force. Pressure exerts a force $-p\hat n\,dA$ on a parcel boundary. Integrating over the boundary and applying the divergence theorem gives a pressure force per unit volume of $-\nabla p$. Dividing by mass density gives the pressure acceleration:

$$\mathbf{a}_{\text{pressure}} = -\frac{1}{\rho}\nabla p$$

The horizontal part is therefore $-(1/\rho)\nabla_z p$.

It remains to understand why this term can be rewritten using geopotential. Define the geopotential

$$\Phi = gz$$

Now consider a surface of constant pressure $p_0$. At fixed time, this surface can be written as

$$z = Z(x,y,p_0)$$

meaning

$$p(x,y,Z(x,y,p_0))=p_0$$

Differentiating this identity with respect to $x$ and $y$ gives

$$p_x + p_z Z_x = 0$$

and

$$p_y + p_z Z_y = 0$$

Therefore,

$$\nabla_p Z = -\frac{1}{p_z}\nabla_z p$$

Using hydrostatic balance,

$$p_z = -\rho g$$

$$\nabla_p Z = \frac{1}{\rho g}\nabla_z p$$

Multiplying by $g$ gives

$$\nabla_p\Phi = \frac{1}{\rho}\nabla_z p$$

Therefore,

$$-\frac{1}{\rho}\nabla_z p = -\nabla_p\Phi$$

Thus the horizontal pressure-gradient acceleration can be written as a horizontal geopotential gradient at fixed pressure. The momentum equation becomes

$$\frac{D\mathbf{v}}{Dt} = -2\Omega\times\mathbf{v} - \nabla_p\Phi + \mathcal{F}$$

For large-scale horizontal weather dynamics, this is often written for the horizontal wind $\mathbf{U}_h$ as

$$\frac{D\mathbf{U}_h}{Dt} = -f\hat r\times\mathbf{U}_h - \nabla_p\Phi + \mathcal{F}_h$$

where $f=2\Omega\cos\phi$ in our colatitude convention.

Pressure-gradient force in sigma coordinates

The pressure-coordinate form is clean:

$$\mathbf{a}_{\text{pressure},h}=-\nabla_p\Phi$$

But sigma coordinates are not exactly pressure coordinates. Since

$$p=\sigma p_s(\theta,\phi,t)$$

holding $\sigma$ fixed is not the same thing as holding $p$ fixed. If $p_s$ changes horizontally, then moving horizontally at fixed $\sigma$ changes pressure.

By the chain rule,

$$\nabla_\sigma\Phi = \nabla_p\Phi + \frac{\partial\Phi}{\partial p}\nabla_\sigma p$$

At fixed $\sigma$,

$$\nabla_\sigma p = \sigma\nabla_h p_s = p\nabla_h\log p_s$$

Using hydrostatic balance in pressure coordinates,

$$\frac{\partial\Phi}{\partial p} = -\frac{R_dT_v}{p}$$

$$\frac{\partial\Phi}{\partial p}\nabla_\sigma p = -R_dT_v\nabla_h\log p_s$$

Therefore,

$$\nabla_\sigma\Phi = \nabla_p\Phi - R_dT_v\nabla_h\log p_s$$

Rearranging,

$$\nabla_p\Phi = \nabla_\sigma\Phi + R_dT_v\nabla_h\log p_s$$

and hence

$$-\nabla_p\Phi = -\nabla_\sigma\Phi - R_dT_v\nabla_h\log p_s$$

This correction appears because fixed sigma surfaces move with surface pressure. This is one of the main places where sigma coordinates affect the primitive equations.

Thermodynamics

Temperature is not just another passive tracer. It is transported by the flow, but it also changes when air expands, compresses, radiatively heats or cools, and exchanges latent heat with water species.

A schematic thermodynamic equation is

$$\frac{DT}{Dt} = \text{adiabatic heating/cooling} + \text{diabatic heating/cooling}$$

A common dry pressure-coordinate form is

$$\frac{DT}{Dt} - \frac{\kappa T}{p}\frac{Dp}{Dt} = \frac{Q}{c_p}$$

where

$$\kappa=\frac{R_d}{c_p}$$

The term involving $Dp/Dt$ represents adiabatic heating and cooling. If air rises, it usually moves to lower pressure, expands, and cools. If air sinks, it moves to higher pressure, compresses, and warms.

The right-hand side $Q/c_p$ represents diabatic heating and cooling. This includes radiation, latent heating from condensation, evaporative cooling, and other unresolved physical processes.

Thus temperature affects the atmosphere in two ways. First, it is carried by the flow. Second, it changes density and geopotential through the equation of state and hydrostatic balance.

Moisture and tracer equations

For atmospheric water species, the model usually stores mixing ratios. A mixing ratio is mass of tracer per mass of air. For example,

$$q_v = \frac{\text{mass of water vapor}}{\text{mass of air}}$$

Similarly, we can define cloud liquid and cloud ice mixing ratios:

$$q_l = \frac{\text{mass of cloud liquid}}{\text{mass of air}}$$$$q_i = \frac{\text{mass of cloud ice}}{\text{mass of air}}$$

For a passive tracer mixing ratio, the equation is simply

$$\frac{Dq}{Dt}=0$$

meaning the parcel carries the tracer ratio with it. For active water species, we include source and sink terms:

$$\frac{Dq_v}{Dt}=S_v$$$$\frac{Dq_l}{Dt}=S_l$$$$\frac{Dq_i}{Dt}=S_i$$

The source terms exchange mass between phases. Condensation decreases water vapor and increases cloud liquid. Evaporation does the reverse. Freezing moves liquid into ice. Melting moves ice into liquid. Precipitation removes condensate from the modeled air column.

Latent heating couples the water equations back to temperature. When vapor condenses, latent heat is released and temperature increases. When liquid evaporates, heat is absorbed and temperature decreases. Thus moisture tracers are not purely passive. They affect thermodynamics and, through virtual temperature, pressure-gradient forces.

Vorticity and divergence

Although the horizontal momentum equation is naturally written in terms of the horizontal wind $\mathbf{U}_h$, these spectral primitive-equation models often do not store raw wind components directly. Instead, they store vorticity and divergence:

$$\zeta = \hat r\cdot(\nabla_h\times\mathbf{U}_h)$$

and

$$\delta = \nabla_h\cdot\mathbf{U}_h$$

This is valid because of the Helmholtz decomposition on the sphere. A sufficiently smooth horizontal vector field can be written as

$$\mathbf{U}_h = \nabla_h\chi + \hat r\times\nabla_h\psi$$

where $\chi$ is a velocity potential and $\psi$ is a streamfunction. Taking divergence gives

$$\delta = \Delta_h\chi$$

and taking the radial component of curl gives

$$\zeta = \Delta_h\psi$$

up to sign conventions. Thus, knowing $\delta$ and $\zeta$ lets us solve two Poisson equations for $\chi$ and $\psi$, and then reconstruct $\mathbf{U}_h$.

This is especially convenient in spherical harmonics, because the spherical harmonics diagonalize the spherical Laplacian:

$$\Delta_hY_l^m = -\frac{l(l+1)}{R^2}Y_l^m$$

So solving $\Delta_h\chi=\delta$ or $\Delta_h\psi=\zeta$ is just division by an eigenvalue in spectral space.

The Coriolis term also naturally combines with vorticity. We define the absolute vorticity

$$\eta = \zeta + f$$

where $f=2\Omega\cos\phi$ in the colatitude convention. Here $\zeta$ is spin relative to Earth, while $f$ is planetary spin from Earth’s rotation.