Numerical Weather Prediction: How Computers Forecast the Weather

1. From equations to a forecast

Given a complete description of the atmosphere's state right now — every wind speed, temperature and humidity value everywhere — the primitive equations from our Navier–Stokes article tell you, in principle, exactly how that state evolves forward in time. Numerical Weather Prediction (NWP) is the practice of turning that "in principle" into an actual number on a screen: representing the atmosphere as billions of grid values, and marching them forward with a computer using discretised versions of those equations.

Modern global models — ECMWF's Integrated Forecast System (IFS), NOAA's Global Forecast System (GFS), the UK Met Office's Unified Model — run at ~9-13 km horizontal grid spacing with 100+ vertical levels, producing a 10-day forecast in under two hours on dedicated supercomputers with tens of thousands of CPU cores.

2. Data assimilation: building the initial state

A forecast is only as good as its starting point. Data assimilation combines millions of heterogeneous observations — satellite radiances, aircraft reports, weather balloons, ships, buoys, ground stations, radar — with the model's own short-range forecast from six hours earlier (the "background") to produce the best possible estimate of the current atmospheric state.

x_b is the background forecast, B its error covariance, y_t the observations at time t, H the "observation operator" mapping model state to what an instrument would measure, and R the observation error covariance. Four-Dimensional Variational assimilation (4D-Var) finds the state x that best fits both the background and every observation across a 6-12 hour window simultaneously — an enormous optimisation problem in millions of dimensions, itself run on supercomputers.

3. Grids, resolution and parameterization

No grid can resolve everything. Anything smaller than roughly twice the grid spacing — individual thunderstorm cells, small cumulus clouds, turbulent boundary-layer eddies — must be represented statistically through parameterization schemes rather than explicitly simulated:

• Cumulus convection: represents the net heating/moistening effect of unresolved thunderstorms using a statistical ensemble of idealised updraft/downdraft plumes.
• Radiation: solves simplified radiative transfer for shortwave (solar) and longwave (infrared) fluxes through cloud and gas layers, updated every 1-3 hours to save compute.
• Microphysics: tracks bulk properties (mixing ratio, number concentration) of cloud water, rain, ice, snow and graupel rather than individual droplets.
• Boundary layer: mixes momentum, heat and moisture vertically in the lowest ~1-2 km using eddy-diffusivity or TKE-based turbulence closures.

Higher resolution reduces the burden on parameterization — km-scale "convection-permitting" models explicitly resolve individual thunderstorm updrafts instead of parameterizing them, at the cost of roughly 1000× more computation than a 10 km global model.

4. Time integration and the CFL condition

Explicit time-stepping schemes are constrained by the Courant–Friedrichs–Lewy (CFL) condition: information must not travel more than one grid cell per time-step, or the scheme becomes numerically unstable and blows up.

For a 10 km grid and 300 km/h jet-stream winds, a naive explicit scheme would need Δt ≈ 2 minutes. Operational models instead use semi-Lagrangian advection (tracking where a parcel came from, rather than pushing values through fixed cells) combined with a semi-implicit treatment of the fast gravity/sound waves, allowing Δt of 10-15 minutes without violating stability — a roughly 5-8× speed-up that made operational global forecasting computationally feasible in the 1990s.

5. Ensemble forecasting and chaos

The atmosphere is a chaotic system: Edward Lorenz's famous 1963 discovery showed that tiny differences in initial conditions — his "butterfly effect" — grow exponentially and eventually produce completely different forecasts. Since the true initial state is never known perfectly (observation errors, assimilation approximations), a single deterministic forecast run ("the model run") cannot tell you how confident to be.

Ensemble forecasting solves this by running the model 20-50 times from slightly perturbed initial conditions (and sometimes perturbed physics parameterizations too), producing a spread of possible outcomes instead of one number. A tight ensemble spread means high confidence; a wide spread — common beyond 7-10 days — means the forecast is genuinely uncertain, and probabilistic products (e.g. "70% chance of rain") come directly from counting how many ensemble members produce rain at that point.

6. The two-week predictability limit

Lorenz's work established that even with a perfect model and near-perfect initial conditions, chaotic error growth imposes a hard theoretical limit on deterministic weather predictability of roughly two weeks — beyond that, forecast skill converges to climatology (long-term average conditions) regardless of computing power. This is fundamentally different from climate prediction, which forecasts statistical distributions decades ahead rather than the exact state on a specific day.

7. Pseudocode: a forecast cycle

function runForecastCycle(observations, previousBackground):

  // 1. Data assimilation: blend background + observations
  analysis = fourDVar(previousBackground, observations)

  // 2. Generate ensemble of perturbed initial states
  ensemble = generatePerturbations(analysis, members=50)

  for each member in ensemble:
    for step in 1..N_STEPS:
      // 3. Advance the primitive equations one time-step
      member = stepAtmosphere(member, dt=900) // 15 min

      // 4. Apply sub-grid parameterizations
      member = applyCumulusScheme(member)
      member = applyRadiation(member)
      member = applyMicrophysics(member)

  // 5. Produce probabilistic forecast from ensemble spread
  return summarizeEnsemble(ensemble)