Part V: Weather Analysis & Forecasting

Observations, data assimilation, numerical prediction, ensemble methods, verification, and synoptic pattern recognition

1. Weather Observations

The global observing system ingests approximately 100 million observations per day from surface stations, upper-air soundings, satellites, radar, aircraft, ships, buoys, and emerging platforms. These observations are the foundation upon which all weather analysis and forecasting rests.

1.1 Surface Observations

Surface weather stations form the backbone of the global observing system, measuring key meteorological variables at standard heights established by the World Meteorological Organization (WMO):

Standard Measurements

- Temperature (2 m above ground, in Stevenson screen)
- Humidity (dewpoint temperature, relative humidity, mixing ratio)
- Pressure (station pressure and reduced to mean sea level)
- Wind (10 m height: sustained speed, direction, gusts)
- Precipitation (tipping-bucket, weighing gauge, optical)
- Visibility and present weather phenomena
- Cloud cover (oktas), cloud base height, cloud type
- Solar radiation (direct, diffuse, global)

Observation Formats & Networks

- SYNOP: Synoptic observations encoded in FM-12 format; reported every 3 or 6 hours at standard synoptic times (00, 06, 12, 18 UTC)
- METAR: Aviation routine weather reports; issued hourly or half-hourly at airports; includes visibility, ceiling, runway visual range
- SPECI: Special aviation reports triggered by rapid weather changes
- ASOS: Automated Surface Observing System (USA); 900+ stations providing 1-minute observations
- AWOS: Automated Weather Observing System; simpler than ASOS
- Buoys: Moored and drifting marine platforms; measure SST, wave height, pressure
- Ships (VOS): Voluntary Observing Ships contribute marine surface data

Sea-Level Pressure Reduction

Station pressure must be reduced to mean sea level for surface analysis. The standard reduction uses the hypsometric equation:

$p_0 = p_s \exp\!\left(\frac{g\, z_s}{R\,\bar{T}_v}\right)$

where $p_s$ is station pressure, $z_s$ is station elevation, $\bar{T}_v$ is the mean virtual temperature of the fictitious air column below the station, and $R$ is the specific gas constant for dry air.

1.2 Upper-Air Soundings

Vertical profiles of the atmosphere are essential for initializing forecast models and understanding the three-dimensional structure of weather systems.

Radiosondes

Approximately 800 stations worldwide launch balloon-borne radiosondes twice daily at 00 UTC and 12 UTC. Each instrument package measures temperature (thermistor or capacitive sensor), relative humidity (thin-film capacitor), and pressure (silicon capacitive sensor). Wind speed and direction are derived from GPS tracking of the balloon position. Ascent rate is approximately 5 m/s, and the balloon typically bursts near 30-35 km altitude (~10 hPa). Data are transmitted in real time via 400-406 MHz radio link. Vertical resolution is approximately 5-10 m when using modern GPS sondes.

Common radiosonde types include Vaisala RS41 (widely considered the gold standard), Meisei iMS-100, and Graw DFM-17. Temperature accuracy is typically 0.3-0.5 K, with known biases from solar radiation heating that require correction.

Dropsondes

Deployed from aircraft (e.g., NOAA Hurricane Hunters, research flights) into weather systems. Fall via parachute at approximately 10-12 m/s, transmitting the same variables as radiosondes but from the top down. Critical for observing tropical cyclones, polar lows, and other systems over data-sparse oceanic regions. The NCAR RD-41 dropsonde is the current standard, providing GPS winds with 0.5 m/s accuracy.

Pilot Balloons (PIBALs) and Wind Profilers

Pilot balloons are tracked optically or by theodolite to obtain wind profiles without thermodynamic data. Doppler wind profilers (915 MHz boundary-layer or 404 MHz tropospheric) use radar to continuously measure vertical profiles of horizontal wind, providing 6-minute temporal resolution. The NOAA Profiler Network operates over 30 sites across the central United States, filling temporal gaps between radiosonde launches.

1.3 Remote Sensing: Satellites

Meteorological satellites provide the only truly global observing capability, delivering imagery and sounding data continuously over both populated and remote regions.

Geostationary Satellites (GEO)

Orbiting at ~36,000 km altitude in the equatorial plane, geostationary satellites maintain a fixed position relative to Earth, enabling continuous monitoring. Key systems include:

- GOES-16/18 (USA): Advanced Baseline Imager (ABI) with 16 channels, 0.5-2 km resolution, full-disk scan every 10 minutes, mesoscale sector every 1 minute
- Meteosat Third Generation (Europe): Flexible Combined Imager (FCI) and Lightning Imager (LI)
- Himawari-8/9 (Japan): 16-channel AHI, 10-minute full-disk scans over western Pacific
- FY-4 series (China): AGRI imager and GIIRS hyperspectral sounder

Key Spectral Channels:

- Visible (VIS, ~0.6 um): Reflected sunlight; shows cloud thickness, texture; daytime only
- Near-IR (NIR, 1.6 um): Ice vs. water cloud discrimination (ice absorbs, water reflects)
- Shortwave IR (SWIR, 3.9 um): Fog detection (day/night), fire hot spots
- Water vapor (WV, 6.2/6.9/7.3 um): Upper and mid-level moisture; reveals jet stream patterns and dry intrusions even in cloud-free regions
- Infrared window (IR, 10.3 um): Cloud-top temperature; available day and night; colder tops = higher clouds
- CO2 absorption (13.3 um): Height assignment for semi-transparent clouds

Polar-Orbiting Satellites (LEO)

Sun-synchronous orbits at ~830 km altitude provide global coverage every 12 hours with higher resolution and additional sounding capabilities:

- NOAA-20/21 (JPSS, USA): VIIRS imager (375 m resolution), CrIS cross-track IR sounder (2,211 channels), ATMS microwave sounder
- MetOp-B/C (EUMETSAT): IASI hyperspectral IR sounder (8,461 channels), AMSU/MHS microwave
- Suomi NPP: VIIRS day-night band enables nighttime visible imagery
- FY-3 series (China): Multiple imagers and sounders for global NWP

Satellite sounders retrieve vertical temperature and humidity profiles by measuring upwelling radiance at multiple wavelengths. Each channel peaks at a different altitude, and mathematical inversion (radiative transfer) recovers the profile. These are the single most impactful observation type for NWP skill in the Southern Hemisphere.

Satellite-Derived Atmospheric Motion Vectors (AMVs)

