← Back to Part 7

7.3 Numerical Weather Prediction

Numerical Weather Prediction (NWP) solves the governing fluid dynamical and thermodynamic equations on computational grids, integrating them forward in time to produce deterministic and probabilistic forecasts. From Richardson's hand-computed forecast in 1922 to today's global ensemble systems on petascale supercomputers, NWP represents one of the great triumphs of applied physics and computational science. Modern operational models achieve useful skill out to 10-14 days for synoptic-scale features, with convection-permitting models providing hourly guidance at 1-3 km resolution.

Historical Development

The idea of forecasting weather by solving equations of motion dates to Bjerknes (1904). Key milestones:

Richardson (1922)

First numerical forecast attempt by hand -- wildly wrong (145 hPa/6hr vs ~1 hPa observed) due to unbalanced initial conditions. Envisioned a "forecast factory" with 64,000 human computers.

Charney, Fjortoft, von Neumann (1950)

First successful computer forecast on ENIAC using the barotropic vorticity equation. Produced a 24-hour 500 hPa forecast for North America. Launched the era of operational NWP.

The Primitive Equations in Detail

The primitive equations are the fundamental governing equations of atmospheric motion under the hydrostatic approximation, expressed in pressure (or sigma/hybrid) coordinates on a rotating sphere:

Horizontal Momentum

$$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + \omega\frac{\partial u}{\partial p} = fv - \frac{\partial \Phi}{\partial x} + F_x$$

$$\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + \omega\frac{\partial v}{\partial p} = -fu - \frac{\partial \Phi}{\partial y} + F_y$$

where $f = 2\Omega\sin\phi$, $\Phi = gz$ is geopotential, $\omega = dp/dt$ is vertical velocity in pressure coordinates

Hydrostatic, Continuity, Thermodynamic, and Moisture Equations

$$\frac{\partial \Phi}{\partial p} = -\frac{R_d T_v}{p}, \quad \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial \omega}{\partial p} = 0$$

$$\frac{\partial T}{\partial t} + \mathbf{v}\cdot\nabla T + \left(\frac{R_d T}{c_p p}\right)\omega = \frac{\dot{Q}}{c_p}$$

$$\frac{\partial q}{\partial t} + \mathbf{v}\cdot\nabla q + \omega\frac{\partial q}{\partial p} = S_q$$

Modern models carry multiple moisture species: $q_v, q_c, q_r, q_i, q_s, q_g$ (vapor, cloud water, rain, ice, snow, graupel)

Sigma Coordinates: $\sigma = p/p_s$ follows terrain. Hybrid: $p(\eta) = A(\eta)p_0 + B(\eta)p_s$ transitions from terrain-following near surface to pure pressure aloft. ECMWF IFS uses 137 hybrid levels to 0.01 hPa (~80 km).

Spectral Methods and Grid Structures

Global spectral models expand fields in spherical harmonics, the eigenfunctions of the Laplacian on the sphere:

$$\psi(\lambda, \phi) = \sum_{n=0}^{N} \sum_{m=-n}^{n} \hat{\psi}_n^m \, P_n^{|m|}(\sin\phi)\, e^{im\lambda}$$

Triangular truncation $|m| \leq n \leq N$. T1279 = ~9 km (ECMWF HRES). Grid spacing $\approx 40000/(2N)$ km.

Gaussian Grids

Transform between spectral and grid-point space via FFT (longitude) and Legendre transforms (latitude). Reduced Gaussian grids have fewer points near poles. Cubic truncation improves efficiency.

Cubed-Sphere (FV3)

6 faces of gnomonic projection. Avoids polar singularity. Locally conservative finite-volume. Used by GFS/GFDL. Quasi-uniform resolution globally.

Icosahedral (ICON, MPAS)

Triangular (ICON) or Voronoi (MPAS) mesh on icosahedron. Supports variable resolution through mesh refinement. Non-hydrostatic dynamics. C-grid staggering.

