The Navier–Stokes Equations in the Atmosphere
Every gust of wind, every rotating cyclone and every thunderstorm updraft obeys the same equations that govern water flowing through a pipe. Here's how the Navier–Stokes equations are adapted, simplified and discretised to describe the motion of the air around us.
1. The atmosphere as a fluid
The air around us is a compressible, viscous, rotating fluid roughly 10 km thick sitting on a sphere spinning once every 24 hours. Despite this complexity, it obeys exactly the same conservation laws as any other Newtonian fluid: conservation of momentum (Newton's second law applied to a fluid parcel), conservation of mass (continuity), and conservation of energy. Together these are the Navier–Stokes equations.
What makes atmospheric flow special is scale. Weather systems span from millimetre-scale turbulent eddies to continent-sized jet streams, and the Earth's rotation makes the Coriolis effect dominant at synoptic scales (hundreds to thousands of kilometres) — something you never need to consider when simulating water in a cup.
Air is made of discrete molecules, but at the scales meteorology cares about (millimetres and up) it behaves as a continuous fluid — the same assumption used in CFD for aircraft wings. This is what allows us to write differential equations for velocity, pressure and density fields instead of tracking individual molecules.
2. The full Navier–Stokes system
For a compressible fluid, the momentum equation (per unit mass) reads:
where u is the 3D velocity vector, ρ is air density, p is pressure, g is gravitational acceleration and μ is the (molecular) dynamic viscosity. The left-hand side is the material derivative — acceleration following a moving parcel of air, made of a local (∂u/∂t) and an advective (u·∇u) term.
Alongside momentum, we track mass and thermal energy:
Thermodynamic energy cp (∂T/∂t + u·∇T) = (1/ρ)(∂p/∂t + u·∇p) + Q
Q represents diabatic heating: solar radiation absorption, latent heat release from
condensation, infrared cooling. This coupling between temperature and pressure through the
ideal gas law p = ρRT is what makes the atmosphere a compressible,
thermodynamically active fluid rather than a simple incompressible liquid.
3. A rotating frame: Coriolis and centrifugal terms
Weather models are written in a reference frame that co-rotates with the Earth, which is non-inertial. Newton's laws only hold in inertial frames, so two fictitious "apparent" forces appear when we rewrite the momentum equation for an Earth-fixed observer:
The term −2Ω × u is the Coriolis acceleration — it deflects moving air to the right in the Northern Hemisphere and to the left in the Southern Hemisphere (the subject of our companion article on the Coriolis effect). The term −Ω × (Ω × r) is the centrifugal acceleration, which is usually absorbed into an "effective gravity" g since it doesn't change direction with wind speed.
Ω is Earth's angular velocity, 7.292 × 10⁻⁵ rad/s (one revolution per sidereal day). Its local vertical component, f = 2Ω sin φ (φ = latitude), is called the Coriolis parameter and appears throughout large-scale dynamics.
4. Hydrostatic balance and the primitive equations
Solving the full 3D compressible Navier–Stokes equations for the whole atmosphere would require resolving sound waves travelling at ~340 m/s, forcing time-steps of a fraction of a second even for continental-scale models. In practice, operational forecast models use the hydrostatic approximation: the vertical momentum equation is replaced by a simple balance between the pressure-gradient force and gravity.
This filters out vertically propagating sound waves while retaining everything relevant to weather — pressure systems, fronts, jet streams. Combined with continuity and thermodynamics, this simplified set is known as the primitive equations, the backbone of models like ECMWF's IFS and NOAA's GFS. High-resolution convection-permitting models (grid spacing below ~4 km) instead solve the full non-hydrostatic equations, because at that scale vertical accelerations inside thunderstorm updrafts are no longer negligible.
Far from the equator and away from strong curvature, the Coriolis force nearly balances the pressure-gradient force: fΩv ≈ −(1/ρ)∇p. The resulting geostrophic wind blows parallel to isobars rather than from high to low pressure — the reason weather-map wind arrows follow the contour lines instead of crossing them.
5. Turbulence closure
Atmospheric flow is turbulent across an enormous range of scales, and no computer can resolve every eddy down to the millimetre dissipation scale. Instead, models split variables into a resolved mean and an unresolved fluctuation, u = ū + u′, and average the Navier–Stokes equations (Reynolds averaging). This produces the Reynolds-Averaged Navier–Stokes (RANS) equations, with an extra unclosed term — the Reynolds stress tensor:
Weather and climate models parameterize τij using boundary-layer schemes (K-theory eddy diffusivity, TKE closures) or, for convection, entirely separate "cumulus parameterization" schemes that represent sub-grid thunderstorms statistically rather than resolving individual updrafts.
6. From continuous equations to a grid
Numerical Weather Prediction (NWP, covered in depth in our dedicated article) discretises the primitive equations onto a 3D grid — either a latitude/longitude/pressure grid, an icosahedral grid, or in spectral models, a truncated series of spherical harmonics. Time is advanced with semi-implicit semi-Lagrangian schemes that allow time-steps of 10-15 minutes even at ~10 km resolution, because the fastest (sound and gravity) waves are treated implicitly while advection is treated along back-trajectories.
7. Pseudocode: one time-step
A drastically simplified single time-step of a shallow, hydrostatic atmospheric model:
function stepAtmosphere(grid, dt):
// 1. Pressure-gradient force from geopotential height
computePressureGradient(grid)
// 2. Coriolis acceleration per grid cell
for each cell (i, j):
f = 2 * OMEGA * sin(grid.lat[i,j])
a_coriolis.u = f * grid.v[i,j]
a_coriolis.v = -f * grid.u[i,j]
// 3. Sub-grid turbulence + convection parameterization
applyBoundaryLayerMixing(grid)
applyCumulusScheme(grid)
// 4. Semi-Lagrangian advection along back-trajectories
advectSemiLagrangian(grid, dt)
// 5. Semi-implicit solve for pressure/geopotential
solveImplicitGravityWaves(grid, dt)
// 6. Update thermodynamic fields (T, humidity)
updateThermodynamics(grid, dt)🌪️ See rotating fluid dynamics in action
The Tornado simulation shows a rotating vortex of particles obeying the same momentum equation described above — minus the Coriolis term, since a tornado is far too small and short-lived for Earth's rotation to matter.
Open simulation →