Wind vectors are derived by tracking cloud or water vapor features in successive satellite images. Height assignment (using IR brightness temperature or CO2-slicing) is the primary source of error. AMVs provide roughly 300,000 wind observations per day, critical over oceans where conventional data are sparse.

1.4 Weather Radar

Weather radar is the primary tool for detecting and monitoring precipitation, severe storms, and mesoscale phenomena in real time.

The Weather Radar Equation

The received power from meteorological targets (distributed scatterers) is given by the Probert-Jones radar equation:

$P_r = \frac{\pi^3 P_t G^2 \theta^2 h |K|^2}{1024\ln(2)\lambda^2} \frac{Z}{r^2}$

- $P_r$: received power (W)

- $P_t$: transmitted power (typically 250-750 kW peak)

- $G$: antenna gain

- $\theta$: beamwidth (typically ~1 degree for weather radar)

- $h$: pulse length (spatial extent $= c\tau/2$)

- $|K|^2$: dielectric factor (0.93 for water, 0.197 for ice)

- $\lambda$: wavelength (S-band ~10 cm, C-band ~5 cm, X-band ~3 cm)

- $Z$: radar reflectivity factor ($\text{mm}^6\,\text{m}^{-3}$)

- $r$: range to target

Reflectivity and the Z-R Relationship

Reflectivity Z is related to the drop size distribution (DSD). For a distribution of spherical drops:

$Z = \int_0^\infty D^6 N(D)\, dD$

Reflectivity is commonly expressed in dBZ: $\text{dBZ} = 10\log_{10}(Z)$. The Marshall-Palmer Z-R relationship converts reflectivity to rainfall rate:

$Z = 200\, R^{1.6}$

where R is in mm/hr. Different Z-R relations exist for tropical rain ($Z = 250R^{1.2}$), snow ($Z = 2000R^{2.0}$), and other precipitation types.

Doppler Velocity

Doppler radar measures the radial velocity of scatterers using the frequency shift between transmitted and received pulses. The Nyquist velocity (maximum unambiguous velocity) is:

$V_{\max} = \frac{\lambda\, \text{PRF}}{4}$

where PRF is the pulse repetition frequency. Velocities beyond $\pm V_{\max}$ are aliased (folded). Doppler data reveal mesocyclones, tornado vortex signatures (TVS), convergence/divergence, and wind shear.

Dual-Polarization Variables

Modern radars transmit both horizontally and vertically polarized pulses, yielding additional variables:

- $Z_{DR}$ (Differential Reflectivity): Ratio of horizontal to vertical reflectivity; indicates particle shape (oblate raindrops: $Z_{DR} > 0$; tumbling hail: $Z_{DR} \approx 0$)
- $\rho_{HV}$ (Correlation Coefficient): Correlation between H and V signals; high for uniform hydrometeors (~0.99), low for mixed-phase or non-meteorological targets (<0.90 for debris)
- $K_{DP}$ (Specific Differential Phase): Rate of change of differential phase; proportional to liquid water content, immune to attenuation and hail contamination
- $\Phi_{DP}$ (Differential Phase): Cumulative phase difference; used to correct for attenuation at C-band and X-band

1.5 Aircraft Observations (AMDAR)

The Aircraft Meteorological Data Relay (AMDAR) program collects meteorological data from commercial aircraft during all phases of flight. Over 300,000 reports per day from 40+ airlines provide temperature (accuracy ~0.3 K), wind speed and direction (accuracy ~1-2 m/s), and increasingly humidity (via WVSS-II sensor). Aircraft profiles during ascent and descent provide excellent high-resolution vertical soundings near major airports, complementing the twice-daily radiosonde network. AMDAR data have been shown to have significant positive impact on NWP, particularly for short-range forecasts in data-rich areas.

TAMDAR (Tropospheric Airborne Meteorological Data Reporting) extends coverage to regional aircraft that fly at lower altitudes, sampling the boundary layer and lower troposphere more frequently.

1.6 GPS Radio Occultation

GPS radio occultation (GPS-RO) measures the bending of GPS signals as they pass tangentially through the atmosphere from a GPS satellite to a receiver on a low-Earth-orbit satellite. The bending angle is related to the vertical gradient of refractivity:

$N = 77.6\frac{p}{T} + 3.73 \times 10^5 \frac{e}{T^2}$

where N is refractivity, p is pressure (hPa), T is temperature (K), and e is water vapor pressure (hPa). Through Abel inversion, the bending angle profile yields a refractivity profile, from which temperature and humidity can be retrieved. GPS-RO provides:

- All-weather capability (unaffected by clouds or precipitation)
- High vertical resolution (~100-200 m in the upper troposphere/stratosphere)
- Self-calibrating: no instrument drift, absolute accuracy ~0.1 K
- Global coverage including data-sparse oceanic and polar regions
- Constellations: COSMIC-2 (6 satellites), Spire, PlanetiQ provide ~10,000 profiles/day

Video: How Weather Satellites Work

NOAA explanation of satellite observing systems and their role in weather forecasting

2. Data Assimilation

Data assimilation is the mathematical and computational process of combining heterogeneous, irregularly distributed observations with a prior estimate (background) from a model forecast to produce the best possible estimate of the atmospheric state -- the analysis. This analysis serves as the initial condition for the next forecast cycle.

2.1 The Analysis Problem

Observations are sparse (covering only a fraction of model grid points), irregular in space and time, heterogeneous (different instruments measuring different quantities), and contain measurement errors. Model forecasts provide complete, dynamically consistent coverage but accumulate errors over time. Data assimilation seeks the analysis $\mathbf{x}^a$ that optimally weights both:

Analysis Equation (Kalman Update)

$\mathbf{x}^a = \mathbf{x}^b + \mathbf{K}(\mathbf{y} - H(\mathbf{x}^b))$

- $\mathbf{x}^a$: Analysis state vector (best estimate of the atmosphere)

- $\mathbf{x}^b$: Background state (short-range forecast, typically 3-6 hours)

- $\mathbf{y}$: Observation vector

- $H$: Observation operator (maps model state to observation space, may include radiative transfer models)

- $\mathbf{y} - H(\mathbf{x}^b)$: Innovation vector (observation minus background in observation space)

- $\mathbf{K}$: Kalman gain matrix determines the weight given to observations

Kalman Gain and Error Covariances

$\mathbf{K} = \mathbf{B}\mathbf{H}^T(\mathbf{H}\mathbf{B}\mathbf{H}^T + \mathbf{R})^{-1}$

- $\mathbf{B}$: Background error covariance matrix (describes uncertainty in $\mathbf{x}^b$)