Arakawa C-Grid: Mass at cell centers, u on east/west faces, v on north/south faces. Best for inertia-gravity waves. Used by WRF, MPAS, ICON, and most modern NWP models.

Time Integration and CFL Condition

$$\Delta t \leq \frac{\Delta x}{c_{max}} \quad \text{(CFL condition)}$$

Acoustic $c \approx 340$ m/s; gravity $c \approx 300$ m/s; advective $c \approx 100$ m/s

Semi-Implicit Semi-Lagrangian

Dominant approach for global models (ECMWF IFS). Semi-implicit treats fast gravity waves implicitly (Helmholtz equation). Semi-Lagrangian traces parcels backward, removing advective CFL. Time steps 10-30x larger than explicit. Robert-Asselin filter: $\bar{\phi}^n = \phi^n + \nu(\phi^{n+1} - 2\phi^n + \phi^{n-1})/2$.

Runge-Kutta (WRF, ICON)

RK3 (Wicker-Skamarock): 3 stages, third-order accuracy. No computational mode. Split-explicit for acoustic modes (small time steps within RK stages). Suited for non-hydrostatic convection-permitting models with $\Delta x \sim 1$-3 km.

Physics Parameterizations

Sub-grid processes are the largest source of NWP model uncertainty. Key parameterization categories:

Planetary Boundary Layer

Turbulent vertical mixing in lowest 1-2 km. Local closure (Mellor-Yamada, TKE-based) or nonlocal (YSU, MYNN). EDMF (eddy-diffusivity mass-flux) schemes represent both diffusive mixing and convective plumes. Controls surface fluxes, diurnal cycle, fog, low-level jets.

Cumulus Convection

Mass-flux schemes: Arakawa-Schubert (quasi-equilibrium), Tiedtke/Bechtold (ECMWF), SAS (GFS). Not needed when $\Delta x < 2$ km. "Gray zone" (2-10 km) is problematic: convection partially resolved, partially parameterized.

Radiation (RRTMG)

16 LW bands, 14 SW bands. Correlated-k method. Accounts for H$_2$O, CO$_2$, O$_3$, clouds, aerosols. ~30% of compute; called every 15-60 min. Shortwave includes Rayleigh scattering, cloud and aerosol scattering via delta-Eddington two-stream.

Land Surface Model

Soil temperature and moisture (4-6 layers), snow cover, vegetation transpiration, runoff. Noah-MP, HTESSEL, CLM. Drives surface sensible and latent heat fluxes. Critical for forecasting 2m temperature, precipitation, and drought.

Operational NWP Models

ECMWF IFS

Spectral, semi-implicit semi-Lagrangian, T1279 (~9 km), 137 hybrid levels to 0.01 hPa. 4D-Var with 12-hour windows. 51-member ensemble at ~18 km (T639). Consistently the most skillful global model. 10-day forecasts 4x daily. Extended-range to Day 46.

GFS (NOAA/NCEP)

Finite-volume cubed-sphere (FV3), ~13 km, 127 hybrid levels. 4x daily to 16 days. Hybrid EnKF-3DVar assimilation. GEFS ensemble: 31 members at ~25 km. Freely available data. Backbone of US weather forecasting. Rapid Refresh (RAP/HRRR) at 3 km for CONUS.

ICON (DWD)

Icosahedral nonhydrostatic. Unstructured triangular grid. Global 13 km, nesting to 6.5 km (Europe) and 2 km (Germany). Non-hydrostatic dynamics allow convection-permitting nesting. Also used as climate model ICON-ESM.

ARPEGE / AROME

ARPEGE: global spectral with variable resolution (stretched grid ~10 km over France, ~60 km antipodal). AROME: non-hydrostatic at 1.3 km, convection-permitting over France. Shared IFS/ARPEGE physics heritage.

Vertical Coordinate Systems

