Numerical weather modeling

How numerical weather prediction (NWP) turns observations into forecasts: the equations, grids, time stepping, physics parameterisations, data assimilation, ensembles and verification, with the operational systems as of September 2026. Basic mechanics and fluids are in physics fundamentals; constants in constants, units and conversions.

What NWP is

NWP solves an initial-value problem: estimate the current state of the atmosphere, then step the fluid equations forward on a computer. Global centers run the whole cycle every 6 hours (00, 06, 12, 18 UTC); rapid-update regional systems run it every hour.

observations (satellite, sonde, aircraft, radar, surface)
    |  quality control, bias correction
    v
data assimilation  <---  background: previous short forecast
    |  analysis (best estimate of the state now)
    v
model integration: dynamics + physics, often an ensemble
    |  fields on the model grid
    v
post-processing: interpolate, bias-correct, calibrate
    v
products: GRIB files, charts, probabilities, warnings

Forecast error has two sources: the initial state (never known exactly) and the model (discretisation and parameterisation errors). Chaos makes the first grow; the second adds bias.

Governing equations

The primitive equations: Newton's second law on a rotating sphere, mass conservation, the first law of thermodynamics, water conservation and the ideal gas law, with the hydrostatic approximation replacing the vertical momentum equation. With the material derivative $\dfrac{D}{Dt} = \dfrac{\partial}{\partial t} + \vec v \cdot \nabla$:

In pressure coordinates the same set is tidier, because density disappears and continuity becomes diagnostic ($\omega = Dp/Dt$ is the vertical velocity in Pa/s):

Balances and scales

The hydrostatic approximation filters sound waves traveling vertically and is accurate while horizontal scales far exceed vertical ones. Below roughly 10 km grid spacing, and certainly at convection-permitting 1–4 km, models solve the full non-hydrostatic equations (ICON, UM, FV3, WRF). IFS at 9 km still runs hydrostatic.

Vertical coordinates

IFS uses 137 hybrid sigma-pressure levels from the surface to 0.01 hPa (about 80 km); GFS uses 127. Levels are packed near the surface to resolve the boundary layer.

Grids and discretisation

Centered second-order differences, with error $O(\Delta x^2)$:

Horizontal grids

Arakawa staggering

Where variables sit on the grid changes how well gravity waves and geostrophic adjustment are represented. $h$ is a mass or pressure point.

Arakawa C: one cell
 
      +------ v ------+
      |               |
      u       h       u
      |               |
      +------ v ------+

Time stepping and the CFL condition

The Courant–Friedrichs–Lewy condition: an explicit scheme is unstable if information crosses more than a fixed fraction of a cell per step.

$C_\text{max}$ is 1 for first-order upwind in 1D and smaller in 2D and 3D. $c_\text{max}$ is the fastest signal the scheme treats explicitly: wind (up to about 100 m/s in jets), gravity waves ($\sqrt{gH}$, about 300 m/s for the external mode) or sound (about 340 m/s).

Leapfrog needs a Robert–Asselin filter to damp its computational mode. First-order upwind is stable for $0 \lt C \le 1$ but diffusive: sharp features smear out.

Toy: 1D linear advection, upwind

A Gaussian blob carried once around a periodic 1000 km domain at 20 m/s. The scheme is exact in mass but loses 38% of the peak to numerical diffusion.

import numpy as np
from numpy.typing import NDArray
 
Field = NDArray[np.float64]
 
def upwind(q: Field, c: float) -> Field:
    # first-order upwind for u > 0, periodic domain
    return q - c * (q - np.roll(q, 1))
 
L, nx, u = 1000e3, 100, 20.0        # m, cells, m/s
dx = L / nx                          # 10 km
dt = 400.0                           # s
c = u * dt / dx                      # Courant number
assert 0 < c <= 1, f"CFL violated: C = {c:.2f}"
 
x = np.arange(nx) * dx
q = np.exp(-(((x - 250e3) / 50e3) ** 2))  # blob
m0 = q.sum()
steps = round(L / (u * dt))          # one full lap
for _ in range(steps):
    q = upwind(q, c)
 