- $\mathbf{R}$: Observation error covariance matrix (describes uncertainty in $\mathbf{y}$)

- If $\mathbf{B} \gg \mathbf{R}$: analysis trusts observations more ($\mathbf{K} \to \mathbf{H}^{-1}$)

- If $\mathbf{R} \gg \mathbf{B}$: analysis trusts background more ($\mathbf{K} \to 0$)

Analysis Increment and Innovation

The analysis increment $\delta\mathbf{x}^a = \mathbf{x}^a - \mathbf{x}^b$ represents the correction made by observations. The innovation $\mathbf{d} = \mathbf{y} - H(\mathbf{x}^b)$ measures the mismatch between observations and the background. Innovation statistics are monitored in real time:

- Mean innovation should be near zero (no systematic bias)
- Innovation variance should match $\mathbf{HBH}^T + \mathbf{R}$ (consistent error statistics)
- Departures from these indicate observation or model biases requiring correction

2.2 Optimal Interpolation (OI)

The classical method of data assimilation, OI computes the analysis at each grid point using a local selection of nearby observations. Background error correlations are prescribed (e.g., Gaussian functions of distance) and assumed stationary and isotropic:

$B_{ij} = \sigma_b^2 \exp\!\left(-\frac{d_{ij}^2}{2L^2}\right)$

where $L$ is the correlation length scale and $d_{ij}$ is the distance between points. OI is computationally efficient but cannot account for flow-dependent error structures. It was the operational standard from the 1970s to 1990s and is still used for some specialized analyses (e.g., sea surface temperature, soil moisture).

2.3 Variational Methods

3D-Var (3-Dimensional Variational Assimilation)

3D-Var finds the analysis $\mathbf{x}^a$ that minimizes the cost function:

$J(\mathbf{x}) = \frac{1}{2}(\mathbf{x}-\mathbf{x}_b)^T\mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}_b) + \frac{1}{2}(\mathbf{y}-H(\mathbf{x}))^T\mathbf{R}^{-1}(\mathbf{y}-H(\mathbf{x}))$

The first term ($J_b$) penalizes departure from the background; the second term ($J_o$) penalizes departure from observations. Setting $\nabla J = 0$ yields the analysis equation. Minimization is performed iteratively using conjugate gradient or quasi-Newton methods (e.g., L-BFGS).

The background error covariance matrix B is the most critical and difficult component. For a global model with $\sim 10^9$ state variables, B is $10^9 \times 10^9$ and cannot be stored explicitly. Instead, B is modeled using a sequence of operators: variable transforms (to decouple balanced/unbalanced components), spectral transforms, and vertical correlations.

4D-Var (4-Dimensional Variational Assimilation)

4D-Var extends the 3D-Var cost function to include observations distributed over a time window (typically 6 or 12 hours):

$J(\mathbf{x}_0) = J_b + \sum_{i=0}^{N} \frac{1}{2}(\mathbf{y}_i - H_i(\mathbf{x}_i))^T \mathbf{R}_i^{-1}(\mathbf{y}_i - H_i(\mathbf{x}_i))$

where $\mathbf{x}_i = M_{0 \to i}(\mathbf{x}_0)$ is the model state propagated from the initial time$t_0$ to observation time $t_i$ using the forecast model M. The gradient of J requires the adjoint model $\mathbf{M}^T$, making development and maintenance expensive.

4D-Var implicitly evolves the background error covariance B using the model dynamics (flow-dependent covariances), enables assimilation of time-distributed observations at their actual time, and uses the model as a strong constraint for dynamical consistency. It is the operational method at ECMWF, Meteo-France, and JMA.

2.4 Ensemble Kalman Filter (EnKF)

The EnKF uses an ensemble of model forecasts to estimate flow-dependent background error covariances:

$\mathbf{P}^b \approx \frac{1}{N-1}\sum_{k=1}^{N}(\mathbf{x}_k^b - \bar{\mathbf{x}}^b)(\mathbf{x}_k^b - \bar{\mathbf{x}}^b)^T$

where $N$ is the ensemble size (typically 20-80 members) and $\bar{\mathbf{x}}^b$ is the ensemble mean. Each member is updated independently using the same observations but with perturbed observations (stochastic EnKF) or a deterministic square-root filter (e.g., LETKF -- Local Ensemble Transform Kalman Filter).

Key advantages of EnKF:

- Flow-dependent covariances automatically (sharp fronts get anisotropic error structures)
- No adjoint model required (easier implementation)
- Naturally provides ensemble for uncertainty estimation
- Parallelizable across ensemble members

Challenges include sampling noise from finite ensemble size (mitigated by localization, which tapers correlations beyond a cutoff distance) and filter divergence (mitigated by multiplicative or additive inflation that increases ensemble spread). Operational at Environment and Climate Change Canada and NOAA (hybrid with 3D-Var).

2.5 Background Error Covariance Modeling

The B matrix encodes how analysis increments spread information from observation locations to surrounding grid points. Proper modeling requires:

- Multivariate balance: A wind observation should adjust both wind and mass fields consistently (via geostrophic or gradient wind balance): $\delta\Phi = f\,\psi$ where $\psi$ is the streamfunction increment
- Horizontal correlations: Typically modeled using recursive filters or spectral methods; length scales vary with variable, level, and latitude (~100-300 km for mass, shorter for wind)
- Vertical correlations: Empirical orthogonal functions (EOFs) of forecast differences capture dominant vertical structures
- Hybrid methods: Modern systems blend static B (from climatological statistics) with flow-dependent ensemble covariances: $\mathbf{B}_{\text{hybrid}} = \beta_1 \mathbf{B}_{\text{static}} + \beta_2 \mathbf{B}_{\text{ensemble}}$

3. Numerical Weather Prediction (NWP)

Numerical weather prediction solves the governing equations of the atmosphere as an initial-value problem: given the current state (from data assimilation), integrate forward in time to produce a forecast. Modern NWP represents one of the greatest achievements in computational science.

3.1 The Primitive Equations

NWP models solve the primitive equations, which are the fundamental conservation laws for atmospheric motion under the hydrostatic approximation:

Vector Form of the Momentum Equation

$\frac{\partial \vec{V}}{\partial t} + (\vec{V}\cdot\nabla)\vec{V} + f\hat{k}\times\vec{V} = -\nabla\Phi - R T\nabla\ln p$

Horizontal momentum (component form):

$\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$

Hydrostatic equation:

$\frac{\partial \Phi}{\partial p} = -\frac{RT}{p} = -\alpha$