Sigma ($\sigma$)

$\sigma = p/p_s$. Terrain-following at all levels. Simple lower boundary ($\sigma=1$). Problem: large pressure gradient errors over steep terrain due to tilted coordinate surfaces at upper levels. Historical use only.

Hybrid ($\eta$)

$p(\eta) = A(\eta)p_0 + B(\eta)p_s$. Near surface $B \to 1$ (terrain-following); aloft $B \to 0$ (pure pressure). Minimizes terrain errors in stratosphere. Standard: ECMWF 137 levels, GFS 127 levels.

Isentropic ($\theta$)

Potential temperature as vertical coordinate. Adiabatic flow is quasi-horizontal on $\theta$ surfaces. Excellent for stratospheric transport. Difficulty: $\theta$ surfaces intersect ground. Used in some research models and isentropic analysis.

Ensemble Forecasting and Predictability

The atmosphere is chaotic: small initial errors grow exponentially. Lorenz (1963) demonstrated sensitive dependence with his three-variable system:

$$\frac{dx}{dt} = \sigma(y-x), \quad \frac{dy}{dt} = x(\rho-z)-y, \quad \frac{dz}{dt} = xy-\beta z$$

$\sigma=10, \rho=28, \beta=8/3$. Lyapunov exponent $\lambda_1 \approx 0.9$; error doubling time $\approx \ln 2/\lambda_1 \approx 0.8$ time units.

Perturbation Methods

Singular vectors (ECMWF): fastest-growing perturbations over optimization period. Bred vectors (NCEP): rescaled nonlinear perturbations. Ensemble Transform. EnKF perturbations provide flow-dependent spread. Stochastic physics (SPPT, SKEB) samples model uncertainty.

Predictability Limits

Synoptic (1000 km): ~10-14 days. Mesoscale (10-100 km): 1-3 days. Convective (1-10 km): 1-6 hours. Error cascade: small-scale errors contaminate large scales (Lorenz 1969). Current skill: Day-5 quality matches 1980-era Day-3.

Spread-Skill Relationship: A well-calibrated ensemble has spread proportional to forecast error. When spread is small, the forecast is confident and usually accurate. Large spread indicates uncertain situations. Reliability diagrams and rank histograms verify calibration.

Forecast Verification

RMSE and Anomaly Correlation

RMSE measures average error magnitude. Anomaly correlation coefficient (ACC) measures pattern similarity against climatology: $ACC = \frac{\sum (f' \cdot a')}{\sqrt{\sum f'^2 \sum a'^2}}$. ACC $> 0.6$ is considered useful skill. 500 hPa ACC at Day 5 now exceeds 0.9 for leading models.

Brier Skill Score

$BSS = 1 - BS/BS_{clim}$ where $BS = \frac{1}{N}\sum(f_k - o_k)^2$. Measures probabilistic forecast quality relative to climatology. Reliability diagrams decompose BS into reliability, resolution, and uncertainty components.

Fortran: 1D Advection with Multiple Numerical Schemes

Solves the 1D advection equation $\partial\phi/\partial t + c\,\partial\phi/\partial x = 0$ using upwind, Lax-Wendroff, and leapfrog schemes, demonstrating stability, diffusion, and dispersion properties fundamental to NWP numerics.