print(f"dx={dx / 1e3:.0f} km  dt={dt:.0f} s  C={c:.2f}")
print(f"steps={steps}  peak 1.000 -> {q.max():.3f}")
print(f"mass ratio={q.sum() / m0:.6f}")

uv run --with numpy python advect.py

dx=10 km  dt=400 s  C=0.80
steps=125  peak 1.000 -> 0.620
mass ratio=1.000000

With dt = 600.0 the check fails with AssertionError: CFL violated: C = 1.20; remove the assert and the solution blows up within a few laps.

Resolution and cost

Halving $\Delta x$ quadruples the horizontal points and, through CFL, halves $\Delta t$:

Effective resolution is 4–8 $\Delta x$: features smaller than that are damped by the numerics. Data volume grows as $\Delta x^{-2}$ per level per output step.

Parameterisations

Processes smaller than the grid, or too costly to simulate, are represented by schemes that feed tendencies of $T$, $q$, $\vec v$ and condensate back into the dynamics.

Stochastic physics (for example SPPT, stochastically perturbed parameterisation tendencies) represents model uncertainty inside ensembles.

Data assimilation

Find the analysis $\mathbf{x}_a$ that best fits both the background $\mathbf{x}_b$ (a short forecast) and the observations $\mathbf{y}$, weighted by their error covariances.

3D-Var minimizes one cost function at a single analysis time:

4D-Var fits all observations in a window (6–12 h) with the model as a constraint, adjusting the initial state $\mathbf{x}_0$:

For linear $H$ the minimum is the Kalman (best linear unbiased) update:

$\mathbf{P}_e = \dfrac{1}{N-1}\mathbf{X}'\mathbf{X}'^\mathsf{T}$ is the ensemble covariance from perturbations $\mathbf{X}'$; $\mathbf{L}$ localizes it (Schur product $\circ$) to kill spurious long-range correlations from a small ensemble.

Ensembles and chaos

Small initial errors grow until forecasts are no better than climatology. Lorenz (1963) showed it with three equations ($\sigma = 10$, $\rho = 28$, $\beta = 8/3$):

Synoptic-scale errors double in about 1–2 days; small-scale (convective) errors in hours. The practical limit for day-to-day weather is around two weeks.

An ensemble runs the model many times from perturbed initial states (singular vectors, ensemble data assimilation members) with perturbed physics (stochastic schemes), and samples the forecast probability distribution.

A well-calibrated ensemble has spread about equal to the RMSE of its mean, and a flat rank histogram. Too little spread (common) means overconfidence.

Global and regional models

Boundary data enter one way from the parent model, so regional forecasts inherit the global model's large-scale errors after a day or two. Nesting steps down in stages (for example ICON 13 km → 6.5 km Europe nest). A new run needs spin-up time for clouds and precipitation.

Operational models

Status on 25 September 2026.

Machine-learning models are trained on reanalysis (ERA5) and fine-tuned on operational analyses. They run in minutes on GPUs, but still depend on a physics-based data assimilation system for their starting state and tend to smooth extremes at longer lead times.

Verification

Compare forecasts $f_i$ with observations or analyses $o_i$, $N$ pairs, climatology $c_i$.

High-resolution precipitation suffers the double penalty (a storm slightly displaced counts as a miss and a false alarm), so use neighborhood scores such as the fractions skill score. Probability forecasts also need reliability diagrams and ROC curves. WeatherBench 2 standardizes these scores for comparing ML and physics models.

Worked examples

Coriolis parameter at a latitude

$f = 2\Omega\sin\varphi$ with $\Omega = 7.292 \times 10^{-5}$ rad/s. At $45^\circ$N:

$f = 0$ at the equator, which is why geostrophic balance fails there.

Geostrophic wind from a pressure gradient

Isobars 4 hPa apart every 300 km at $45^\circ$N, $\rho = 1.2$ kg/m³:

On a 500 hPa chart, height contours 60 m apart over 500 km: $V_g = \dfrac{g}{f}\dfrac{\Delta z}{\Delta n} = \dfrac{9.81 \times 60}{(1.031 \times 10^{-4})(5 \times 10^5)} = 11.4$ m/s, blowing parallel to the contours with low pressure on the left (Northern Hemisphere).

CFL time step for a 9 km grid

Jet-stream wind 100 m/s, $\Delta x = 9$ km, $C_\text{max} = 1$:

Explicit treatment of sound waves would force 20 s steps. Semi-implicit treatment removes the fast waves and semi-Lagrangian advection removes the wind limit, so IFS runs at 450 s: an advective Courant number of $100 \times 450 / 9000 = 5$.

Cost of higher resolution

Going from 9 km to 4.5 km: $2^2$ more columns $\times$ 2 more steps $= 8\times$ the cost ($16\times$ with twice the levels). From 9 km to 1.4 km: $(9/1.4)^3 \approx 266\times$.

Columns at 9 km: $\dfrac{4\pi(6.371 \times 10^6)^2}{(9000)^2} \approx 6.3 \times 10^6$; with 137 levels about $8.6 \times 10^8$ grid points, each carrying about ten prognostic variables.

Height of the 500 hPa surface

Isothermal layer at 250 K: $H = R_dT/g = 287 \times 250 / 9.81 = 7.3$ km. From a 1013 hPa surface:

so the 500 hPa surface sits near 5.5 km, higher in warm air and lower in cold air.

RMSE, bias and a skill score

Four 2 m temperature forecasts $f = (12, 15, 9, 20)$ °C against observations $o = (10, 16, 9, 17)$ °C, errors $(2, -1, 0, 3)$:

If persistence (tomorrow = today) has $MSE_\text{ref} = 7.0$ K², the skill score is $1 - 3.5/7.0 = 0.50$: half the persistence error variance removed.

References

• ECMWF: IFS documentation (opens in a new tab): dynamics, physics and data assimilation of the IFS, cycle by cycle
• ECMWF: IFS Cycle 50r1 and AIFS v2 go live (opens in a new tab): the May 2026 upgrade
• ECMWF: Implementation of IFS Cycle 49r1 (opens in a new tab): TCo1279, 137 levels and the HRES/ENS control unification
• ECMWF: ensemble AI forecasts become operational (opens in a new tab): AIFS ENS
• ECMWF: AIFS Single 2.0 model card (opens in a new tab): architecture, resolution and training data of the ML model
• ECMWF: Learning (opens in a new tab): training courses and e-learning on NWP, data assimilation and ensembles
• UCAR COMET MetEd (opens in a new tab): free modules on NWP, model physics and ensembles
• NOAA EMC: Global Forecast System (opens in a new tab): GFS configuration and version history
• NWS Public Information Statement 26-29 (opens in a new tab): the proposed GFS v17 upgrade
• NWS Service Change Notice 26-48 (opens in a new tab): RRFS and REFS implementation
• NOAA GSL: HRRR (opens in a new tab): the High-Resolution Rapid Refresh
• DWD: ICON model (opens in a new tab): icosahedral grid and model design
• Met Office: Unified Model (opens in a new tab): UM configurations
• Met Office: LFRic (opens in a new tab): the cubed-sphere successor to the UM
• Google: WeatherNext models (opens in a new tab): GraphCast, GenCast and WeatherNext lineage
• WeatherBench 2 (opens in a new tab): benchmark scores for ML and physics models
• Lorenz (1963): Deterministic nonperiodic flow (opens in a new tab): the paper behind chaos in weather
• Kalnay: Atmospheric Modeling, Data Assimilation and Predictability (opens in a new tab): the standard NWP and data assimilation text
• Warner: Numerical Weather and Climate Prediction (opens in a new tab): practical guide to using models, verification and ensembles
• Durran: Numerical Methods for Fluid Dynamics (opens in a new tab): finite difference, spectral and semi-Lagrangian methods in depth