Continuity equation (pressure coordinates):

$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial \omega}{\partial p} = 0$

Thermodynamic energy equation:

$\frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} + \omega\left(\frac{\partial T}{\partial p} - \frac{\alpha}{c_p}\right) = \frac{\dot{Q}}{c_p}$

Moisture equation:

$\frac{\partial q}{\partial t} + u\frac{\partial q}{\partial x} + v\frac{\partial q}{\partial y} + \omega\frac{\partial q}{\partial p} = S_q$

Equation of state:

$p = \rho R_d T_v, \quad T_v = T(1 + 0.608q)$

3.2 Spectral Methods vs. Grid-Point Methods

Spectral Transform Method

Fields are represented as truncated series of spherical harmonics $Y_n^m(\lambda, \phi)$:

$\psi(\lambda,\phi) = \sum_{n=0}^{N}\sum_{m=-n}^{n} \hat{\psi}_n^m Y_n^m(\lambda,\phi)$

Linear terms (Coriolis, pressure gradient) are computed in spectral space; nonlinear terms are computed on a Gaussian grid and transformed back. Spectral methods provide excellent accuracy for global models (no pole problems) but are expensive at high resolution due to Legendre transforms. Used by ECMWF IFS (up to T1279 ~ 9 km) and GFS (T1534 ~ 13 km).

Grid-Point and Finite-Volume Methods

Discretize equations directly on a computational grid. Modern approaches include:

- Finite difference: Approximate derivatives using neighboring grid point values; simple but diffusive
- Finite volume: Conserve fluxes across cell boundaries; excellent conservation properties; used by MPAS, FV3 (GFS dynamical core)
- Spectral element / discontinuous Galerkin: High-order accuracy on unstructured grids; used by ICON (DWD/MPI), NUMA
- Cubed-sphere grids: Avoid pole singularities by mapping six cube faces to the sphere; used by FV3, LFRic
- Icosahedral grids: Quasi-uniform triangular or hexagonal meshes; used by ICON, NICAM, MPAS

3.3 Semi-Lagrangian Advection

The semi-Lagrangian (SL) method traces trajectories backward from grid points to find departure points, then interpolates the field at the departure point:

$\psi(\mathbf{x}_a, t+\Delta t) = \psi(\mathbf{x}_d, t) + \Delta t\, S$

where $\mathbf{x}_d$ is the departure point found by solving $\mathbf{x}_a - \mathbf{x}_d = \int_t^{t+\Delta t} \mathbf{V}\, dt'$iteratively. The key advantage is that SL is unconditionally stable for advection, allowing much larger time steps than Eulerian methods. This relaxes the CFL condition, enabling time steps 5-10 times larger than explicit methods. Used by ECMWF IFS, GEM (Canada), and Unified Model (UKMO).

3.4 CFL Condition and Time Stepping

Courant-Friedrichs-Lewy (CFL) Condition

For explicit numerical schemes, stability requires that the numerical domain of dependence contain the physical domain of dependence:

$\Delta t \leq \frac{\Delta x}{c_{\max}}$

where $c_{\max}$ is the fastest wave speed in the system. For atmospheric models, the fastest modes are acoustic waves (~330 m/s) and gravity waves (~100 m/s). At $\Delta x = 10$ km with explicit acoustics, this would require $\Delta t \leq 30$ s. Strategies to allow larger time steps:

- Semi-implicit methods: Treat fast waves implicitly (solve elliptic equation each step); $\Delta t$ limited by advection CFL (~600 s at 10 km)
- Split-explicit: Integrate acoustic modes with small substeps; used in WRF
- Semi-Lagrangian + semi-implicit: No advective CFL; $\Delta t$ limited only by accuracy (~900 s at 10 km); used by ECMWF, UKMO

3.5 Physics Parameterization Schemes

Subgrid-scale processes that cannot be explicitly resolved must be parameterized -- their bulk effects on resolved scales are represented by simplified models:

Cumulus Convection

Represents vertical transport of heat, moisture, and momentum by convective clouds that are smaller than the grid spacing. Approaches include mass-flux schemes (e.g., Arakawa-Schubert, Tiedtke, Zhang-McFarlane) that model an ensemble of cloud plumes, and adjustment schemes (Betts-Miller-Janjic) that relax profiles toward observed post-convective states. At convection-permitting resolution ($\Delta x < 4$ km), deep convection is increasingly resolved explicitly.

Radiation

Computes heating/cooling rates from shortwave (solar, 0.2-4 um) and longwave (terrestrial, 4-100 um) radiation transfer through the atmosphere. Must account for absorption/emission by gases (H2O, CO2, O3), scattering by clouds and aerosols, and surface reflection. Computationally expensive; often called every 15-60 minutes. Methods: RRTMG (correlated-k), Fu-Liou, ecRad (ECMWF).

Planetary Boundary Layer (PBL)

Parameterizes turbulent mixing in the lowest 1-2 km. Local schemes (e.g., Mellor-Yamada) relate fluxes to local gradients via eddy diffusivity $K_h$. Nonlocal schemes (e.g., YSU, ACM2) add countergradient terms for convective boundary layers where large eddies transport against the local gradient. The surface layer uses Monin-Obukhov similarity theory.

Cloud Microphysics

Represents phase changes and precipitation processes. Bulk schemes predict mass mixing ratios of hydrometeor categories (cloud water, rain, ice, snow, graupel, hail). Bin/spectral schemespredict the full drop/ice size distribution (much more accurate but 100x more expensive). Key processes: autoconversion, accretion, evaporation, melting, riming, aggregation, ice nucleation.

Land Surface Models

Compute surface energy and water budgets: $R_n = H + LE + G$ where $R_n$ is net radiation, H sensible heat flux, LE latent heat flux, G ground heat flux. Multi-layer soil temperature and moisture, snowpack, vegetation transpiration, canopy interception. Models: Noah-MP, JULES, CLM, HTESSEL.

Gravity Wave Drag

Parameterizes momentum deposition from orographic gravity waves (generated by flow over mountains) and non-orographic gravity waves (from convection, fronts, jet imbalance). Critical for correct stratospheric circulation, QBO simulation, and surface pressure patterns. The drag force decelerates the flow where waves break or reach critical levels.

3.6 Major Operational Models

ECMWF IFS (Integrated Forecasting System)