program advection_schemes
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: nx = 200, nt = 300
  real(dp), parameter :: Lx = 10.0_dp, c = 1.0_dp, CFL = 0.8_dp
  real(dp) :: dx, dt, mu
  real(dp) :: phi_up(nx), phi_lw(nx), phi_lf(nx), phi_lfm1(nx)
  real(dp) :: phi0(nx), phi_exact(nx), phi_tmp(nx)
  real(dp) :: x, rmse_up, rmse_lw, rmse_lf
  integer  :: i, n, ip, im

  dx = Lx / real(nx, dp); dt = CFL * dx / c; mu = c * dt / dx
  write(*,'(A,F8.5,A,F8.5,A,F8.4)') 'dx=', dx, ' dt=', dt, ' CFL=', mu

  ! Initial condition: Gaussian pulse
  do i = 1, nx
    x = (i - 1) * dx
    phi0(i) = exp(-((x - Lx/4.0_dp)**2) / 0.2_dp)
  end do
  phi_up = phi0; phi_lw = phi0; phi_lf = phi0; phi_lfm1 = phi0

  ! First step for leapfrog: use upwind (Euler forward)
  do i = 1, nx
    im = mod(i - 2 + nx, nx) + 1
    phi_tmp(i) = phi0(i) - mu * (phi0(i) - phi0(im))
  end do
  phi_lf = phi_tmp

  ! Time integration
  do n = 2, nt
    ! --- Upwind scheme ---
    phi_tmp = phi_up
    do i = 1, nx
      im = mod(i - 2 + nx, nx) + 1
      phi_up(i) = phi_tmp(i) - mu * (phi_tmp(i) - phi_tmp(im))
    end do

    ! --- Lax-Wendroff scheme ---
    phi_tmp = phi_lw
    do i = 1, nx
      ip = mod(i, nx) + 1; im = mod(i - 2 + nx, nx) + 1
      phi_lw(i) = phi_tmp(i) &
        - 0.5_dp*mu*(phi_tmp(ip) - phi_tmp(im)) &
        + 0.5_dp*mu**2*(phi_tmp(ip) - 2.0_dp*phi_tmp(i) + phi_tmp(im))
    end do

    ! --- Leapfrog scheme ---
    phi_tmp = phi_lf
    do i = 1, nx
      ip = mod(i, nx) + 1; im = mod(i - 2 + nx, nx) + 1
      phi_tmp(i) = phi_lfm1(i) - mu*(phi_lf(ip) - phi_lf(im))
    end do
    phi_lfm1 = phi_lf; phi_lf = phi_tmp
  end do

  ! Exact solution
  do i = 1, nx
    x = (i - 1) * dx
    x = mod(x - c * nt * dt, Lx)
    if (x < 0.0_dp) x = x + Lx
    phi_exact(i) = exp(-((x - Lx/4.0_dp)**2) / 0.2_dp)
  end do

  ! RMSE
  rmse_up = 0.0_dp; rmse_lw = 0.0_dp; rmse_lf = 0.0_dp
  do i = 1, nx
    rmse_up = rmse_up + (phi_up(i) - phi_exact(i))**2
    rmse_lw = rmse_lw + (phi_lw(i) - phi_exact(i))**2
    rmse_lf = rmse_lf + (phi_lf(i) - phi_exact(i))**2
  end do
  rmse_up = sqrt(rmse_up/nx); rmse_lw = sqrt(rmse_lw/nx)
  rmse_lf = sqrt(rmse_lf/nx)

  write(*,'(A,ES12.5)') 'RMSE Upwind:       ', rmse_up
  write(*,'(A,ES12.5)') 'RMSE Lax-Wendroff: ', rmse_lw
  write(*,'(A,ES12.5)') 'RMSE Leapfrog:     ', rmse_lf

  write(*,'(/,A)') '  i    x      Exact     Upwind    Lax-W     Leapfrog'
  do i = 1, nx, 10
    x = (i-1)*dx
    write(*,'(I4,F7.3,4F10.5)') i, x, phi_exact(i), &
         phi_up(i), phi_lw(i), phi_lf(i)
  end do
end program advection_schemes

Interactive Simulation: 1D Advection Equation Solver

Python

Solve the 1D advection equation using three numerical schemes: FTCS (unstable), Upwind (diffusive), and Lax-Wendroff (dispersive). Demonstrates how a Gaussian pulse evolves and compares numerical diffusion and dispersion errors.

advection_solver.py84 lines

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← 7.2 Weather Stations7.4 Data Assimilation →