Intro to weather modeling
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, we just need the net flux through each pair of opposite faces. The flux in through the face at $\theta$ is $V_\theta(\phi, \theta) \cdot R \Delta \phi$ and the flux out at $\theta + \Delta\theta$ is the same thing evaluated at $\theta + \Delta\theta$, so the net flux through the $\theta$-faces is:
$$(V_\theta(\phi, \theta + \Delta \theta) - V_\theta(\phi, \theta)) \cdot R \Delta \phi \approx \frac{\partial V_\theta}{\partial \theta} R \Delta \phi \Delta \theta$$The $\phi$-faces work the same way, except the face lengths $R \sin\phi \, \Delta\theta$ themselves depend on $\phi$, so the $\sin \phi$ rides along inside the derivative:
$$(V_\phi(\phi + \Delta \phi, \theta)\sin(\phi + \Delta \phi) - V_\phi(\phi, \theta)\sin(\phi)) \cdot R \Delta \theta \approx R\frac{\partial \sin\phi V_\phi}{\partial \phi} \Delta \phi \Delta \theta$$Adding the two and dividing by the cell area $R^2 \sin \phi \, \Delta \theta \, \Delta \phi$ gives the divergence:
$$\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 try the separable ansatz $f(\phi, \theta) = g(\phi) h(\theta)$. The partial derivatives factor nicely, and substituting them in:
$$\frac{1}{\sin^2 \phi} g(\phi)h''(\theta) + \frac{\cos \phi}{\sin \phi} g'(\phi) h(\theta) + g''(\phi) h(\theta) = c R^2 g(\phi) h(\theta)$$Dividing by $g(\phi) h(\theta)$ and rearranging:
$$\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 solutions to this associated Legendre equation are the associated Legendre polynomials. Putting the two pieces back together, the spherical harmonics satisfy:
$$\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)$$so the velocity 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'$$In other words, pressure at a given height is just the weight per unit area of the air sitting above it. 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. In terms of the specific volume $\alpha = 1/\rho$, this reads:
$$\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, \quad \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.
Let’s rewrite hydrostatic balance in pressure coordinates. Since $d\Phi = g\,dz$ and hydrostatic balance says $dp = -\rho g\,dz$, we have $dp = -\rho\,d\Phi$. Flipping this around with $\alpha = 1/\rho$ and the equation of state:
$$\frac{\partial \Phi}{\partial p} = -\alpha = -\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$$So 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$$This is how mountains enter the dynamics: they change the lower boundary geopotential, which shifts the whole geopotential field above.
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$ and matching components against $\mathbf{v} = u e_\theta + v_m e_\phi + w\hat r$:
$$\dot\theta = \frac{u}{R\sin\phi}, \quad \dot\phi = \frac{v_m}{R}, \quad \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$$Here $\rho\mathbf{v}$ is the mass flux (mass per area per time). Applying the divergence theorem to the right side, and noting that the identity must hold for every volume we could have picked, the integrand itself has to vanish:
$$\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0$$Expanding the divergence with the product rule and recognizing the material derivative in the first two terms:
$$\frac{D\rho}{Dt}+\rho\nabla\cdot\mathbf{v}=0$$The interpretation is nice: a moving parcel’s density only changes if the parcel expands ($\nabla\cdot\mathbf{v}>0$, density drops) or compresses ($\nabla\cdot\mathbf{v}<0$, density rises).
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.
Inside a fixed column, $p_s$ doesn’t change as we move vertically, so $dp = p_s\,d\sigma$. Combining with hydrostatic balance $dp=-\rho g\,dz$:
$$d\sigma = -\frac{\rho g}{p_s}\,dz$$which is just the Jacobian between physical height and sigma.
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}$$So in sigma coordinates, mass continuity does two jobs at once. It evolves the column mass (equivalently $p_s$), and it diagnoses the vertical coordinate velocity $\dot\sigma$. Schematically:
$$\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$$Reading left to right: local change of column mass, horizontal advection of surface pressure, horizontal convergence or divergence of the wind, and vertical mass flux through sigma levels.
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.
So the model very much has vertical motion. There is just no prognostic vertical momentum equation, and vertical transport happens through $\dot\sigma$ instead.
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$.
The neat trick is that this term can be rewritten using 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$, which packaged together says:
$$\nabla_p Z = -\frac{1}{p_z}\nabla_z p = \frac{1}{\rho g}\nabla_z p$$where the last step used hydrostatic balance $p_z = -\rho g$. Multiplying by $g$:
$$-\frac{1}{\rho}\nabla_z p = -\nabla_p\Phi$$So the horizontal pressure-gradient acceleration is just a 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$ we have $\nabla_\sigma p = \sigma\nabla_h p_s = p\nabla_h\log p_s$, and hydrostatic balance gives $\partial\Phi/\partial p = -R_dT_v/p$. Putting the pieces together:
$$-\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.
So temperature pulls double duty. It gets carried around by the flow, and it feeds back into the winds by changing 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, \quad \frac{Dq_l}{Dt}=S_l, \quad \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$$Why is this enough to recover the wind? The Helmholtz decomposition on the sphere says any sufficiently smooth horizontal vector field splits 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 the divergence and the radial component of the curl (up to sign conventions):
$$\delta = \Delta_h\chi, \quad \zeta = \Delta_h\psi$$So given $\delta$ and $\zeta$, we solve two Poisson equations and reconstruct $\mathbf{U}_h$. This is where all our spherical harmonics work pays off, since they diagonalize the Laplacian:
$$\Delta_hY_l^m = -\frac{l(l+1)}{R^2}Y_l^m$$Solving $\Delta_h\chi=\delta$ or $\Delta_h\psi=\zeta$ in spectral space is literally just dividing by an eigenvalue.
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.
Comments