European Centre for Medium-Range Weather Forecasts. Spectral dynamical core (cubic octahedral reduced Gaussian grid). HRES: TCo1279 (~9 km), 137 vertical levels to 0.01 hPa, 4D-Var data assimilation. ENS: 51 members at TCo639 (~18 km). Extended range to 46 days. Generally ranked as the most skillful global model. Transitioning to a finite-volume dynamical core on an icosahedral grid. Computational cost: ~100 million core-hours per year.

GFS (Global Forecast System) -- NOAA/NCEP

FV3 dynamical core on cubed-sphere grid. Resolution: C768 (~13 km), 127 vertical levels. Hybrid 4D-EnVar data assimilation (combines ensemble covariances with static B). GEFS ensemble: 31 members at C384 (~25 km). Forecasts to 16 days (deterministic) and 35 days (ensemble). Freely available model output has made GFS the world's most widely used global model for downstream applications.

ICON (Icosahedral Nonhydrostatic Model) -- DWD/MPI

Unstructured triangular grid on the icosahedron. Nonhydrostatic dynamical core enabling seamless application from global (13 km) to convection-permitting (2.2 km for Europe, 1 km nests). Two-moment microphysics, RRTMG radiation. Designed for GPU computing. Operational at DWD since 2015.

WRF (Weather Research and Forecasting Model)

Community mesoscale model developed by NCAR, NOAA, and partners. Compressible nonhydrostatic equations, Arakawa-C staggered grid, Runge-Kutta time integration with split-explicit acoustic steps. Highly configurable with dozens of physics options. Supports one-way and two-way nesting from ~30 km to ~100 m (large-eddy simulation). Used operationally and in research worldwide.

HRRR (High-Resolution Rapid Refresh) -- NOAA

Convection-permitting model covering CONUS at 3 km resolution, 18-hour forecasts updated hourly. WRF-based with extensive radar data assimilation (3D reflectivity, radial velocity). Excels at predicting convective initiation, storm evolution, aviation weather, and renewable energy forecasts. RRFS (Rapid Refresh Forecast System) on FV3 is the next-generation replacement.

3.7 Resolution and Computational Cost

Doubling horizontal resolution increases computational cost by approximately a factor of 8-10 in 3D (factor 4 from grid points, factor 2 from smaller time step per CFL):

$\text{Cost} \propto \frac{N_x \cdot N_y \cdot N_z}{\Delta t} \propto \Delta x^{-4}$

A single 10-day ECMWF HRES forecast at 9 km uses approximately 50,000 core-hours. The full operational suite (HRES + ENS + extended range + reanalysis) requires one of the world's largest supercomputers (~100 PFLOPS). Moving to 1 km global resolution (the "digital twin" vision) would require exascale computing (~1000 PFLOPS), currently feasible only for limited-area or short forecasts.

Video: How Weather Models Work

Met Office explanation of numerical weather prediction fundamentals

Video: Numerical Weather Prediction

Comprehensive overview of NWP: from equations to operational forecasting

4. Ensemble Forecasting

4.1 Initial Condition Uncertainty and Chaos

The atmosphere is a chaotic dynamical system in the sense of Lorenz (1963): small errors in initial conditions grow exponentially, eventually saturating at the scale of natural variability. This fundamentally limits deterministic predictability:

Error Growth in Chaotic Systems

$|\delta \mathbf{x}(t)| \sim |\delta \mathbf{x}_0|\, e^{\lambda t}$

- $\delta\mathbf{x}_0$: Initial error (analysis uncertainty, typically ~1 K, ~1 m/s)

- $\lambda$: Leading Lyapunov exponent ($\approx 0.35-0.5$ day$^{-1}$ for the atmosphere)

- Error doubling time: $t_d = \ln(2)/\lambda \approx 1.4-2.0$ days

- Saturation occurs when errors reach climatological variance

- Practical limit of deterministic forecasting: approximately 10-14 days

Since a single deterministic forecast cannot capture this inherent uncertainty, ensemble forecasting runs multiple forecasts (members) with perturbed initial conditions and/or model physics to sample the probability distribution of possible outcomes.

4.2 Ensemble Generation Methods

Bred Vectors (BVs)

Developed by Toth and Kalnay (1993). A random perturbation is added to the analysis, integrated forward, and periodically rescaled to prevent nonlinear saturation. After several breeding cycles (typically 24-48 hours), the perturbation aligns with the fastest-growing instabilities of the current flow. Simple to implement but tends to collapse onto a small subspace. Used historically by NCEP GEFS.

Singular Vectors (SVs)

Developed by Buizza and Palmer (1995) at ECMWF. Singular vectors are the perturbations that grow most rapidly over a finite optimization time (typically 48 hours) in a linear sense:

$\max_{\delta\mathbf{x}_0} \frac{||\mathbf{M}\,\delta\mathbf{x}_0||^2}{||\delta\mathbf{x}_0||^2}$

where M is the tangent-linear model propagator. SVs are computed using the adjoint model and Lanczos iteration. They identify baroclinically sensitive regions and are combined to produce initial perturbations that optimally span the space of growing errors. Used operationally at ECMWF for the ENS.

Ensemble Transform Kalman Filter (ETKF) Perturbations

Uses the analysis-cycle ensemble to generate perturbations that are consistent with the analysis uncertainty. The ETKF computes a transformation matrix that converts the forecast perturbation matrix into analysis perturbations while preserving the ensemble variance. This approach naturally accounts for observation network density (smaller perturbations in well-observed regions). Used by NCEP GEFS (since 2019), UKMO MOGREPS.

Model Uncertainty Representation

Initial condition perturbations alone underestimate forecast uncertainty because they neglect model error. Approaches to represent model uncertainty:

- SPPT (Stochastically Perturbed Parameterization Tendencies): Multiply physics tendencies by a random field $1 + r$ where $r$ is spatially and temporally correlated noise
- SKEB (Stochastic Kinetic Energy Backscatter): Injects energy at scales where it is spuriously dissipated by the dynamical core
- SPP (Stochastic Parameter Perturbation): Randomly varies uncertain parameters within physics schemes
- Multi-model ensembles: Combine forecasts from different NWP systems (TIGGE, NAEFS)

4.3 Ensemble Mean and Spread

Ensemble Mean

$\bar{x} = \frac{1}{N}\sum_{i=1}^{N} x_i$

The ensemble mean filters out unpredictable small-scale features, retaining only the signal common to all members. Beyond day 3-4, the ensemble mean typically has lower RMSE than any individual member (including the unperturbed control) because averaging reduces random error variance by a factor of$1/N$ for independent errors.

Ensemble Spread (Standard Deviation)

$\sigma^2 = \frac{1}{N-1}\sum_{i=1}^{N} (x_i - \bar{x})^2$

The ensemble spread quantifies forecast uncertainty at each grid point and time. A well-calibrated ensemble has spread that, on average, equals the RMSE of the ensemble mean. If spread is too small (underdispersive), the ensemble is overconfident; if too large, it is underconfident. The spread-skill relationship is a key diagnostic of ensemble quality.

4.4 Probability Forecasts and Products

Exceedance probabilities: The probability that a threshold is exceeded is estimated as the fraction of members exceeding it:$P(x > \tau) = \frac{1}{N}\sum_{i=1}^{N} \mathbf{1}(x_i > \tau)$. Example: P(rainfall > 25 mm/24h) = 15/51 = 29%.

Spaghetti plots: Overlay a single contour (e.g., 500 hPa 5400 m height) from all ensemble members on one map. Tight clustering = high confidence; large spread = low confidence. Particularly useful for identifying regime uncertainty (e.g., trough vs. ridge).

Probability maps: Color-coded maps showing the probability of events (precipitation exceeding a threshold, temperatures below freezing, wind gusts above warning criteria). Most useful communication tool for end users and emergency managers.

Plume (meteogram) diagrams: Time series of all ensemble members for a specific location, showing the evolution of uncertainty over the forecast period. Often displayed as a box-and-whisker plot or fan chart.

Ensemble Prediction Intervals (EPI): The 10th-90th percentile range provides a probabilistic "cone of uncertainty" (analogous to tropical cyclone track cones).

4.5 Multi-Model Ensembles

Combining ensembles from multiple NWP centers further improves probabilistic forecast skill by sampling model structural uncertainty:

- TIGGE (THORPEX Interactive Grand Global Ensemble): Archives ensemble forecasts from 10+ global centers; widely used in research
- NAEFS (North American Ensemble Forecast System): Combines GEFS (USA) and GEPS (Canada) for improved probabilistic guidance
- EPS-grams: Multi-model meteograms available from ECMWF combining multiple ensemble systems
- Bayesian Model Averaging (BMA): Weights models by their recent performance to produce calibrated probability distributions

Skill vs. Lead Time:

- Days 1-3: Deterministic (HRES) forecast often best; ensemble spread still small

- Days 4-7: Ensemble mean beats HRES; probability forecasts add clear value

- Days 7-14: Only probabilistic forecasts are skillful; regime-dependent skill

- Weeks 3-4: Subseasonal prediction; skill from MJO, soil moisture, sea ice

- Beyond 4 weeks: Seasonal forecasts; skill from ENSO, boundary conditions

5. Forecast Verification

Forecast verification is the objective assessment of forecast quality against observations or analyses. It is essential for model development, user confidence, and understanding predictability limits. Verification should address multiple attributes: accuracy, skill, reliability, resolution, and sharpness.

5.1 Deterministic Verification Scores

Mean Absolute Error (MAE)

$\text{MAE} = \frac{1}{N}\sum_{i=1}^N |f_i - o_i|$

Average absolute difference between forecast (f) and observation (o). Easy to interpret in physical units (e.g., MAE = 1.5 K for temperature). Less sensitive to outliers than RMSE.

Root Mean Square Error (RMSE)

$\text{RMSE} = \sqrt{\frac{1}{N}\sum_{i=1}^N (f_i - o_i)^2}$

Penalizes large errors more heavily than MAE due to the squaring. Standard metric for NWP verification. Can be decomposed into bias and random error components:$\text{RMSE}^2 = \text{Bias}^2 + \text{Variance}$.

Anomaly Correlation Coefficient (ACC)

$\text{ACC} = \frac{\sum_i (f_i' - c_i')(o_i' - c_i')}{\sqrt{\sum_i (f_i' - c_i')^2 \sum_i (o_i' - c_i')^2}}$

Correlation between forecast and observed anomalies from climatology (c). Measures pattern accuracy independent of bias. ACC = 1 is perfect; ACC > 0.6 is traditionally considered the threshold of useful skill. The day at which 500 hPa ACC drops to 0.6 has improved from ~5 days (1980) to ~8-9 days (2024) for ECMWF.

5.2 Categorical (Contingency Table) Scores

For yes/no forecasts (e.g., precipitation exceeding a threshold), the 2x2 contingency table defines: Hits (a), False Alarms (b), Misses (c), Correct Negatives (d). From these:

Probability of Detection (POD / Hit Rate): $\text{POD} = a/(a+c)$

False Alarm Ratio (FAR): $\text{FAR} = b/(a+b)$

Critical Success Index (CSI / Threat Score): $\text{CSI} = a/(a+b+c)$

Equitable Threat Score (ETS / Gilbert Skill Score):

$\text{ETS} = \frac{a - a_r}{a + b + c - a_r}, \quad a_r = \frac{(a+b)(a+c)}{a+b+c+d}$

where $a_r$ is the expected number of hits by random chance. ETS = 0 for no skill, ETS = 1 perfect.

Heidke Skill Score (HSS):

$\text{HSS} = \frac{2(ad - bc)}{(a+c)(c+d) + (a+b)(b+d)}$

Measures fractional improvement over random chance. HSS = 0 for random, HSS = 1 for perfect.

5.3 Probabilistic Verification Scores

Brier Score (BS)

$\text{BS} = \frac{1}{N}\sum_{i=1}^N (f_i - o_i)^2$

where $f_i$ is the forecast probability (0 to 1) and $o_i$ is the binary observation (0 or 1). BS = 0 is perfect; BS = 0.25 for constant 50% forecasts of a 50% event. The Brier Score decomposes into three illuminating components:

$\text{BS} = \underbrace{\text{Reliability}}_{\text{calibration}} - \underbrace{\text{Resolution}}_{\text{separation}} + \underbrace{\text{Uncertainty}}_{\text{climatology}}$

Reliability: Measures how well forecast probabilities match observed frequencies (smaller = better). Resolution: Measures how much forecast probabilities differ from climatological frequency (larger = better). Uncertainty: Depends only on event frequency (constant for a given sample).

Brier Skill Score (BSS)

$\text{BSS} = 1 - \frac{\text{BS}}{\text{BS}_{\text{ref}}}$

Compares the Brier Score to a reference forecast (typically climatological frequency). BSS = 1 perfect, BSS = 0 no better than climatology, BSS < 0 worse than climatology.

Continuous Ranked Probability Score (CRPS)

$\text{CRPS} = \int_{-\infty}^{\infty} [F(x) - F_o(x)]^2\, dx$

where $F(x)$ is the forecast cumulative distribution function and $F_o(x)$ is the observation CDF (a Heaviside step function at the observed value). CRPS generalizes MAE to probabilistic forecasts: for a deterministic forecast, CRPS reduces to MAE. It has units of the forecast variable, making it intuitive. Lower is better. The CRPS is the standard metric for evaluating ensemble and statistical post-processing systems.

Reliability Diagrams

Plot observed relative frequency (y-axis) against forecast probability (x-axis) for binned probability categories. A perfectly reliable forecast lies on the 1:1 diagonal. Points above the diagonal indicate underforecasting (events occur more often than forecast); points below indicate overforecasting. Sharpness histograms (showing how often each probability is issued) complement reliability: a system that always forecasts climatology has perfect reliability but no resolution.

ROC Curves and AUC

The Relative Operating Characteristic (ROC) curve plots hit rate (POD) versus false alarm rate (POFD = b/(b+d)) at different probability thresholds. The Area Under the Curve (AUC) measures discrimination ability:

- AUC = 1.0: perfect discrimination (events and non-events perfectly separated)
- AUC = 0.5: no discrimination (equivalent to random forecasting)
- AUC < 0.5: worse than random (should invert the forecast)

ROC measures potential skill (can be improved by recalibration), while reliability measures actual calibration. A forecast system should have both good ROC (discrimination) and good reliability.

5.4 Predictability Limits: Lorenz (1963)

Edward Lorenz demonstrated in his 1963 paper that deterministic systems can exhibit sensitive dependence on initial conditions, famously illustrated by his three-variable convection model:

$\dot{x} = \sigma(y - x)$

$\dot{y} = x(\rho - z) - y$

$\dot{z} = xy - \beta z$

For $\sigma = 10$, $\rho = 28$, $\beta = 8/3$, the system exhibits a strange attractor with a leading Lyapunov exponent $\lambda_1 \approx 0.91$, meaning nearby trajectories diverge exponentially:

$\delta(t) = \delta_0\, e^{\lambda t}$

For the real atmosphere, Lorenz (1969) further showed that even small-scale turbulent errors can propagate upscale, contaminating larger scales. This establishes a finite-time horizon for deterministic predictability of roughly 2-3 weeks, though the practical limit (useful skill) is reached at approximately 10 days for synoptic-scale weather patterns. This is why probabilistic (ensemble) forecasting is fundamentally necessary beyond a few days.

6. Synoptic Analysis & Pattern Recognition

Despite the power of numerical models, human pattern recognition and conceptual models remain essential for interpreting model output, communicating hazards, and catching model errors. Synoptic-scale analysis integrates observations, satellite imagery, radar data, and model fields into a coherent picture of the atmosphere.

6.1 Frontal Analysis

Fronts are boundaries between air masses with significantly different temperature and moisture characteristics. Frontal analysis identifies these boundaries and their type, motion, and intensity:

Cold Fronts

The leading edge of an advancing cold air mass, displacing warmer air. Characterized by:

- Sharp temperature drop (5-15 K) behind the front
- Wind shift (typically SW to NW in Northern Hemisphere)
- Pressure trough (minimum at frontal passage, then rapid rise)
- Steep slope (~1:50 to 1:100); narrow precipitation band
- Ana-cold fronts: Active, with ascending warm air above the front; heavy precipitation
- Kata-cold fronts: Subsiding warm air above; clearing skies, weak precipitation
- Frontal lifting produces narrow but intense convective bands; severe weather possible

Warm Fronts

The leading edge of advancing warm air, overriding cooler air ahead. Characterized by:

- Gradual temperature increase behind the front
- Gentle slope (~1:100 to 1:300); wide precipitation zone (200-400 km)
- Cloud sequence: Ci, Cs, As, Ns approaching from the west
- Steady precipitation (stratiform), possible embedded convection
- Fog and low stratus common in the warm sector
- Wind shift (typically SE to SW in NH)
- Warm-frontal inversions trap pollution and create poor visibility

Occluded Fronts

Form when the cold front overtakes the warm front, lifting the warm sector air off the surface:

- Cold occlusion: Air behind the cold front is colder than air ahead of the warm front; behaves like a cold front (common in continental climates)
- Warm occlusion: Air behind the cold front is warmer than air ahead of the warm front; behaves like a warm front (common over oceans)
- Often marks the mature-to-decaying phase of the cyclone life cycle
- Trowal (trough of warm air aloft) produces significant precipitation

Stationary Fronts

Fronts with negligible movement (surface wind component perpendicular to front < 5 kt on both sides). Persistent ascent along the boundary can produce extended periods of precipitation. Quasi-stationary fronts, particularly the Mei-Yu/Baiu front in East Asia, are responsible for prolonged heavy rainfall events and flooding. Frontogenesis function $F = -\frac{1}{|\nabla\theta|}\frac{d|\nabla\theta|}{dt}$ diagnoses the rate at which fronts are being strengthened or weakened.

6.2 Conveyor Belt Model of Extratropical Cyclones

The conveyor belt model (Carlson, 1980; Browning, 1990) describes the three-dimensional airflow through extratropical cyclones using three coherent airstreams:

Warm Conveyor Belt (WCB)

A broad, slantwise ascending airstream originating in the warm sector boundary layer near the surface warm front. Air ascends from ~1000 hPa to ~300 hPa over 1-2 days, producing the main cloud shield and precipitation of the cyclone. The WCB rises along and above the warm front, turns anticyclonically at upper levels, and forms the cirrus outflow. Latent heat release in the WCB is the primary energy source for rapid cyclogenesis. WCB trajectories transport moisture poleward and are the dominant precipitation-producing mechanism in extratropical cyclones.

Cold Conveyor Belt (CCB)

A low-level easterly flow north of the warm front, running beneath and opposite to the WCB. Originates in the cool air ahead of the cyclone and wraps cyclonically around the low center. The CCB produces low clouds and drizzle in the cold sector. Some branches ascend near the cyclone center, contributing to the "comma head" cloud pattern visible in satellite imagery. Strong CCB flow can produce damaging surface winds on the poleward side of the cyclone (sting jet).

Dry Intrusion (DI)

A descending airstream originating in the upper troposphere/lower stratosphere behind the cold front. Characterized by very low humidity, high potential vorticity, and high ozone concentration. The dry intrusion wraps cyclonically, creating the "dry slot" visible in water vapor imagery as a dark band curving around the cyclone. It enhances frontogenesis, sharpens the cold front, and contributes to the development of the cyclone's characteristic comma shape. Dry intrusions can bring stratospheric air to near-surface levels in the post-cold-frontal region.

6.3 Norwegian Cyclone Model

The Bergen School model (Bjerknes and Solberg, 1922) describes the life cycle of extratropical cyclones in terms of frontal wave development:

Stage 1 - Initial perturbation: A small wave develops on a stationary front (polar front) between polar and tropical air masses, triggered by upper-level forcing or terrain.

Stage 2 - Open wave: The wave amplifies with a warm front extending northeast and a cold front extending southwest from the developing low center. The warm sector narrows. Pressure falls rapidly at the center.

Stage 3 - Mature cyclone: The cold front, moving faster, begins overtaking the warm front. Central pressure reaches its minimum. Maximum precipitation and wind intensity.

Stage 4 - Occlusion: The warm sector is lifted off the surface. The occluded front wraps around the low center. The system becomes vertically stacked (cold core).

Stage 5 - Dissipation: The cyclone weakens as it becomes cut off from the baroclinic energy source. Fronts weaken and the system fills. May transition to a barotropic vortex.

While simplified, the Norwegian model remains the conceptual foundation for surface analysis and frontal placement. Modern extensions include the Shapiro-Keyser model (T-bone frontal structure common over oceans) and the inclusion of upper-level dynamics.

6.4 Jet Streak Dynamics (Four-Quadrant Model)

A jet streak is a localized wind maximum embedded within the jet stream. The four-quadrant model describes the ageostrophic circulation associated with jet streaks and its effect on vertical motion:

Right entrance region: Ageostrophic flow directed from the anticyclonic (equatorward) to the cyclonic (poleward) side. Upper-level convergence and subsidence. Surface high pressure tendency.

Left entrance region: Upper-level divergence and ascending motion. Favorable for cyclogenesis and precipitation. Surface low pressure tendency.

Right exit region: Upper-level divergence and ascending motion. Favorable for cyclogenesis. Often coupled with surface fronts and precipitation.

Left exit region: Upper-level convergence and subsidence. Suppressed weather. Surface high pressure tendency.

The ageostrophic wind in the entrance region is described by the isallobaric wind approximation:

$\vec{V}_{ag} = -\frac{1}{f}\hat{k} \times \frac{\partial \vec{V}_g}{\partial t} \approx -\frac{1}{f}\hat{k} \times (\vec{V}_g \cdot \nabla)\vec{V}_g$

This transverse ageostrophic circulation drives a thermally direct (entrance) or indirect (exit) secondary circulation that links upper-level and surface dynamics. Forecasters routinely overlay jet streak positions with surface features to diagnose regions of forcing for ascent.

6.5 Surface Analysis Techniques

Operational surface analysis integrates multiple data sources to produce a coherent depiction of the surface weather pattern:

Isobar analysis: Contours of sea-level pressure at 4 hPa intervals. Isobars should be smooth, parallel to the wind (slightly crossing toward low pressure due to friction), and packed tightly in regions of strong wind.

Frontal placement: Identify fronts using: (1) temperature gradient maxima; (2) wind shift lines; (3) pressure troughs; (4) dewpoint discontinuities; (5) cloud and precipitation patterns; (6) satellite and radar imagery. All criteria should be consistent.

Station model decoding: The station model plots temperature, dewpoint, wind (barbs), pressure, pressure tendency, weather, cloud cover, and cloud type in a standardized format around the station location.

Thickness analysis: The 1000-500 hPa thickness is proportional to the mean virtual temperature of the layer. The 540 dam (5400 m) line approximates the rain-snow boundary. Thickness advection patterns diagnose temperature advection.

Isentropic analysis: Analysis on constant potential temperature surfaces reveals the three-dimensional moisture transport (since air flows along isentropic surfaces in adiabatic flow). Moisture plumes ascending along isentropic surfaces indicate regions of precipitation.

Potential vorticity (PV) analysis: PV is conserved for adiabatic, frictionless flow. Upper-level PV anomalies (tropopause folds) induce surface cyclogenesis through the Bretherton (1966) PV inversion framework. The dynamic tropopause (2 PVU surface) is a powerful diagnostic tool.

Q-Vector Diagnosis of Vertical Motion

Q-vectors provide a rigorous diagnostic of quasigeostrophic vertical motion that avoids the cancellation problems of separate vorticity and temperature advection terms:

$\vec{Q} = -\frac{R}{p}\left(\frac{\partial \vec{V}_g}{\partial x}\cdot\nabla T,\; \frac{\partial \vec{V}_g}{\partial y}\cdot\nabla T\right)$

$\nabla^2\omega \propto -2\nabla\cdot\vec{Q}$

Convergence of Q-vectors ($\nabla \cdot \vec{Q} < 0$) implies forcing for ascent (upward vertical motion, $\omega < 0$). Divergence of Q-vectors implies forcing for descent. Q-vectors point toward the warm air in frontogenetic regions and are perpendicular to the thermal gradient in regions of pure vorticity advection.

Summary

Part V has covered the complete weather forecasting chain, from observations through analysis, prediction, uncertainty quantification, and verification:

1. The global observing system ingests ~100 million observations daily from surface stations, radiosondes, satellites, radar, aircraft, and GPS-RO, each with distinct strengths and error characteristics
2. Data assimilation (3D-Var, 4D-Var, EnKF, hybrid methods) optimally combines observations with model backgrounds, accounting for error covariances in both
3. NWP models solve the primitive equations using spectral, finite-volume, or spectral-element methods, with parameterizations for subgrid physics including convection, radiation, microphysics, and boundary layer turbulence
4. Ensemble forecasting quantifies uncertainty from initial condition and model errors using perturbed forecasts, providing probabilistic guidance essential beyond day 3-4
5. Forecast verification using scores (RMSE, ACC, Brier Score, CRPS, ROC) and reliability diagrams objectively assesses forecast quality and guides model development
6. Synoptic analysis using conceptual models (Norwegian cyclone model, conveyor belts, jet streak dynamics, Q-vectors) remains essential for interpreting model output and communicating weather hazards

The predictability horizon continues to advance: 500 hPa anomaly correlation useful skill has improved from ~5 days in 1980 to ~9 days today, driven by better observations, improved data assimilation, higher model resolution, and enhanced physics parameterizations. Machine learning is emerging as a complementary approach, with models like GraphCast and Pangu-Weather showing competitive skill at a fraction of the computational cost.

The following parts will explore Earth's climate history (paleoclimatology), extreme weather events, and the science of climate change.