39 Generative Models for Weather Prediction

“The atmosphere is a chaotic symphony of ${10}^{44}$ molecules, yet we dare predict its song days in advance.”

On August 29, 2005, Hurricane Katrina made landfall near Buras, Louisiana, as a Category 3 storm with sustained winds of 205 kilometres per hour. The storm surge, peaking at nine metres along the Mississippi coast, overwhelmed the levée system protecting New Orleans. Roughly 1,800 people lost their lives; economic damage exceeded $125 billion, making Katrina the costliest natural disaster in United States history up to that point. Yet the storm itself was not a surprise. Five days before landfall, the European Centre for Medium-Range Weather Forecasts (ECMWF) model had already shown a coherent Gulf of Mexico track. The US National Hurricane Center issued watches and warnings 36 hours ahead-well within the window for organised evacuation. What failed was not weather science; it was the translation of probabilistic forecast information into decisive action.

Nearly two decades later, in September 2024, Hurricane Hel\-ene swept across Florida's Big Bend coast, then tracked hundreds of kilometres inland before devastating the Appalachian communities of western North Carolina. Asheville, far from any coastline and nestled between mountain ridges, received over half a metre of rain in 48 hours; the French Broad River crested at a record 7.3 metres above flood stage. More than 200 people died. Yet the operational forecast guidance-from both traditional Numerical Weather Prediction (NWP) models and the new generation of machine-learning weather models, including GraphCast and Pangu-Weather-had correctly identified the extreme rainfall potential days in advance. The fundamental scientific challenge of predicting these events had, in large measure, been solved.

The story of Katrina and Helene illustrates the astonishing progress of weather forecasting over the past half-century, and simultaneously reveals where the remaining gaps lie. Forecasting the track of a major hurricane is no longer the primary challenge; forecasting the distribution of possible impacts-the tail of the rainfall distribution, the compound flooding risk, the probability of any given neighbourhood seeing sustained winds above 50 metres per second-remains difficult, consequential, and, as we shall see, precisely the problem for which generative models offer transformative potential.

This chapter tells the story of how generative artificial intelligence entered one of the most practically important and scientifically demanding prediction problems humanity faces: the state of the atmosphere. We develop the mathematical foundations, survey the landmark models, analyse the theoretical guarantees and known failure modes, and situate weather prediction within the broader landscape of this book. The chapter connects directly to the physics-constrained PDE modelling of Chapter 35 and draws on the interpretability tools of Chapter 38; it will also prepare the ground for subsequent chapters on geophysical and climate modelling.

The Challenge of Atmospheric Prediction

Why Weather Forecasting Matters

Weather is not merely a backdrop to human activity; it is one of the principal drivers of the global economy. A landmark 2009 analysis by the United States Weather Research Programme estimated that weather-sensitive economic sectors account for roughly one third of US GDP, approximately $3 trillion annually. Scaling globally, the weather-sensitive share of the world economy exceeds $5 trillion per year.

The channels through which weather affects economic value are diverse and sometimes surprising:

• Agriculture. Approximately 80% of the variability in US crop yields is attributable to weather, principally temperature and precipitation anomalies during sensitive phenological windows. A three-day forecast of a late frost can allow a Florida citrus grower to activate smudge pots and save an entire harvest. Sub-seasonal to seasonal forecasts of El Niño–Southern Oscillation (ENSO) events allow commodity traders, insurance underwriters, and government food security programmes to hedge risks months in advance.

• Aviation. The US Federal Aviation Administration estimates that weather causes approximately 70% of all National Airspace System delays, costing the airline industry over $4 billion per year in direct operating costs. Turbulence-related injuries cause more than 30 severe injuries annually in the US alone; clear-air turbulence, which leaves no visual signature and cannot be detected by airborne radar, is a particular concern as a warming atmosphere alters upper-level wind shear patterns.

• Energy. The rapid expansion of wind and solar generation has made weather forecasting central to grid operation. Electricity demand peaks with temperature extremes; solar output varies with cloud cover on timescales of minutes; offshore wind turbines must be shut down in extreme conditions to prevent mechanical damage. A 1% improvement in day-ahead wind power forecast accuracy across the European grid has been estimated to save $\sim$$100 million annually in balancing costs.

• Disaster preparedness. The World Meteorological Organisation estimates that improved early warning systems for weather-related hazards could save 2 billion euros per year in disaster losses across Europe alone. Globally, the mortality from tropical cyclones has decreased dramatically since the 1970s-from roughly 3 deaths per 100,000 population to under 0.5-primarily due to improved track forecasts and warning lead times. Economic losses, however, continue to grow as population and infrastructure expand into vulnerable coastal zones.

A Philosophical Hook: Predicting Chaos

The atmosphere contains approximately ${10}^{44}$ molecules. Their collective behaviour is governed by the laws of classical mechanics and thermodynamics, combined with the transfer of radiation, the microphysics of cloud droplet nucleation, and the turbulent mixing of the planetary boundary layer. No computer will ever simulate this system at the molecular scale; even a description at the kilometre scale requires resolving $\sim {10}^9$ grid cells, each carrying dozens of physical variables, over thousands of time steps.

And yet, we predict. Not perfectly-the atmosphere is chaotic, and perfect prediction is provably impossible beyond a finite horizon-but usefully. A modern medium-range forecast for surface temperature is accurate to within about 2 degrees Celsius at seven days lead time, a level of skill that would have seemed miraculous to the meteorologists of the 1960s.

The philosophical tension is acute: the atmosphere obeys deterministic physical laws (the Navier–Stokes equations, thermodynamics, Maxwell's equations for radiation), yet its long-term behaviour is unpredictable in any practical sense. How do we navigate between the rock of determinism and the hard place of chaos?

The answer, as Edward Lorenz showed in 1963, is that determinism and unpredictability are not mutually exclusive. Tiny differences in initial conditions-differences smaller than any measurement can detect-grow exponentially in time, eventually overwhelming the signal. This is the celebrated butterfly effect: the metaphorical flap of a butterfly's wings in Brazil can, in principle, alter whether a tornado forms in Texas three weeks hence. The practical implication is that there exists a fundamental predictability horizon beyond which deterministic single-valued forecasts become useless, and the only honest statement is a probability distribution over possible atmospheric states.

This insight, which we will formalise below, is the philosophical foundation for everything that follows in this chapter. Weather forecasting is, at its core, a problem of learning and sampling from a conditional probability distribution-precisely the kind of problem generative models are designed to solve.

Historical Development: From Manual Computation to Artificial Intelligence

The timeline of these developments is shown in Figure fig:weather:timeline.

Lorenz and the Mathematics of Chaos

In 1963, Edward Lorenz published what would become one of the most cited papers in all of science: Deterministic Nonperiodic Flow, in the Journal of the Atmospheric Sciences. Working on a simplified model of atmospheric convection, Lorenz stumbled upon behaviour that would transform our understanding of determinism.

His model consisted of three coupled ordinary differential equations: (Lorenz1)$$\begin{array}{cccc}& \frac{dx}{dt}& =\sigma (y-x),& \text{<span class="tml-eqn"></span>}\\ & \frac{dy}{dt}& =x(\rho -z)-y,& \text{<span class="tml-eqn"></span>}\\ & \frac{dz}{dt}& =xy-\beta z,& \text{<span class="tml-eqn"></span>}\end{array}$$ with canonical parameters $\sigma =10$, $\rho =28$, $\beta =8/3$. These equations exhibit what Lorenz called sensitive dependence on initial conditions: two trajectories starting arbitrarily close together diverge exponentially, with a characteristic Lyapunov exponent ${\lambda }_1\approx 0.9$ per nondimensional time unit of the system.

The practical implication for weather forecasting is profound. If the leading Lyapunov exponent of the atmosphere corresponds to an error-doubling time of two days, and if we cannot measure initial conditions with infinite precision, then beyond some finite horizon the forecast error will be comparable to the climatological variance-at which point the forecast has no skill.

Formal Definitions: The Atmospheric State and Forecast Operator

We now formalise the mathematical framework for weather prediction, establishing notation that will be used throughout the chapter.

The Operational Bottleneck: Compute Cost and the NWP Scaling Problem

Modern operational NWP represents one of the most computationally demanding regular tasks in all of science. ECMWF's Integrated Forecasting System (IFS), which produces the world's leading medium-range global forecasts, runs on a dedicated supercomputer cluster with a peak performance of approximately ${10}^{17}$ floating-point operations per second (100 petaFLOPS). A single ten-day global forecast at $9$ km resolution requires roughly ${10}^{20}$ floating-point operations-approximately ${10}^7$ times more than a single forward pass through a large language model.

The computational cost of NWP scales approximately as: (NWP COST)$$\begin{array}{ccc}& C_{\mathrm{NWP}}\sim N_x\cdot N_y\cdot N_z\cdot N_t,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $N_x\times N_y$ is the horizontal grid resolution, $N_z$ is the number of vertical levels, and $N_t$ is the number of time steps. Halving the horizontal grid spacing increases cost by a factor of $2^3=8$ (two factors in horizontal resolution, one in time step required for numerical stability by the Courant–Friedrichs–Lewy condition).

The Quiet Revolution: Skill Improvement in Traditional NWP

The sixty-year history of operational NWP reveals a remarkable and sustained trend in forecast skill: roughly one day of useful forecast range gained per decade. This improvement is most clearly measured by the anomaly correlation coefficient (ACC) of the 500 hPa geopotential height field-a standard metric that compares the forecast anomaly (deviation from climatology) with the verifying analysis. An ACC above 0.6 is conventionally considered the threshold of useful skill.

In 1980, ECMWF forecasts crossed the ACC,$>0.6$ threshold at approximately 5 days lead time. By 2020, the threshold had extended to approximately 8.5 days. This represents a gain of 3.5 days over four decades. Equivalently, the quality of today's seven-day forecast is roughly comparable to the quality of a three-day forecast in 1980. The improvement reflects advances in:

• Data assimilation: Moving from optimal interpolation to four-dimensional variational assimilation (4D-Var), which optimally blends observations with model trajectories over a 6–12 hour assimilation window.

• Model resolution: From 150 km in 1979 to 9 km in 2023, enabling better representation of orographic and mesoscale processes.

• Physical parameterisations: Improved schemes for convection, boundary layer turbulence, radiation, and cloud microphysics.

• Observations: The global observing network has grown from $\sim {10}^4$ surface stations in 1980 to $\sim {10}^7$ observations per day from the current satellite constellation.

Exercises for Section The Challenge of Atmospheric Prediction

Why Generative AI for Weather?

The previous section established that weather forecasting is, at its mathematical core, a problem of learning and sampling from a high-dimensional conditional probability distribution. It also established that the classical NWP paradigm, while enormously successful, faces a cost wall: operational ensemble forecasting with 50 members at high resolution consumes resources available only to the largest national meteorological centres. In this section we examine more precisely why generative models offer a compelling solution to these challenges, and how they relate to the broader taxonomy of machine learning approaches to weather prediction.

The Three Pillars: Cost, Ensembles, and Extremes

The case for generative models in weather prediction rests on three mutually reinforcing pillars.

Pillar 1: Computational Cost

As established in Proposition Proposition 1, the cost of a single NWP forecast is enormous, and the cost of an ensemble forecast is this number multiplied by the ensemble size. A 50-member ensemble at $9$ km global resolution, run for a 10-day lead time, requires on the order of ${10}^{22}$ floating-point operations-enough to occupy the ECMWF supercomputer for several hours. Only a handful of national centres worldwide have the resources to run such systems operationally.

By contrast, a trained generative model produces ensemble members by forward sampling, which requires a single forward pass through a neural network (or, for diffusion models, a sequence of $T_{\mathrm{diff}}\approx 20$–100 such passes). The total inference cost is: (GEN COST)$$\begin{array}{ccc}& C_{\mathrm{gen}}=N_{\mathrm{ens}}\times T_{\mathrm{diff}}\times C_{\mathrm{fwd}},& \text{<span class="tml-eqn"></span>}\end{array}$$ where $C_{\mathrm{fwd}}$ is the cost of a single neural network forward pass. For a model like GenCast or SEEDS, this amounts to roughly ${10}^{15}$ operations per 50-member ensemble-seven orders of magnitude less than the NWP equivalent.

Pillar 2: Ensemble Generation and Calibration

Operational ensemble prediction systems generate spread by perturbing initial conditions (using singular vectors or ensemble Kalman filter perturbations) and sometimes by stochastic model perturbations. These perturbations are designed to sample the leading error modes of the system, but they are necessarily limited in their ability to represent the full uncertainty distribution.

A well-trained generative model, by contrast, learns the full conditional distribution $p({𝒙}_{t+\mathrm{\Delta }t}|{𝒙}_t)$ from data. If the training data is sufficiently rich and diverse, the model can in principle capture uncertainty sources that are difficult to parameterise explicitly in physical models: mesoscale convective organisation, land-surface feedback uncertainty, orographic precipitation spread, and more. This offers the prospect of better-calibrated ensembles-ensembles whose spread more faithfully represents the true forecast uncertainty.

Generative models offer a systematic path to calibration improvement through post-processing: rather than running the expensive NWP ensemble and then applying bias correction, one can train a conditional generative model to map a deterministic (or low-resolution) NWP forecast directly to a calibrated ensemble distribution. This is the ensemble emulation paradigm explored in GenCast, NowcastNet, and similar systems.

Pillar 3: Extreme Events

Extreme weather events-heat waves, atmospheric rivers, tropical cyclones, flash floods-are precisely the events for which forecast skill matters most, both in terms of economic impact and loss of life. They are also, almost by definition, the hardest to forecast.

Traditional ML weather models trained with mean squared error (MSE) loss face a fundamental problem with extremes: MSE minimisation encourages predictions that are close to the conditional mean, which for multimodal distributions is not a typical sample. A model predicting whether a storm will track north or south of a coastal city, each with 50% probability, should predict one or the other-not an average track that lies in the middle of the coast and is unlikely to verify. MSE-optimal forecasts are necessarily smooth, tending to regress toward climatology for uncertain events. This is the blurring problem.

Generative models, by learning to sample from the full conditional distribution, can in principle produce sharp forecasts of extreme events. A diffusion model conditioned on a predicted atmospheric state can generate individual sample trajectories that include coherent, physically realistic extreme events, even if those events appear in only a fraction of the ensemble members. The probability of the event is estimated from the ensemble frequency.

Deterministic MLWP: Limitations and Where Generative Models Help

The first generation of modern ML weather models-FourCastNet (2022), Pangu-Weather (2023), GraphCast (2023)-are deterministic systems that directly predict the atmospheric state at the next time step. They have demonstrated remarkable skill, matching or exceeding the ECMWF IFS for many variables at 7–10 day lead times while running at a fraction of the computational cost.

However, their fundamental architecture-a neural network trained to minimise MSE between predicted and verifying states-imposes a systematic bias: the regression to the mean problem. Consider the formal decomposition: (MSE Solution)$$\begin{array}{cccc}& {\widehat{𝒙}}_{t+\mathrm{\Delta }t}^{\mathrm{MSE}}& =𝔼[{𝒙}_{t+\mathrm{\Delta }t}|{𝒙}_t]& \text{<span class="tml-eqn"></span>}\\ & & =\int {𝒙}^{\prime }p({𝒙}^{\prime }|{𝒙}_t)d{𝒙}^{\prime },& \text{<span class="tml-eqn"></span>}\end{array}$$ where the equality in holds because MSE minimisation yields the conditional expectation. For multimodal distributions $p({𝒙}^{\prime }|{𝒙}_t)$-which are common in the atmosphere, especially at longer lead times when the flow can bifurcate into qualitatively different regimes-the conditional mean is not a typical sample. It lies between the modes, in a region of low probability density.

Figure fig:weather:det_vs_gen illustrates the architectural difference between deterministic and generative MLWP. The key structural insight is that the conditional mean-the optimal deterministic forecast under MSE loss-is not the same as the most likely forecast, and for multimodal distributions it may lie in a region of essentially zero probability.

Ensemble Prediction Systems and the Cost of NWP Ensembles

The state-of-the-art in operational probabilistic forecasting is the ensemble prediction system (EPS). The flagship systems are:

• ECMWF ENS: 51 members (1 control + 50 perturbed) at $\sim 18$ km resolution, 15-day lead time. Initial perturbations generated by singular vectors and ensemble data assimilation; model uncertainty represented by stochastic parameterisation schemes.

• NOAA GEFS (Global Ensemble Forecast System): 31 members at $\sim 25$ km, 16-day lead time; being upgraded to higher resolution under the Next Generation Global Prediction System (NGGPS) programme.

• UK Met Office MOGREPS-G: 18 members at $\sim 20$ km; uses the ensemble transform Kalman filter (ETKF) for perturbation generation.

• Multi-model ensembles: The THORPEX Interactive Grand Global Ensemble (TIGGE) archive collects outputs from all major operational EPS systems, enabling multi-model ensemble combinations that typically outperform any individual system.

The computational cost of running 50 ensemble members is, approximately, 50 times the cost of a single deterministic forecast. For the ECMWF system in 2023, this amounts to approximately 30 megawatt-hours of electricity per day-comparable to the annual electricity consumption of 3,000 European households. The energy and financial cost of ensemble NWP is a fundamental barrier to democratising access to probabilistic weather services: most national meteorological services in developing countries cannot afford to run their own ensemble systems and rely on the ECMWF or NOAA outputs.

Formal Framework: Weather Forecasting as Conditional Generative Modelling

We now place weather forecasting squarely within the conditional generative modelling framework developed in earlier chapters of this book.

This framework immediately suggests the applicability of the main generative model families covered in this book:

• Diffusion models (Chapter 19): Learn $p_{\mathrm{\Delta }t}$ by training a denoising score network $s_{𝜽}(\cdot ,\sigma )\approx {\nabla }_{𝒙}\log p_{\mathrm{\Delta }t}(𝒙|{𝒙}_t)$. Sampling proceeds via Langevin dynamics or DDIM. Key models: GenCast, CorrDiff, EarthFormer-Diffusion.

• Normalising flows (Chapter 18): Model $p_{\mathrm{\Delta }t}$ as a pushforward of a standard Gaussian through an invertible map $f_{𝜽}$. Exact likelihoods enable principled calibration assessment. Key models: RainFall (Gong et al., 2023), NF-based post-processing (Groenke et al., 2022).

• Variational Autoencoders (Chapter 15): A latent-variable model for atmospheric state that enables compressed representations and fast sampling. Key models: WeatherVAE, GenWP.

• GANs (Chapter 9): Generative Adversarial Networks for downscaling and super-resolution of weather fields. Key models: DeepSD (Stengel et al., 2020), DeepMind's downscaling GAN, ESRGAN-based approaches.

A Taxonomy of ML Weather Prediction Approaches

The landscape of ML weather prediction has grown dramatically since 2020. To navigate it, we introduce a systematic taxonomy along four axes: (1) model type (deterministic vs. generative), (2) spatial scope (global vs. regional/limited-area), (3) temporal architecture (autoregressive vs. direct multi-step), and (4) target variable (full atmosphere vs. single-variable).

Table Table 1 provides a structured overview of the most influential models along these axes.

Model	Type	Scope	Temporal	Res.	Key Feature
Deterministic Global Models
FourCastNet	Det.	Global	AR	${0.25}^{\circ }$	Spherical FNO
Pangu-Weather	Det.	Global	Direct	${0.25}^{\circ }$	3D Earth transformer
GraphCast	Det.	Global	AR	${0.25}^{\circ }$	Graph neural network
FengWu	Det.	Global	AR	${0.25}^{\circ }$	Cross-modal fusion
Aardvark	Det.	Global	AR	${0.25}^{\circ }$	Sparse observations in
Generative Global Models
GenCast	Diffusion	Global	AR	${0.25}^{\circ }$	Calibrated ensembles
SEEDS	Diffusion	Global	Emulator	$1^{\circ }$	EPS emulation
PreDiff	Diffusion	Regional	AR	$4$ km	Nowcasting
CorrDiff	Diffusion	Regional	AR	$3$ km	Downscaling
EarthFormer	Diff+Attn	Regional	AR	$1^{\circ }$	Earth attention
Nowcasting and Short-Range
NowcastNet	GAN	Regional	Direct	$1$ km	Extrapolation
Skilful Prec.	GAN	Regional	Direct	$1$ km	DeepMind GAN
DGM-R	NF	Regional	AR	$1^{\circ }$	Radar nowcast
Downscaling / Super-Resolution
DeepSD	GAN	Regional	Direct	$3$ km	Statistical SD
CorrDiff	Diffusion	Regional	Direct	$3$ km	Score-based SD
GenRain	Flow	Regional	Direct	$1$ km	Probabilistic

Chapter Roadmap and Dependencies

This chapter is structured to develop the material from foundational principles through state-of-the-art models and open problems. Figure fig:weather:roadmap shows the dependency structure between sections.

Exercises for Section Why Generative AI for Weather?

The Governing Equations of the Atmosphere

“The atmosphere is a nonlinear fluid on a rotating sphere, forced by solar radiation and interacting with the ocean, land, and biosphere. Its equations are simple to write, impossible to solve exactly, and extraordinarily difficult to approximate well.” - Edward N. Lorenz (paraphrased)

Weather prediction is, at its mathematical core, an initial value problem for a system of nonlinear partial differential equations. The atmosphere is a thin shell of gas-roughly 100,km thick over a planet of radius $\mathrm{6,371}$,km-set in rotation at one revolution per day, heated unevenly from below by sunlight, and stirred by the irregular topography of continents and mountain ranges. Despite this complexity, the fundamental physical laws governing atmospheric motion are well known. They form the primitive equations, a closed system of PDEs derived from Newtonian mechanics, thermodynamics, and the ideal gas law.

In this section we derive and discuss these equations in detail, connect them to the general Navier–Stokes framework introduced in Chapter 35 (Section The Navier--Stokes equations: beauty and terror), and explain why, despite being fully known, they resist exact solution and force us toward numerical approximation.

Coordinate systems and notation

The Earth's atmosphere is most naturally described in spherical coordinates $(\lambda ,\varphi ,r)$, where $\lambda \in [0,2\pi )$ is longitude, $\varphi \in [-\pi /2,\pi /2]$ is latitude, and $r=a+z$ is the radial distance from Earth's centre with mean radius $a\approx \mathrm{6,371}\text{km}$ and height above the surface $z\ge 0$. In practice, NWP models use a hybrid pressure–terrain-following coordinate $\eta$ or $\sigma$, which simplifies the treatment of the lower boundary.

We denote the three-dimensional velocity field as $𝐮=(u,v,w)$, where $u$ is the zonal (east–west) component, $v$ is the meridional (north–south) component, and $w=\mathrm{d}z/\mathrm{d}t$ is the vertical velocity. For notational simplicity in the following subsections, we work in local Cartesian coordinates on the sphere (the tangent plane approximation), which is exact in the limit of small horizontal scales and standard in textbook derivations.

The material derivative (Lagrangian rate of change following a fluid parcel) is (Material)$$\begin{array}{ccc}& \frac{D}{Dt}\equiv \frac{\partial }{\partial t}+𝐮\cdot \nabla =\frac{\partial }{\partial t}+u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}+w\frac{\partial }{\partial z},& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\nabla$ denotes the three-dimensional gradient operator. The distinction between Eulerian ($\partial /\partial t$ at a fixed point) and Lagrangian ($D/Dt$ following a parcel) perspectives is essential throughout atmospheric dynamics.

The momentum equation: Navier–Stokes on a rotating sphere

Newton's second law for a fluid parcel in an inertial frame gives the Navier–Stokes equations (Chapter 35, Section The Navier--Stokes equations: beauty and terror). Transforming to a reference frame rotating with the Earth at angular velocity $𝛀=\mathrm{\Omega }\widehat{𝐤}$ (where $\mathrm{\Omega }=7.292\times {10}^{-5}{\mathrm{s}}^{-1}$ and $\widehat{𝐤}$ points toward the North Pole) introduces two pseudo-forces: the Coriolis acceleration $-2𝛀\times 𝐮$ and the centrifugal acceleration $-𝛀\times (𝛀\times 𝐫)$. The centrifugal term is absorbed into a modified gravitational acceleration $𝐠$ (pointing effectively downward, of magnitude $g\approx 9.81{\mathrm{m}\mathrm{s}}^{-2}$). The result is the atmospheric momentum equation: (Momentum)$$\begin{array}{ccc}& \frac{D𝐮}{Dt}+2𝛀\times 𝐮=-\frac{1}{\rho }\nabla p+𝐠+𝐅,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $p(𝐱,t)$ is pressure, $\rho (𝐱,t)$ is air density, $𝐠=-g\widehat{𝐤}$ is effective gravity, and $𝐅$ represents frictional and turbulent drag forces (which act primarily in the planetary boundary layer and through parameterised gravity-wave drag).

The term $2𝛀\times 𝐮$ is the Coriolis acceleration. In local coordinates at latitude $\varphi$, defining the Coriolis parameter $f=2\mathrm{\Omega }\sin \varphi$, the horizontal momentum equations take the familiar form: (Momentum U)$$\begin{array}{cccc}& \frac{Du}{Dt}-fv& =-\frac{1}{\rho }\frac{\partial p}{\partial x}+F_x,& \text{<span class="tml-eqn"></span>}\\ & \frac{Dv}{Dt}+fu& =-\frac{1}{\rho }\frac{\partial p}{\partial y}+F_y,& \text{<span class="tml-eqn"></span>}\\ & \frac{Dw}{Dt}& =-\frac{1}{\rho }\frac{\partial p}{\partial z}-g+F_z.& \text{<span class="tml-eqn"></span>}\end{array}$$ At mid-latitudes ($\varphi \approx {45}^{\circ }$), $f\approx {10}^{-4}{\mathrm{s}}^{-1}$, so the Coriolis force is significant for synoptic-scale motions (cyclones, anticyclones) but negligible at small scales (thunderstorms, boundary layer turbulence).

The continuity equation: conservation of mass

Mass conservation for a compressible fluid yields the continuity equation: (Continuity)$$\begin{array}{ccc}& \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho 𝐮)=0,& \text{<span class="tml-eqn"></span>}\end{array}$$ or equivalently, using the material derivative: (Continuity Lagrange)$$\begin{array}{ccc}& \frac{D\rho }{Dt}+\rho \nabla \cdot 𝐮=0.& \text{<span class="tml-eqn"></span>}\end{array}$$ This states that the density of an air parcel changes only through divergence of the velocity field. A convergent flow ($\nabla \cdot 𝐮<0$) compresses air parcels and increases their density.

In the Boussinesq approximation, valid for small density variations (useful for ocean dynamics but less so for the full atmosphere), $\rho$ is treated as constant and one recovers the incompressibility condition $\nabla \cdot 𝐮=0$. The full atmosphere requires the compressible form.

The thermodynamic energy equation

The first law of thermodynamics applied to a fluid parcel gives the thermodynamic energy equation, governing temperature evolution: (Thermo)$$\begin{array}{ccc}& c_p\frac{DT}{Dt}-\frac{1}{\rho }\frac{Dp}{Dt}=Q,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $T(𝐱,t)$ is temperature, $c_p\approx 1005\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{K}}^{-1}$ is the specific heat at constant pressure, and $Q$ is the diabatic heating rate per unit mass (units $\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{s}}^{-1}$ $=\mathrm{W}\mathrm{k}{\mathrm{g}}^{-1}$). The diabatic heating $Q$ encompasses:

• Radiative heating: absorption of solar shortwave radiation and emission/absorption of terrestrial longwave radiation;

• Latent heat release: phase changes of water (condensation releasing heat, evaporation absorbing heat);

• Sensible heat fluxes: exchange of heat between the atmosphere and the underlying surface;

• Diffusion and turbulent mixing: small-scale transport of heat by unresolved eddies.

When $Q=0$, the process is adiabatic: temperature changes only through compression and expansion as the parcel moves vertically in the pressure gradient.

Using the ideal gas law (Section The equation of state: the ideal gas law) to eliminate $Dp/Dt$, the thermodynamic equation can be rewritten in terms of potential temperature $\theta =T{(p_0/p)}^{R_d/c_p}$, where $p_0=1000\mathrm{hPa}$ is a reference pressure and $R_d=287\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{K}}^{-1}$ is the gas constant for dry air: (Pottemp)$$\begin{array}{ccc}& \frac{D\theta }{Dt}=\frac{\theta }{c_{p}T}Q.& \text{<span class="tml-eqn"></span>}\end{array}$$ Potential temperature $\theta$ is conserved following an adiabatic parcel and is the natural prognostic variable for the thermodynamic equation in NWP models.

The moisture equation

Water in the atmosphere exists in three phases (vapour, liquid, ice) and plays a crucial role in weather through cloud formation, precipitation, and latent heat release. The prognostic equation for specific humidity $q$ (mass of water vapour per unit mass of moist air, $\mathrm{k}\mathrm{g}\mathrm{k}{\mathrm{g}}^{-1}$) is: (Moisture)$$\begin{array}{ccc}& \frac{Dq}{Dt}=S_q,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $S_q$ is the net source/sink term including:

• condensation and evaporation from cloud microphysics;

• precipitation (removing liquid water from the column);

• surface evaporation (adding moisture at the lower boundary);

• turbulent transport within the boundary layer.

In full models, separate equations are carried for cloud liquid water $q_{\ell }$, cloud ice $q_i$, rain $q_r$, and snow $q_s$, each with its own source/sink parameterisation. The latent heat release when water vapour condenses directly feeds back into the thermodynamic equation via $Q=-L_v{S}_q$, where $L_v\approx 2.5\times {10}^6\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}$ is the latent heat of vaporisation.

The equation of state: the ideal gas law

Air behaves to excellent approximation as an ideal gas. The equation of state closes the system by relating pressure, density, and temperature: (EOS)$$\begin{array}{ccc}& p=\rho R_d{T}_v,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $R_d=287\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{K}}^{-1}$ is the gas constant for dry air and $T_v=T(1+(R_v/R_d-1)q)\approx T(1+0.608q)$ is the virtual temperature, which accounts for the lower molecular weight of water vapour ($R_v=461\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{K}}^{-1}$). For dry air ($q\approx 0$), $T_v\approx T$ and the familiar ideal gas law $p=\rho R_{d}T$ applies.

The primitive equations as a formal system

We can now state the complete governing system.

Conservation properties

The primitive equations conserve fundamental physical quantities in the absence of diabatic forcing and friction. This is a cornerstone of both the theoretical analysis and the numerical implementation of atmospheric models.

The hydrostatic approximation

The atmosphere's extreme aspect ratio-horizontal scales of thousands of kilometres versus vertical scales of tens of kilometres-implies that vertical accelerations $Dw/Dt$ are negligibly small compared to the dominant balance between gravity and the vertical pressure gradient. This yields the hydrostatic approximation, obtained by setting $Dw/Dt=0$ in Eq. : (Hydrostatic)$$\begin{array}{ccc}& \frac{\partial p}{\partial z}=-\rho g.& \text{<span class="tml-eqn"></span>}\end{array}$$ This single equation replaces the vertical momentum equation, reducing the prognostic system by one equation. Hydrostatic models are appropriate for horizontal scales larger than $\sim 10$,km; at smaller scales (cloud-resolving models, large-eddy simulations), the full non-hydrostatic equations are required.

Combined with the shallow atmosphere approximation-treating the atmosphere as a thin shell so that $r\approx a$ throughout-one arrives at the equations used in virtually all global operational NWP models. Departures from the shallow atmosphere assumption become non-negligible only at very high altitudes.

Boundary conditions

The primitive equations require boundary conditions at two surfaces:

Lower boundary (Earth's surface).

At the surface $z=z_s(x,y)$ (terrain height), the kinematic boundary condition requires that no air passes through the surface: (BC Lower)$$\begin{array}{ccc}& w=\frac{\partial z_s}{\partial t}+u\frac{\partial z_s}{\partial x}+v\frac{\partial z_s}{\partial y}\text{at }z=z_s.& \text{<span class="tml-eqn"></span>}\end{array}$$ For a stationary surface ($\partial z_s/\partial t=0$), this says that vertical velocity equals the slope of the terrain times the horizontal wind-air is deflected over mountains. Thermodynamic and moisture boundary conditions prescribe surface fluxes of heat and water vapour, computed from land-surface and ocean models.

Upper boundary (top of atmosphere).

At the model top $z=z_{\mathrm{top}}$, the pressure is set to a small value ($p_{\mathrm{top}}\approx 0.1$–$1\mathrm{hPa}$, corresponding to the stratosphere or mesosphere). A sponge layer absorbs spurious wave reflections near the model top.

Atmospheric structure

The primitive equations govern motion throughout the atmosphere, but different vertical layers have distinct dynamical regimes. Figure fig:weather:layers shows the standard atmospheric stratification.

Why direct numerical simulation is impossible

Parameterisation of sub-grid processes

Because a global NWP model operates on a grid with typical spacing $\mathrm{\Delta }x\sim 5$–$50\mathrm{km}$, all physical processes occurring at smaller scales must be parameterised: represented through equations relating their aggregate effect to the resolved grid-scale variables. The main parameterisation schemes in a state-of-the-art NWP model are:

Radiation.

Atmospheric radiation involves absorption and emission across millions of spectral lines of ${\mathrm{H}}_2\mathrm{O}$, $\mathrm{C}{\mathrm{O}}_2$, ${\mathrm{O}}_3$, and other gases. Full line-by-line calculations are computationally prohibitive for NWP; instead, two-stream approximations or correlated-$k$ methods provide accurate broadband fluxes at manageable cost. Radiation parameterisation is among the most computationally expensive components of a GCM.

Convection.

Thunderstorms and moist convection organise at scales of 1–10,km, far below the grid scale of global models. Mass flux or entraining plume schemes represent the effects of deep and shallow convection on the column temperature and moisture profiles. Errors in convective parameterisation are among the leading sources of forecast error in the tropics.

Cloud microphysics.

The formation and growth of cloud droplets, ice crystals, and precipitation particles involves processes from nanometre to millimetre scales. Bulk microphysics schemes represent clouds through moment equations for the distribution of hydrometeor sizes; spectral bin schemes are more accurate but far more expensive.

Planetary boundary layer.

The lowest 1–2,km of the atmosphere is the planetary boundary layer (PBL), where turbulent eddies transfer momentum, heat, and moisture between the surface and the free atmosphere. PBL schemes typically use mixing-length or turbulence kinetic energy (TKE) closures to represent this transport.

Gravity wave drag.

Sub-grid-scale orographic and convective gravity waves transport momentum vertically, decelerating the large-scale flow when they break. This drag is critical for the correct simulation of the stratospheric circulation and must be parameterised in global models.

The parameterisation problem is not merely a technical inconvenience-it is a fundamental scientific challenge. Errors in parameterisation account for a large fraction of model bias and uncertainty, particularly in climate projections. Machine learning approaches to learned parameterisation (replacing hand-tuned schemes with data-driven approximations) are an active frontier, discussed in later sections.

Historical background: from Richardson to ENIAC

The idea that weather could be predicted by solving the equations of fluid motion numerically dates to Lewis Fry Richardson (1881–1953). In his landmark 1922 book Weather Prediction by Numerical Process, Richardson described in extraordinary detail a scheme for solving the primitive equations on a grid, using a network of 64,000 human computers working in parallel in a vast “forecast factory”. His own attempt at a 6-hour forecast for a point in central Europe-carried out by hand during World War I-took approximately six weeks to compute and produced a catastrophically wrong answer (a pressure change of 145,hPa, when the observed change was near zero) due to an inadequate initialisation procedure.

Richardson's vision proved prophetic. The missing ingredient was data assimilation: a method for initialising the model state consistently with the underlying equations. Without it, the initial conditions contain fast-moving acoustic and gravity waves that swamp the slower weather-bearing Rossby waves, causing the forecast to explode.

The breakthrough came in 1950, when Jule Charney, Ragnar Fjortoft, and John von Neumann produced the first successful numerical weather forecast on the ENIAC computer at the Aberdeen Proving Ground in Maryland (the team was based at the Institute for Advanced Study in Princeton). Working with a simplified barotropic vorticity equation (a single-level filter that eliminates the problematic fast waves), they ran 24-hour forecasts for the 500,hPa geopotential height field over North America. The computation took 24 hours of ENIAC time to produce a 24-hour forecast-a ratio of 1:1. Today's operational models run at ratios of order ${10}^3$:1 (about ${10}^3$ seconds of forecast produced per second of wall-clock time), reflecting seven decades of algorithmic and hardware advances.

Exercises for Section The Governing Equations of the Atmosphere

Numerical Weather Prediction Foundations

“NWP is the crown jewel of computational science: a discipline where every factor of two in resolution, every percentage point improvement in initialisation, every better parameterisation scheme translates directly into lives saved and economic losses prevented.” - Attributed to Peter Lynch

The primitive equations of Section The Governing Equations of the Atmosphere form a continuous system of PDEs. To solve them on a computer, we must discretise them: approximate the continuous fields on a finite set of grid points or basis functions, and advance the solution forward in time using a discrete time-stepping scheme. This is the domain of Numerical Weather Prediction (NWP).

NWP is not simply “computational weather”. It is a sophisticated interplay between numerical analysis (stability, accuracy, convergence), dynamical meteorology (wave propagation, balanced flow, geostrophy), computational science (parallelisation, I/O at scale), and data science (assimilation, ensemble generation, probabilistic post-processing). Modern NWP systems represent one of the most complex software engineering endeavours in science-the ECMWF Integrated Forecasting System (IFS) comprises more than two million lines of Fortran code, accumulated over 40 years of continuous development.

Discretisation strategies

Three main families of spatial discretisation are used in NWP:

Finite differences

The simplest approach: replace spatial derivatives by algebraic differences between neighbouring grid points. For a uniform grid with spacing $\mathrm{\Delta }x$, the second-order centred difference for the first derivative is (FD)$$\begin{array}{ccc}& \frac{\partial u}{\partial x}{|}_i\approx \frac{u_{i+1}-u_{i-1}}{2\mathrm{\Delta }x}+O(\mathrm{\Delta }{x}^2).& \text{<span class="tml-eqn"></span>}\end{array}$$ Finite differences are simple to implement and easily adapted to stretched or variable grids, but achieving high accuracy requires large stencils (many neighbouring points), and formulating them on the sphere requires care near the poles (the “pole problem”).

Spectral methods

In spectral methods, the prognostic variables are expanded in a truncated series of global basis functions. For the sphere, the natural basis is spherical harmonics: (Spectral)$$\begin{array}{ccc}& u(\lambda ,\varphi ,t)=\sum _{n=0}^N\sum _{m=-n}^n{\hat{u}}_n^m(t)Y_n^m(\lambda ,\varphi ),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $Y_n^m$ are spherical harmonics of wavenumber $n$ and zonal wavenumber $m$, and ${\hat{u}}_n^m(t)$ are spectral coefficients. A triangular truncation at total wavenumber $N$ (denoted T$N$) resolves features with a minimum wavelength of approximately $\pi a/N$ where $a$ is Earth's radius. ECMWF's operational IFS uses cubic-octahedral grids: TCo1279 ($\approx 9\mathrm{km}$) for high-resolution deterministic forecasting (HRES) and TCo639 ($\approx 18\mathrm{km}$) for the ensemble.

Spectral methods achieve spectral accuracy (exponentially fast convergence for smooth fields) and automatically avoid the pole problem. They dominated global NWP from the 1970s through the 2000s. Their main limitation is that non-linear terms must be evaluated in grid space and transformed (via fast Legendre and Fourier transforms), adding cost, and that they are poorly suited to massively parallel architectures with local memory-motivating the shift to finite-volume and finite-element approaches.

Semi-Lagrangian schemes

In the semi-Lagrangian approach, the advection step is solved by tracking trajectories backward in time: to find the value of a field $\psi$ at grid point ${𝐱}_i$ at time $t+\mathrm{\Delta }t$, one solves for the departure point ${𝐱}^{\ast }$ such that a parcel currently at ${𝐱}^{\ast }$ will reach ${𝐱}_i$ by time $t+\mathrm{\Delta }t$, then interpolates $\psi ({𝐱}^{\ast },t)$ from surrounding grid points: (SL)$$\begin{array}{ccc}& \psi ({𝐱}_i,t+\mathrm{\Delta }t)\approx \psi ({𝐱}^{\ast },t)+\mathrm{\Delta }t\tilde{S}({𝐱}^{\ast },t),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\tilde{S}$ collects all non-advective source terms. The key advantage: unlike Eulerian advection, the semi-Lagrangian scheme is unconditionally stable with respect to the Courant number $c\mathrm{\Delta }t/\mathrm{\Delta }x$, allowing time steps 5–10 times larger than explicit Eulerian schemes. Combined with semi-implicit treatment of gravity waves (below), this gives the dominant time-stepping strategy for global NWP.

Grid types

The choice of horizontal grid profoundly affects the performance and scalability of an NWP model:

Gaussian grid.

The standard grid for spectral models: latitudinal grid points are placed at the zeros of Legendre polynomials (Gaussian latitudes), with a regular spacing in longitude. Efficient for Legendre transforms but has a very high point density near the poles (leading to extremely small effective grid spacings and very short stable time steps for explicit schemes).

Reduced Gaussian grid.

A practical modification: the number of longitudinal points is reduced near the poles so that the grid spacing remains roughly uniform in physical distance. ECMWF's IFS uses the octahedral reduced Gaussian grid (O$N$), where the number of points near the poles decreases following a specific octahedral tessellation. This reduces the total number of grid points by $\sim 20$% compared to the full Gaussian grid and improves computational efficiency.

Cubed sphere.

The sphere is divided into six faces of a cube, each discretised with a regular Cartesian grid. The cubed sphere avoids the pole singularity and maps well to parallel architectures (six independent patches). It is used in GFDL's FV3 dynamical core (which underlies NOAA's GFS and GEFS after 2019), the NASA GEOS model, and others.

Icosahedral (triangular and hexagonal) grids.

The sphere is tiled with triangles or hexagons derived from the geodesic refinement of an icosahedron. These grids achieve near-uniform resolution across the globe and are well-suited to adaptive mesh refinement. Germany's ICON model (Icosahedral Non-hydrostatic) uses a triangular-cell icosahedral grid; Japan's Nonhydrostatic ICosahedral Atmospheric Model (NICAM) uses a hexagonal formulation.

Time integration: semi-implicit and semi-Lagrangian schemes

Atmospheric dynamics contain waves of very different speeds. Acoustic waves travel at $\sim 300{\mathrm{m}\mathrm{s}}^{-1}$; gravity waves at $\sim 100$–$300{\mathrm{m}\mathrm{s}}^{-1}$; Rossby (weather) waves at $\sim 10$–$30{\mathrm{m}\mathrm{s}}^{-1}$. An explicit time-stepping scheme must satisfy the Courant–Friedrichs–Lewy (CFL) condition $\mathrm{\Delta }t\lesssim \mathrm{\Delta }x/c_{\max }$ with respect to the fastest wave speed $c_{\max }$, making time steps impractically small.

The semi-implicit approach treats fast linear waves (gravity and acoustic) implicitly and slow nonlinear terms explicitly. Writing the governing equations as (Semiimplicit Abstract)$$\begin{array}{ccc}& \frac{\partial 𝐗}{\partial t}=L𝐗+N(𝐗),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $L$ contains the linear (fast-wave) terms and $N$ the nonlinear (slow) terms, the semi-implicit time step reads: (Semiimplicit)$$\begin{array}{ccc}& \frac{{𝐗}^{n+1}-{𝐗}^n}{\mathrm{\Delta }t}=L\left(\frac{{𝐗}^{n+1}+{𝐗}^n}{2}\right)+N({𝐗}^n).& \text{<span class="tml-eqn"></span>}\end{array}$$ This is the two-time-level semi-implicit scheme. Solving for ${𝐗}^{n+1}$ requires inverting a linear system $(I-\frac{\mathrm{\Delta }t}{2}L){𝐗}^{n+1}=\mathrm{RHS}$, which in spectral space decouples into independent systems for each wavenumber-a crucial efficiency. Combined with semi-Lagrangian advection, this approach allows time steps of $\mathrm{\Delta }t\sim 1800\mathrm{s}$ (30 minutes) at $\mathrm{\Delta }x\sim 25\mathrm{km}$ resolution, compared to $\mathrm{\Delta }t\sim 60$,s required for explicit Eulerian schemes.

Data assimilation fundamentals

The most accurate numerical model is useless without good initial conditions. The process of constructing the best estimate of the atmospheric state at a given time, combining information from a short-range model forecast (the background) with all available observations, is called data assimilation.

The variational framework: 3D-Var and 4D-Var

Data assimilation can be formulated as a Bayesian inference problem. Let $𝐱\in {ℝ}^n$ denote the atmospheric state vector (all model variables on all grid points; $n\sim {10}^9$ for modern models) and let $𝐲\in {ℝ}^m$ denote the observation vector ($m\sim {10}^7$ observations per 6-hour window for ECMWF). The background state ${𝐱}_b$ is a short-range forecast from the previous cycle, with error covariance $𝐁\in {ℝ}^{n\times n}$. Each observation $y_i\approx h_i(𝐱)+{ϵ}_i$ is related to the state through a (possibly nonlinear) observation operator $h_i$, with observational error ${ϵ}_i$ having covariance $𝐑\in {ℝ}^{m\times m}$.

The maximum a posteriori (MAP) estimate minimises the cost function: (3dvar)$$\begin{array}{ccc}& 𝒥(𝐱)=\frac{1}{2}{(𝐱-{𝐱}_b)}^{\top }{𝐁}^{-1}(𝐱-{𝐱}_b)+\frac{1}{2}{(𝐲-H(𝐱))}^{\top }{𝐑}^{-1}(𝐲-H(𝐱)),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $H:{ℝ}^n\to {ℝ}^m$ is the observation operator mapping the model state to observation space.

• 3D-Var (three-dimensional variational assimilation): minimise $𝒥(𝐱)$ over a single time level, treating all observations as simultaneous. This is computationally cheap but ignores the temporal structure of observation coverage.

• 4D-Var (four-dimensional variational assimilation): minimise a time-extended cost function (4dvar)$$\begin{array}{ccc}& 𝒥({𝐱}_0)=\frac{1}{2}{({𝐱}_0-{𝐱}_b)}^{\top }{𝐁}^{-1}({𝐱}_0-{𝐱}_b)+\frac{1}{2}\sum _{k=0}^K{({𝐲}_k-H_k({𝐱}_k))}^{\top }{𝐑}_k^{-1}({𝐲}_k-H_k({𝐱}_k)),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝐱}_k={ℳ}_{0\to k}({𝐱}_0)$ is the model state at time $t_k$ obtained by running the forecast model $ℳ$ forward from the initial state ${𝐱}_0$. The gradient ${\nabla }_{{𝐱}_0}𝒥$ is computed using the adjoint model-the transpose of the linearised forward model's tangent linear. 4D-Var is used operationally by ECMWF and is the gold standard for NWP initialisation, but requires maintaining the adjoint code of the full forecast model.

The ensemble Kalman filter

The Ensemble Kalman Filter (EnKF), introduced by Evensen (1994) and widely used in NWP since the 2000s, approximates the background error covariance $𝐁$ by the sample covariance of a finite ensemble of model states.

An ensemble of $K$ members ${\{{𝐱}^{(k)}\}}_{k=1}^K$ (typically $K=20$–$100$) is evolved forward from the previous analysis: (ENKF Forecast)$$\begin{array}{ccc}& {𝐱}_k^{(k)}=ℳ({𝐱}_{k-1}^{(k)})+{𝜼}^{(k)},k=1,\dots ,K,& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝜼}^{(k)}$ is a model error perturbation. The sample background covariance is (ENKF B)$$\begin{array}{ccc}& 𝐁\approx \widehat{𝐁}=\frac{1}{K-1}\sum _{k=1}^K({𝐱}^{(k)}-\overline{𝐱}){({𝐱}^{(k)}-\overline{𝐱})}^{\top },& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\overline{𝐱}=K^{-1}{\sum }_k{𝐱}^{(k)}$ is the ensemble mean.

The Kalman gain matrix (Kalman GAIN)$$\begin{array}{ccc}& 𝐊=\widehat{𝐁}{H}^{\top }(H\widehat{𝐁}{H}^{\top }+𝐑{)}^{-1}\in {ℝ}^{n\times m}& \text{<span class="tml-eqn"></span>}\end{array}$$ provides the analysis update (analysis step) for each ensemble member: (ENKF Analysis)$$\begin{array}{cccc}& {𝐱}_a^{(k)}& ={𝐱}_b^{(k)}+𝐊(𝐲+{𝝐}^{(k)}-H({𝐱}_b^{(k)})),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝝐}^{(k)}\sim 𝒩(\mathrm{𝟎},𝐑)$ are perturbed observations that prevent ensemble collapse. The ensemble mean of $\{{𝐱}_a^{(k)}\}$ is the analysis state used to initialise the forecast.

Reanalysis datasets

The three most widely used reanalysis products in current AI weather research are:

ERA5 (ECMWF Reanalysis v5).

ERA5 covers 1940 to present at ${0.25}^{\circ }$ ($\approx 28\mathrm{km}$) horizontal resolution with 137 vertical pressure levels and hourly output. It assimilates more than ${10}^9$ observations per year using ECMWF's 4D-Var system. ERA5 is the primary training dataset for virtually all published AI weather models (FourCastNet, Pangu-Weather, GraphCast, GenCast, AIFS, and others). The full ERA5 archive (1940–present, all variables, all levels, hourly) occupies $\sim 5$,PB.

MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2).

NASA's MERRA-2 reanalysis covers 1980–present at ${0.5}^{\circ }\times {0.625}^{\circ }$ resolution with 72 vertical levels. It is particularly valued for its treatment of aerosols and its long-term consistency in the hydrological cycle, and is widely used in climate research and machine learning applications requiring aerosol and composition data.

JRA-55 (Japanese 55-year Reanalysis).

Produced by the Japan Meteorological Agency (JMA), JRA-55 covers 1958–present at $\approx 55\mathrm{km}$ resolution. Its extended temporal coverage (starting before the satellite era) makes it valuable for long-term trend analysis and for training models that must generalise across a wider range of climate states.

The NWP pipeline

Figure fig:weather:pipeline shows the complete pipeline of a modern operational NWP system.

Operational NWP systems

The world's major operational NWP centres maintain global models that run multiple times per day, 365 days per year, and provide the backbone of public weather forecasts. Table Table 2 summarises the key operational systems as of 2024.

System	Centre	Det. res.	Ens. members	Top level	Dynamical core
IFS (HRES)	ECMWF	$\sim 9\mathrm{km}$	1 (det.)	$0.01\mathrm{hPa}$	Spectral + SL
IFS (ENS)	ECMWF	$\sim 18\mathrm{km}$	50	$0.01\mathrm{hPa}$	Spectral + SL
GFS	NOAA/NCEP	$\sim 13\mathrm{km}$	1 (det.)	$0.002\mathrm{hPa}$	FV3 (cubed sphere)
GEFS	NOAA/NCEP	$\sim 25\mathrm{km}$	30	$0.2\mathrm{hPa}$	FV3 (cubed sphere)
UM (MOGREPS-G)	UKMO	$\sim 10\mathrm{km}$	18	$\sim 80\mathrm{km}$	ENDGame (semi-impl.)
ICON	DWD/MPI-M	$\sim 13\mathrm{km}$	40 (ICON-EPS)	$0.01\mathrm{hPa}$	Icosahedral FV
GSM	JMA	$\sim 20\mathrm{km}$	27 (GEPS)	$0.01\mathrm{hPa}$	Spectral + SL
SL = semi-Lagrangian; FV = finite volume; det. = deterministic; ens. = ensemble; res. = resolution.
ECMWF's AIFS (AI model) runs alongside IFS but is not listed as it uses a different paradigm.

ECMWF IFS.

The European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System is widely considered the world's most accurate global NWP model, consistently outperforming competing systems on objective verification metrics over the 10-day range. It uses a spectral dynamical core with TCo1279 ($\sim 9\mathrm{km}$) for deterministic and TCo639 ($\sim 18\mathrm{km}$) for ensemble forecasting, 4D-Var data assimilation, and runs on the ECMWF supercomputer in Bologna (currently $\sim 6$,petaflop/s sustained).

NOAA GFS.

The Global Forecast System is the primary operational model of the US National Centers for Environmental Prediction (NCEP), updated to the FV3 dynamical core on the cubed sphere in 2019. GFS data are freely available in near-real time, making it the primary training target for many ML weather benchmarks.

UKMO UM.

The UK Met Office Unified Model uses a non-hydrostatic, semi-implicit, semi-Lagrangian dynamical core (ENDGame) and the same model code from global scales down to storm-resolving nest domains-the “unified” philosophy. The probabilistic MOGREPS-G ensemble provides calibrated uncertainty estimates.

ICON.

The ICosahedral Nonhydrostatic model, developed jointly by the German Weather Service (DWD) and the Max Planck Institute for Meteorology, is the primary example of a production icosahedral-grid model. ICON is fully non-hydrostatic and can be run from global scales down to cloud-resolving resolution without changing the dynamical core.

JMA GSM.

The Japan Meteorological Agency's Global Spectral Model is the primary NWP model for East Asia and the western Pacific, where the JMA provides authoritative guidance on typhoon track and intensity.

Ensemble forecasting and uncertainty quantification

The atmosphere is a chaotic system: tiny differences in initial conditions amplify exponentially, so that after $\sim 14$ days even a perfect model loses all predictive skill. This predictability limit (the Lorenz limit) means that a single deterministic forecast is incomplete: it provides no information about forecast uncertainty.

Ensemble forecasting addresses this by running $K$ forecasts from slightly perturbed initial conditions (and, in modern systems, stochastic perturbations to the model physics): (Ensemble)$$\begin{array}{ccc}& {\{{𝐱}_T^{(k)}\}}_{k=1}^K={\{{ℳ}_{0\to T}({𝐱}_0^{(k)})\}}_{k=1}^K,{𝐱}_0^{(k)}={𝐱}_0+{𝜹}_k,& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝜹}_k$ are carefully constructed initial perturbations (singular vectors, bred vectors, or EnKF perturbations) and $K=50$–$100$ members. The ensemble mean approximates the expectation, and the ensemble spread provides a flow-dependent estimate of forecast uncertainty. Ensemble forecasting has transformed operational meteorology from deterministic prediction to probabilistic weather forecasting, enabling products such as probability of precipitation, percentile temperature forecasts, and severe weather risk maps.

Exercises for Section Numerical Weather Prediction Foundations

Model Output Statistics and Early Post-Processing

By the early 1970s, numerical weather prediction had achieved a remarkable technical milestone: operational atmospheric models were running daily at major meteorological centres around the world, producing gridded forecasts of temperature, wind, humidity, and precipitation on scales of hundreds of kilometres. Yet operational forecasters faced a persistent and embarrassing problem. The model output was systematically wrong in ways that were entirely predictable. Temperature forecasts at Denver, Colorado were routinely too warm because the coarse model grid smoothed the sharp elevation gradients of the Rocky Mountain Front Range. Precipitation probabilities at coastal stations were biased because the model could not resolve mesoscale orographic lifting. Wind speeds at airport sites were consistently underestimated because the aerodynamic roughness of urban terrain was absent from the land-surface parameterisation.

The operational forecaster could see these errors accumulating. They knew, from years of experience, that the model always ran cold in the San Joaquin Valley in summer, always ran dry over the Appalachian windward slopes in winter, always underforecast the morning low at coastal stations. But there was no systematic framework for correcting these biases – only institutional memory, forecaster experience, and manual adjustment.

It was into this gap that Harry Glahn and Dale Lowry stepped in 1972 with a paper that would reshape operational meteorology for the next five decades.

The Glahn–Lowry Paper and the MOS Framework

The fundamental insight of Glahn and Lowry was to treat the relationship between NWP model output and observed weather at a specific surface station as a statistical estimation problem. Rather than attempting to correct the model's physics – an enormously difficult task – one could instead learn a statistical mapping from model-predicted quantities to locally observed quantities, using historical pairs of model forecasts and verifying observations.

This idea rested on three prior observations. First, model errors are not random; they are systematic, location-specific, and seasonally varying. Second, the model's large-scale fields contain genuine skill – correctly capturing synoptic-scale patterns such as frontal passages, cyclone tracks, and upper-level ridges – even when local surface variables are poorly forecast. Third, multiple model output fields (upper-level temperatures, geopotential heights, absolute vorticity, precipitable water) contain complementary information about local weather conditions that no single predictor captures alone.

These three observations collectively motivate the use of multiple linear regression applied to a carefully selected subset of model output fields as predictors of a surface meteorological variable of interest. The method Glahn and Lowry developed for selecting this subset – screening regression – became the methodological core of operational MOS.

The definition captures several important degrees of freedom that distinguish MOS from naive bias correction. First, the predictor vector $𝐱$ may contain not only the model's direct forecast of the predictand (e.g., the model's 2-metre temperature forecast) but also indirect predictors: upper-level temperatures, thickness values, vorticity, relative humidity at multiple pressure levels, and even squared or product terms formed from primary predictors. Second, the mapping is station-specific: a separate regression equation is fitted for each station and each forecast lead time. Third, the mapping implicitly performs both bias correction (removing systematic mean errors) and downscaling (bridging the spatial resolution gap between the coarse model grid and the point observation).

Sequential Screening Regression

The central algorithmic challenge in MOS equation development is predictor selection. A modern NWP model may produce forecasts of hundreds of variables at dozens of pressure levels and multiple lead times, yielding a predictor pool of potentially thousands of candidates for a single predictand. Including all predictors in the regression would produce severe overfitting on the developmental sample and numerically unstable coefficient estimates.

The solution adopted by Glahn and Lowry – and still used in operational MOS development today – is sequential forward screening, a greedy procedure that builds the regression equation one predictor at a time.

The screening algorithm is conservative by design. The entry and retention criteria prevent the inclusion of predictors whose apparent skill is an artefact of the finite developmental sample rather than genuine meteorological signal. The check-retention step (step 4) is particularly important: as new predictors enter the equation, previously entered predictors may become redundant due to multicollinearity, and their removal improves the equation's out-of-sample performance.

The result of the screening procedure for a specific station, predictand, and lead time is a small set of regression coefficients $\widehat{𝜷}=({\hat{\beta }}_0,{\hat{\beta }}_1,\dots ,{\hat{\beta }}_k)$ that defines the MOS equation. These coefficients are stored and applied in real time: when the NWP model produces a new forecast, the predictor values are extracted, fed into the equation, and the local station forecast is produced within seconds.

Predictands: Temperature, Precipitation, and Wind

Operational MOS systems develop separate equation sets for each predictand. The three predictand families of greatest economic and social importance are surface temperature, precipitation occurrence and amount, and wind.

Temperature predictands.

The most important temperature predictands are maximum daytime temperature ($T_{\max }$) and minimum overnight temperature ($T_{\min }$) at the 2-metre standard observation height. These quantities drive energy demand, agricultural frost risk, human health impacts, and transportation planning. Classical MOS predictors for $T_{\max }$ include the model's 850-hPa temperature (which is largely free of surface boundary-layer errors), the model's 2-metre temperature (which carries the direct low-level information despite its biases), 1000–500-hPa thickness (a proxy for the mean tropospheric temperature), relative humidity, and the sine and cosine of the day of year (to capture seasonal residual biases after the climatological mean is removed).

Precipitation predictands.

The primary precipitation predictand in classical MOS is the probability of precipitation (PoP), defined as the probability that measurable precipitation ($\ge 0.01$ inches) occurs within a specified time window (typically 12 or 24 hours) at or near the forecast station. Because PoP is a probability rather than a continuous quantity, the linear regression output must be transformed; the classical approach uses a linear regression for a conditional expectation and then truncates predictions to $[0,1]$, though in practice the regression coefficients are chosen so that the fitted values naturally fall within a physically plausible range. More sophisticated extensions include the use of logistic regression or probit regression applied to a dichotomous (occurrence/non-occurrence) observation.

Wind predictands.

Wind is the most challenging classical MOS predictand because it is inherently two-dimensional and exhibits strong non-Gaussianity due to the directional wrapping of the wind rose. Operational MOS systems typically forecast wind speed and wind direction separately, using the model's 10-metre wind components $u_{10}$ and $v_{10}$ as primary predictors but also including terrain-adjusted wind exposure indices, mixing-layer depth, and roughness-length corrections specific to each station's exposure category.

Bias Correction and Downscaling

The two most important functions performed by MOS regression equations – beyond pure statistical prediction – are systematic bias correction and spatial downscaling. These are conceptually distinct but are accomplished simultaneously by the same regression coefficients.

Systematic bias correction.

Let ${\mu }_y$ denote the climatological mean of the predictand at the station and ${\mu }_{{\tilde{x}}_1}$ denote the long-run mean of the primary model predictor. If the model has a systematic warm bias $b={\mu }_{{\tilde{x}}_1}-{\mu }_y>0$, the intercept term ${\hat{\beta }}_0$ will adjust downward to compensate. More precisely, taking expectations of the MOS equation over the developmental sample: (Unbiasedness)$$\begin{array}{ccc}& 𝔼[\hat{y}]={\hat{\beta }}_0+\sum _{j=1}^k{\hat{\beta }}_j𝔼[x_j]=\bar{y},& \text{<span class="tml-eqn"></span>}\end{array}$$ so the MOS forecast is unbiased over the developmental period by construction. This is the climatological unbiasedness property of ordinary least-squares regression, and it guarantees that MOS forecasts will not exhibit the systematic warm or cold biases present in the raw NWP output.

Spatial downscaling.

The NWP model produces forecasts on a grid with spacing of order 10–100 km in operational systems. A surface station may sit at an elevation, land-surface type, or coastal proximity that differs substantially from the grid-cell average. The MOS regression coefficients encode this discrepancy implicitly: if the station sits in a valley that is systematically 3${}^{\circ }$C cooler than the surrounding plateau that dominates the model grid cell, that constant offset is absorbed by the regression intercept ${\hat{\beta }}_0$; a slope different from unity instead encodes an amplified or damped sensitivity of the station to the large-scale thermal signal. No explicit terrain correction is applied; the statistics of the developmental sample speak for themselves.

Example: MOS Equations for Maximum Temperature

To make the abstract framework concrete, consider a hypothetical MOS equation for maximum temperature ($T_{\max }$ in ${}^{\circ }$F) at a mid-latitude inland station, valid 24–36 hours after model initialisation. After screening from a pool of 312 candidate predictors, suppose the screening algorithm selects the following five predictors:

Variable	Description	Pressure level	$r$
$x_1$	2-metre temperature	Surface	$0.89$
$x_2$	Temperature	850 hPa	$0.82$
$x_3$	1000–500 hPa thickness	–	$0.74$
$x_4$	Relative humidity	700 hPa	$-0.31$
$x_5$	$\cos (2\pi d/365)$	–	$0.27$

The fitted MOS equation takes the form: (TMAX Example)$$\begin{array}{ccc}& {\hat{T}}_{\max }=\underset{{\hat{\beta }}_0}{\underbrace{4.7}}+\underset{{\hat{\beta }}_1}{\underbrace{0.61}}{x}_1+\underset{{\hat{\beta }}_2}{\underbrace{0.28}}{x}_2+\underset{{\hat{\beta }}_3}{\underbrace{0.019}}{x}_3-\underset{{\hat{\beta }}_4}{\underbrace{0.15}}{x}_4+\underset{{\hat{\beta }}_5}{\underbrace{1.8}}{x}_5,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $x_3$ is expressed in decametres and $x_4$ in per cent.

Several features of equation are worth noting. The dominant predictor is $x_1$, the model's own 2-metre temperature forecast, with a coefficient less than unity: this is regression attenuation, since the noisy predictor makes the least-squares fit shrink its slope toward the mean, so the regression pulls the model's temperature slightly toward the climatological mean. The 850-hPa temperature $x_2$ provides an independent estimate of the synoptic thermal environment, adding skill beyond $x_1$ alone because the 850-hPa level is largely above the boundary layer and hence less affected by surface parameterisation errors. The negative coefficient on $x_4$ (mid-tropospheric relative humidity) captures the physical relationship between cloud cover (which suppresses maximum temperatures) and moisture. The seasonal harmonic $x_5$ corrects for residual seasonal biases that persist after the large-scale predictors are included.

The MOS Pipeline: From NWP Output to Local Forecast

fig:weather:mos:pipeline illustrates the complete operational MOS workflow, from NWP model output through predictor extraction, regression application, and local forecast dissemination.

SAPEMOS and Regional Variants

The original Glahn–Lowry MOS framework was developed for the United States and relied on the dense observational network and long archive of NWP model runs that the National Weather Service had accumulated. Applying MOS in regions with sparse observational networks, short model archives, or less consistent station records required significant methodological adaptation.

SAPEMOS (South African Probabilistic Ensemble MOS) is a regional adaptation developed to provide probabilistic surface weather forecasts for southern Africa, where the observational network is far less dense than in North America or Europe. SAPEMOS addresses the data-sparsity problem through several innovations. First, it uses ensemble NWP output (rather than a single deterministic forecast) as the predictor set, so that forecast spread is available as an additional predictor alongside the ensemble mean; this allows the regression to encode the relationship between model uncertainty and forecast skill. Second, it employs spatial pooling across nearby stations to augment the effective developmental sample size, using station-specific dummy variables and elevation-adjusted corrections to retain station-specific information while borrowing statistical strength from the regional pool. Third, it outputs full probability distributions over the predictand rather than point forecasts, matching the needs of decision-makers in sectors such as agriculture, water management, and disaster preparedness.

Classical Statistical Ensemble Post-Processing

The shift from single deterministic NWP runs to ensemble NWP systems – multiple model runs with perturbed initial conditions and model parameters – opened a new chapter in statistical post-processing. An ensemble is, in principle, a probabilistic forecast: the spread of ensemble members should encode forecast uncertainty. In practice, however, raw ensembles are both underdispersive (their spread is too narrow, underestimating uncertainty in smooth flows) and biased (systematic errors in initial perturbations and model parameterisations shift the ensemble mean away from the true state). Classical statistical ensemble post-processing methods were developed to correct these deficiencies.

Bayesian Model Averaging (BMA).

Raftery and colleagues introduced Bayesian Model Averaging (BMA) for weather forecasts in a landmark 2005 paper. The BMA approach treats each ensemble member as generating a conditional predictive distribution for the observation, and the overall predictive distribution is a weighted mixture: (BMA)$$\begin{array}{ccc}& p(y|{𝐱}_1,\dots ,{𝐱}_K)=\sum _{k=1}^K{w}_k{p}_k(y|{𝐱}_k),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $w_k\ge 0$, ${\sum }_k{w}_k=1$ are mixture weights fitted by maximum likelihood (via the EM algorithm) over a training window, and $p_k(y|{𝐱}_k)$ is a Gaussian predictive density centred at a bias-corrected version of ensemble member $k$.

Ensemble Model Output Statistics (EMOS).

Gneiting and colleagues proposed EMOS (also called Non-homogeneous Gaussian Regression, NGR) as a more parsimonious alternative to BMA. Rather than fitting a mixture distribution, EMOS fits a single Gaussian whose mean and variance are affine functions of the ensemble mean and variance: (EMOS)$$\begin{array}{ccc}& y|𝐱\sim \mathrm{Normal}(a+b\bar{x},c+d{s}^2),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\bar{x}$ is the ensemble mean, $s^2$ is the ensemble variance, and $(a,b,c,d)$ are fitted parameters. The parameters are estimated by minimising the continuous ranked probability score (CRPS) over a rolling training window.

Non-homogeneous Gaussian Regression (NGR).

NGR is essentially a synonym for EMOS, emphasising the statistical interpretation: the predictive variance is allowed to be a function of the ensemble spread (hence “non-homogeneous”), unlike classical homoscedastic regression where the residual variance is constant. This flexibility is critical for ensemble post-processing because the model's predictability varies substantially by weather regime: zonal, predictable flows yield small ensemble spread and tight predictive distributions, while amplified blocked patterns yield large spread and broad distributions.

Limitations of Linear MOS

For all its practical success, classical linear MOS carries fundamental limitations that became increasingly apparent as operational forecasters pushed into higher-resolution, higher-stakes prediction problems.

Linearity.

The central assumption of MOS is that the relationship between model output predictors and local observations is well approximated by a linear combination. This assumption is reasonable for large-scale temperature and thickness predictors, where the relationship is genuinely quasi-linear over synoptic scales. But it fails for a broad class of meteorologically important phenomena. Convective initiation – the transition from a stable, cloudless atmosphere to active convective precipitation – is intrinsically nonlinear: small changes in lower-tropospheric moisture and stability can flip the atmosphere from suppressed to convectively active within hours, a bifurcation that no linear model can represent.

Interactions.

A linear regression cannot capture the interaction between predictors. Severe wind gusts at a coastal station depend not on wind speed or atmospheric stability alone, but on their joint configuration: a marine layer that is both deep and strongly stable traps coastal jets more effectively than either condition alone. The product term $x_{\text{depth}}\cdot x_{\text{stability}}$ could in principle be added as a new predictor, but the curse of dimensionality quickly limits the number of interaction terms that can be included without overfitting on the developmental sample.

Spatial structure.

Classical MOS equations are fitted independently for each station, ignoring the spatial correlation structure of the atmosphere. Two stations 50 km apart share large fractions of their forecast error, and a method that models them jointly could exploit this correlation to improve both forecasts. Station-independent regression discards this information.

Regime dependence.

The statistical relationship between model output and observations can change dramatically between weather regimes. A regression equation fitted on climatological averages of all regime types will perform poorly in any specific regime, particularly extreme events where the atmosphere is furthest from its climatological state and where forecast accuracy matters most.

Historical Note: Harry Glahn at NOAA-MDL

Exercises

Early Machine Learning in Meteorology

The closing years of the twentieth century brought a surge of interest in applying artificial neural networks to meteorological problems. The appeal was obvious: neural networks are universal approximators, capable in principle of learning arbitrary nonlinear mappings from input to output, and the limitations of linear MOS that we documented in the previous section were precisely the limitations that nonlinear function approximation promised to cure. Convective initiation, regime-dependent forecast skill, complex multiway predictor interactions – all of these phenomena involve nonlinearities that a sufficiently expressive neural network could in principle capture.

The reality, as we shall see, was more complicated. The neural networks of the 1990s and early 2000s were shallow – almost universally single-hidden-layer multilayer perceptrons trained by backpropagation. They were trained on modest datasets that were large by the standards of their era but tiny by those of modern deep learning. They ran on hardware orders of magnitude slower than contemporary GPUs. And they were applied to a forecasting problem whose spatiotemporal complexity would not be revealed until convolutional and recurrent architectures, trained on continental- scale radar archives, arrived two decades later.

Nevertheless, the early ML literature in meteorology is far from insignificant. It established the conceptual vocabulary of learned weather prediction, documented the first systematic experiments with nonlinear post-processing, and laid intellectual groundwork on which deep learning methods would later build. Understanding why these early methods were ultimately limited – and what was needed to overcome those limits – is essential context for appreciating the deep learning revolution in weather forecasting.

The First Wave: 1990s ANN Applications

The first serious applications of artificial neural networks to weather problems appeared in the early 1990s, driven by the renewed interest in connectionism that followed the rediscovery of backpropagation by Rumelhart, Hinton, and Williams in 1986. Meteorologists and environmental scientists were attracted by the promise of a nonlinear regression substitute that required no explicit specification of the functional form.

Chow and Leung (1996): Weather-Compensated Load Forecasting

An influential early application bridging meteorology and electrical engineering came from Chow and Leung, who in 1996 applied neural networks to short-term load forecasting – predicting electricity demand for grid management – with explicit weather-variable compensation.

The problem is meteorologically interesting because electricity demand depends nonlinearly on temperature: demand rises sharply for cooling at high temperatures and for heating at low temperatures, with a minimum in the mild temperature range. This U-shaped or V-shaped temperature–load relationship is precisely the kind of nonlinearity that linear regression handles poorly and that neural networks can in principle represent.

Chow and Leung used a two-stage network architecture. The first stage was a weather compensation module: a small MLP trained to predict the load residual attributable to temperature and humidity deviations from climatological norms, given the forecast temperature and humidity as inputs. The second stage was a base-load predictor that modelled the calendar-driven component of demand (day of week, hour of day, holiday flags). The outputs were combined linearly to produce the final load forecast.

The key contribution was demonstrating that the weather compensation could be learned automatically from historical data rather than parameterised analytically, and that the resulting nonlinear mapping outperformed linear regression-based weather adjustment on a dataset of Hong Kong hourly load observations. The authors reported mean absolute percentage errors that were reduced by 15–20% relative to the linear baseline in summer months, when the nonlinear cooling load dominated.

Hall, Brooks, and Doswell (1999): ANN Review for Meteorology

By the late 1990s, enough neural network studies had appeared in the meteorological literature to warrant a systematic review. Hall, Brooks, and Doswell published an influential survey in 1999 in the Weather and Forecasting journal that assessed the state of ANN applications across multiple meteorological problem domains: severe weather prediction, precipitation estimation, tropical cyclone intensity, and synoptic pattern recognition.

Several consistent findings emerged from the review. First, neural networks consistently outperformed linear regression on tasks with demonstrably nonlinear input–output relationships, with the largest improvements in heavy precipitation prediction and severe convective event discrimination. Second, the performance of neural networks was highly sensitive to the training sample size and composition: networks trained on climatologically balanced samples generalised poorly to rare extreme events because the extreme cases were underrepresented in the training data. Third, the “black box” nature of the single-hidden-layer MLP was a significant practical barrier to operational deployment: operational forecasters were reluctant to use tools whose internal reasoning they could not interrogate or verify.

The last finding anticipates a tension that remains unresolved in modern deep learning for weather: the trade-off between model expressiveness and interpretability. A linear MOS equation can be examined term by term; a forecaster can understand why the model predicted a warm day by looking at which predictors contributed most to the forecast. A neural network provides no such transparency.

Fuzzy and Bayesian Approaches

Chen and Hwang (2000): Fuzzy Time Series for Temperature

Alongside neural network approaches, the late 1990s and early 2000s saw a parallel tradition of fuzzy time series methods applied to meteorological forecasting. Chen and Hwang extended Song and Chissom's fuzzy time series framework to the problem of temperature prediction, introducing a method that represented historical observations as fuzzy sets and used fuzzy relational equations to model the temporal dynamics.

The approach partitioned the historical temperature range into a set of fuzzy linguistic variables (e.g., “very cold”, “cold”, “mild”, “warm”, “very hot”), defined membership functions for each linguistic variable, and then estimated a fuzzy transition matrix encoding the conditional probability that a temperature classified as “mild” today transitions to “warm” tomorrow. Forecasting consisted of applying the fuzzy matrix to the fuzzified current observation and defuzzifying the result to obtain a crisp temperature prediction.

The method's attraction was its interpretability: the fuzzy rules had a natural linguistic interpretation that could in principle be examined and validated by domain experts. Its limitation was its inability to use information from more than one or two previous time steps without an exponential explosion in the size of the fuzzy relational matrix.

Cofino et al. (2002): Bayesian Networks for Probabilistic Weather

A more principled probabilistic approach appeared in 2002 with the work of Cofino and colleagues, who applied Bayesian networks – directed acyclic graphs encoding conditional independence relationships between random variables – to the prediction of precipitation occurrence and amount at Spanish meteorological stations.

The Bayesian network framework represented the joint distribution of NWP predictors and surface observations as a factored product: (Bayesnet)$$\begin{array}{ccc}& P(𝐲,𝐱)=\prod _{i=1}^{n}P(y_i|\mathrm{Pa}(y_i)),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\mathrm{Pa}(y_i)$ denotes the parent nodes of $y_i$ in the directed graph and the factorisation structure encodes the meteorologist's prior beliefs about which NWP fields causally influence which surface observations.

The key advantage over MOS was the ability to represent complex joint distributions over multiple simultaneous predictands and to propagate probabilistic uncertainty through the network. If the 500-hPa geopotential height was uncertain (as in an ensemble forecast), that uncertainty could be propagated through the network structure to yield a predictive distribution for precipitation that reflected both the model uncertainty and the statistical relationship between the large-scale predictor and the local observation.

Schizas, Baboo, and Backpropagation Networks

The application of backpropagation networks to the specific problem of minimum temperature forecasting represents one of the most thoroughly documented threads in early ML meteorology. Schizas, in a series of studies conducted at the University of Cyprus, and Baboo and colleagues in subsequent independent work, systematically compared MLP architectures of varying depth and width on the problem of predicting next-day minimum temperature from NWP predictors and surface observations.

The findings were consistent across different geographic settings and data periods. A single-hidden-layer MLP with between 10 and 20 hidden units, trained by batch backpropagation with momentum, outperformed linear MOS on temperature predictands when the training period included multiple weather regimes. The improvement was largest for minimum temperature – the predictand most sensitive to boundary-layer processes that the NWP model parameterised crudely – and smallest for maximum temperature, where the large-scale synoptic predictors (850-hPa temperature, 1000–500-hPa thickness) captured most of the variance.

The networks were typically trained with 5–10 years of daily data, yielding training sets of order $N\approx \mathrm{2,000}$–$\mathrm{4,000}$ samples. Generalisation to independent test periods was assessed by root-mean-squared error and mean absolute error. While the improvements over MOS were statistically significant at many stations, they were modest in absolute terms – typically 0.1–0.3 ${}^{\circ }$C reductions in RMSE.

Hsieh (2009): Nonlinear Methods for Environmental Science

The most comprehensive attempt to synthesise the scattered early ML literature in atmospheric and oceanic science came with the publication of William Hsieh's 2009 textbook Machine Learning Methods in the Environmental Sciences. Published by Cambridge University Press, the book provided the first systematic treatment of neural networks, kernel methods, self-organising maps, and related techniques from the perspective of a practicing atmospheric scientist.

Hsieh's central contribution was to give the disparate ML applications in meteorology a unified statistical foundation. He showed that MOS could be seen as a special case of the broader framework of nonlinear regression, that Bayesian methods provided a principled approach to uncertainty quantification, and that the overfitting problems observed empirically in early neural network studies could be understood through the lens of the bias–variance trade-off.

The Limits of Shallow Networks

To understand why the early ML approaches in meteorology remained marginal – improving on MOS but never replacing it operationally – we must examine the specific limitations of shallow networks when applied to atmospheric data.

Vanishing gradients.

The standard activation function in 1990s neural networks was the logistic sigmoid $\sigma (z)={(1+e^{-z})}^{-1}$ or the hyperbolic tangent. Both functions saturate: for $|z|\gg 1$, the gradient ${\sigma }^{\prime }(z)\approx 0$. In networks with more than one or two hidden layers, the backpropagation gradient signal diminishes exponentially as it propagates from the output layer toward the input layer, making the early layers effectively untrainable. The practical consequence was that deep networks were inferior to shallow ones in the 1990s, not because depth was theoretically less powerful, but because the training algorithm failed to exploit it.

The vanishing gradient problem meant that early ML practitioners confined themselves to single-hidden-layer architectures. The universal approximation theorem guarantees that a single-hidden-layer network with sufficiently many units can approximate any continuous function, but the number of units required may be exponentially large in the input dimension. A single hidden layer of 20 units cannot capture the hierarchical spatial structure of a weather system – the synoptic scale pattern of ridges and troughs that organises the mesoscale convective structure that organises the microscale turbulent boundary layer – without an enormous number of parameters that the small training datasets of the era could not support.

No exploitation of spatial structure.

The single-hidden-layer MLP treats all inputs as exchangeable real numbers, with no structure imposed by the spatial proximity of adjacent grid points. A meteorological field – a 2-D array of temperature or geopotential values on a regular grid – contains rich spatial structure: nearby points are correlated, fronts appear as coherent oriented discontinuities, vorticity filaments have characteristic length scales. A fully connected MLP sees this field as a flat vector and must learn spatial relationships entirely from data, requiring far larger training sets than were available in the 1990s.

No capture of temporal dynamics.

Similarly, the early MLP architectures did not exploit the temporal structure of meteorological sequences. A weather state at time $t$ is not statistically independent of the states at times $t-1,t-2,\dots$; the atmosphere is a dynamical system, and its future state depends on its recent history. MLP architectures that use a fixed window of lagged inputs cannot generalise to longer-range dependencies and cannot represent the recurrent dynamics of developing weather systems.

Data scarcity.

Digital archives of NWP model output paired with surface observations were limited in the 1990s. Reanalysis products (ERA-Interim did not begin until 1979, and was not available until the late 2000s; NCEP/NCAR reanalysis extended to 1948 but at low resolution) provided modest training datasets by modern standards. The consequence was systematic undertraining of models with more than a few hundred free parameters.

Hardware constraints.

Training a neural network on 4,000 samples with 50 hidden units by backpropagation required hours on a workstation in 1996. Experimenting with multiple architectures, regularisation strengths, and predictor sets – the hyperparameter search that is routine today – was prohibitively expensive. Practitioners typically trained one or two architectures and reported results on the architecture that appeared to work best, a practice that inflated apparent gains and made inter-study comparisons unreliable.

The ConvLSTM Breakthrough: Shi et al. (2015)

The transition from shallow post-processing networks to deep spatiotemporal architectures for weather forecasting has a remarkably precise historical date: 2015, the year of the ConvLSTM paper by Shi and colleagues at the Hong Kong University of Science and Technology.

The problem they addressed was radar nowcasting: predicting the spatial structure of radar-observed precipitation at lead times of 0–2 hours. This is a qualitatively different problem from the MOS and ANN applications surveyed above. It is not a post-processing problem (there is no NWP model whose output is being corrected); it is a direct spatiotemporal prediction problem, mapping a sequence of past radar images to a sequence of future radar images.

The key insight is that spatial convolution replaces matrix multiplication at every gate, so the hidden state ${𝖧}_t$ and cell state ${𝖢}_t$ remain spatial tensors rather than collapsing to vectors. This preserves the spatial structure of the input throughout the temporal unrolling, avoiding the information destruction that occurs when a fully connected LSTM flattens a spatial field to a vector.

Shi et al. trained a stacked ConvLSTM on the HKO-7 dataset, a 7-year archive of Hong Kong Observatory radar images at 6-minute intervals. The architecture used three ConvLSTM layers, each with $64$ filters of size $3\times 3$, followed by a convolutional output layer predicting the next 20 radar frames (2 hours at 6-minute intervals).

The results were decisive. The ConvLSTM architecture outperformed ROVER (the optical-flow nowcasting baseline), LSTM (without convolution), and other deep learning baselines on all three evaluation metrics: CSI (Critical Success Index), FAR (False Alarm Rate), and POD (Probability of Detection). The improvement over the LSTM without convolution was particularly informative: it demonstrated that spatial structure was the critical missing ingredient, not just temporal depth.

Shi et al. (2017): TrajGRU and the HKO-7 Benchmark

In a 2017 follow-up, Shi and colleagues identified a limitation of the ConvLSTM: its convolutional connections are location-invariant (the same kernel is applied everywhere), which means the network implicitly assumes the same flow structure at every spatial location. Weather patterns, however, are inherently location-variant: a cyclone generates counterclockwise flow in the northern hemisphere, so the contribution of the state to the north differs systematically from the contribution of the state to the south.

The Trajectory GRU (TrajGRU) architecture they introduced replaced the fixed recurrent convolution with a learnable flow field: at each time step, the network predicted an optical-flow-like warp field ${𝐔}_t\in {\mathrm{ℝ}}^{H\times W\times 2}$ that specified how to sample the previous hidden state when computing the current state. This sampling was implemented via bilinear interpolation, making the operation differentiable and compatible with standard backpropagation.

(WARP)$$\begin{array}{ccc}& {\widetilde{𝖧}}_{t-1}=\mathrm{Warp}({𝖧}_{t-1},{𝐔}_t),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝐔}_t$ is predicted from ${𝖷}_t$ and ${𝖧}_{t-1}$ by a small convolutional sub-network, and $\mathrm{Warp}$ denotes bilinear sampling at the positions displaced by ${𝐔}_t$. The TrajGRU cell then computes standard GRU transitions using ${\widetilde{𝖧}}_{t-1}$ in place of ${𝖧}_{t-1}$.

The same paper introduced a comprehensive benchmark evaluation framework for the HKO-7 dataset, defining standard train/validation/ test splits, evaluation metrics, and baselines (optical flow, ROVER, ConvLSTM). This benchmark infrastructure was an important contribution independent of the TrajGRU architecture itself: it provided a reproducible comparison point that subsequent methods could use, an example that the weather forecasting deep learning community would follow in establishing later benchmarks such as WeatherBench and ERA5-based evaluation suites.

MetNet and MetNet-2: Toward Continental-Scale Nowcasting

The conceptual arc from MOS through ConvLSTM reached a significant milestone with Google's MetNet (2020) and MetNet-2 (2021), presented in the paper by Sonderby and colleagues. These models extended the spatiotemporal deep learning paradigm to the full problem of precipitation probability forecasting over the continental United States at 1–8 km resolution, with lead times up to 8 hours (MetNet) and 12 hours (MetNet-2).

MetNet departed from the ConvLSTM paradigm in several important respects. First, it used a multi-scale convolutional encoder (based on a deep ResNet-style architecture with axial self-attention) rather than a recurrent architecture, extracting spatial context at multiple resolutions from a large spatial patch of radar and satellite observations. Second, it conditioned on the forecast lead time explicitly, using a learnable lead-time embedding that allowed a single model to produce forecasts at all lead times from 1 minute to 8 hours. Third, it was trained on a very large archive of MRMS (Multi-Radar Multi-Sensor) observations over the continental US, with years of training data covering diverse weather regimes.

The MetNet-2 paper reported that at lead times up to 6 hours, the model outperformed HRRR (the High-Resolution Rapid Refresh NWP model) in terms of CSI and other standard metrics. This was a landmark result: a learned model, trained entirely from observations without any explicit atmospheric physics, had matched or exceeded a state-of-the-art operational NWP model at mesoscale precipitation nowcasting.

The Evolution of Architectures: From MOS to Transformers

fig:weather:earlyml:evolution charts the conceptual evolution from classical MOS through shallow ANNs, ConvLSTM, and transformer-based foundation models.

Comparison of Early ML Methods

Table 4 summarises the key characteristics and limitations of the early ML approaches surveyed in this section.

Method	Spatial	Temporal	Probabilistic	Primary limitation
Linear MOS	No	No	Partial	Linearity; no interactions
Logistic/NGR	No	No	Yes	Parametric family assumption
BMA	No	No	Yes	Gaussian members; computational cost
Fuzzy time series	No	Partial	No	Exponential state space
Bayesian network	No	No	Yes	Manual graph elicitation
Shallow ANN (MLP)	No	No	No	Vanishing gradient; no spatial
ConvLSTM	Yes	Yes	No	Short lead time; no physics
TrajGRU	Yes	Yes	No	Short lead time; no physics
MetNet/MetNet-2	Yes	Yes	Yes	$\le$ 12-hour horizon

The Transition from Shallow to Deep

The history surveyed in this section invites a natural question: why did the transition from shallow post-processing networks to deep spatiotemporal architectures take twenty years, and what eventually made it possible?

The forty-three years between Glahn and Lowry's 1972 paper and the 2015 ConvLSTM paper were not wasted. The early ML community established the problem formulation: that weather forecasting could be cast as a statistical learning problem with NWP output and observations as inputs and outputs. They documented the failure modes of linear and shallow nonlinear models, motivating the specific architectural innovations that would arrive later. And they built the evaluation infrastructure – standard datasets, metrics, and comparison baselines – on which the deep learning era would build.

What they could not provide was the essential technical substrate: the hardware, the data, and the architecture. These arrived in a compressed window between 2012 and 2017, and the field of weather forecasting has not been the same since.

Exercises

FourCastNet and the Spectral Revolution

The history of numerical weather prediction is in part a history of choosing the right basis. Richardson's finite-difference scheme (1922) tiled the atmosphere onto a rectilinear grid and advanced each cell forward in time by finite differences. It was correct in principle but required ten thousand human computers to advance a six-hour forecast for the northern hemisphere. Charney's spectral model (1949) transformed the governing equations into Fourier space, expressing every field as a superposition of spherical harmonics. The spectral formulation did not change the physics; it changed the algebra, converting the differential operators of the primitive equations into diagonal multiplications in frequency space. The ECMWF's Integrated Forecasting System (IFS) still relies on a spectral core for its dynamical equations.

The first generation of deep learning weather models-U-Net architectures trained to minimise mean-squared error-largely ignored this spectral heritage. They processed data on the physical grid, applying local convolution kernels that could only communicate information one receptive-field width at a time. The atmosphere is global: a blocking anticyclone over the Urals changes the phase of the North Atlantic Oscillation, which modulates precipitation over Britain, which alters the sensible-heat flux that feeds back into large-scale dynamics. Capturing these teleconnections requires global receptive fields, and achieving them with local convolutions requires either very deep networks or explicitly non-local attention.

Transformers offer global attention: every token communicates with every other token at each layer. But the cost of dense self-attention scales as $𝒪(N^2)$ in the number of tokens $N$. At the resolution used by ERA5-${0.25}^{\circ }\times {0.25}^{\circ }$, corresponding to a $720\times 1440$ longitude-latitude grid-a weather state comprises $N=\mathrm{1,036,800}$ spatial locations. Standard self-attention on this grid is computationally impossible. FourCastNet, introduced by Pathak et al. from NVIDIA in 2022, solved this problem by moving the token-mixing operation from physical space into Fourier (spectral) space, where global communication can be achieved in $𝒪(N\log N)$ operations via the Fast Fourier Transform.

Complexity and the Case for Spectral Mixing

To understand why spectral processing is attractive for global atmospheric fields, consider the three main approaches to spatial mixing and their asymptotic costs on a two-dimensional grid of $H\times W$ points, with $C$ channels.

Local convolution.

A $k\times k$ convolution over a single channel pair requires $k^{2}HW$ multiply-add operations. The receptive field after $L$ layers grows as ${(L(k-1)+1)}^2$; to achieve a global receptive field on a $720\times 1440$ grid requires hundreds of layers. Communication is local; teleconnections emerge only indirectly.

Dense self-attention.

Treating each grid cell as a token, the attention matrix $𝐀\in {\mathrm{ℝ}}^{N\times N}$ (where $N=HW$) has $N^2$ entries. Computing the matrix product $\mathrm{Attention}(𝐐,𝐊,𝐕)=\mathrm{softmax}\left(\frac{𝐐{𝐊}^{\top }}{\sqrt{d}}\right)𝐕$ costs $𝒪(N^{2}d)$ operations and $𝒪(N^2)$ memory. For $N\approx {10}^6$ this is prohibitive.

Spectral mixing (AFNO).

The Discrete Fourier Transform (DFT) of a signal of length $N$ can be computed in $𝒪(N\log N)$ operations via the FFT. In Fourier space, convolution becomes pointwise multiplication; hence a global convolution-one whose kernel spans the full domain-costs only $𝒪(N\log N)$. The Adaptive Fourier Neural Operator augments this with a learned, frequency-wise mixing step, achieving global communication at $𝒪(N\log N)$ cost, five orders of magnitude cheaper than dense attention for ERA5 resolution.

Fourier Neural Operators: Mathematical Foundations

The Fourier Neural Operator (FNO), introduced by Li et al. (2021), provides the mathematical basis for spectral-domain learning. Let $\mathrm{\Omega }\subset {\mathrm{ℝ}}^d$ be a bounded domain (in weather applications, $\mathrm{\Omega }$ is the sphere ${𝕊}^2$ discretised as a regular longitude-latitude grid). An FNO layer maps a function $v:\mathrm{\Omega }\to {\mathrm{ℝ}}^{d_v}$ to a function $u:\mathrm{\Omega }\to {\mathrm{ℝ}}^{d_u}$ via

(FNO Layer)$$\begin{array}{ccc}& u(x)=\sigma (Wv(x)+(𝒦(v))(x)),& \text{<span class="tml-eqn"></span>}\end{array}$$

where $W\in {\mathrm{ℝ}}^{d_u\times d_v}$ is a pointwise linear map (the residual branch), $\sigma$ is a nonlinearity applied elementwise, and $𝒦$ is the integral kernel operator

(Integral Kernel)$$\begin{array}{ccc}& (𝒦(v))(x)={\int }_{\mathrm{\Omega }}\kappa (x-y)v(y)\mathrm{d}y.& \text{<span class="tml-eqn"></span>}\end{array}$$

By the convolution theorem, equation becomes, in Fourier space,

(CONV Theorem)$$\begin{array}{ccc}& \widehat{𝒦(v)}(\xi )=\hat{\kappa }(\xi )\cdot \hat{v}(\xi ),& \text{<span class="tml-eqn"></span>}\end{array}$$

where $\hat{\cdot }$ denotes the Fourier transform and $\xi$ is a frequency index. The FNO parametrises $\hat{\kappa }(\xi )$ as a learnable complex-valued tensor $R_{\xi }\in {ℂ}^{d_u\times d_v}$ for each retained frequency $\xi$. Only the lowest $k_{\max }$ modes are learned; high-frequency components are zeroed (a soft truncation that acts as a spectral filter). The full FNO layer in discrete form is thus

(FNO Discrete)$$\begin{array}{ccc}& (𝒦(v){)}_n={ℱ}^{-1}(\sum _{|\xi |\le k_{\max }}{R}_{\xi }\cdot (ℱv{)}_{\xi }{)}_n,& \text{<span class="tml-eqn"></span>}\end{array}$$

where $ℱ$ denotes the (fast) Discrete Fourier Transform and ${ℱ}^{-1}$ its inverse.

The key properties of the FNO layer that make it attractive for weather modelling are:

• Global receptive field: every grid point communicates with every other via the Fourier transform.

• Permutation equivariance: the layer is equivariant to translations of the input field.

• Resolution invariance: the spectral representation is basis-independent; a model trained at one resolution can be evaluated at another by adjusting the number of retained modes.

• Computational efficiency: the FFT costs $𝒪(N\log N)$, and the learned mixing costs $𝒪(N{k}_{\max }d)$, much less than $𝒪(N^{2}d)$ for dense attention.

Adaptive Fourier Neural Operators

FourCastNet does not use FNO directly; it uses the Adaptive Fourier Neural Operator (AFNO), introduced by Guibas et al. (2021) as part of the FNet/AFNO line of work. The AFNO modifies the FNO in two key ways:

1. Token mixing via FFT: Rather than treating each channel independently, the AFNO applies the 2D FFT across spatial dimensions, treating the $C$-dimensional feature vector at each frequency as a “token” and applying a shared MLP in frequency space.

2. Sparsity via soft-thresholding: The AFNO applies a learnable soft-thresholding operator to the Fourier coefficients, adaptively pruning frequency components that carry little information for the current input. This is the “adaptive” part.

Formally, let $𝒙\in {\mathrm{ℝ}}^{H\times W\times C}$ be the input feature map. The AFNO layer proceeds as follows:

1. Forward FFT: Compute $\widehat{𝒙}={ℱ}_{2D}(𝒙)\in {ℂ}^{H\times W\times C}$.

2. Frequency-wise MLP: For each frequency bin $(m,n)\in \{0,\dots ,H-1\}\times \{0,\dots ,W-1\}$, apply a two-layer MLP with hidden dimension $r$ (typically $r=C/4$ for efficiency): (AFNO MLP)$$\begin{array}{ccc}& {\widehat{𝒚}}_{m,n}={\mathrm{MLP}}_{\theta }({\widehat{𝒙}}_{m,n}),& \text{<span class="tml-eqn"></span>}\end{array}$$ where the MLP operates on the complex-valued channel vector at that frequency. The weights $\theta$ are shared across all frequency bins, giving translation equivariance.

3. Soft-thresholding (sparsification): Apply a learnable element-wise shrinkage: (SOFT Threshold)$$\begin{array}{ccc}& {\tilde{y}}_{m,n,c}=\mathrm{sign}({\hat{y}}_{m,n,c})\cdot \max (|{\hat{y}}_{m,n,c}|-{\lambda }_c,0),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\lambda }_c\ge 0$ is a learnable threshold per channel. This promotes sparsity in frequency space, concentrating representational energy in the most informative modes.

4. Inverse FFT: Compute $𝒚={ℱ}_{2D}^{-1}(\widetilde{𝒚})\in {\mathrm{ℝ}}^{H\times W\times C}$.

5. Residual connection: Output $𝒛=𝒙+𝒚$.

FourCastNet Architecture

FourCastNet (Pathak et al., 2022) is a vision-transformer-like architecture that replaces every attention layer with an AFNO layer. The full architecture proceeds in four stages.

Stage 1: Patch embedding.

Each atmospheric state is represented as a multi-channel field $𝑿\in {\mathrm{ℝ}}^{H\times W\times C_{\text{in}}}$ where $H=720$, $W=1440$, and $C_{\text{in}}$ is the number of input variables (20 in FourCastNet-v1, including surface variables and pressure-level fields). The field is first partitioned into non-overlapping $p\times p$ patches (patch size $p=4$ in the original model), and each patch is linearly projected to an embedding dimension $d_e$:

(Patch Embed)$$\begin{array}{ccc}& {𝑬}_{i,j}={𝑾}_e\mathrm{vec}({𝑿}_{ip:(i+1)p,jp:(j+1)p,:})+{𝒃}_e+{𝑷}_{i,j},& \text{<span class="tml-eqn"></span>}\end{array}$$

where ${𝑬}_{i,j}\in {\mathrm{ℝ}}^{d_e}$ is the embedding of patch $(i,j)$, ${𝑾}_e\in {\mathrm{ℝ}}^{d_e\times p^2{C}_{\text{in}}}$ and ${𝒃}_e\in {\mathrm{ℝ}}^{d_e}$ are the projection weights and biases, and ${𝑷}_{i,j}$ is a learnable positional embedding. After embedding, the field becomes a 2D array of tokens $𝑬\in {\mathrm{ℝ}}^{(H/p)\times (W/p)\times d_e}$.

Stage 2: AFNO blocks.

The embedded tokens are processed by $L$ sequential AFNO blocks, each consisting of:

1. Layer normalisation (pre-norm);

2. An AFNO layer performing spectral token mixing;

3. A residual connection;

4. Layer normalisation;

5. A two-layer feed-forward MLP with expansion ratio $r$ and GELU nonlinearity;

6. A residual connection.

Letting ${𝑯}^{(\ell )}\in {\mathrm{ℝ}}^{(H/p)\times (W/p)\times d_e}$ denote the hidden state after block $\ell$, the recurrence is

(AFNO Block)$$\begin{array}{cccc}& {𝑯}^{\prime (\ell )}& =\mathrm{AFNO}(\mathrm{LN}({𝑯}^{(\ell -1)}))+{𝑯}^{(\ell -1)},& \\ & {𝑯}^{(\ell )}& =\mathrm{MLP}(\mathrm{LN}({𝑯}^{\prime (\ell )}))+{𝑯}^{\prime (\ell )},& \text{<span class="tml-eqn"></span>}\end{array}$$

for $\ell =1,\dots ,L$. FourCastNet-v1 uses $L=12$ AFNO blocks with embedding dimension $d_e=768$.

Stage 3: Decoding to grid.

The final hidden state ${𝑯}^{(L)}\in {\mathrm{ℝ}}^{(H/p)\times (W/p)\times d_e}$ is projected back to the physical grid using a learned patch-expansion layer:

(Patch Expand)$$\begin{array}{ccc}& \widehat{𝑿}=\mathrm{PatchExpand}\left({𝑯}^{(L)}\right)\in {\mathrm{ℝ}}^{H\times W\times C_{\text{out}}},& \text{<span class="tml-eqn"></span>}\end{array}$$

where $C_{\text{out}}$ is the number of predicted output variables. The patch-expansion uses a transposed linear projection followed by a pixel-shuffle (sub-pixel convolution) to upsample from patch resolution to full grid resolution.

Stage 4: Autoregressive rollout.

FourCastNet is trained to predict the atmospheric state ${𝑿}_{t+\mathrm{\Delta }t}$ from ${𝑿}_t$ for a fixed time step $\mathrm{\Delta }t=6\text{h}$. For longer forecast horizons, the model is rolled out autoregressively: ${\widehat{𝑿}}_{t+k\mathrm{\Delta }t}=f_{\theta }({\widehat{𝑿}}_{t+(k-1)\mathrm{\Delta }t})$ for $k=1,2,\dots$

Training Data and Objective

Dataset.

FourCastNet is trained on the ERA5 reanalysis dataset produced by ECMWF. ERA5 provides hourly estimates of the global atmospheric state at ${0.25}^{\circ }$ horizontal resolution (approximately 28,km at the equator), with 37 pressure levels from 1000,hPa to 1,hPa. The training period spans 1979–2015, the validation set covers 2016–2017, and the test set covers 2018 onward.

The model uses 20 atmospheric channels, a sparse subset of variables and levels (not a full variable-by-level product):

• Surface fields (5): 10-metre wind components $U_{10}$, $V_{10}$, 2-metre temperature $T_{2m}$, surface pressure, and mean sea-level pressure.

• 1000,hPa (3): $U$, $V$, geopotential $Z$.

• 850,hPa (5): $T$, $U$, $V$, $Z$, relative humidity.

• 500,hPa (5): $U$, $V$, $Z$, $T$, relative humidity.

• 50,hPa (1): $Z$; and total column water vapour $TCWV$ (1).

(The full $5\times 13+5=70$ level set listed for other models is Pangu-Weather's, not FourCastNet-v1's.)

Each channel is normalised per-level to zero mean and unit variance using statistics computed from the training split, ensuring equal gradient contributions from fields with very different physical magnitudes.

Loss function.

The training objective is the weighted mean-squared error

(Fourcastnet LOSS)$$\begin{array}{ccc}& ℒ(\theta )=\frac{1}{|T_{\text{train}}|}\sum _{t\in T_{\text{train}}}\sum _{c=1}^C\sum _{i=1}^H\sum _{j=1}^W{w}_c{a}_i({\hat{X}}_{t+\mathrm{\Delta }t,i,j}^{(c)}-X_{t+\mathrm{\Delta }t,i,j}^{(c)}{)}^2,& \text{<span class="tml-eqn"></span>}\end{array}$$

where $w_c>0$ is a channel-specific weight that up-weights physically important variables (e.g., geopotential $Z_{500}$ is given higher weight than humidity), and $a_i=\cos ({\varphi }_i)/{\sum }_i\cos ({\varphi }_i)$ is the area weighting for latitude ${\varphi }_i$, compensating for the convergence of meridians toward the poles where grid cells have smaller physical area.

The area weighting is crucial: without it, gradient contributions from high-latitude cells (which are small in physical area but numerous in the regular grid) would dominate training, causing the model to over-optimise polar regions at the expense of the tropics.

Multi-step fine-tuning.

FourCastNet is first trained with the single-step loss . It is then fine-tuned with a multi-step loss that includes up to 8 autoregressive steps, backpropagating through the rollout:

(Multistep LOSS)$$\begin{array}{ccc}& {ℒ}_{\text{multi}}(\theta )=\sum _{k=1}^K{\gamma }^k{ℒ}_k(\theta ,{\widehat{𝑿}}_{t+(k-1)\mathrm{\Delta }t}),& \text{<span class="tml-eqn"></span>}\end{array}$$

where $K=8$ steps and $\gamma \in (\mathrm{0,1})$ is a discount factor that down-weights distant predictions. This fine-tuning phase is essential: without it, errors compound quickly during autoregressive rollout because the model was only trained to map ERA5 states to ERA5 states, not its own imperfect predictions to ERA5 states.

Worked Example: Seven-Day Surface Wind Forecast

Performance Against Operational Forecasting

Pathak et al. (2022) evaluated FourCastNet against the ECMWF IFS high-resolution deterministic forecast (HRES) using latitude-weighted RMSE and anomaly correlation coefficient (ACC) on the 2018 test year. Key findings:

• Geopotential $Z_{500}$: FourCastNet matches IFS HRES at lead times up to 5 days; IFS is superior at 7 days (FourCastNet shows $\sim$10% higher RMSE).

• Surface wind $U_{10}$, $V_{10}$: FourCastNet matches IFS at 1–3 days; slight degradation at 5+ days.

• Precipitation: FourCastNet performs substantially worse than IFS for total precipitation, because precipitation is poorly constrained in ERA5 (it is derived, not directly observed) and involves nonlinear convective processes that the deterministic regression objective cannot capture.

• Speed: FourCastNet inference on a single A100 GPU produces a 7-day global forecast in under 5 seconds; IFS HRES requires approximately 150,000 CPU-hours on ECMWF's Atos supercomputer. The speed-up is five to six orders of magnitude.

FourCastNet-v2 and Scalable Geometric ML

The original FourCastNet established the AFNO paradigm. Two subsequent iterations significantly extended its capability.

FourCastNet-v2 (Bonev et al., 2023).

FourCastNet-v2 replaces the AFNO token mixer with a spherical Fourier neural operator (SFNO) that operates natively on the sphere ${𝕊}^2$ using spherical harmonics as the Fourier basis. The motivation is geometric: the regular latitude-longitude grid treats the sphere as a flat rectangle, introducing distortions near the poles. Spherical harmonics $Y_{\ell }^m(\varphi ,\lambda )$ (where $\ell$ is the degree, $m$ the order, $\varphi$ the latitude, and $\lambda$ the longitude) are the natural basis for functions on the sphere, satisfying ${\nabla }_{{𝕊}^2}^2{Y}_{\ell }^m=-\ell (\ell +1)Y_{\ell }^m$.

The SFNO layer replaces the 2D FFT with the Spherical Harmonic Transform (SHT) and its inverse, and replaces the frequency-wise MLP with a degree-wise (rotationally equivariant) mixing:

(SFNO MIX)$$\begin{array}{ccc}& {\hat{y}}_{\ell }^m=\sum _{{\ell }^{\prime }}{W}_{\ell ,{\ell }^{\prime }}{\hat{x}}_{{\ell }^{\prime }}^m,& \text{<span class="tml-eqn"></span>}\end{array}$$

where $W_{\ell ,{\ell }^{\prime }}\in {\mathrm{ℝ}}^{d_e\times d_e}$ is a learnable coupling matrix between spherical harmonic degrees $\ell$ and ${\ell }^{\prime }$. This formulation is equivariant to $SO(3)$ rotations of the sphere, a property the original Cartesian AFNO lacks. FourCastNet-v2 achieves improved accuracy at polar latitudes and better long-range teleconnection representation.

FourCastNet 3: Geometric Foundation Models.

Subsequent work (2024–2025) extended FourCastNet into a geometric foundation model framework, in which a single large model pre-trained on diverse atmospheric data (ERA5, MERRA-2, JRA-55, satellite radiances) can be fine-tuned for downstream tasks including:

• Medium-range deterministic forecasting (days 1–10);

• Subseasonal-to-seasonal prediction (weeks 2–8);

• Downscaling to high-resolution regional grids via learned geometric interpolation operators;

• Data assimilation, generating analysis states from sparse observations.

The architecture scales from $\sim$70M parameters (FourCastNet-v1) to $\sim$2B parameters (FourCastNet 3), with performance improvements following approximate power laws in parameter count, mirroring scaling behaviour observed in large language models.

Exercises

Pangu-Weather and 3D Transformer Architectures

In 2023, a paper published in Nature reported that a neural network had, for the first time, beaten the world's best operational numerical weather prediction system-ECMWF's Integrated Forecasting System-on over 90% of variables at all tested lead times. The model was called Pangu-Weather, developed by researchers at Huawei Cloud Research and Huawei Noah's Ark Laboratory, led by Kaifeng Bi and colleagues. The result was not a marginal improvement; at 5-day lead times, Pangu-Weather reduced the RMSE for 500-hPa geopotential by approximately 3% over IFS HRES while running ten thousand times faster.

The architectural leap that enabled this result was three-dimensionality. Where FourCastNet processed the atmosphere as a 2D field of horizontal patches-treating vertical levels as separate channels-Pangu-Weather introduced a genuine three-dimensional transformer that attends simultaneously over latitude, longitude, and pressure level. This apparently minor change has profound consequences: it allows the model to learn the vertical coupling that governs phenomena as diverse as tropical cyclone development, jet-stream variability, and the propagation of planetary waves from the stratosphere into the troposphere.

From Vision Transformers to Weather Transformers

Pangu-Weather's architecture has a clear lineage in computer vision.

Vision Transformer (ViT).

Dosovitskiy et al. (2021) introduced the ViT by treating an image as a sequence of non-overlapping patches and applying standard transformer self-attention. The key insight was that the inductive biases of convolutional networks (translation equivariance, local connectivity) could be abandoned in favour of global attention, with the model free to learn its own spatial inductive biases from data. ViTs achieve state-of-the-art performance on ImageNet with sufficient data, outperforming CNNs of comparable parameter count.

Swin Transformer.

Liu et al. (2021) observed that global attention in ViT scales as $𝒪(N^2)$. The Swin Transformer replaces global attention with shifted window (SW) attention: each token attends only to tokens within a local $w\times w$ window, reducing cost to $𝒪(N{w}^2/w^2)=𝒪(N)$ for fixed window size $w$. A “shift” mechanism introduces cross-window communication by shifting the window partition between consecutive layers. Swin won the ImageNet-22K top-1 accuracy benchmark in 2021 and became the dominant backbone for dense prediction tasks (detection, segmentation).

3D Swin Transformer.

The natural extension to video and volumetric data replaces 2D windows with 3D windows spanning height $\times$ width $\times$ time (or depth). Liu et al. (2022) introduced Video Swin for video recognition. Pangu-Weather adapts this idea to atmospheric data, where the three dimensions are latitude $\times$ longitude $\times$ pressure level.

The 3D Earth-Specific Transformer (3DEST)

The core of Pangu-Weather is the 3D Earth-Specific Transformer (3DEST), which differs from a standard 3D Swin Transformer in two critical respects: (1) it uses an Earth-specific positional bias that accounts for the non-uniform physical spacing of pressure levels and the spherical geometry of the planet; and (2) it processes upper-air and surface variables in separate streams that interact only at the decoding stage.

Input representation.

The full atmospheric state at a given time is represented as two tensors:

• Upper-air tensor: ${𝖷}_{\text{upper}}\in {\mathrm{ℝ}}^{13\times H\times W\times C_{\text{upper}}}$, where 13 is the number of pressure levels (1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 50,hPa), $H=721$, $W=1440$, and $C_{\text{upper}}=5$ channels (geopotential $Z$, temperature $T$, zonal wind $U$, meridional wind $V$, specific humidity $Q$).

• Surface tensor: ${𝖷}_{\text{surf}}\in {\mathrm{ℝ}}^{H\times W\times C_{\text{surf}}}$, where $C_{\text{surf}}=4$ channels (mean sea-level pressure $P_{\text{msl}}$, 10m wind $U_{10}$, $V_{10}$, and 2m temperature $T_{2m}$).

Patch embedding.

The upper-air tensor is partitioned into $2\times 4\times 4$ patches (2 pressure levels $\times$ 4 latitude points $\times$ 4 longitude points), giving $(13/2)\times (721/4)\times (1440/4)\approx 6\times 180\times 360=\mathrm{388,800}$ tokens. Each token is a $2\times 4\times 4\times 5=160$-dimensional vector, linearly projected to embedding dimension $d_e$. The surface tensor uses $1\times 4\times 4$ patches.

3D window attention.

The 3DEST divides the token grid into non-overlapping 3D windows of size $W_z\times W_y\times W_x$ (with typical values $W_z=2$, $W_y=6$, $W_x=12$). Standard self-attention is applied within each window:

(3D ATTN)$$\begin{array}{ccc}& \mathrm{Attn}(𝐐,𝐊,𝐕)=\mathrm{softmax}(\frac{𝐐{𝐊}^{\top }}{\sqrt{d}}+B)𝐕,& \text{<span class="tml-eqn"></span>}\end{array}$$

where $𝐐,𝐊,𝐕\in {\mathrm{ℝ}}^{n_w\times d}$ are query, key, and value matrices for the $n_w=W_z{W}_y{W}_x$ tokens in the window, and $B$ is the Earth-specific relative position bias defined below.

The cost of 3D window attention is $𝒪(N{n}_w)$ rather than $𝒪(N^2)$; for fixed window size $n_w$ this is linear in $N$. Cross-window communication is achieved by shifting the window partition by $(W_z/2,W_y/2,W_x/2)$ in alternating layers.

Earth-Specific Positional Bias

The Earth-specific positional bias $B$ in equation is the central innovation of Pangu-Weather. It encodes three geometric facts about the atmosphere that generic positional encodings ignore:

1. Non-uniform pressure level spacing: The 13 pressure levels are not uniformly spaced in altitude. In geometric altitude, the spacing between adjacent levels varies from $\sim$0.6,km (1000 to 925,hPa) to $\sim$5.2,km (100 to 50,hPa, in the lower stratosphere). Pairs of tokens at the same horizontal distance but different vertical positions should have different attention biases reflecting their actual physical proximity.

2. Latitude-dependent longitude spacing: At latitude $\varphi$, one degree of longitude corresponds to physical distance $111\text{km}\cdot \cos (\varphi )$. Two tokens separated by the same number of grid points in the east-west direction are physically much closer near the poles than near the equator.

3. Rotational symmetry: The atmosphere is approximately rotationally symmetric about the Earth's axis (in the climatological mean), so biases should depend on latitude but not longitude.

Let $(z_i,y_i,x_i)$ and $(z_j,y_j,x_j)$ denote the vertical, meridional, and zonal indices of tokens $i$ and $j$ within a window. The relative position is $(\delta z,\delta y,\delta x)=(z_i-z_j,y_i-y_j,x_i-x_j)$. The Earth-specific bias is a learnable function

(ESPB)$$\begin{array}{ccc}& B_{ij}=b(\delta z,\delta y,\delta x\cdot \cos ({\varphi }_{y_i}),{\varphi }_{y_i}),& \text{<span class="tml-eqn"></span>}\end{array}$$

where ${\varphi }_{y_i}$ is the latitude of the window row containing token $i$ and the longitude displacement is scaled by $\cos ({\varphi }_{y_i})$ to account for the physical narrowing of longitude spacing at higher latitudes. The function $b$ is parametrised as a lookup table over quantised offsets, learnable during training.

Lead-Time-Specific Models

The most architecturally distinctive decision in Pangu-Weather is the use of four separate models, each optimised for a different forecast lead time:

• $f_1$: trained to predict the state 1 hour ahead;

• $f_3$: trained to predict the state 3 hours ahead;

• $f_6$: trained to predict the state 6 hours ahead;

• $f_{24}$: trained to predict the state 24 hours ahead.

For an arbitrary forecast horizon $\tau$, the models are combined hierarchically: $\tau$ is decomposed into a sum of 1, 3, 6, and 24-hour increments (analogous to expressing $\tau$ in a “time-base-24” numeral system), and the corresponding models are applied in sequence. For example, a 72-hour forecast is produced as $f_{24}\circ f_{24}\circ f_{24}$ (three applications of $f_{24}$), a 7-hour forecast as $f_6\circ f_1$, and a 37-hour forecast as $f_{24}\circ f_6\circ f_6\circ f_1$.

Why multiple models?

The key insight is that autoregressive error accumulation is dominated by the number of model applications, not the total lead time. A single 6-hour model applied 16 times to produce a 4-day forecast accumulates error through 16 composition steps. The same 4-day forecast from $f_{24}$ requires only 4 composition steps. Since each model application introduces systematic biases (the model output differs from ERA5 in systematic ways), reducing the number of compositions dramatically reduces error accumulation.

More formally, let $ϵ$ denote the per-step error standard deviation and suppose errors accumulate linearly (a conservative upper bound). The RMSE after $K$ steps is at most $Kϵ$. A single 6-hour model applied $K=16$ times achieves RMSE $\le 16ϵ$. The 24-hour model applied $K=4$ times achieves RMSE $\le 4{ϵ}^{\prime }$, where ${ϵ}^{\prime }$ may be larger than $ϵ$ (harder single-step prediction) but the reduced $K$ typically dominates.

Hierarchical time aggregation.

The four-model hierarchy can be viewed as a mixed-radix number system for time. Let $\tau$ be the desired forecast horizon in hours. Write

(Pangu Decompose)$$\begin{array}{ccc}& \tau =24n_{24}+6n_6+3n_3+1n_1,& \text{<span class="tml-eqn"></span>}\end{array}$$

where $n_{24},n_6,n_3,n_1\ge 0$ are non-negative integers chosen to minimise $n_{24}+n_6+n_3+n_1$ subject to satisfying the equation. The greedy algorithm (fill with the largest unit first) achieves this minimum. The forecast is then

(Pangu Compose)$$\begin{array}{ccc}& {\widehat{𝑿}}_{t+\tau }=\underset{\text{1-h steps}}{\underbrace{f_1^{\circ n_1}}}\circ \underset{\text{3-h steps}}{\underbrace{f_3^{\circ n_3}}}\circ \underset{\text{6-h steps}}{\underbrace{f_6^{\circ n_6}}}\circ \underset{\text{24-h steps}}{\underbrace{f_{24}^{\circ n_{24}}}}({𝑿}_t).& \text{<span class="tml-eqn"></span>}\end{array}$$

Training Objective and Data

Each of the four Pangu-Weather models is trained independently on the ERA5 reanalysis using a weighted mean-absolute error (MAE) loss. Unlike FourCastNet, which uses MSE, Pangu-Weather adopts MAE because it is less sensitive to the large outliers introduced by extreme weather events, produces sharper predictions of wind and pressure features, and aligns better with the meteorological practice of evaluating forecasts by absolute error in physically meaningful units.

Training loss.

For model $f_{\mathrm{\Delta }}$ with lead time $\mathrm{\Delta }$:

(Pangu LOSS)$$\begin{array}{ccc}& {ℒ}_{\mathrm{\Delta }}(\theta )=\frac{1}{|T_{\text{train}}|}\sum _{t\in T_{\text{train}}}\sum _{c=1}^C\sum _{l=1}^{13}\sum _{i=1}^H\sum _{j=1}^W{w}_{c,l}{a}_i|{\hat{X}}_{t+\mathrm{\Delta },i,j}^{(c,l)}-X_{t+\mathrm{\Delta },i,j}^{(c,l)}|,& \text{<span class="tml-eqn"></span>}\end{array}$$

where $w_{c,l}$ is a variable-and-level-specific weight (empirically set to balance gradient magnitudes across physically diverse channels), $a_i=\cos ({\varphi }_i)$ is the latitude area weight, the outer sum runs over all $C_{\text{upper}}\times 13+C_{\text{surf}}=69$ channels (5 upper-air $\times$ 13 levels + 4 surface), and the inner sums run over the $721\times 1440$ grid.

Data split.

ERA5 (1979–2017) is used for training with training on 1979–2017, validation on 2019, and testing on 2018. A consistent split is used for all four models to ensure comparability.

Optimisation.

Each model is trained with AdamW (${\beta }_1=0.9$, ${\beta }_2=0.999$, weight decay $=0.05$) for 200 epochs with a cosine annealing learning rate schedule, starting from a peak learning rate of $5\times {10}^{-4}$. Each model has approximately 200 million parameters. Training each model on 192 NVIDIA A100 GPUs requires approximately 16 days. Despite this training cost, inference for a 7-day forecast takes approximately 1.4 seconds on a single A100.

Pangu-Weather Architecture: TikZ Illustration

Performance: Beating IFS on 90% of Variables

Bi et al. (2023) evaluated Pangu-Weather against IFS HRES on the 2018 ERA5 test set using both latitude-weighted RMSE and ACC. The headline result is remarkable:

• Pangu-Weather outperforms IFS HRES on over 90% of all variable-level combinations at lead times from 1 day to 7 days.

• For the critical metric $Z_{500}$ (500-hPa geopotential, a primary measure of large-scale weather pattern skill), Pangu-Weather achieves approximately 3% lower RMSE than IFS at 5-day lead time.

• Pangu-Weather's ACC (anomaly correlation coefficient) for $Z_{500}$ exceeds 0.6 at day 7, the traditional “useful forecast” threshold, whereas IFS falls below 0.6 at approximately day 6.5 in the test year.

• Inference time for a 7-day global forecast: approximately 1.4 seconds on a single A100 GPU. IFS HRES requires approximately 14,000 CPU-hours on ECMWF's Atos supercomputer. The speed-up factor is $\approx {10}^4$ (ten thousand times faster).

Variable-by-variable breakdown.

The improvement is not uniform. Pangu-Weather most dramatically outperforms IFS for:

• Upper-tropospheric temperature (200,hPa $T$): RMSE reduced by $\sim$8% versus IFS at 5 days.

• Geopotential at most pressure levels (500, 850, 1000,hPa $Z$).

• Meridional wind at multiple levels.

IFS retains advantages (or ties) for:

• Precipitation (Pangu does not directly predict precipitation).

• Boundary-layer variables in complex terrain (Pangu is on a regular grid; IFS uses terrain-following levels).

• Near-surface wind in regions of strong orographic forcing.

Case Study: Hurricane Lee (2023)

The 2023 Atlantic hurricane season provided a high-profile operational test of AI weather models. Hurricane Lee, which formed on 5 September 2023 in the tropical Atlantic, intensified rapidly and tracked across the Caribbean toward Atlantic Canada.

The Smoothing Caveat: Extreme Events and ERA5 Inheritance

Pangu-Weather inherits a fundamental limitation from its training data: ERA5 is a reanalysis product, not a direct observation. ERA5 is produced by 4D-Var data assimilation that blends a prior from the ECMWF IFS model with available observations. This process is regularising: it smooths sharp features in the analysis, penalises departures from the model climatology, and represents phenomena below the assimilation resolution imperfectly.

The smoothing issue motivates two active research directions: (1) post-processing Pangu outputs with diffusion models or quantile regression to recover physically plausible sub-grid distributions of extreme values; and (2) training weather AI models directly on observation data (station networks, radar, satellite) using data-assimilation-inspired loss functions that account for observation operators and error covariances.

Connection to Vision Transformers: Transfer and Adaptation

The architectural similarity between Pangu-Weather and Swin-based vision models raises the question: can ImageNet-pretrained weights help weather AI? The answer is nuanced.

What transfers.

Low-level features in ViTs (first few attention layers) detect oriented edges and patches of coherent colour or intensity-analogous to the horizontal gradients and vorticity filaments in atmospheric fields. Pre-training on large vision datasets can provide a useful initialisation for the attention weights, reducing the number of gradient steps needed to reach good weather performance.

What does not transfer.

The Earth-specific positional bias is domain-specific and has no ImageNet equivalent. The physical semantics of channels (geopotential, temperature, humidity) differ fundamentally from image RGB channels. The temporal structure of weather (quasi-periodic diurnal and annual cycles, teleconnections) has no direct analogue in static image datasets.

In practice, Pangu-Weather is trained from scratch rather than initialised from ViT weights. The training set (ERA5, $\sim$40 years $\times$ 1460 analysis times per year = $\sim$58,000 snapshots) is small by deep learning standards, but the physical constraints imposed by the loss function and the ESPB provide sufficient regularisation.

Cross-architecture lessons.

The Pangu-Weather success validates a design principle from computer vision: hierarchical, windowed self-attention with local-to-global spatial coupling (the Swin paradigm) generalises well beyond natural images. The critical adaptations required for atmospheric science-3D tokenisation, Earth-specific bias, lead-time-specific models-are relatively small changes to the Swin blueprint, suggesting that future advances in vision transformers will continue to benefit weather AI, and vice versa.

Historical Note: Huawei's Entry into Weather AI

Exercises

GraphCast and Graph Neural Networks

Weather prediction is, at its core, a problem of propagating information across a sphere. Temperature anomalies in the Pacific influence European precipitation weeks later. A tropical cyclone's track depends on hemispheric-scale pressure gradients. Jet-stream meanders couple the stratosphere to the surface boundary layer. The fundamental challenge for any data-driven forecasting architecture is therefore how to represent and process spatially coherent signals on a curved, planetary domain.

The approaches studied in earlier sections-Pangu-Weather (Section Pangu-Weather and 3D Transformer Architectures) and FourCastNet (Section FourCastNet and the Spectral Revolution)-answered this question with variants of the Transformer operating on regular latitude-longitude (lat-lon) grids, borrowing the vision-transformer paradigm of dividing the globe into patches and computing attention between them. This choice is natural: ERA5 data is natively delivered on lat-lon grids, and the rectangular-grid paradigm inherits decades of software infrastructure.

Yet lat-lon grids carry a hidden flaw that becomes catastrophic at the poles. On a regular ${0.25}^{\circ }\times {0.25}^{\circ }$ grid, meridional grid lines converge toward the poles, creating cells whose east-west spacing shrinks from roughly 28 km at the Equator to less than 1 km at ${89}^{\circ }$N. Numerical instabilities accumulate at these singularities. Convolutional and attention operations that assume uniform cell areas must apply ad-hoc area-weighting corrections that do not compose cleanly with the rest of the computational graph. The geometry is, in a precise sense, wrong.

In 2023, Remi Lam and colleagues at Google DeepMind published GraphCast: Learning Skillful Medium-Range Global Weather Forecasting in Science, demonstrating that an entirely different approach could outperform the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) HRES operational model on the vast majority of standard verification targets. The key insight was deceptively simple: replace the lat-lon grid with a multi-mesh icosahedron, a hierarchy of increasingly refined triangulations of the sphere that tile the globe with nearly-equal-area cells, and process information on this mesh using a graph neural network (GNN). The result was the first data-driven system to systematically surpass HRES, doing so in a fraction of the computational cost, in roughly 60 seconds on a single Google TPU v4 chip.

Why Graphs? The Failure Mode of Lat-Lon Grids

To appreciate why a graph formulation helps, consider the computational representation of a function $f:{𝕊}^2\to {ℝ}^d$ on the sphere. The standard meteorological choice is a regular lat-lon grid with resolution $\mathrm{\Delta }\lambda \times \mathrm{\Delta }\varphi$ (longitude $\times$ latitude), producing $N_{\lambda }\times N_{\varphi }$ grid cells. For ERA5 at ${0.25}^{\circ }$ resolution, this gives $1440\times 721\approx {10}^6$ cells.

The area of cell $(i,j)$ at latitude ${\varphi }_j$ is (Latlon AREA)$$\begin{array}{ccc}& A_{ij}=R_{\oplus ︎}^2\mathrm{\Delta }\lambda \mathrm{\Delta }\varphi \cos {\varphi }_j,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $R_{\oplus ︎}\approx 6371\mathrm{km}$ is Earth's radius. As $|{\varphi }_j|\to {90}^{\circ }$ we have $\cos {\varphi }_j\to 0$ and the cell area collapses to zero, while the number of cells per latitude circle remains fixed. A $3\times 3$ convolution kernel applied at ${89.75}^{\circ }$N integrates information from cells spanning radically different physical distances: neighbouring cells in the east-west direction are only $\sim 100\mathrm{m}$ apart, while cells in the north-south direction are $\sim 28\mathrm{km}$ apart. The implicit assumption of uniform spatial scale is violated by three orders of magnitude.

Gauss-Legendre grids and reduced Gaussian grids (used by ECMWF) partially address this by halving the number of longitude points at high latitudes, but they introduce irregular connectivity that complicates convolutional operations. Spherical harmonic representations avoid the singularity entirely but sacrifice the locality of spatial operations.

An icosahedral tessellation takes a different approach. Starting from a regular icosahedron inscribed in the unit sphere (20 equilateral triangular faces, 12 vertices, 30 edges), each triangular face is subdivided into $4^k$ smaller triangles by connecting the midpoints of edges and projecting them onto the sphere. After $k$ refinements, the mesh has $20\times 4^k$ faces, $10\times 4^k+2$ vertices, and $30\times 4^k$ edges. Crucially, all triangles have nearly equal areas: the ratio of the largest to smallest face area is bounded by a constant independent of refinement level, typically around 1.03 for standard icosahedral subdivision.

This near-uniformity has profound implications. A graph convolution or message-passing operation applied at any node integrates information from geometrically similar neighbourhoods regardless of latitude. There are no poles, no singularities, no ad-hoc area-weighting corrections. The mesh is topologically a graph $𝒢=(𝒱,ℰ)$ with vertex set $𝒱$ and edge set $ℰ$, and graph neural network machinery applies directly.

Multi-Mesh Icosahedron: Six Refinement Levels

GraphCast employs not one but six icosahedral meshes at different resolution levels, labelled ${ℳ}_0$ through ${ℳ}_5$. Mesh ${ℳ}_k$ is obtained by applying $k$ rounds of icosahedral subdivision. The number of vertices and mean edge lengths are summarised in Table Table 5.

Level $k$	Subdivisions	$N_V$	$N_E$	Mean edge (km)
0	0	12	30	7,054
1	1	42	120	3,765
2	2	162	480	1,914
3	3	642	1,920	961
4	4	2,562	7,680	481
5	5	10,242	30,720	241

The finest mesh ${ℳ}_5$ has a mean inter-node spacing of roughly 241 km, coarser than the ${0.25}^{\circ }$ ERA5 grid resolution at mid-latitudes (approximately 28 km); the encoder still ingests the full ERA5 grid, so fine-scale information is not discarded. However, the mesh nodes are not placed at ERA5 grid points; instead, the encoder maps from the lat-lon ERA5 grid to the icosahedral mesh, and the decoder maps back.

The multi-mesh structure means that the processor graph ${𝒢}_{\mathrm{proc}}$ is the union of all six icosahedral meshes: (Multimesh Union)$$\begin{array}{ccc}& {𝒢}_{\mathrm{proc}}=\bigcup _{k=0}^5{ℳ}_k,& \text{<span class="tml-eqn"></span>}\end{array}$$ with nodes ${𝒱}_{\mathrm{proc}}={𝒱}_5$ (the finest level's vertices, which subsume all coarser levels' vertices) and edges ${ℰ}_{\mathrm{proc}}={\bigcup }_{k=0}^5{ℰ}_k$. This is the crucial design: by including edges from all six refinement levels in the processor graph, each message-passing step can propagate information over a range of spatial scales simultaneously. A single message-passing step on an edge from ${ℳ}_0$ connects vertices 2,559 km apart, while edges from ${ℳ}_5$ provide fine-grained local refinement. This multi-scale connectivity replaces the explicit down-sampling and up-sampling of U-Net-style architectures with a more organic, graph-theoretic analogue.

Encoder-Processor-Decoder Architecture

GraphCast adopts a three-stage architecture that cleanly separates the concerns of grid-to-mesh mapping, atmospheric dynamics processing, and mesh-to-grid reconstruction. Figure fig:weather:graphcast:arch illustrates the pipeline.

Each of the three stages is a graph neural network with a different bipartite or homogeneous graph structure.

Encoder (Grid-to-Mesh).

The encoder maps features from ERA5 grid nodes $v_g\in {𝒱}_g$ to mesh nodes $v_m\in {𝒱}_5$. It is a single-layer bipartite GNN defined by a bipartite graph ${𝒢}_{g\to m}$ whose edges connect each grid node to the nearest $K=4$ mesh nodes. Let ${𝐡}_g^{(0)}$ denote the initial feature embedding of grid node $g$ (including atmospheric variables, static fields like orography and land-sea mask, and sinusoidal positional encodings of latitude, longitude, and time of year). The mesh node features after encoding are (Encoder)$$\begin{array}{ccc}& {𝐡}_m^{\mathrm{enc}}={\varphi }_m({𝐡}_m^{(0)},\underset{g:(g,m)\in {ℰ}_{g\to m}}{⨁}{\varphi }_e({𝐡}_g^{(0)},{𝐡}_m^{(0)},{𝐚}_{gm})),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\varphi }_m$ and ${\varphi }_e$ are learned MLPs, ${𝐚}_{gm}$ are edge features (3D unit vector from $g$ to $m$, arc length, etc.), and $⨁$ denotes sum aggregation.

Processor (Mesh-to-Mesh).

The processor is the heart of GraphCast. It applies $T_{\mathrm{proc}}=16$ rounds of message passing on the multi-mesh graph ${𝒢}_{\mathrm{proc}}$. Each round implements a standard edge-update then node-update GNN step (detailed in Section Graph Neural Network for Weather Forecasting).

Decoder (Mesh-to-Grid).

The decoder is the mirror image of the encoder: a single-layer bipartite GNN on a graph ${𝒢}_{m\to g}$ whose edges connect each mesh node $m$ to the nearest $K^{\prime }$ grid nodes. It maps processed mesh features back to ERA5 grid locations, predicting the residual (change in state) rather than the absolute next state: (Decoder Residual)$$\begin{array}{ccc}& {\widehat{𝐱}}_{t+6}={𝐱}_t+\mathrm{Decoder}\left({𝐡}^{(T_{\mathrm{proc}})}\right).& \text{<span class="tml-eqn"></span>}\end{array}$$ Predicting residuals rather than absolute states improves training stability: the model only needs to learn the change over 6 hours, which is small relative to the typical magnitudes of atmospheric variables.

Graph Neural Network for Weather Forecasting

We now give the formal definition of the graph neural network used within GraphCast's processor.

The residual connections in and are critical: without them, 16 layers of message passing would suffer from over-smoothing and vanishing gradients. GraphCast uses $d_n=d_e=512$ and two-hidden-layer MLPs ($512\to 512\to 512$) with ReLU activations followed by layer normalisation.

Training on ERA5: Autoregressive Rollout and Loss

Data: ERA5 at ${0.25}^{\circ }$, 37 Years

GraphCast is trained on ERA5 reanalysis data spanning 1979 to 2017 (approximately 37 years), held out: 2018 for validation, 2019 and 2020 for testing. The model operates at ${0.25}^{\circ }$ horizontal resolution with 37 pressure levels (from 1 hPa to 1000 hPa) and six surface variables (2-metre temperature, 10-metre winds $u$ and $v$, mean sea-level pressure, total precipitation, and total column water vapour; four of these enter the state vector). The full input-output dimensionality is (DIM)$$\begin{array}{ccc}& d=(37\times 5+4)\times 1440\times 721\approx 1.98\times {10}^8\text{ values per time step.}& \text{<span class="tml-eqn"></span>}\end{array}$$

Two-State Input

A fundamental design choice is to provide two consecutive atmospheric states as input: ${𝐱}_{t-6\mathrm{h}}$ and ${𝐱}_t$. This gives the model access to the instantaneous state and its recent tendency (the 6-hour difference ${𝐱}_t-{𝐱}_{t-6\mathrm{h}}$), mimicking the “two time levels” scheme of finite-difference NWP schemes that must approximate first-order time derivatives. Without this second time level, the model would need to infer the effective “velocity” of the atmospheric state from a single snapshot, a harder problem.

Autoregressive Rollout Training

Rather than training only on one-step predictions, GraphCast performs autoregressive rollout training: the model is unrolled for $K$ steps, and the loss is computed as the sum over all rollout steps. For a rollout of length $K$ starting at time $t$: (Rollout)$$\begin{array}{ccc}& {\widehat{𝐱}}_{t+k\cdot 6\mathrm{h}}=f_{\theta }({\widehat{𝐱}}_{t+(k-1)\cdot 6\mathrm{h}},{\widehat{𝐱}}_{t+(k-2)\cdot 6\mathrm{h}}),k=1,\dots ,K,& \text{<span class="tml-eqn"></span>}\end{array}$$ with ${\widehat{𝐱}}_{t-6\mathrm{h}}={𝐱}_{t-6\mathrm{h}}$ and ${\widehat{𝐱}}_t={𝐱}_t$ (ground truth initialisation for the first step). The published GraphCast model was trained with increasing rollout length: beginning with single-step ($K=1$) pre-training and progressing to $K=12$ fine-tuning (three days of autoregressive rollout), one step at a time.

Normalised MSE Loss

The loss function must treat all variables comparably despite enormous differences in magnitude and variance (e.g., temperature in Kelvin vs. specific humidity in kg/kg). GraphCast uses a normalised mean squared error weighted by pressure level and variable: (LOSS)$$\begin{array}{ccc}& ℒ=\frac{1}{K}\sum _{k=1}^K\frac{1}{|{𝒱}_g|}\sum _{v\in {𝒱}_g}\sum _{\ell \in 𝒫}\sum _{c\in 𝒞}{w}_{\ell }{w}_c{\left(\frac{{\hat{x}}_{v,\ell ,c}^{(k)}-x_{v,\ell ,c}^{(k)}}{{\sigma }_{\ell ,c}}\right)}^2,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝒫$ is the set of pressure levels, $𝒞$ the set of predicted channels, $w_{\ell }$ is the pressure-level weight (proportional to the pressure thickness of the layer), $w_c$ is the per-variable weight, and ${\sigma }_{\ell ,c}$ is the temporal standard deviation of channel $c$ at level $\ell$ computed over the training set. The normalisation by ${\sigma }_{\ell ,c}$ ensures that variables with small natural variability (e.g., total column ozone) are not overwhelmed by variables with large variability (e.g., near-surface wind speed).

An area-weighting factor $a_v\propto \cos ({\varphi }_v)$ (where ${\varphi }_v$ is the latitude of ERA5 grid point $v$) is also applied to correct for the unequal areas of lat-lon cells, so polar regions do not dominate the loss despite having more grid points per unit area.

Operational Results: Surpassing HRES

The headline result of the Lam et al. (2023) paper is striking: GraphCast outperformed ECMWF HRES on 90.4% of the 1380 verification targets used in the standard WMO evaluation framework. The 1380 targets factor as $69\times 20$:

• 5 atmospheric variables (geopotential, temperature, $u$-wind, $v$-wind, specific humidity) at 13 pressure levels (50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000 hPa), that is 65,

• plus 4 surface variables, giving 69 variable-level combinations,

• at 20 lead times from $12\mathrm{h}$ to $10\times 24\mathrm{h}=240\mathrm{h}$ in 12-hour steps.

GraphCast showed the largest skill improvements at medium range (3–7 days), particularly for upper-tropospheric geopotential (Z500) and wind, which are most relevant for synoptic-scale pattern forecasting. At short range (0–48 h), high-resolution NWP still held advantages on fine-scale convective precipitation.

Lead time range	Fraction (%) where GraphCast $>$ HRES
12–48 h (short range)	81
60–120 h (medium range)	95
132–240 h (extended)	88
All 1380 targets (reported)	90.3

A particularly noteworthy achievement was GraphCast's skill at predicting the track of tropical cyclones. Tropical cyclone track prediction requires accurate representation of the large-scale steering flow over multiple days, a task where medium-range skill is paramount. GraphCast demonstrated competitive cyclone track skill even before any task-specific fine-tuning, a strong validation of the general representational capacity of the model.

Open Source Release and Google DeepMind's Vision

In November 2023, Google DeepMind released GraphCast as open source, publishing model weights, code, and evaluation benchmarks on GitHub. This release was significant: it was the first time a state-of-the-art MLWP system was made freely available for research and operational experimentation, catalysing a wave of derivative work that adapted GraphCast to regional forecasting, ensemble generation, and climate downscaling.

The open release reflected a deliberate strategic choice by Google DeepMind. The AlphaFold release in 2021 had demonstrated that open publication of model weights (coupled with the AlphaFold Protein Structure Database) could catalyse an entire scientific field. The GraphCast open release carried a similar ambition: to accelerate weather science by providing the community with a production-quality forecast tool that could be fine-tuned on regional or application-specific data without the computational cost of training from scratch.

DeepMind has articulated a broader vision of AI for science in which deep learning serves as a universal approximator for complex physical systems, replacing or augmenting expensive first-principles simulations. GraphCast, AlphaFold, AlphaTensor (for discovering fast matrix multiplication algorithms), and AlphaCode represent different facets of this programme: applying large-scale deep learning and reinforcement learning to problems that have traditionally been the domain of physics, biology, and mathematics.

Exercises

Extending Horizons: FuXi, FengWu, NeuralGCM

The publication of GraphCast in late 2022 (preprint) and late 2023 (Science) (Science) ignited a global race. Within months, teams from Fudan University, Shanghai AI Lab, and Google Research had published competing architectures that matched or exceeded GraphCast on different aspects of the evaluation benchmark. This section examines three of the most scientifically influential: FuXi from Fudan University, FengWu from the Shanghai AI Lab, and NeuralGCM from Google Research, each of which contributes a distinct conceptual advance.

FuXi addresses the fundamental tension in extended-range forecasting: a single model optimised for short-range accuracy tends to become oversmoothed at long range, while a model optimised for long-range plausibility may sacrifice short-range detail. FuXi resolves this through a cascaded architecture of specialised models with explicit handoffs. FengWu addresses the catastrophic forgetting problem that emerges when autoregressive models are trained with very long rollouts, introducing a replay buffer and uncertainty-weighted multi-task loss. NeuralGCM takes the most radical departure: it is not a pure ML model but a hybrid system that combines a differentiable atmospheric dynamical core with learned parameterisations, bridging the gap between MLWP and climate modelling.

FuXi: Cascaded Architecture for 15-Day Forecasts

FuXi (Chen et al., 2023, Fudan University) takes its name from the legendary Chinese deity associated with cosmic order and the origins of civilisation-a fitting metaphor for a model that seeks to impose structure on the chaotic atmosphere over unprecedented lead times. The key technical contribution is a cascaded architecture consisting of three specialised models:

• FuXi-Short: lead times 0–5 days (0–120 h),

• FuXi-Medium: lead times 5–10 days (120–240 h),

• FuXi-Long: lead times 10–15 days (240–360 h).

Each model is a U-Transformer built on SwinTransformer V2 blocks (described below), trained on ERA5 at ${0.25}^{\circ }$ resolution with 37 pressure levels.

U-Transformer with SwinTransformer V2 Blocks

The base model for each of the three FuXi components is a U-Transformer: a hierarchical encoder-decoder architecture in the style of U-Net, but with SwinTransformer V2 blocks replacing convolutional layers.

SwinTransformer V2 recap.

SwinTransformer V2 (Liu et al., 2022) computes shifted-window self-attention over local $w\times w$ windows of patches. Given an input feature map of shape $H\times W\times C$ (height, width, channels), the map is partitioned into non-overlapping windows of $w\times w$ patches and self-attention is computed independently within each window: (WMSA)$$\begin{array}{ccc}& \text{W-MSA}(𝐐,𝐊,𝐕)=\mathrm{softmax}(\frac{𝐐{𝐊}^{\top }}{\sqrt{d}}+𝐁)𝐕,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝐁$ is a learned relative position bias and $d$ is the head dimension. By alternating between regular and shifted-window partitions, consecutive SwinTransformer blocks achieve cross-window information exchange without computing global attention.

U-Transformer design.

The U-Transformer encoder reduces the spatial resolution through patch merging (analogous to strided convolution) while increasing channel depth, and the decoder up-samples through patch splitting (analogous to transposed convolution) with skip connections. For FuXi at ERA5 ${0.25}^{\circ }$ resolution, the encoder operates on an input of $721\times 1440\times C_0$ and produces representations at multiple scales ($721\times 1440$, $361\times 720$, $181\times 360$), all processed by SwinTransformer V2 blocks.

Cascaded Handoffs

The cascade works as follows. FuXi-Short is applied first, rolling out in 6-hour steps from the initial state to produce forecasts at 6, 12, 18, , 120 h (20 rollout steps). At the 120 h mark, the predicted state ${\widehat{𝐱}}_{120\mathrm{h}}$ is passed to FuXi-Medium as its initial condition. FuXi-Medium then rolls out from 120 h to 240 h (another 20 steps). Finally, FuXi-Long takes the ${\widehat{𝐱}}_{240\mathrm{h}}$ state and extends to 360 h (15 days).

The rationale for this decomposition is rooted in the predictability horizon of the atmosphere. Lorenz (1969) showed that the predictability of synoptic-scale weather is limited to roughly two weeks by the intrinsic chaotic divergence of nearby trajectories. Different timescales are dominated by different dynamical processes: the first five days are primarily governed by baroclinic instability and the detailed synoptic-scale flow; days 5–10 are increasingly influenced by the large-scale circulation patterns and teleconnections; beyond 10 days, predictability depends on slowly varying boundary conditions (sea surface temperatures, soil moisture) and the Madden-Julian Oscillation. A single model trained to excel across all three regimes must compromise; the cascade avoids this by specialising.

Training and Achievements

Each FuXi sub-model is trained independently on ERA5 data from 1979 to 2015, with 2016–2017 for validation and 2018 for testing. Training uses 6-hourly data with autoregressive rollout of $K$ steps: FuXi-Short with $K=20$ (120 h); FuXi-Medium initialised from FuXi-Short predictions and trained for $K=20$ additional steps; FuXi-Long similarly.

The published results demonstrate skillful forecasts out to 15 days, measured by the latitude-weighted RMSE and anomaly correlation coefficient (ACC) for Z500 and T850. FuXi achieves positive ACC for Z500 at 15-day lead time in the Northern Hemisphere extra-tropics-a threshold traditionally associated with the “predictability barrier” of NWP. In head-to-head comparisons on the WeatherBench 2 benchmark (Rasp et al., 2024), FuXi-Medium outperforms GraphCast and Pangu-Weather at lead times beyond 7 days, while FuXi-Short is competitive with GraphCast at 1–5 day lead times.

FuXi-ENS: Ensemble Extension

Building on the deterministic FuXi system, the Fudan team subsequently developed FuXi-ENS: an ensemble extension that generates probabilistic forecasts by applying stochastic perturbations to the latent representations of the U-Transformer at each rollout step. Rather than training separate ensemble members from scratch, FuXi-ENS learns a perturbation network that produces spatially coherent additive noise fields in the latent space, conditioned on the current atmospheric state. This approach avoids the combinatorial explosion of training many ensemble members independently while producing spread-skill relationships that improve on deterministic single forecasts for verification of extreme events.

FengWu: Replay Buffer and Multi-Task Loss

FengWu (Chen et al., 2023, Shanghai AI Lab) is named after the ancient Chinese wind god. Its architecture is a multi-modal encoder-decoder that processes atmospheric variables in groups based on their physical characteristics, feeding separately encoded representations into a shared Transformer processor.

Multi-Modal Encoder-Decoder

The key design insight of FengWu is that the 70+ atmospheric variables in ERA5 span very different physical scales and statistical distributions. Rather than treating all variables identically (as GraphCast and Pangu-Weather do, modulo normalisation), FengWu divides variables into modality groups:

1. Upper-air dynamics: geopotential, temperature, wind components ($u$, $v$), specific humidity at 13 pressure levels (5 variables $\times$ 13 levels = 65 channels);

2. Surface fields: 2-metre temperature, 10-metre wind, mean sea-level pressure, total precipitation (4 channels);

3. Static fields: orography, land-sea mask, soil type (time-invariant, encoded once and broadcast).

Each modality group is processed by a dedicated encoder: a convolutional stem followed by cross-attention to a learned set of “modality tokens.” The modality tokens serve as a compressed representation of each group; they are then concatenated and processed by a shared PerceiverIO-style Transformer. The decoder inverts this process, mapping shared tokens back to per-modality spatial fields.

This architecture has two advantages. First, it allows the model to allocate representational capacity differentially: the upper-air group, which dominates forecast error at medium range, receives more parameters. Second, the modality tokens provide a bottleneck that enforces cross-variable communication: the upper-air and surface encoders must coordinate through the shared Transformer to produce a consistent forecast, rather than evolving independently.

Replay Buffer for Long Rollouts

A fundamental challenge in training MLWP models with long autoregressive rollouts is catastrophic forgetting. During training, as the model learns to minimise the loss at long lead times (e.g., 7–10 days), the gradient signal from those late steps may interfere with the representations learned for short lead times (1–3 days). This is precisely the catastrophic forgetting phenomenon from continual learning applied to the temporal axis of a single forecasting model.

FengWu addresses this with a replay buffer: a store of previously generated forecast trajectories (both from the model's own outputs and from ERA5) that is mixed into each training batch. At each gradient step, the model sees not only fresh ERA5 sequences but also “replayed” trajectories from earlier training phases, preserving the gradient signal that maintains short-range skill.

Formally, let $ℬ=\{({𝐱}_0^{(i)},{𝐱}_{k_i}^{(i)})\}$ be the replay buffer, where $k_i$ is the lead time (in 6-hour steps) of trajectory $i$. The training objective mixes a standard one-step MSE loss with a replay loss: (Replay LOSS)$$\begin{array}{ccc}& {ℒ}_{\mathrm{FengWu}}=\underset{\text{one-step loss}}{\underbrace{{𝔼}_{({𝐱}_t,{𝐱}_{t+1})\sim 𝒟}\left[{\ell }_{\mathrm{MSE}}(f_{\theta }({𝐱}_t),{𝐱}_{t+1})\right]}}+{\lambda }_r\underset{\text{replay loss}}{\underbrace{{𝔼}_{({𝐱}_0,{𝐱}_k)\sim ℬ}\left[{\ell }_{\mathrm{MSE}}(f_{\theta }^k({𝐱}_0),{𝐱}_k)\right]}},& \text{<span class="tml-eqn"></span>}\end{array}$$ where $f_{\theta }^k$ denotes the model rolled out $k$ steps, and ${\lambda }_r$ is a balancing coefficient (typically $0.1$–$0.3$).

Uncertainty-Based Multi-Task Loss Balancing

FengWu predicts 69 atmospheric variable-level fields (65 upper-air plus 4 surface). Naively summing the MSE across all variables in equal proportion gives disproportionate weight to variables with high variance (e.g., wind speed at 10 m) and effectively ignores low-variance variables that may be meteorologically significant (e.g., stratospheric temperatures).

FengWu adopts uncertainty-based task weighting (Kendall et al., 2018), treating each variable as a separate task in a multi-task learning framework. For task $c$, the loss contribution is (MTL)$$\begin{array}{ccc}& {ℒ}_c=\frac{1}{2s_c^2}{\Vert {\widehat{𝐱}}_c-{𝐱}_c\Vert }^2+\log s_c,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $s_c>0$ is a learnable uncertainty (standard-deviation) parameter for variable $c$. The total loss is $ℒ={\sum }_c{ℒ}_c$. Minimising over $s_c$ for a fixed prediction error is equivalent to maximising the Gaussian log-likelihood with variance $s_c^2$, and the regularisation term $\log s_c$ prevents $s_c\to \infty$ (which would collapse the loss to zero). Variables for which the model achieves low error naturally have small $s_c$ (low uncertainty), amplifying their gradient signal; variables for which error is irreducible (high atmospheric chaos) have large $s_c$, down-weighting their contribution.

Energy Efficiency

One of FengWu's headline claims is exceptional energy efficiency: a complete 10-day forecast (40 autoregressive steps at 6-hourly resolution) requires approximately $12\mathrm{kJ}$ of energy, compared to roughly $\mathrm{10,000}\mathrm{kJ}$ for an equivalent ECMWF HRES run. This three-orders-of-magnitude reduction in energy cost is a consequence of two factors: (i) the ML forward pass replaces millions of finite-difference operations on a petascale supercomputer with a single neural network evaluation on a modern GPU, and (ii) the smaller number of parameters in FengWu ($\sim$270 million) versus the ECMWF model (which implements a full dynamical core with many parameterisation schemes).

The energy efficiency advantage has profound implications for operational meteorology in resource-limited settings. A national weather service in a developing country that cannot afford supercomputer access can run FengWu on a single GPU server and produce medium-range forecasts of quality that would previously have required bilateral data-sharing agreements with European or North American centres.

Comparison: Pangu vs. GraphCast vs. FuXi vs. FengWu

Table Table 7 summarises the key properties of the four major MLWP systems studied in this chapter.

System	Architecture	Resolution	Horizon	10-day cost	Open
ECMWF HRES${}^{†}$	NWP (spectral)	9 km	10 d	$\sim$3.6 GJ	–
Pangu-Weather	3D + 2D ViT	0.25${}^{\circ }$	7 d	$\sim$12 kJ	Yes
GraphCast	GNN (icosahedral)	0.25${}^{\circ }$	10 d	$\sim$15 kJ	Yes
FuXi	U-Trans. cascade	0.25${}^{\circ }$	15 d	$\sim$15 kJ	Yes
FengWu	Multi-modal Enc-Dec	0.25${}^{\circ }$	10 d	$\sim$12 kJ	Yes

The table reveals several striking patterns. First, all MLWP systems are multiple orders of magnitude more energy-efficient than ECMWF HRES. Second, maximum horizon varies considerably: FuXi uniquely pushes to 15 days through its cascaded design. Third, resolution remains constant at ${0.25}^{\circ }$ across MLWP systems, not because finer resolution is undesirable, but because ERA5 training data at sub-${0.25}^{\circ }$ resolution is unavailable at the scale needed for training. Fourth, all major MLWP systems have been open-sourced, a remarkable departure from the proprietary culture of traditional NWP.

NeuralGCM: Differentiable Physics Meets Learned Parameterisations

NeuralGCM (Kochkov et al., 2024, Google Research) represents the most conceptually radical departure from the transformer-stacking approach. Rather than replacing the entire NWP pipeline with a neural network, NeuralGCM hybridises a differentiable atmospheric dynamical core with learned neural parameterisations. This hybrid design was motivated by a key observation: pure MLWP systems excel at weather forecasting (initial-value problems on timescales of days) but fail at climate simulation (boundary-value problems on timescales of years to centuries) because they do not conserve mass, energy, and momentum, and they cannot generalise to atmospheric states outside their training distribution.

Differentiable Atmospheric Dynamical Core

The foundation of NeuralGCM is a differentiable implementation of the primitive equations of atmospheric dynamics: (Momentum)$$\begin{array}{cccc}& \frac{D𝐮}{Dt}& =-\frac{1}{\rho }\nabla p-f\hat{k}\times 𝐮+{𝐅}_{\mathrm{phys}},& \text{<span class="tml-eqn"></span>}\\ & \frac{\partial T}{\partial t}& =-𝐮\cdot \nabla T+\frac{\kappa T\omega }{p}+Q_{\mathrm{phys}},& \text{<span class="tml-eqn"></span>}\\ & \frac{\partial \ln p_s}{\partial t}& =-{\int }_0^1\nabla \cdot 𝐮\mathrm{d}\eta +S_{\mathrm{phys}},& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝐮$ is the horizontal wind vector, $\rho$ is air density, $p$ is pressure, $f$ is the Coriolis parameter, $T$ is temperature, $\omega =Dp/Dt$ is vertical velocity in pressure coordinates, $\kappa =R_d/c_p$ is the ratio of the gas constant to the specific heat at constant pressure, $p_s$ is surface pressure, and $\eta$ is the terrain-following vertical coordinate.

The terms ${𝐅}_{\mathrm{phys}}$, $Q_{\mathrm{phys}}$, and $S_{\mathrm{phys}}$ are physical parameterisations: approximations of sub-grid-scale processes (convection, clouds, radiation, turbulence, land-surface exchange) that cannot be explicitly resolved at the model's grid resolution. In traditional GCMs, these parameterisations are elaborate hand-coded schemes with dozens of empirical parameters. NeuralGCM replaces them with learned neural networks.

The dynamical core is implemented in JAX, making it end-to-end differentiable with respect to all parameters (including the neural parameterisations), so the entire hybrid system can be trained jointly using automatic differentiation. This is a departure from the traditional “offline training” approach in which ML parameterisations are fitted to column-level inputs and outputs from the GCM without backpropagating through the dynamical core.

Learned Parameterisations

The neural parameterisations in NeuralGCM are column-local neural networks: small MLPs that take as input the vertical profile of atmospheric state at a grid column (temperature, humidity, wind, pressure) and produce tendencies (time derivatives) for the parameterised processes. Formally, for grid column $(i,j)$: (Param)$$\begin{array}{ccc}& {({𝐅}_{\mathrm{phys}},Q_{\mathrm{phys}},S_{\mathrm{phys}})}_{ij}=g_{\varphi }\left({\{T,q,u,v,p\}}_{ij,\ell =1}^L\right),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\varphi$ are the neural network parameters, $q$ is specific humidity, and $L$ is the number of vertical levels. The network $g_{\varphi }$ is shared across all horizontal grid columns-it is not location-specific-so it must learn universal parameterisations that apply equally in the tropics, mid-latitudes, and poles.

This column-locality is a deliberate inductive bias: physical parameterisations are largely local in the horizontal (a convective event in one column does not directly affect parameterised tendencies in a column 100 km away), so imposing this structure reduces the parameter count dramatically compared to a global attention mechanism.

Weather and Climate Capability

A crucial test of NeuralGCM's hybrid design is whether it can operate stably for climate-length integrations (years to decades). Pure MLWP models like GraphCast and Pangu-Weather accumulate systematic errors that cause the model state to drift off the atmospheric manifold after a few weeks; the predicted fields develop physically implausible structures (negative humidity, spurious temperature inversions) that make long-run stability impossible.

NeuralGCM avoids this pathology because its dynamical core explicitly conserves mass and angular momentum through the discretised primitive equations. The neural parameterisations can produce physically implausible tendencies locally, but the dynamical core constrains the global evolution to remain on the manifold of physically consistent atmospheric states. Kochkov et al. (2024) demonstrated NeuralGCM's ability to:

• Match or exceed GraphCast on medium-range weather forecasting (1–10 day skill);

• Simulate multiple years of stable atmospheric evolution without state-space divergence;

• Reproduce climatological statistics (mean temperature, precipitation patterns, variability indices) that drift catastrophically in pure MLWP models.

The ability to run in climate mode unlocks a new set of applications: century-scale climate projections at weather-resolving resolutions, rapid generation of large ensembles for uncertainty quantification, and “storyline” simulations that ask how specific weather events would change under different greenhouse gas concentrations.

Example: FuXi's 15-Day Forecast of an Atmospheric River

The Race for Weather AI Supremacy

Exercises

The Probabilistic Foundations of Weather Forecasting

The history of weather forecasting is, at its deepest level, the history of humanity coming to terms with a fundamental philosophical shift: the atmosphere does not have a single, deterministic future. It has a distribution of futures, and any honest forecast must communicate that distribution. This insight, which now seems obvious, was radical when Edward Epstein first articulated it in 1969. The intervening half-century has transformed it from a philosophical curiosity into the operational backbone of every major national weather service on the planet.

In this section we develop the mathematical foundations of probabilistic weather prediction from first principles. We trace the intellectual lineage from Epstein's stochastic-dynamic equations through Leith's theoretical analysis of Monte Carlo skill, to the operational ensemble systems introduced by Toth and Kalnay at the National Meteorological Center and by the European Centre for Medium-Range Weather Forecasts. We introduce the formal vocabulary - ensemble prediction systems, initial-condition uncertainty, model uncertainty, boundary uncertainty - and we derive the primary verification metric for probabilistic forecasts, the Continuous Ranked Probability Score. The section concludes with a spread-skill theorem that characterises when an ensemble is well-calibrated, and with a discussion of the immense computational cost that makes this framework ripe for machine-learning acceleration.

Epstein's Stochastic-Dynamic Prediction

In 1969, Edward Epstein published a short but transformative paper in the Journal of Applied Meteorology that asked a question no operational forecaster had dared to formalise: if we cannot know the initial state of the atmosphere exactly, what is the correct mathematical object to forecast? His answer was unambiguous. The correct object is not a single state trajectory, but the probability density function over all possible atmospheric trajectories consistent with our uncertain knowledge of the initial conditions.

Epstein observed that the atmosphere evolves according to a system of nonlinear partial differential equations - the primitive equations of fluid dynamics. Let $𝒙(t)\in {\mathrm{ℝ}}^n$ denote the state vector of the atmosphere at time $t$, where $n$ is the number of degrees of freedom (in practice, of order ${10}^8$ for a global model). The dynamics are governed by (Primitive)$$\begin{array}{ccc}& \frac{d𝒙}{dt}=𝒇(𝒙,t),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝒇:{\mathrm{ℝ}}^n\times \mathrm{ℝ}\to {\mathrm{ℝ}}^n$ encodes the physical laws of fluid dynamics, thermodynamics, and radiative transfer. In the deterministic paradigm, one specifies an initial condition $𝒙(0)={𝒙}_0$ and integrates forward to produce a unique forecast trajectory.

Epstein's insight was to replace the point initial condition with a probability distribution $p_0(𝒙)$ over possible initial states and to ask: how does this distribution evolve under the dynamics? If $p(𝒙,t)$ denotes the probability density at time $t$, then by the continuity equation (the Liouville equation), (Liouville)$$\begin{array}{ccc}& \frac{\partial p}{\partial t}+{\nabla }_{𝒙}\cdot (p𝒇)=0.& \text{<span class="tml-eqn"></span>}\end{array}$$ This is the stochastic-dynamic equation: a partial differential equation for the evolution of the full probability distribution over atmospheric states. A forecast is no longer a single trajectory but the entire time-evolving density $p(𝒙,t)$.

Leith's Theoretical Analysis of Monte Carlo Skill

Epstein's stochastic-dynamic formulation was theoretically clean but computationally intractable. Solving the Liouville equation in the full state space of a global weather model is impossible; the dimensionality is far too large. The natural alternative is Monte Carlo sampling: instead of propagating the distribution analytically, one propagates a finite ensemble of $N$ sample trajectories and uses the ensemble to approximate the distribution.

Cecil Leith (1974) provided the first rigorous theoretical analysis of how the skill of a Monte Carlo ensemble scales with ensemble size $N$ and with the accuracy of the underlying deterministic model. Leith's central result is that the mean-square error of an ensemble mean forecast decomposes into a variance term and a bias term, and that the variance term decreases as $1/N$.

Let $x_t$ be a scalar forecast variable (a simplification for exposition; the full theory applies componentwise to vector states). Suppose the true future state is $x_t^{\ast }$ and we have an ensemble ${\{x_t^{(k)}\}}_{k=1}^N$ of $N$ forecast members drawn from the forecast distribution $p(x_t|{𝒙}_0)$. Define the ensemble mean (Ensemble MEAN)$$\begin{array}{ccc}& {\bar{x}}_t=\frac{1}{N}\sum _{k=1}^N{x}_t^{(k)}.& \text{<span class="tml-eqn"></span>}\end{array}$$ Leith showed that under mild regularity conditions, (Leith Decomp)$$\begin{array}{ccc}& 𝔼[(x_t^{\ast }-{\bar{x}}_t{)}^2]=\underset{\text{irreducible error}}{\underbrace{𝔼[(x_t^{\ast }-{\mu }_t{)}^2]}}+\underset{\text{sampling variance}}{\underbrace{\frac{{\sigma }_t^2}{N}}},& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\mu }_t=𝔼[x_t]$ and ${\sigma }_t^2=\mathrm{𝖵ar}[x_t]$ are the mean and variance of the forecast distribution. The first term is the irreducible forecast error due to chaotic growth of uncertainties; it cannot be reduced by running more ensemble members. The second term is the finite-ensemble sampling error; it vanishes as $N\to \infty$.

Operational Ensemble Forecasting at NMC

The theoretical case for ensemble forecasting was well-established by the early 1990s, but operational implementation required both the computational resources and the methodological innovations to make it practical. Toth and Kalnay (1993) introduced the breeding method for generating initial perturbations and demonstrated the first operational ensemble forecast system at the National Meteorological Center (NMC, later renamed the National Centers for Environmental Prediction, NCEP).

The core challenge in designing an ensemble is the choice of initial perturbations. One could, in principle, draw perturbations at random from the analysis error covariance. But Toth and Kalnay observed that random perturbations project mostly onto stable or neutral directions in phase space - directions in which errors decay or remain bounded. What matters for forecast skill is the growing error directions: the modes in which small initial errors amplify most rapidly.

The breeding method constructs perturbations by iteratively running the model forward and then normalising the evolved perturbation back to the initial amplitude. Formally, let ${𝒙}_0$ be the analysis state. Choose a random initial perturbation $\delta {𝒙}_0$ with small amplitude $ϵ$. Integrate both the perturbed and unperturbed forecasts forward by a breeding cycle of length $\mathrm{\Delta }t$: (Breeding Cycle)$$\begin{array}{ccc}& \delta {𝒙}_{\mathrm{\Delta }t}={\mathrm{\Phi }}_{\mathrm{\Delta }t}({𝒙}_0+\delta {𝒙}_0)-{\mathrm{\Phi }}_{\mathrm{\Delta }t}({𝒙}_0),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\mathrm{\Phi }}_{\mathrm{\Delta }t}$ is the model's forward-integration operator. The rescaled perturbation (Rescale)$$\begin{array}{ccc}& \delta {𝒙}_0^{\text{new}}=ϵ\cdot \frac{\delta {𝒙}_{\mathrm{\Delta }t}}{\Vert \delta {𝒙}_{\mathrm{\Delta }t}\Vert }& \text{<span class="tml-eqn"></span>}\end{array}$$ is then used as the starting point for the next breeding cycle. After many cycles, the perturbation aligns with the fastest-growing error mode - the bred vector - and captures the directions in phase space where the ensemble spread is most meaningful.

The NCEP operational ensemble, as described by Toth and Kalnay, used five pairs of bred vectors (each pair consisting of a positive and a negative perturbation) to create ten initial states offset from the analysis. Together with the unperturbed control run, this produced an eleven-member ensemble at a fraction of the cost of running eleven independent forecasts.

Stochastic Weather Generators

Ensemble prediction systems target the spread of forecast trajectories at medium range (up to about two weeks). A complementary tradition, reviewed comprehensively by Wilks and Wilby (1999), addresses a different but related problem: generating synthetic sequences of weather variables that are statistically indistinguishable from observations over long time horizons. These stochastic weather generators (SWGs) are used for climate impact assessment, hydrological modelling, agricultural planning, and risk analysis.

The classic SWG for daily precipitation models wet-day occurrence as a first-order Markov chain: the probability of rain on day $t$ depends only on whether day $t-1$ was wet or dry. Let $W_t\in \{\mathrm{0,1}\}$ be a binary indicator of wet days. The Markov model is specified by two transition probabilities, (Markov Precip)$$\begin{array}{ccc}& p_{01}=P(W_t=1|W_{t-1}=0),p_{11}=P(W_t=1|W_{t-1}=1),& \text{<span class="tml-eqn"></span>}\end{array}$$ estimated from historical records. Given $W_t=1$, the precipitation amount $R_t>0$ is drawn from a fitted distribution (often gamma or mixed exponential).

The stochastic weather generator tradition is important for our purposes because it represents an early, if simple, form of generative modelling for weather: the goal is to learn a distribution $p({𝒛}_{1:T})$ over weather sequences and sample from it, rather than to predict the next state deterministically. The limitations of classical SWGs - their reliance on prescribed parametric distributions, their inability to capture complex spatial structures, their difficulty with rare extremes - are precisely the limitations that modern deep generative models are designed to overcome.

The Conceptual Framework for Ensemble Prediction

The most influential modern synthesis of why ensembles matter is due to Leutbecher and Palmer (2008), who provide a comprehensive conceptual and practical guide to ensemble prediction. Their framework rests on three interacting sources of forecast uncertainty that together determine the spread of an ensemble.

Leutbecher and Palmer emphasise that a well-designed ensemble must represent all three sources simultaneously. An ensemble that perturbs only the initial conditions will be under-dispersive at long lead times, because model errors grow and accumulate in a way that the ensemble cannot capture. Conversely, an ensemble that represents model uncertainty but ignores initial condition uncertainty will produce misleading probability estimates in the first few days when initial condition uncertainty dominates.

Singular Vector Perturbations at ECMWF

The European Centre for Medium-Range Weather Forecasts (ECMWF) developed singular vector (SV) perturbations as the primary mechanism for capturing initial condition uncertainty in its ensemble prediction system (EPS). The mathematical foundation of SVs rests on the theory of optimal linear perturbation growth.

Let ${\mathrm{\Phi }}_{\tau }:{\mathrm{ℝ}}^n\to {\mathrm{ℝ}}^n$ denote the nonlinear forward model integrated over time interval $\tau$, and let ${𝐋}_{\tau }=D_{𝒙}{\mathrm{\Phi }}_{\tau }({𝒙}_0)$ be its linearisation (the tangent linear model) about the analysis trajectory. Define an energy inner product ${⟨𝒖,𝒗⟩}_E={𝒖}^{\top }𝐄𝒗$ on the state space, where $𝐄$ is a positive definite energy-metric matrix. The singular vectors are the solutions to the constrained optimisation problem (SV Optim)$$\begin{array}{ccc}& {𝒖}_j={\mathrm{arg\hspace{0.17em}max}}_{𝒖:\Vert 𝒖{\Vert }_E=1}\Vert {𝐋}_{\tau }𝒖{\Vert }_E,j=1,\dots ,M,& \text{<span class="tml-eqn"></span>}\end{array}$$ i.e., the initial perturbations of unit energy that experience the greatest energy growth over the interval $\tau$. These are the right singular vectors of ${𝐄}^{1/2}{𝐋}_{\tau }{𝐄}^{-1/2}$ ordered by decreasing singular value ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0$. The $j$-th singular value ${\sigma }_j$ quantifies the amplification factor: an initial perturbation aligned with ${𝒖}_j$ will grow by a factor of ${\sigma }_j$ in energy over time $\tau$.

ECMWF's operational implementation uses $M=50$ singular vectors per hemisphere, computed for a target time of $\tau =48$ hours. The leading SVs are those that grow most rapidly, and they correspond to the meteorologically active regions where forecast uncertainty is expected to be largest (typically jet stream regions, baroclinic zones, and cyclone development areas).

Stochastic Physics Parameterisations

Initial-condition perturbations address the first source of uncertainty in Leutbecher and Palmer's framework, but they cannot capture model errors that accumulate during integration. To represent model uncertainty, modern ensemble systems augment their physics parameterisations with stochastic components.

The most widely used approach is the Stochastically Perturbed Parameterisation Tendencies (SPPT) scheme, introduced at ECMWF in the early 2000s. In SPPT, the parameterised tendency ${\dot{𝒙}}^{\text{param}}$ from all physics schemes (convection, boundary layer, radiation, etc.) is multiplied by a spatially and temporally correlated random field $r(𝒙,t)$: (SPPT)$$\begin{array}{ccc}& \frac{d𝒙}{dt}={\dot{𝒙}}^{\text{dyn}}+r\cdot {\dot{𝒙}}^{\text{param}},r(𝒙,t)\sim 𝒢({\rho }_s,{\rho }_t,{\sigma }_r),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\dot{𝒙}}^{\text{dyn}}$ is the dynamical (resolved-scale) tendency and $𝒢({\rho }_s,{\rho }_t,{\sigma }_r)$ denotes a Gaussian random field with spatial correlation length ${\rho }_s$, temporal decorrelation time ${\rho }_t$, and standard deviation ${\sigma }_r$. Typical operational values are ${\rho }_s\approx 500$,km, ${\rho }_t\approx 6$ hours, and ${\sigma }_r\approx 0.5$.

The Continuous Ranked Probability Score

Evaluating a probabilistic forecast requires a verification metric that rewards both accuracy and reliability - a metric that assigns a high score to a forecast distribution that places most of its mass near the verifying observation, and a low score to a distribution that is either biased or over-/under-confident. The standard metric for this purpose is the Continuous Ranked Probability Score (CRPS).

Several properties of the CRPS make it the preferred metric for ensemble verification. It is a proper scoring rule: the expected CRPS is minimised if and only if the predictive distribution $F$ equals the true data-generating distribution. This means a forecaster cannot artificially improve the expected CRPS score by hedging or by reporting a distribution that differs from their true belief.

Decomposition of CRPS

The Spread-Skill Relationship

A well-calibrated ensemble should have a spread that equals its skill: if the ensemble members are drawn from the forecast distribution, the root-mean-square error (RMSE) of the ensemble mean should equal the ensemble standard deviation (spread), on average.

In practice, most operational ensemble systems are moderately under-dispersive, particularly at longer lead times, because model errors introduce systematic biases that the ensemble perturbations do not fully capture. This is one reason why post-processing techniques - statistical methods that inflate ensemble spread to match observed error statistics - are routinely applied before ensemble forecasts are issued to users.

The Computational Cost of Ensemble NWP

The theoretical and practical case for ensemble prediction is compelling, but the cost is formidable. A state-of-the-art global ensemble such as ECMWF's ENS runs 50 ensemble members plus one unperturbed control forecast at the same resolution as the members (the separate high-resolution deterministic run is HRES). Each member is itself a NWP model integration at approximately 9,km (since IFS Cycle 48r1, June 2023; 18,km before it) horizontal resolution (as of 2024), with roughly ${10}^9$ grid-point degrees of freedom and $\sim {10}^3$ vertical levels after accounting for the spectral truncation.

The high cost of ensemble NWP motivates two complementary machine-learning strategies. First, one can train a deterministic ML model to approximate a single NWP integration, reducing the cost of each member while maintaining ensemble diversity through perturbations - the approach taken by GraphCast, Pangu-Weather, and others. Second, one can train a generative model that directly samples from the forecast distribution without running $N$ separate integrations - the approach of GenCast, AIFS-CRPS, and related systems. The next section examines why the second approach requires generative modelling and cannot be achieved by deterministic MLmodels alone.

Exercises

Why Generative Models? The Blurring Problem

The previous section established that weather forecasting is inherently probabilistic: the atmosphere is a chaotic system, and an honest forecast must be a distribution, not a point. Machine learning has produced a generation of deterministic weather models - GraphCast, Pangu-Weather, FourCastNet, and their successors - that achieve remarkable skill at short to medium lead times, matching or exceeding operational NWP at a fraction of the compute cost. But there is a deep mathematical reason why these deterministic models cannot replace ensemble forecasting, no matter how large they are trained or how much data they see. This section develops that argument rigorously.

The argument proceeds in four steps. First, we show that any model trained to minimise mean squared error converges to the conditional expectation, which is a smooth average over all possible futures. Second, we show that for chaotic systems, this conditional expectation loses the sharp spatial features of any individual realisation. Third, we document empirically that GraphCast and Pangu-Weather exhibit exactly this behaviour at long lead times. Fourth, we argue that the only principled remedy is to train a generative model that explicitly targets the conditional distribution rather than its mean.

The MSE Minimisation Trap

This proposition is a classical result in estimation theory, but its implications for weather forecasting are severe. A neural network trained with an MSE loss on a dataset of input–output pairs $({𝒙}_t,{𝒙}_{t+1})$ will converge, given sufficient capacity, to the conditional mean $𝔼[{𝒙}_{t+1}|{𝒙}_t]$. This is not a sample from the conditional distribution; it is an average over all possible next states.

Why the Conditional Mean Is Blurry

The conclusion of Proposition Proposition 4 is stark: a deterministic ML model trained to minimise MSE will produce forecasts with systematically damped spatial variability compared to the truth. At short lead times, when the conditional distribution is tightly concentrated around a single trajectory, the conditional mean closely approximates a typical realisation. But as lead time increases and the forecast distribution spreads over multiple distinct weather patterns, the conditional mean becomes an average of incompatible structures - a field that looks like no particular weather realisation, with blurred fronts, smoothed cyclones, and damped precipitation maxima.

Empirical Evidence: GraphCast and Pangu-Weather

Proposition Proposition 4 is not merely a theoretical concern; it is observed in practice. Both GraphCast and Pangu-Weather - the two most widely benchmarked deterministic MLWP models - exhibit increasing spatial blurring at lead times beyond approximately 5–7 days.

The spectral signature of this blurring was systematically documented by Kim et al. (2026), who computed the kinetic energy spectrum of GraphCast and GenCast forecasts as a function of lead time and compared them with the ERA5 reanalysis power spectrum. At lead times of 1–3 days, both deterministic and generative models track the ERA5 spectrum closely across all resolved wavenumbers. At lead times of 7–10 days, GraphCast shows a pronounced spectral depression at synoptic and mesoscales (wavenumbers 10–100): power is systematically deficient by 10–30% relative to ERA5. GenCast, trained as a generative diffusion model, maintains spectral fidelity across all lead times.

Brenowitz et al. (2025) document a related paradox in MLWP benchmarking: when deterministic MLWP models such as GraphCast and Pangu are evaluated using CRPS in a lagged-ensemble configuration - using multiple initialisation times to construct an ensemble - their CRPS scores are substantially better than their deterministic metrics suggest. In other words, these models produce diverse forecasts across different runs (capturing climatological variability), but any individual run is smooth (failing to capture the spatial sharpness of a single realisation). This creates the paradox that a model's ensemble CRPS can be competitive while its individual forecast is visually blurred.

The Physical Consequences of Blurring

The blurring of deterministic MLWP forecasts has two particularly damaging physical consequences.

Underprediction of extreme events. Extreme weather events - heavy precipitation, strong winds, deep cyclones - are precisely the features that users care most about. They correspond to the tails of the conditional distribution, regions far from the conditional mean. An MSE-trained model produces forecasts near the mean and therefore systematically underestimates the intensity of extremes. Operationally, this means that a deterministic ML forecast of a category-4 hurricane at $\tau =10$ days may appear as a broad, weak low-pressure system with no clear intensification signal, because the conditional mean averages over trajectories that do and do not intensify.

Damped energy spectrum. The atmosphere has a power-law kinetic energy spectrum $E(k)\sim k^{-3}$ in the synoptic range (wavelengths above roughly 500 to 1000,km), transitioning to $E(k)\sim k^{-5/3}$ at mesoscales (below roughly 500,km). A blurred forecast lacks energy at small scales, violating the physical scaling law. This has practical consequences for downstream applications: wind energy estimation, convection-permitting downscaling, and any application that depends on mesoscale structure will be degraded by blurred input fields.

The Generative Weather Model

The preceding analysis makes clear what a weather model must do to avoid the blurring problem: it must learn the full conditional distribution $p({𝒙}_{t+1}|{𝒙}_t)$, not merely its mean.

The phrase “learns $p({𝒙}_{t+1}|{𝒙}_t)$ not $𝔼[{𝒙}_{t+1}|{𝒙}_t]$” is the essential distinction. A generative model, when sampled, produces sharp, physically realistic fields that look like actual weather. When sampled many times, the collection of samples forms an ensemble that correctly quantifies forecast uncertainty.

Physical Consistency: Spectral Tests

Kim et al. (2026) provide a comprehensive evaluation of the physical consistency of GenCast, the diffusion-based generative ensemble model of Price et al. (2025). Their analysis introduces a suite of spectral diagnostic tests that go beyond conventional RMSE and CRPS verification.

The central test compares the time-averaged kinetic energy spectrum $E_{\text{model}}(k,\tau )$ of model forecasts to the corresponding spectrum of ERA5 reanalysis fields $E_{\text{ERA5}}(k)$, as a function of lead time $\tau$ and wavenumber $k$. Formally, they define the spectral ratio (Spectral Ratio)$$\begin{array}{ccc}& R(k,\tau )=\frac{E_{\text{model}}(k,\tau )}{E_{\text{ERA5}}(k)},& \text{<span class="tml-eqn"></span>}\end{array}$$ with $R(k,\tau )=1$ indicating perfect spectral consistency. Values $R<1$ indicate spectral damping (blurring); values $R>1$ indicate spectral over-amplification (noise).

Their key finding is that GenCast maintains $R(k,\tau )\approx 1$ across all wavenumbers $k$ and lead times $\tau \le 15$ days, whereas GraphCast shows $R(k,\tau )<0.85$ for $k>15$ and $\tau >7$ days. The spectral deficiency grows approximately linearly with lead time, consistent with the blurring mechanism described in Proposition Proposition 4.

Kim et al. also test for physical consistency in the precipitation field, using the distribution of daily accumulations and their spatial correlation structure. Generative models correctly reproduce the heavy-tailed marginal distribution of precipitation, while deterministic models produce over-smooth precipitation fields with too little probability in the tail beyond the 95th percentile.

The AI Weather Benchmarks Paradox

Brenowitz et al. (2025) identify what they call the AI weather benchmarks paradox: the metrics most commonly used to evaluate deterministic MLWP models (RMSE of geopotential height at 500,hPa) may have led the field to optimise for the wrong objective.

Their argument proceeds as follows. When a deterministic model such as GraphCast is evaluated at $\tau =10$ days on RMSE, it achieves a score comparable to the ECMWF HRES operational deterministic forecast. This has been widely reported as evidence that ML models have “matched” NWP. But when the same models are evaluated on CRPS, comparing each deterministic forecast against the ECMWF ENS ensemble, a striking pattern emerges: GraphCast's CRPS at $\tau =10$ days is actually better than its RMSE would suggest, because the model is implicitly diverse across initialisation times.

However, the paradox cuts both ways. The CRPS improvement of deterministic models is largely due to their climatological diversity (they are sensitive to the initial condition), not to their ability to represent the forecast distribution at a fixed initial time. A model that is evaluated on CRPS through lagged ensembles is rewarded for correlating strongly with observed variability across time, not for producing calibrated distributional forecasts at a fixed lead time. Brenowitz et al. show that when models are evaluated on strictly synchronous CRPS (fixed initialisation time, $N$ samples from the model vs $N$ ensemble members from ECMWF ENS), deterministic models perform much worse than their lagged-ensemble evaluation suggests.

Visualising the Blurring Problem

The Two-Way Relationship

The argument so far might suggest that generative models entirely supersede deterministic ones. This is too strong. In practice, the two paradigms are complementary and each benefits from the other.

Generative models condition on recent analysis states. The leading operational generative weather models - GenCast, AIFS-CRPS, and related systems - condition the generator on the most recent analysis states rather than perturbing a separate deterministic model's output. GenCast, for instance, is a pure conditional diffusion model that takes the two most recent analysis states as input and denoises directly into the next state. The conditioning states provide the large-scale flow structure (which is predictable at medium range); the generative component adds the small-scale variability that accounts for forecast uncertainty. Without the deterministic backbone, the generative model would need to learn all of the large-scale dynamics from scratch, which is far more difficult.

Deterministic models need generative corrections. Conversely, a deterministic model run at multiple initialisations - with small perturbations applied at each run - can be used to approximate an ensemble. The quality of this lagged ensemble depends on how sensitive the model is to initial perturbations, which is itself a function of how well the model has learned the chaotic dynamics. A generative model that captures the conditional distribution provides, in principle, a better characterisation of what perturbations are physically plausible, which can be used to design better initialisation perturbations for deterministic models.

Training Directly on CRPS: AIFS-CRPS

A particularly direct approach to the blurring problem is to replace the MSE training objective with the CRPS itself. If the CRPS is a proper scoring rule that rewards distributional accuracy, why not train a model to minimise it directly?

Lang et al. (2026) describe AIFS-CRPS, the ECMWF AI Integrated Forecasting System trained with a CRPS loss. The architecture is a graph transformer operating on a multi-scale icosahedral mesh, following the AIFS architecture of Ben Bouallègue et al. (2024). The key modification is the training objective: instead of minimising the MSE between predicted and verifying fields, the model minimises the energy score loss, which is equivalent to the multivariate CRPS: (CRPS LOSS)$$\begin{array}{ccc}& {ℒ}_{\mathrm{CRPS}}(\theta )=𝔼[\frac{1}{M}\sum _{k=1}^M\Vert {𝒙}_{t+1}^{(k)}-{𝒙}_{t+1}^{\ast }\Vert -\frac{1}{2M^2}\sum _{k,j}\Vert {𝒙}_{t+1}^{(k)}-{𝒙}_{t+1}^{(j)}\Vert ],& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\{{𝒙}_{t+1}^{(k)}\}}_{k=1}^M$ are $M$ samples drawn from the model and ${𝒙}_{t+1}^{\ast }$ is the verifying state. Minimising ((CRPS LOSS)) with respect to $\theta$ encourages the model to produce a sample distribution that matches the verifying distribution in expectation.

Lang et al. find that AIFS-CRPS achieves competitive CRPS scores compared to the 50-member ECMWF ENS ensemble, at a fraction of the computational cost (each ensemble member is one AIFS-CRPS sample, which is several orders of magnitude cheaper than an NWP member). Crucially, AIFS-CRPS samples do not exhibit the spectral damping seen in MSE-trained models: the energy spectrum tracks ERA5 closely at lead times up to 10 days.

Comparison: Deterministic vs. Generative MLWP

Dimension	Deterministic MLWP	Generative MLWP
Representative models	GraphCast, Pangu, FourCastNet	GenCast, AIFS-CRPS, NeuralGCM
[2pt] Training objective	MSE / MAE loss	CRPS, Energy Score, ELBO
[2pt] Output type	Single deterministic field	Distribution; sampleable
[2pt] Forecast uncertainty	None (point estimate)	Full distribution
[2pt] Ensemble generation	Requires lagged runs	Native: draw $N$ samples
[2pt] Spatial sharpness	Degrades at $\tau \gtrsim 5$,d	Maintained across lead times
[2pt] Energy spectrum	Damped at short scales	Consistent with ERA5
[2pt] Extreme events	Systematically underpredicted	Correctly represented in tails
[2pt] Compute per forecast	1 model forward pass	$N$ forward passes ($N\sim 50$)
[2pt] Calibration (CRPS)	Poor (not designed for this)	Good (trained for this)
[2pt] Physical consistency	Weakens with lead time	Maintained if trained on CRPS
[2pt] Interpretability	Single trajectory, easy to inspect	Requires ensemble diagnostics
[2pt] Anchor needed?	Self-contained	Often uses deterministic anchor
[2pt] Key limitation	Blurs at long lead times	More expensive; harder to train

Exercises

GenCast and Diffusion Ensemble Forecasting

Deterministic neural weather models such as GraphCast and Pangu-Weather demonstrated that deep learning could surpass the European Centre for Medium-Range Weather Forecasts (ECMWF) deterministic High-Resolution Forecast (HRES) on a wide array of verification metrics. Yet deterministic forecasts are fundamentally incomplete: the atmosphere is a chaotic system, and a single trajectory through state space cannot convey the range of possible futures that operational forecasters need to communicate risk. The state of the art in operational NWP therefore relies on ensemble forecasting: running many perturbed versions of a model to sample the forecast probability distribution. ECMWF's Ensemble Prediction System (ENS) runs 50 members at reduced resolution (0.2 degrees) with carefully constructed initial condition perturbations, at a cost of several hours of supercomputer time.

GenCast (Price et al., 2024/2025, Nature) is a conditional diffusion model trained on the ERA5 reanalysis archive that converts the paradigm of ensemble generation from running many perturbed deterministic models to drawing many independent samples from a learned probabilistic generative model. Each sample is produced by running a reverse diffusion process conditioned on the two most recent atmospheric states, requiring 8 minutes on a single tensor processing unit (TPU)- far less than the compute required to run a single ENS member. In a head-to-head comparison across 1,320 verification targets spanning lead times of 12 to 360 hours, GenCast outperformed ECMWF ENS on 97.4% of those targets, establishing a new state of the art in probabilistic medium-range weather forecasting.

Diffusion Models: A Probabilistic Primer

Before developing GenCast's specific architecture, we establish the mathematical framework of denoising diffusion probabilistic models (DDPMs) in the setting relevant to conditional weather forecasting.

The data distribution.

Let ${𝒙}_0\in {\mathrm{ℝ}}^d$ denote a sample from the target data distribution $q({𝒙}_0)$. In the weather forecasting context, ${𝒙}_0$ represents the atmospheric state at the next forecast step: a high-dimensional vector encoding temperature, wind, humidity, and surface variables on a global grid at 0.25-degree resolution across 13 pressure levels, with $d$ exceeding ${10}^6$ entries.

The forward process.

The DDPM forward process is a fixed Markov chain that gradually corrupts ${𝒙}_0$ with Gaussian noise according to a noise schedule ${\{{\beta }_t\}}_{t=1}^T$, with $0<{\beta }_1<{\beta }_2<\cdots <{\beta }_T<1$: (Forward STEP)$$\begin{array}{ccc}& q({𝒙}_t|{𝒙}_{t-1})=𝒩({𝒙}_t;\sqrt{1-{\beta }_t}{𝒙}_{t-1},{\beta }_t𝐈),t=1,\dots ,T.& \text{<span class="tml-eqn"></span>}\end{array}$$ The schedule $\{{\beta }_t\}$ is chosen so that at $t=T$ the distribution $q({𝒙}_T|{𝒙}_0)\approx 𝒩(\mathrm{𝟎},𝐈)$ is essentially pure Gaussian noise, destroying all information about ${𝒙}_0$.

The reverse process.

The key insight of DDPM is that the reverse conditional distribution $q({𝒙}_{t-1}|{𝒙}_t,{𝒙}_0)$ is also Gaussian, with computable mean and variance: (Reverse Posterior)$$\begin{array}{ccc}& q({𝒙}_{t-1}|{𝒙}_t,{𝒙}_0)=𝒩({𝒙}_{t-1};{\widetilde{𝝁}}_t({𝒙}_t,{𝒙}_0),{\tilde{\beta }}_t𝐈),& \text{<span class="tml-eqn"></span>}\end{array}$$ where (Reverse MEAN)$$\begin{array}{cccc}& {\widetilde{𝝁}}_t({𝒙}_t,{𝒙}_0)& =\frac{\sqrt{{\bar{\alpha }}_{t-1}}{\beta }_t}{1-{\bar{\alpha }}_t}{𝒙}_0+\frac{\sqrt{1-{\beta }_t}(1-{\bar{\alpha }}_{t-1})}{1-{\bar{\alpha }}_t}{𝒙}_t,& \text{<span class="tml-eqn"></span>}\\ & {\tilde{\beta }}_t& =\frac{1-{\bar{\alpha }}_{t-1}}{1-{\bar{\alpha }}_t}{\beta }_t.& \text{<span class="tml-eqn"></span>}\end{array}$$ Because ${𝒙}_0$ is unknown at inference time, we train a neural network to predict it (or equivalently, to predict the noise $𝜺$ added to ${𝒙}_0$ to obtain ${𝒙}_t$).

The score network.

The denoising network ${𝜺}_{𝜽}({𝒙}_t,t,𝒄)$ takes the noisy sample ${𝒙}_t$, the diffusion step index $t$, and a conditioning context $𝒄$ (the recent atmospheric states in GenCast's case), and predicts the noise $𝜺$. The approximate reverse process is then: (Learned Reverse)$$\begin{array}{ccc}& p_{𝜽}({𝒙}_{t-1}|{𝒙}_t,𝒄)=𝒩({𝒙}_{t-1};{𝝁}_{𝜽}({𝒙}_t,t,𝒄),{\tilde{\beta }}_t𝐈),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝝁}_{𝜽}$ is obtained by substituting the network prediction ${\widehat{𝒙}}_0=({𝒙}_t-\sqrt{1-{\bar{\alpha }}_t}{𝜺}_{𝜽}({𝒙}_t,t,𝒄))/\sqrt{{\bar{\alpha }}_t}$ in place of ${𝒙}_0$ into equation .

Conditional Diffusion for Weather: Formal Definition

GenCast models the conditional distribution of the next atmospheric state ${𝒙}_{t+1}$ given the current and previous states $𝒄=({𝒙}_t,{𝒙}_{t-1})$.

Forward SDE formulation.

GenCast adopts the continuous-time SDE framework. The forward process corresponds to an Itô SDE: (Forward SDE)$$\begin{array}{ccc}& d𝑿=-\frac{1}{2}\beta (t)𝑿dt+\sqrt{\beta (t)}d{𝑾}_t,t\in [0,T],& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝑾}_t$ is a standard Wiener process and $\beta (t)$ is a continuous noise schedule. The marginal at time $t$ satisfies $q_t(𝒙|{𝒙}_0)=𝒩(𝒙;\alpha (t){𝒙}_0,{\sigma }^2(t)𝐈)$ with $\alpha (t)=e^{-\frac{1}{2}{\int }_0^t\beta (s)ds}$ and ${\sigma }^2(t)=1-e^{-{\int }_0^t\beta (s)ds}$.

Reverse SDE.

By Anderson's (1982) time-reversal theorem, the reverse-time SDE is: (Reverse SDE)$$\begin{array}{ccc}& d𝑿=[-\frac{1}{2}\beta (t)𝑿-\beta (t){\nabla }_{𝑿}\log q_t(𝑿|𝒄)]dt+\sqrt{\beta (t)}d{\overline{𝑾}}_t,& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\overline{𝑾}}_t$ is a reverse-time Wiener process and ${\nabla }_{𝑿}\log q_t(𝑿|𝒄)$ is the conditional score function. GenCast approximates this score with the neural network $s_{𝜽}(𝑿,t,𝒄)\approx {\nabla }_{𝑿}\log q_t(𝑿|𝒄)$.

Training Objective

The GenCast training objective is the noise-prediction loss (epsilon-matching), which minimises: (LOSS)$$\begin{array}{ccc}& ℒ(𝜽)={𝔼}_{t\sim 𝒰[1,T],{𝒙}_0\sim q(\cdot |𝒄),𝜺\sim 𝒩(\mathrm{𝟎},𝐈)}[w(t)\Vert 𝜺-{𝜺}_{𝜽}(\sqrt{{\bar{\alpha }}_t}{𝒙}_0+\sqrt{1-{\bar{\alpha }}_t}𝜺,t,𝒄\left){\Vert }^2\right],& \text{<span class="tml-eqn"></span>}\end{array}$$ where $w(t)$ is a loss weighting function that may up-weight certain noise levels to improve sample quality. This objective is a variational lower bound (ELBO) on the negative log-likelihood of the data under the model (Ho et al., 2020), weighted by the SNR schedule.

Variable weighting in weather.

In the atmospheric context, a naive ${\ell }_2$ loss in state space gives undue weight to variables with large physical magnitude (e.g., surface pressure in Pascals vs. specific humidity in kg/kg). GenCast normalises each variable and pressure level by its climatological standard deviation, computed from the ERA5 training set. Additionally, grid-cell area weighting ensures that equatorial cells (with larger area on the sphere) receive proportional weight: (Weighted LOSS)$$\begin{array}{ccc}& {ℒ}_{\mathrm{GenCast}}(𝜽)={𝔼}_{t,{𝒙}_0,𝜺}\sum _{v,p,\ell }\frac{A_{\ell }}{{\sigma }_{v,p}^2}({\epsilon }_{v,p,\ell }-[{𝜺}_{𝜽}(\cdot ){]}_{v,p,\ell }{)}^2,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $v$ indexes variables, $p$ indexes pressure levels, $\ell$ indexes grid cells, $A_{\ell }$ is the grid-cell area, and ${\sigma }_{v,p}^2$ is the climatological variance of variable $v$ at level $p$.

Architecture: Graph-Based Diffusion

GenCast inherits the encode–process–decode graph neural network architecture of GraphCast, augmenting it with the noise-conditioning and time-step embedding required for diffusion.

Spatial discretisation.

The global atmosphere at 0.25-degree resolution is represented on an icosahedral–octahedral (Ico-O) reduced Gaussian grid with approximately 1.1 million grid points. GenCast covers the full globe and models 80+ atmospheric variables across 13 pressure levels (50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000 hPa) plus several surface and single-level variables.

Encoder.

The encoder maps the conditioning context $𝒄=({𝒙}_t,{𝒙}_{t-1})$ and the noisy target ${𝒙}_s^{(\tau )}$ (at diffusion step $\tau$) from the physical grid to a multi-mesh graph latent space. The physical grid nodes are connected to a coarser multi-resolution icosahedral mesh via learned bipartite message-passing layers: (Encoder)$$\begin{array}{ccc}& {𝒉}_i^{\mathrm{lat}}=f_{\mathrm{enc}}(\{{𝒄}_j,{𝒙}_{s,j}^{(\tau )}{\}}_{j\in 𝒩(i)};{𝜽}_{\mathrm{enc}}),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝒩(i)$ denotes the physical grid neighbours of mesh node $i$.

Diffusion step embedding.

The diffusion time step $\tau$ is embedded via a sinusoidal positional encoding followed by a two-layer MLP, producing a vector ${𝒆}_{\tau }\in {\mathrm{ℝ}}^{d_{\mathrm{emb}}}$. This embedding is injected into every graph network layer through feature-wise linear modulation (FiLM): (FILM)$$\begin{array}{ccc}& {\widetilde{𝒉}}_i=(\mathrm{𝟏}+\gamma ({𝒆}_{\tau }))\odot {𝒉}_i+\delta ({𝒆}_{\tau }),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\gamma ,\delta :{\mathrm{ℝ}}^{d_{\mathrm{emb}}}\to {\mathrm{ℝ}}^{d_{\mathrm{lat}}}$ are learned linear projections and $\odot$ denotes element-wise multiplication.

Processor.

The processor applies $L$ rounds of message-passing on the multi-mesh graph, enabling global information exchange: (MSG)$$\begin{array}{cccc}& {𝒎}_{ij}^{(l)}& ={\varphi }_e^{(l)}({𝒉}_i^{(l)},{𝒉}_j^{(l)},e_{ij};{𝜽}_e^{(l)}),& \text{<span class="tml-eqn"></span>}\\ & {𝒉}_i^{(l+1)}& ={\varphi }_v^{(l)}({𝒉}_i^{(l)},{\sum }_{j\in 𝒩(i)}{𝒎}_{ij}^{(l)};{𝜽}_v^{(l)}),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $e_{ij}$ are edge features encoding relative position and distance in the mesh, and ${\varphi }_e^{(l)},{\varphi }_v^{(l)}$ are MLPs with residual connections. The diffusion embedding ${𝒆}_{\tau }$ is injected via FiLM in each layer.

Decoder.

The decoder maps the updated mesh representations back to the physical grid to produce the predicted noise: (Decoder)$$\begin{array}{ccc}& {\widehat{𝜺}}_j=f_{\mathrm{dec}}(\{{𝒉}_i^{(L)}{\}}_{i\in ℳ(j)};{𝜽}_{\mathrm{dec}}),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $ℳ(j)$ denotes the mesh neighbours of physical grid point $j$. The output $\widehat{𝜺}\in {\mathrm{ℝ}}^d$ has the same dimension as the atmospheric state.

12-Hour Autoregressive Steps and Rollout

GenCast generates extended forecasts (up to 15 days, 360 hours) by autoregressively applying the conditional diffusion model with 12-hour time steps. At each step, a fresh sample is drawn from $p_{𝜽}(\cdot |𝒄)$ using the learned reverse process.

Why 12-hour steps?

The choice of 12-hour autoregressive steps is deliberate and reflects a trade-off between error accumulation and computational cost. Shorter steps (e.g., 6 hours) require more autoregressive iterations and accumulate more sampling variance; longer steps (e.g., 24 hours) reduce the conditioning resolution and increase per-step prediction error. The 12-hour step matches ECMWF ENS's output frequency and was found empirically to minimise rollout drift on the ERA5 validation period (2019–2021).

Independent ensemble members.

At each autoregressive step, a new Gaussian noise sample ${𝑿}_T\sim 𝒩(\mathrm{𝟎},𝐈)$ is drawn independently for each ensemble member. The diffusion reverse process then denoises this noise, conditioned on the current ensemble member's atmospheric state. Crucially, unlike NWP ensemble methods that construct initial perturbations through singular vectors or ensemble Kalman filters, GenCast's ensemble spread arises entirely from the stochasticity of the reverse diffusion sampling. This is both its strength-no perturbation design required-and a source of potential miscalibration if the learned distribution differs from the true climatological spread.

Fast Sampling: DDIM-Style Inference

Running the full $T=1000$ step reverse diffusion at inference time is impractical for operational forecasting. GenCast employs a denoising diffusion implicit model (DDIM; Song et al., 2020) sampler that achieves high-quality samples in a small number of denoising steps $S\ll T$ by using a non-Markovian reverse process.

Operational speed.

Using $S=20$ DDIM steps, GenCast generates a single ensemble member in 8 minutes on a single TPU v4. By comparison, a single ECMWF ENS member at operational resolution requires approximately 30–60 minutes of dedicated high-performance computing. This means an ensemble of 50 members-comparable to the operational ENS-can be generated in 50 $\times$ 8 = 400 minutes on a single TPU, or roughly 8 minutes on a TPU pod of 50 chips.

TikZ Figure: Forward Diffusion and Reverse Denoising

Performance: Verification Against ECMWF ENS

GenCast was evaluated against ECMWF ENS on an independent test set covering 2019–2023, using the standard WMO verification framework. The comparison used 50-member ensembles from both systems.

Headline result.

Across 1,320 verification targets-consisting of combinations of variables (geopotential, temperature, wind components, precipitation, mean sea-level pressure), pressure levels, and lead times from 12 to 360 hours-GenCast outperformed ECMWF ENS on 97.4% of targets as measured by the Continuous Ranked Probability Score (CRPS).

Spread–skill ratio.

A well-calibrated ensemble should have spread (ensemble standard deviation) approximately equal to skill (root-mean-square error of the ensemble mean). GenCast achieves spread–skill ratios close to 1.0 across lead times, slightly outperforming ECMWF ENS at medium range (days 7–15), where NWP models tend to be overdispersive.

Rank histogram calibration.

For a perfectly calibrated ensemble of $M$ members, the observation should be equally likely to fall in any of the $M+1$ rank positions defined by the sorted ensemble. GenCast's rank histograms are approximately flat across variables and lead times, indicating good calibration. Minor u-shaped departures at very short lead times ($<$48 hours) reflect the challenge that the diffusion model was not specifically tuned for deterministic day-1 forecasting.

Extreme Event Forecasting

Extreme weather events-tropical cyclones, heatwaves, extratropical storms, wind power extremes-are precisely the situations where ensemble forecasting matters most. A deterministic forecast conveys no information about low-probability, high-impact outcomes.

Hurricane Beryl (2024).

Hurricane Beryl (2024) was the earliest Category 5 Atlantic hurricane on record, reaching peak intensity in late June. GenCast ensembles initialised 5 days prior to landfall captured the storm's track and intensification with calibrated uncertainty. The ensemble spread correctly indicated high track uncertainty at day 4 but converged to tight agreement on landfall location at day 2, providing operationally useful guidance.

Wind power tails.

GenCast's probabilistic forecasts of 10-metre wind speed are evaluated using the Energy Score for multivariate distributions and the Threshold-Weighted CRPS (twCRPS) for upper-tail events: (Twcrps)$$\begin{array}{ccc}& \mathrm{twCRPS}(F,y;\tau )={\int }_{\tau }^{\infty }{[F(z)-\mathrm{𝟏}(z\ge y)]}^{2}dz,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\tau$ is a threshold (e.g., the 90th or 99th climatological percentile). GenCast outperforms ECMWF ENS on twCRPS for wind speed at the 95th percentile threshold at most lead times, indicating better characterisation of wind energy extremes.

Heatwave skill.

For 2-metre temperature extremes, GenCast is evaluated using the Brier Score for the probability of exceeding the 90th climatological percentile: (Brier)$$\begin{array}{ccc}& \mathrm{BS}(p,y)={(p-{\mathrm{𝟏}}_{y>\tau })}^2,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $p$ is the forecast probability and $y$ is the observation. At 240-hour (10-day) lead time, GenCast achieves Brier Skill Scores (BSS, relative to climatology) that exceed ECMWF ENS by 3–8 percentage points over Europe during the 2021–2023 test period.

Key Idea: Diffusion as Ensemble Generator

Exercises

SEEDS, DEFfusion, and the Diffusion Ecosystem

GenCast demonstrated that a single conditional diffusion model trained end-to-end on reanalysis data could produce a state-of-the-art probabilistic medium-range forecast. The broader research community quickly recognised that the diffusion framework is not a monolithic solution but rather a family of tools, each suited to a different part of the atmospheric prediction pipeline: ensemble generation, data assimilation, statistical downscaling, nowcasting, and superensemble combination. This section surveys the emerging ecosystem of diffusion-based weather models, with particular attention to SEEDS (Li et al., 2024), DEFfusion (Andrae et al., 2025), GenEPS (Nai et al., 2025), CasCast (Gong et al., 2024), and DiffDA (Huang et al., 2024).

SEEDS: Scalable Ensemble Envelope Diffusion Sampler

SEEDS (Scalable Ensemble Envelope Diffusion Sampler; Li et al., 2024, Google Research) takes a different philosophy from GenCast. Rather than learning the atmospheric distribution from scratch, SEEDS trains a diffusion model to emulate an existing operational ensemble system-in this case, the ECMWF ENS-using a small number of operational ensemble members as conditioning information. The goal is to generate thousands of statistically consistent additional members at negligible marginal cost, thereby mapping the tail of the forecast distribution with far higher resolution than the operational 50-member ensemble allows.

Motivation: tail coverage.

An operational ensemble of $M=50$ members can estimate the $p$-th percentile of the forecast distribution with standard error $\sqrt{p(1-p)/M}$. For $p=0.99$ (99th percentile, relevant for extreme events), this gives a standard error of approximately 0.014 in percentile units-or equivalently, about $\pm 0.7$ members. To estimate the 99.9th percentile with comparable relative accuracy requires $M\approx \mathrm{10,000}$ members. SEEDS is designed to generate exactly this order of magnitude.

SEEDS architecture.

SEEDS is a latent diffusion model operating in a compressed embedding space of the atmospheric state. The forward process corrupts a target ensemble member ${𝒙}^{\ast }$, and the reverse process is conditioned on a small set of operational ensemble members $𝒄=\{{𝒙}^{(1)},\dots ,{𝒙}^{(K)}\}$ (typically $K=2$ or $K=4$): (Latent)$$\begin{array}{ccc}& p_{𝜽}^{\mathrm{SEEDS}}({𝒙}^{\ast }|𝒄)=\int p_{𝜽}({𝒙}^{\ast }|{𝒛}^{\ast },𝒄)p_{𝜽}({𝒛}^{\ast }|𝒄)d{𝒛}^{\ast },& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒛}^{\ast }$ is a latent variable obtained by encoding ${𝒙}^{\ast }$ into a lower-dimensional space via a variational autoencoder (VAE). The latent diffusion model $p_{𝜽}({𝒛}^{\ast }|𝒄)$ operates in this compressed space, dramatically reducing the cost of the diffusion forward and reverse passes.

Key capability: 10,000+ members.

Because the latent diffusion operates in a compressed space and the conditioning context is small (2–4 operational members), SEEDS can generate a single additional ensemble member in seconds. Li et al. (2024) demonstrate generation of 10,000 ensemble members for operational forecast periods, enabling high-resolution mapping of the 99.9th percentile of the forecast distribution across the globe.

Training procedure.

SEEDS is trained on historical pairs $(𝒄,{𝒙}^{\ast })$ where $𝒄$ are operational ENS members available at forecast initialisation and ${𝒙}^{\ast }$ is a held-out ENS member treated as the target. The training objective is the standard noise-prediction loss applied in the latent space: (LOSS)$$\begin{array}{ccc}& {ℒ}^{\mathrm{SEEDS}}(𝜽)={𝔼}_{t,{𝒛}_0^{\ast },𝜺}\Vert 𝜺-{𝜺}_{𝜽}^{\mathrm{SEEDS}}(\sqrt{{\bar{\alpha }}_t}{𝒛}_0^{\ast }+\sqrt{1-{\bar{\alpha }}_t}𝜺,t,f_{\mathrm{enc}}(𝒄)){\Vert }^2,& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒛}_0^{\ast }=f_{\mathrm{enc}}({𝒙}^{\ast })$ and $f_{\mathrm{enc}}$ is a frozen VAE encoder. At inference time, the latent code is decoded via $f_{\mathrm{dec}}$ to recover a full-resolution ensemble member.

DEFfusion: Direct Ensemble Forecasting with Diffusion

DEFfusion (Direct Ensemble Forecasting with Diffusion; Andrae et al., 2025) takes yet another approach: rather than emulating an existing NWP ensemble or replacing deterministic models, it trains a diffusion model directly as an autoregressive ensemble forecasting system, using observations (ERA5) rather than NWP ensembles as training targets.

Architecture.

DEFfusion combines a convolutional autoencoder (to compress the global atmospheric state) with a denoising diffusion model operating in the latent space. Unlike GenCast's graph-based architecture, DEFfusion uses a U-Net backbone with attention at multiple resolutions, following the design used in image-generation diffusion models (Ho et al., 2020; Rombach et al., 2022).

(UNET)$$\begin{array}{ccc}& \widehat{𝜺}=𝒰({𝒛}_t^{(\tau )},\tau ,{𝒛}_{𝒄};{𝜽}_{𝒰}),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝒰$ is the U-Net score network, ${𝒛}_t^{(\tau )}$ is the latent noisy state at diffusion step $\tau$, and ${𝒛}_{𝒄}$ is the latent encoding of the conditioning context produced by the autoencoder encoder $f_{\mathrm{enc}}^{\mathrm{DEF}}(𝒄)$.

CRPS and spread–skill improvements.

DEFfusion was evaluated on medium-range forecasting benchmarks including the WeatherBench-2 dataset. On 2-metre temperature at 72-hour lead time, DEFfusion achieves:

• CRPS improvement: 8.3% reduction versus the persistence ensemble baseline.

• Spread–skill ratio: 0.98 versus 0.82 for a deterministic model with Gaussian noise added post-hoc.

• Energy Score: 5.1% improvement for joint verification of temperature and geopotential over Europe.

The improved spread–skill ratio reflects that DEFfusion's ensemble spread is dynamically generated by the diffusion process, intrinsically adapted to the atmospheric state, rather than being a fixed noise level.

Comparison with deterministic baseline.

A natural question is whether the CRPS improvement from DEFfusion is simply due to the ensemble size (more members reduce sampling variance in CRPS estimation) or a genuine improvement in distributional accuracy. Andrae et al. (2025) control for this by comparing against a Gaussian ensemble with the same mean and variance as DEFfusion. DEFfusion outperforms the Gaussian ensemble by 4.2% on CRPS at 72 hours, indicating that the non-Gaussian structure learned by the diffusion model is genuinely informative.

GenEPS: Generative Superensemble Prediction System

GenEPS (Generative Superensemble Prediction System; Nai et al., 2025) exploits a key property of diffusion models: they can serve as unconditional priors over the atmospheric phase space, independently of any specific forecast system.

Unconditional diffusion prior.

Nai et al. train an unconditional diffusion model on ERA5 reanalysis data, learning the climatological distribution $q(𝒙)$ of atmospheric states without any conditioning. The resulting model can evaluate the density $q_{𝜽}(𝒙)$ of any atmospheric state, approximated via the Stein score: (Score)$$\begin{array}{ccc}& {\nabla }_{𝒙}\log q_{𝜽}(𝒙)\approx -\frac{1}{{\sigma }^2}(𝒙-{𝜺}_{𝜽}^{-1}(𝒙,\sigma )),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\sigma$ is a small noise level and ${𝜺}_{𝜽}^{-1}$ denotes the denoised prediction at noise level $\sigma$. States with high climatological probability have large gradient norms pointing toward the high-density region of the atmospheric manifold.

Hybrid superensemble.

The GenEPS superensemble combines output from multiple AI weather models (e.g., GraphCast, Pangu-Weather, FourCastNet) with the NWP ensemble (ECMWF ENS) using the diffusion prior to weight members according to their climatological plausibility: (Weights)$$\begin{array}{ccc}& w^{(m)}\propto \exp (\lambda \log q_{𝜽}({𝒙}^{(m)})),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒙}^{(m)}$ is the $m$-th superensemble member (from any source), $q_{𝜽}$ is the diffusion prior density, and $\lambda >0$ is a temperature parameter. Members with higher climatological plausibility receive higher weight. The superensemble mean is then: (MEAN)$$\begin{array}{ccc}& {\overline{𝒙}}_{\mathrm{GenEPS}}=\frac{{\sum }_m{w}^{(m)}{𝒙}^{(m)}}{{\sum }_m{w}^{(m)}}.& \text{<span class="tml-eqn"></span>}\end{array}$$

Results.

On a superensemble combining 50 ECMWF ENS members, 50 GraphCast trajectory samples, and 50 Pangu-Weather samples, GenEPS achieves CRPS improvements of 6–12% over the unweighted superensemble at lead times of 5–10 days, particularly for rare atmospheric states such as blocking anticyclones and atmospheric rivers.

CasCast: Cascaded Diffusion for Precipitation Nowcasting

Precipitation nowcasting-predicting rainfall at high spatial resolution (1–4 km) and short lead times (0–6 hours)-presents qualitatively different challenges from medium-range forecasting. The precipitation field is highly non-Gaussian, spatially intermittent (sparse with intense local maxima), and temporally rapid-evolving. Classical deep learning nowcasters (e.g., DeepMind's MetNet, DGMRs) have been applied with some success, but tend to produce blurry predictions at longer nowcast lead times.

Cascaded diffusion architecture.

CasCast (Gong et al., 2024) addresses this by using a cascaded diffusion model: a series of diffusion models operating at progressively finer spatial resolutions, each conditioned on the output of the previous (coarser) model. (Cascade)$$\begin{array}{ccc}& p_{𝜽}({𝒙}^{(1)},\dots ,{𝒙}^{(R)}|𝒄)=p_{𝜽}^{(1)}({𝒙}^{(1)}|𝒄)\prod _{r=2}^R{p}_{𝜽}^{(r)}({𝒙}^{(r)}|{𝒙}^{(r-1)},𝒄),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒙}^{(r)}$ denotes the precipitation field at resolution level $r$ ($r=1$ is coarsest, $r=R$ is finest), and each $p_{𝜽}^{(r)}$ is a separate diffusion model conditioned on the coarser-resolution output.

Why cascading works.

The coarse-to-fine cascade exploits the multi-scale structure of precipitation: large-scale organisation (fronts, convective systems) is captured at coarse resolution, while fine-scale structure (individual convective cells, orographic precipitation) is resolved at fine resolution conditioned on the coarse prediction. This avoids the mode collapse that single-resolution diffusion models can suffer when trying to simultaneously capture both scales.

Verification.

On the SEVIR benchmark (Storm EVent ImagRy; Veillette et al., 2020) with 384$\times$384 pixel frames at 1-km resolution, CasCast achieves:

• CSI (Critical Success Index) at 74 dBZ: 0.387 (CasCast) vs. 0.342 (DGMR) vs. 0.318 (EarthFormer).

• MSE at 60-minute lead: 35.2 vs. 41.8 (DGMR).

• FID (Fréchet Inception Distance for sharpness): 4.1 vs. 8.6 (DGMR), indicating significantly sharper samples.

DiffDA: Diffusion for Data Assimilation

Data assimilation (DA) is the process of combining prior model forecasts (the background state) with observations to produce an analysis-an optimal estimate of the true atmospheric state-that serves as the initial condition for the next forecast. Classical DA methods (4D-Var, Ensemble Kalman Filter) are expensive because they require either adjoint model integration or large ensembles to estimate the background error covariance.

Diffusion DA framework.

DiffDA (Huang et al., 2024) formulates DA as a conditional sampling problem: given observations $𝒚$ and a prior (background) state ${𝒙}^b$, draw samples from the posterior $p({𝒙}^a|𝒚,{𝒙}^b)$. By training a diffusion model as a learned prior over analysis states, DiffDA avoids the explicit construction of background error covariance matrices.

The observation operator $\mathscr{H}$ maps from state space to observation space (e.g., extracting temperature at specific radiosonde locations, or simulating satellite radiances). The observation likelihood is: (Likelihood)$$\begin{array}{ccc}& p(𝒚|𝒙)\propto \exp (-\frac{1}{2}{(𝒚-\mathscr{H}(𝒙))}^{𝖳}{𝐑}^{-1}(𝒚-\mathscr{H}(𝒙))),& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝐑$ is the observation error covariance matrix. DiffDA combines this likelihood with the diffusion prior via diffusion posterior sampling (DPS; Chung et al., 2022): at each reverse diffusion step, the gradient of the log-likelihood with respect to the current noisy state is used to guide the denoising toward states consistent with the observations: (DPS)$$\begin{array}{ccc}& {𝒙}_{t-1}\leftarrow {𝒙}_{t-1}^{\mathrm{uncond}}-{\rho }_t{\nabla }_{{𝒙}_t}\Vert 𝒚-\mathscr{H}({\widehat{𝒙}}_0({𝒙}_t)){\Vert }_{{𝐑}^{-1}}^2,& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒙}_{t-1}^{\mathrm{uncond}}$ is the unconditional reverse step, ${\widehat{𝒙}}_0({𝒙}_t)$ is the current denoised estimate, and ${\rho }_t>0$ is a step-size schedule.

Results.

On the ERA5-based DA benchmark (assimilating simulated radiosonde and satellite observations), DiffDA achieves analysis errors within 8–15% of 4D-Var at 10% of the computational cost, making it a promising direction for operational DA-particularly in data-sparse regions where classical DA methods perform poorly due to insufficient observations to constrain the error covariance.

Continuous Ensemble Forecasting

A natural extension of the discrete-time diffusion framework is continuous ensemble forecasting, which models the forecast distribution as the output of a continuous normalising flow or score-based SDE. In this framework, the ensemble is not a finite set of members but a continuous family of trajectories parameterised by the initial noise ${𝑿}_0\sim 𝒩(\mathrm{𝟎},𝐈)$.

Flow matching as an alternative.

Conditional flow matching (Lipman et al., 2022; Liu et al., 2022) offers an alternative continuous framework that trains a vector field ${𝒖}_{𝜽}(𝑿,\tau ,𝒄)$ directly by regressing onto the optimal transport (OT) flow between the prior $𝒩(\mathrm{𝟎},𝐈)$ and the data distribution $q(\cdot |𝒄)$: (FLOW Matching)$$\begin{array}{ccc}& {ℒ}^{\mathrm{FM}}(𝜽)={𝔼}_{\tau ,{𝒙}_0,𝒁}\Vert {𝒖}_{𝜽}(\tau {𝒙}_0+(1-\tau )𝒁,\tau ,𝒄)-({𝒙}_0-𝒁){\Vert }^2.& \text{<span class="tml-eqn"></span>}\end{array}$$ Flow matching converges faster than DDPM training and avoids the numerical instabilities of the reverse SDE integration, making it an active research direction for atmospheric modelling (Watt-Meyer et al., 2024; Nguyen et al., 2024).

Comparison Table

Model	Training data	Sampling cost${}^{†}$	Use case	Key feature
GenCast	ERA5 40yr	8 min/member	Medium-range	End-to-end
(Price et al., 2025)	reanalysis	(single TPU)	ensemble	probabilistic
[2pt] SEEDS	ECMWF ENS	$<$1 min/member	Tail mapping	10k+ members
(Li et al., 2024)	members	(CPU feasible)	risk assessment	cheap
[2pt] DEFfusion	ERA5	12 min/member	Direct	U-Net latent
(Andrae et al., 2025)	reanalysis	(single GPU)	ensemble	diffusion
[2pt] GenEPS	ERA5 (prior)	Negligible	Superensemble	Diffusion
(Nai et al., 2025)	+ AI models	(weighting)	combination	prior weighting
[2pt] CasCast	Radar +	3 min/member	Precipitation	Cascaded
(Gong et al., 2024)	ERA5	(GPU)	nowcasting	multi-scale
[2pt] DiffDA	ERA5	15 min/analysis	Data	Posterior
(Huang et al., 2024)	reanalysis	(GPU)	assimilation	sampling

TikZ Figure: Taxonomy of Diffusion Weather Approaches

Diffusion for Statistical Downscaling

Statistical downscaling-learning a mapping from coarse-resolution model output to fine-resolution fields-is a natural application of image-to-image translation diffusion models. In climate science, downscaling bridges global climate projections (typically 25–100 km resolution) with regional impact assessments that require 1–10 km resolution for hydrology, agriculture, and urban planning.

Formulation.

Let ${𝒙}^{\mathrm{coarse}}\in {\mathrm{ℝ}}^{d_c}$ be a coarse-resolution atmospheric field and ${𝒙}^{\mathrm{fine}}\in {\mathrm{ℝ}}^{d_f}$ with $d_f\gg d_c$ be the corresponding fine-resolution field. A downscaling diffusion model learns: (Model)$$\begin{array}{ccc}& p_{𝜽}^{\mathrm{down}}({𝒙}^{\mathrm{fine}}|{𝒙}^{\mathrm{coarse}})=\int \prod _{s=1}^T{p}_{𝜽}^{(s)}({𝑿}_{s-1}^{\mathrm{fine}}|{𝑿}_s^{\mathrm{fine}},{𝒙}^{\mathrm{coarse}})d{𝑿}_{1:T}^{\mathrm{fine}},& \text{<span class="tml-eqn"></span>}\end{array}$$ where the conditioning on ${𝒙}^{\mathrm{coarse}}$ (typically upsampled and concatenated with the noisy fine-resolution state) steers the diffusion toward fine-resolution fields consistent with the coarse input.

Stochastic downscaling.

A key advantage over deterministic super-resolution methods is that the diffusion model produces a distribution of plausible fine-resolution fields, each consistent with the coarse input but differing in fine-scale detail. This is essential for climate impact studies where the uncertainty in small-scale features (e.g., local extreme precipitation) is as important as the mean.

Insight: Diffusion as a Family of Atmospheric Tools

Warning: Sampling Cost and Speed-Up Research

The Diffusion–NWP Interface

An important architectural question is where in the forecasting pipeline diffusion models are most effectively deployed relative to classical NWP. Three distinct paradigms have emerged:

Replace-NWP paradigm (GenCast, DEFfusion).

The diffusion model replaces the NWP ensemble entirely, generating atmospheric states from scratch conditioned on the analysis. This requires the model to learn all physical dynamics from data, placing heavy demands on training data coverage and model capacity. The risk is physical inconsistency: samples may violate conservation laws or produce unphysical atmospheric structures not penalised by the training loss.

Augment-NWP paradigm (SEEDS, GenEPS).

The diffusion model augments an existing NWP ensemble by generating additional members or reweighting existing ones. This leverages NWP's physical consistency guarantees while adding statistical power. The risk is that the diffusion model inherits biases from the NWP system it is trained to emulate.

Correct-NWP paradigm (post-processing diffusion).

A diffusion model learns the mapping from NWP forecasts to verifying analyses, correcting systematic biases in the NWP output while preserving its large-scale structure. This is the most conservative paradigm and requires the least training data but cannot generate truly independent ensemble members.

Exercises

Functional Generative Networks and Joint Distributions

The probabilistic revolution in weather forecasting has been extraordinarily successful at producing calibrated marginal distributions. Given a single location, a trained ensemble model can output a probability distribution over tomorrow's temperature that correctly assigns, say, 10% probability to values above ${32}^{\circ }$C and passes every standard calibration test. Yet a meteorologist looking at the resulting ensemble members may feel profound unease: something is wrong with the spatial pattern. Individual ensemble members look like random noise fields rather than physically coherent weather maps. A temperature front in one member appears five hundred kilometres north of where it sits in another, even though each member is individually well-calibrated.

This tension between marginal correctness and joint coherence is the central problem of this section. We formalise it, exhibit a sharp counterexample, and then develop the Functional Generative Network (FGN) framework introduced by Alet et al. (2025) as the principled solution.

The Joint Distribution Challenge

Let $\mathrm{\Omega }\subset {ℝ}^2$ be the spatial domain of a numerical weather prediction grid, with $|\mathrm{\Omega }|=N$ grid points. At lead time $\tau$, the atmosphere is described by a state vector ${𝒙}^{(\tau )}\in {ℝ}^{N\times V}$, where $V$ indexes meteorological variables (temperature, wind components, geopotential, humidity, etc.). A probabilistic forecast is a distribution $p({𝒙}^{(\tau )}|{𝒙}^{(0)})$, where ${𝒙}^{(0)}$ is the analysis state at initialisation.

Marginal CRPS optimisation.

The workhorse evaluation metric in NWP is the Continuous Ranked Probability Score (CRPS). For a scalar forecast $F$ (a distribution) and observation $y$, it is defined as:

The marginal CRPS for a gridded forecast averages ((CRPS Ensemble)) over all grid points independently: (MEAN CRPS)$$\begin{array}{ccc}& \overline{\mathrm{CRPS}}=\frac{1}{N}\sum _{i=1}^N{\widehat{\mathrm{CRPS}}}_i,& \text{<span class="tml-eqn"></span>}\end{array}$$ where the index $i$ runs over grid points and ${\widehat{\mathrm{CRPS}}}_i$ is computed using the marginal ensemble at location $i$. The critical observation is that $\overline{\mathrm{CRPS}}$ depends only on marginals: it is completely insensitive to the joint structure of the ensemble.

A Counterexample: Perfect Marginals, Incoherent Storms

We now exhibit a concrete construction showing that attaining the minimum achievable marginal CRPS at every location is entirely compatible with producing storm forecasts of zero physical utility.

The counterexample is not a pathology. It reflects a fundamental mathematical fact: the marginal distributions of a joint distribution carry strictly less information than the joint distribution itself. Optimising a sum of marginal objectives cannot, in general, constrain the joint. What is needed is an architectural or training mechanism that enforces joint coherence without requiring a joint likelihood objective (which would be intractable for high-dimensional forecasts).

Functional Generative Networks: Formal Definition

Alet et al. (2025) introduce Functional Generative Networks (FGN) as an architectural solution to the joint coherence problem. The key insight is deceptively simple: rather than injecting stochasticity into the data space (as diffusion models do), inject it into the function space of the network itself.

The 32-Dimensional Noise Vector

In the WeatherNext 2 implementation, $d=32$. This dimension is not arbitrary. The atmospheric state on a global grid has $𝒪({10}^7)$ degrees of freedom, but the meteorologically relevant stochastic modes-large-scale circulation regimes, blocking events, jet-stream positions, cyclone tracks-lie on a much lower-dimensional manifold. Empirically, 32 dimensions captures the dominant sources of forecast uncertainty without over-parametrising the noise space.

Formally, the modulation at layer $l$ can be written in terms of a scale-and-shift: (Adanorm)$$\begin{array}{ccc}& {\widetilde{𝒉}}^{(l)}(𝒙,𝒛)=(\mathrm{𝟏}+{𝜸}^{(l)}(𝒛))\odot \frac{{𝒉}^{(l)}(𝒙)-{𝝁}^{(l)}}{{𝝈}^{(l)}}+{𝜷}^{(l)}(𝒛),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝜸}^{(l)},{𝜷}^{(l)}:{ℝ}^d\to {ℝ}^{d_l}$ are learned affine maps from the noise vector, and ${𝝁}^{(l)},{𝝈}^{(l)}$ are the layer normalisation statistics. This is analogous to the Adaptive Layer Norm (AdaLN) conditioning used in DiT architectures, but here the conditioning signal is a stochastic latent vector rather than a deterministic class embedding.

Why FGN Works: Coherent Propagation Through Layers

The key theorem underlying FGN is that functional perturbations propagate equivariantly through the network, preserving spatial correlations that data-space noise would break.

The practical implication is profound. Because every output grid point $i$ shares the same $𝒛$, the correlation between ${\hat{y}}_i$ and ${\hat{y}}_j$ is mediated by the network's receptive field-precisely the physical correlations that the deterministic backbone was trained to respect. A large receptive field (as in graph-based models) automatically yields large-scale spatial coherence in the ensemble.

Training with Marginal CRPS

The remarkable property of FGNs is that they can be trained purely on marginal CRPS objectives, yet the architecture guarantees joint coherence in the resulting ensemble. The training objective is:

WeatherNext 2: Operational FGN at Google DeepMind

WeatherNext 2 (Google DeepMind, 2025) represents the first large-scale operational deployment of the FGN paradigm. It inherits the graph-based architecture of GraphCast and GenCast but replaces diffusion-based ensemble generation with functional noise injection.

Architecture Overview

The WeatherNext 2 architecture rests on three components:

1. Multi-mesh graph encoder. The global atmosphere is represented on an icosahedral multi-mesh with $\sim \mathrm{40,000}$ nodes at ${0.25}^{\circ }$ resolution (the operational ERA5 grid). Grid-point features are lifted to a latent graph representation via learned encoder edges, following the GraphCast design.

2. Graph neural network processor. A deep GNN with message-passing layers processes the latent mesh. Each layer applies: (WN2 Message)$$\begin{array}{ccc}& {𝒉}_i^{(l+1)}={\varphi }_v({𝒉}_i^{(l)},\sum _{j\in 𝒩(i)}{\varphi }_e({𝒉}_i^{(l)},{𝒉}_j^{(l)},{𝒆}_{ij})),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${\varphi }_v$ and ${\varphi }_e$ are MLP update and edge functions, and ${𝒆}_{ij}$ encodes geographic edge features. The FGN noise injection modulates both ${\varphi }_v$ and ${\varphi }_e$ through AdaLN conditioning on $𝒛$.

3. Decoder and ensemble generation. The latent mesh representation is decoded back to the physical grid. An ensemble of $K=50$ members is generated by drawing 50 independent ${𝒛}^{(k)}\sim 𝒩(0,{𝐈}_{32})$ and running a single forward pass per member-no iterative sampling required.

Speed and Accuracy

Replacing 20 iterative diffusion steps with a single forward pass per member cuts the network evaluations per member by roughly $20\times$. A 50-member 15-day global ensemble forecast requires approximately 45 seconds on a single TPUv5 chip; a single-chip GenCast ensemble of the same size needs of order hours (Section GenCast and Diffusion Ensemble Forecasting gives about $50\times 8=400$ minutes on one TPU v4), so the wall-clock speedup on a single chip is large.

On the WeatherBench,2 evaluation suite, WeatherNext 2 achieves superior skill to GenCast on 99.9% of all evaluated targets, spanning:

• Deterministic metrics: RMSE and anomaly correlation coefficient (ACC) for all upper-air and surface variables.

• Probabilistic metrics: CRPS, spread/skill ratio, and reliability diagrams at all lead times from 6 hours to 15 days.

• Joint metrics: energy score and variogram score over multi-location subsets of the forecast domain.

The gains are especially large on joint metrics, confirming the theoretical expectation that FGNs improve joint coherence relative to diffusion-based models.

15-Day Forecasts and Joint Coherence

At 15-day lead times-well into the chaotic regime where deterministic skill has vanished-WeatherNext 2 ensemble members exhibit qualitatively correct synoptic-scale patterns. Individual members show coherent frontal systems, realistic cyclone tracks, and physically plausible blocking configurations. This stands in contrast to ensemble members from models that optimise only marginal objectives, which at 15 days often degenerate into spatially incoherent “spaghetti” fields.

The improvement is quantified by the variogram score:

TikZ: FGN vs. Diffusion Noise Injection Comparison

Example: Continental Heatwave with Spatial Coherence

Skillful Joint Probabilistic Forecasting from Marginals

The paper “Skillful Joint Probabilistic Weather Forecasting from Marginals” (Alet et al., 2025) proves a more general result that underpins the FGN framework. The central result can be paraphrased as follows. (Note that requiring zero marginal CRPS would force each marginal to a Dirac mass and hence a deterministic forecast; the substantive statement uses the minimum achievable marginal score, so the coherence must come from the FGN's architecture, not from the marginal objective alone.)

Exercises

Flow Matching, Latent Diffusion, and Data Assimilation

Score-based diffusion models and their variants have proven extraordinarily capable at sampling from complex high-dimensional distributions. Yet for weather applications they carry a practical burden: the iterative reverse-diffusion process requires many network evaluations, and the diffusion path through data space produces trajectories with little physical meaning. Two complementary directions attack this problem.

Flow matching (Lipman et al., 2023) replaces the stochastic SDE with a deterministic ODE whose velocity field is learned directly, enabling faster and smoother trajectories. Latent diffusion compresses the high-dimensional atmospheric state into a compact latent code before applying the generative model, reducing both memory and compute. Together, these techniques open the door to two of the most computationally demanding problems in weather science: subseasonal-to-seasonal (S2S) prediction out to 30 days, and data assimilation-the reconstruction of the full atmospheric state from sparse, heterogeneous observations.

Score-Based Generative Models: Continuous-Time SDEs

Recall that the modern diffusion framework (Song et al., 2021) unifies DDPM and SMLD under a common continuous-time SDE forward process.

The score-matching training objective is: (Score Matching)$$\begin{array}{ccc}& {ℒ}_{\mathrm{SM}}(𝜽)={𝔼}_{t,{𝒙}_0,{𝒙}_t}[\lambda (t)\Vert s_{𝜽}({𝒙}_t,t)-{\nabla }_{{𝒙}_t}\log p({𝒙}_t|{𝒙}_0){\Vert }^2],& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\lambda (t)>0$ is a loss weighting function and $p({𝒙}_t|{𝒙}_0)$ is available in closed form as a Gaussian for linear forward processes.

Flow Matching: Deterministic Continuous-Time Modelling

Flow matching (Lipman et al., 2023) provides an alternative to score-based models that avoids stochastic sampling entirely. Instead of learning a score function, the network learns a velocity field whose integration transports a simple base distribution to the data distribution.

Marchuk: Generative Latent Flow-Matching DiT for S2S

Marchuk (2026) introduces a generative latent flow-matching Diffusion Transformer (DiT) that achieves state-of-the-art subseasonal-to-seasonal (S2S) prediction at lead times up to 30 days. The architecture combines three innovations:

Latent Space Compression with VAE

Raw atmospheric fields are compressed into a latent space via a convolutional variational autoencoder (VAE): (Marchuk VAE)$$\begin{array}{ccc}& {𝒛}_{\mathrm{atm}}={\mathrm{Encoder}}_{\varphi }(𝒙),\widehat{𝒙}={\mathrm{Decoder}}_{\psi }({𝒛}_{\mathrm{atm}}),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒛}_{\mathrm{atm}}\in {ℝ}^{H^{\prime }\times W^{\prime }\times C}$ with $H^{\prime }=H/r$, $W^{\prime }=W/r$ for spatial compression ratio $r$. At $r=4$, the global atmospheric state at ${0.25}^{\circ }$ resolution ($1440\times 721\times V$ variables) is compressed to $360\times 180\times C$ latent codes-a $V/C$-fold reduction in the channel dimension plus a $16\times$ spatial compression.

The VAE is trained with a combined loss: (Marchuk VAE LOSS)$$\begin{array}{ccc}& {ℒ}_{\mathrm{VAE}}={ℒ}_{\mathrm{recon}}+\beta D_{\mathrm{KL}}(q_{\varphi }(𝒛|𝒙)\Vert p(𝒛)),& \text{<span class="tml-eqn"></span>}\end{array}$$ with $\beta ={10}^{-4}$ in the Marchuk implementation, preserving fine-scale atmospheric structure at the cost of minimal additional KL regularisation.

DiT Backbone with Spatial and Temporal Embeddings

The generative model operates on the latent codes ${𝒛}_{\mathrm{atm}}$. The backbone is a Diffusion Transformer (DiT) with:

• 2D spatial embeddings: rotary positional embeddings (RoPE) applied in the latitude and longitude dimensions, enabling the model to generalise across spatial positions.

• 1D temporal RoPE: a separate rotary embedding for the forecast lead time $\tau$, allowing smooth interpolation between lead times and generalisation beyond the training horizon. Formally, for query/key vectors $𝒒,𝒌$ at lead time $\tau$: (ROPE)$$\begin{array}{ccc}& {\tilde{q}}_j=q_j{e}^{i{\omega }_j\tau },{\tilde{k}}_j=k_j{e}^{i{\omega }_j{\tau }^{\prime }},{\omega }_j={\tau }_{\mathrm{base}}^{-2j/d},& \text{<span class="tml-eqn"></span>}\end{array}$$ where the inner product $⟨\widetilde{𝒒},\widetilde{𝒌}⟩$ depends only on $\tau -{\tau }^{\prime }$ (relative lead time), encoding temporal equivariance.

• Adaptive Layer Norm (AdaLN-Zero): following the DiT convention, flow time $t$ conditions each transformer block through learned scale and shift parameters.

Efficiency: 276M Parameters vs. 1.6B Baseline

The latent-space operation reduces the sequence length for the transformer from $\sim {10}^6$ (full resolution) to $\sim 6.5\times {10}^4$ (latent resolution), enabling a 276M parameter model to match or exceed the performance of 1.6B parameter autoregressive baselines that operate in data space. The effective parameter count for the generative component is: (Marchuk Params)$$\begin{array}{ccc}& N_{\mathrm{Marchuk}}\approx N_{\mathrm{DiT}}+N_{\mathrm{VAE}}=220\text{M}+56\text{M}=276\text{M},& \text{<span class="tml-eqn"></span>}\end{array}$$ compared to $1.6\times {10}^9$ for the best non-latent baseline.

S2S Results

Marchuk (2026) evaluates on the ECMWF S2S benchmark, which provides verification data for weeks 3–6 post-initialisation. Key findings:

• At lead times of 20–30 days, Marchuk achieves higher ACC than ECMWF's 51-member ENS-extended on $Z_{500}$, $T_{850}$, and MJO phase.

• The latent flow-matching ensemble provides better reliability than both IFS-ENS and the 1.6B autoregressive baseline on probabilistic skill scores at weeks 3–4.

• Inference for a 50-member ensemble to 30 days requires $\sim 3$ minutes on a single A100 GPU, enabling rapid regeneration of ensembles for initialisation studies.

Latent Diffusion Architecture for Weather: TikZ Figure

Appa: Latent Diffusion for Data Assimilation

Data assimilation (DA) is the process of combining a prior forecast (the background) with new observations to produce an optimal estimate of the current atmospheric state (the analysis). Classical DA methods-3D-Var, 4D-Var, and the Ensemble Kalman Filter (EnKF)-solve a Bayesian inverse problem: (DA Posterior)$$\begin{array}{ccc}& p(𝒙|{𝒚}_{\mathrm{obs}})\propto p({𝒚}_{\mathrm{obs}}|𝒙)p(𝒙),& \text{<span class="tml-eqn"></span>}\end{array}$$ where ${𝒚}_{\mathrm{obs}}\in {ℝ}^M$ is the observation vector ($M\ll D$) and $p(𝒙)$ is the forecast prior.

Appa (Manshausen et al., 2025) is a 565M-parameter implementation of this framework trained on ERA5 reanalysis data. Its key capabilities:

• Global atmospheric trajectories from sparse observations: given only radiosonde profiles, satellite radiances, and surface station data (representing $<0.1\%$ of grid points), Appa generates full 3D atmospheric states consistent with the observations.

• Direct posterior inference: unlike EnKF, which approximates the posterior as Gaussian, Appa samples from the full non-Gaussian posterior, capturing multimodal uncertainty in data-sparse regions.

• Heterogeneous observation types: the observation operator $𝐇$ can encode any differentiable mapping from atmospheric state to observed quantity-temperature profiles, radiances, radar reflectivities, GPS radio occultations-without requiring separate DA systems for each observation type.

FlowDA: Flow Matching for Low-Latency Data Assimilation

While Appa achieves high accuracy, its score-based reverse SDE requires $𝒪(100)$ network evaluations per sample. For operational DA, which must produce analysis increments within minutes to meet NWP model initialisation deadlines, this latency is prohibitive. FlowDA (Cheng and Zhu, 2026) addresses this by replacing the reverse SDE with a conditional flow-matching ODE.

The key advantage of FlowDA over Appa is inference speed: with $N_{\mathrm{step}}=10$ ODE solver steps, FlowDA generates a full global analysis in $\sim 15$ seconds on a single GPU, compared to $\sim 8$ minutes for Appa with 100 diffusion steps. This $32\times$ speedup makes it viable for operational 6-hourly analysis cycles.

SwinVRNN: VAE-Style Stochastic Ensemble Forecasting

SwinVRNN (Hu et al., 2023) is a variational recurrent neural network for weather ensemble forecasting that predates the diffusion era but anticipates several of the ideas in FGN and Marchuk. It combines a Swin Transformer backbone with a variational latent-variable model:

SwinVRNN generates ensemble members by sampling ${𝒛}_t\sim p({𝒛}_t|{𝒙}_{t-1})$ (the prior, not the posterior) at inference time. Each sample produces a different forecast trajectory, with diversity controlled by the learned prior variance.

LO-SDA: Latent Optimisation for Score-Based Atmospheric DA

An alternative to the sequential SDE approach of Appa is to perform latent optimisation (LO) directly in the latent space of a pre-trained generative model.

The key properties of LO-SDA are:

1. No generative sampling required: the analysis is a point estimate in latent space, obtained by gradient descent through the pre-trained decoder. This is $10$–$100\times$ faster than ensemble SDE approaches.

2. Physical consistency guaranteed by the decoder: since $\widehat{𝒙}$ lies in the image of ${\mathrm{Decoder}}_{\psi }$, it is by construction a sample-like atmospheric state from the training distribution.

3. Uncertainty quantification: ensemble spread can be obtained by running LO-SDA from multiple random initialisations ${𝒛}_0^{(k)}\stackrel{\mathrm{iid}}{\sim }𝒩(0,𝐈)$, approximating the posterior mode structure.

Comparison: Flow Matching vs. Diffusion for Weather

Dimension	Score-Based Diffusion	Flow Matching
Generative path	Stochastic SDE; noisy trajectories through data space	Deterministic ODE; smooth straight-line interpolant
Inference steps	20–100 steps (DDPM/DDIM)	5–20 steps (Euler/RK4)
Exact inversion	Approximate (DDIM) or intractable (DDPM)	Exact (integrate ODE backward)
DA integration	Perturbed reverse SDE; gradient through decoder	Conditional ODE with observation cross-attention
Training signal	Denoising score matching $\nabla \log p_t$	Conditional velocity matching ${𝒙}_1-{𝒙}_0$
Flexibility for novel obs. operators	High (no retraining needed)	Lower (conditional training amortised)
Sample quality (weather)	High; established track record (GenCast)	Comparable; fewer evaluations
Preferred use case	Research DA, full posterior; medium-range	Operational DA, S2S; speed-critical

Exercises

Generative Downscaling and Super-Resolution

A global weather model sees the world at a resolution of roughly 25,km per grid cell. At that scale, a city block is invisible, a valley too narrow to represent, and a cold front compressed into a blurry gradient across three or four cells. Yet the questions that matter most to society - will it flood on this street, will this wind farm produce enough power tonight, will this ski resort have snow on Saturday - demand answers at the kilometre scale or below. The gap between what global models supply and what local decision-makers demand is one of the most consequential open problems in computational meteorology.

Downscaling is the family of methods that bridge this gap: given a coarse-resolution representation of the atmosphere, produce a fine-resolution one. For several decades downscaling divided cleanly into two camps - expensive physics-based simulations and cheap but crude statistical regressions. Generative deep learning has opened a third path that achieves the fidelity of physics-based methods at a fraction of the cost, while simultaneously treating fine-scale weather as the stochastic variable that it truly is.

The Downscaling Problem

Let ${𝒙}_{\text{coarse}}\in {\mathrm{ℝ}}^{H_c\times W_c\times C}$ denote a coarse atmospheric state: temperature, pressure, wind, humidity, and possibly dozens of other fields on a global grid of spacing ${\mathrm{\Delta }}_c\approx 25$,km. Let ${𝒙}_{\text{fine}}\in {\mathrm{ℝ}}^{H_f\times W_f\times C}$ denote the corresponding fine-scale state on a local grid of spacing ${\mathrm{\Delta }}_f\approx 1$–$2$,km. The downscaling ratio is $r={\mathrm{\Delta }}_c/{\mathrm{\Delta }}_f$; typical values range from $r=10$ (25,km to 2.5,km) to $r=25$ (25,km to 1,km).

The probabilistic form is not merely a philosophical nicety. At kilometre scales, the atmosphere is genuinely chaotic: two coarse states that are indistinguishable at 25,km resolution may evolve into qualitatively different fine-scale patterns. A cold front that crosses Chicago at 25,km resolution does not determine whether a particular neighbourhood receives an intense convective squall or merely light rain. Any method that collapses the conditional distribution to a single prediction is suppressing real meteorological uncertainty.

Classical Approaches

Dynamical Downscaling

The most physically rigorous approach is dynamical downscaling: run a high-resolution numerical weather prediction (NWP) model over a limited area, driven at its boundaries by the coarse global model. Models such as WRF (Weather Research and Forecasting), COSMO, and HARMONIE are routinely used for this purpose. The fine-scale model integrates the full primitive equations - conservation of momentum, mass, and energy - at a resolution of 1–4,km, resolving convective updrafts, orographic gravity waves, and sea-breeze circulations that the global model cannot see.

The cost is prohibitive. A single day of 2,km WRF simulation over a domain the size of Taiwan requires on the order of ${10}^5$ CPU core-hours. Running ensembles - necessary for probabilistic forecasting - scales this linearly with the ensemble size. Producing retrospective climate datasets at kilometre resolution over entire continents is effectively impossible with dynamical downscaling alone; the computational demand would require centuries of wall-clock time on current hardware.

Statistical Downscaling

Statistical downscaling fits a function $f$ in equation from a training set of paired coarse–fine observations or reanalyses. Classical methods include:

• Perfect prognosis (PP) methods. Regress fine-scale observations (e.g., rain-gauge records) directly on coarse atmospheric predictors such as 500,hPa geopotential height, boundary-layer moisture, and temperature gradients. Multiple linear regression, logistic regression for precipitation occurrence, and analogue methods (find historical days with similar coarse patterns) all belong to this family.

• Model output statistics (MOS). Correct the systematic biases of a specific NWP model by fitting regression equations between model output and observed surface variables. MOS is operationally deployed by national weather services worldwide.

• Quantile mapping. Apply a transfer function to match the marginal distribution of model output to that of observations, variable by variable and grid cell by grid cell. Widely used in climate impact assessment to correct precipitation and temperature biases.

• Lookup-table and climatological scaling. Multiply or add spatially varying correction factors derived from historical climatology. The simplest method: long-term mean ratio between high-resolution observations and coarse reanalysis at each grid cell.

All statistical methods share a fundamental limitation: they learn marginal statistics at each location independently, failing to capture spatial correlations. A statistically downscaled precipitation field may have the correct distribution at every grid point but may look unrealistically noisy - lacking the organised convective clusters, frontal rainbands, and orographic streaks that characterise real kilometre-scale precipitation.

Generative Downscaling: The Stochastic Inverse Perspective

Downscaling is a textbook inverse problem. The coarse field ${𝒙}_{\text{coarse}}$ is a blurry, incomplete observation of the true fine-scale state ${𝒙}_{\text{fine}}$. We seek to invert this “blurring” operator, but the problem is severely under-determined: many fine-scale states are consistent with a single coarse observation. Classical inversion seeks a unique solution; the correct approach is to sample from the posterior $p({𝒙}_{\text{fine}}|{𝒙}_{\text{coarse}},{𝒛}_{\text{topo}})$.

Generative models - GANs, diffusion models, normalising flows, and latent-variable models - are natural tools for this task because they are designed to sample from high-dimensional conditional distributions. The training objective is to match the joint distribution $p({𝒙}_{\text{coarse}},{𝒙}_{\text{fine}},{𝒛}_{\text{topo}})$ on a dataset of paired coarse–fine atmospheric states, then use the learned model to sample ${𝒙}_{\text{fine}}$ conditioned on new coarse inputs at inference time.

CorrDiff: Two-Stage Diffusion Downscaling

The most influential generative downscaling system to date is CorrDiff (Mardani et al., 2025), developed by NVIDIA Research. CorrDiff introduces a two-stage architecture grounded in the Reynolds decomposition of fluid dynamics, a classical tool for separating the mean and fluctuating components of a turbulent flow.

Reynolds Decomposition and the Two-Stage Idea

In Reynolds decomposition, a turbulent field $u(𝒙,t)$ is split into a time-averaged mean $\bar{u}(𝒙)$ and a fluctuating residual $u^{\prime }(𝒙,t)=u(𝒙,t)-\bar{u}(𝒙)$. The mean and residual follow different dynamical equations and different statistical structures.

CorrDiff applies the same decomposition to the spatial dimension of the downscaling problem: (Decomp)$$\begin{array}{ccc}& {𝒙}_{\text{fine}}=\underset{\text{mean state}}{\underbrace{\overline{𝒙}}}+\underset{\text{stochastic residual}}{\underbrace{{𝒙}^{\prime }}},& \text{<span class="tml-eqn"></span>}\end{array}$$ where $\overline{𝒙}$ is the deterministic mean conditioned on the coarse input, and ${𝒙}^{\prime }$ encodes the fine-scale stochastic structures - convective cells, turbulent eddies, precipitation streaks - that are under-determined by the coarse state.

The separation of concerns is crucial. The U-Net is tasked with a well-conditioned regression: it learns to map coarse structure to mean fine-scale structure, a task where MSE is an appropriate loss. The diffusion model is tasked with a generative problem: it learns the distribution of fine-scale residuals, recovering the power-law tail of the spatial spectrum that the U-Net suppresses.

Training Data and Results

CorrDiff is trained on pairs of 25,km ERA5 reanalysis fields (inputs) and 2,km regional NWP analyses (targets). Two domains were studied in the original paper: Taiwan (a challenging orographic domain with typhoon activity) and the contiguous United States (CONUS, a continental domain with diverse weather regimes). The target fields were generated by a regional NWP system at 2,km, yielding an effective downscaling ratio of $r=12.5$.

The training dataset spans several years of hourly snapshots, comprising temperature at 2,m, wind speed and direction at 10,m, relative humidity, surface pressure, and precipitation accumulation. Static predictors include digital elevation model (DEM) data, land-sea mask, and vegetation fraction.

Recovery of Power-Law Spectra

A key diagnostic in the CorrDiff paper is the spatial power spectrum of downscaled fields. Deterministic downscaling (Stage 1 alone) produces spectrally deficient outputs: at high wavenumbers, the predicted spectrum falls well below the reference spectrum of the 2,km NWP analysis. Adding Stage 2 restores the correct spectral slope.

Formally, let ${\hat{S}}_{\text{fine}}(k)$ denote the power spectrum of a CorrDiff sample and $S_{\text{ref}}(k)$ the spectrum of the 2,km NWP target. The deterministic U-Net achieves ${\hat{S}}_{\text{fine}}(k)\ll S_{\text{ref}}(k)$ for $k>k^{\ast }$, where $k^{\ast }$ is the wavenumber corresponding to roughly 10,km. CorrDiff with both stages achieves ${\hat{S}}_{\text{fine}}(k)\approx S_{\text{ref}}(k)$ across all resolved wavenumbers, matching the power-law slope $\beta \approx 3$ of the reference field.

This result is not surprising from first principles: MSE loss is equivalent to minimising KL divergence in a Gaussian model, and minimum-MSE estimators produce the posterior mean, which is smoother than individual realisations. The diffusion model adds the missing variance.

Computational Efficiency

CorrDiff achieves approximately $1000\times$ speedup over dynamical downscaling with WRF at comparable spatial resolution. On NVIDIA H100 GPUs, generating a 2,km output field over Taiwan takes roughly $30$–$60$,seconds, compared to $\sim$8,hours of wall-clock time for an equivalent WRF simulation. Energy efficiency is even more dramatic: the diffusion-based approach is estimated to be roughly $3000\times$ more energy efficient per generated field, primarily because it avoids the iterative numerical integration of the full primitive equations.

RainScaleGAN: Adversarial Rainfall Downscaling

Precipitation is a special case in atmospheric downscaling. Unlike temperature or wind, precipitation is fundamentally discrete and intermittent: at any given time and place, either it is raining or it is not. The spatial distribution of precipitation is characterised by sharp boundaries (convective cells, frontal rainbands), heavy tails (extreme events that dominate mean annual totals), and multi-scale organisation (from sub-kilometre convective cells to synoptic-scale frontal systems spanning thousands of kilometres).

Mean-squared-error losses are poorly suited to these properties. Minimising MSE over precipitation fields produces diffuse, over-smoothed outputs that spread precipitation over too large an area, miss the sharp edges of convective cells, and underestimate extreme intensities. Generative adversarial networks are better aligned with the task: the discriminator can be trained to distinguish realistic from unrealistic precipitation patterns, enforcing the sharp gradients and heavy tails that MSE suppresses.

RainScaleGAN (Iotti et al., 2025) is a conditional GAN designed specifically for precipitation downscaling at ratios of $r=10$–$20$. The generator $G_{\theta }$ takes a coarse precipitation field and terrain predictors as input and produces a fine-scale precipitation map. The discriminator $D_{\varphi }$ is conditioned on both the coarse input and the terrain, ensuring that it evaluates realism within the physically constrained context.

The intensity constraint addresses a known failure mode of unconditional GANs on precipitation: they tend to hallucinate spurious heavy-rain events that inflate the area-mean precipitation far above the coarse input. By anchoring the spatial mean, RainScaleGAN can freely redistribute precipitation into realistic convective structures while remaining physically consistent.

Tail-Preserving Spatial Statistics

A standard diagnostic for precipitation downscaling is the empirical variogram: the expected squared difference between precipitation values at two points as a function of their separation distance $h$: (Variogram)$$\begin{array}{ccc}& \gamma (h)=\frac{1}{2}𝔼\left[{(P(𝒙+h)-P(𝒙))}^2\right].& \text{<span class="tml-eqn"></span>}\end{array}$$ A model that recovers the correct variogram reproduces the spatial correlation structure of precipitation at all scales, from sub-kilometre to mesoscale.

RainScaleGAN is evaluated against classical PP methods and deterministic CNN baselines. The deterministic CNN achieves low RMSE but fails to match the variogram at fine scales ($h<5$,km), where it underestimates spatial variance by a factor of $2$–$3$. RainScaleGAN produces variograms that match the reference within statistical uncertainty bounds across all spatial scales, including the sub-5,km range critical for convective precipitation. The adversarial loss is key: by training the discriminator to distinguish realistic from unrealistic spatial patterns at multiple scales, the GAN implicitly enforces the variogram.

SpateGAN-ERA5 and Global Precipitation Downscaling

SpateGAN-ERA5 (Glawion et al., 2025) extends adversarial precipitation downscaling to the global scale. The key challenge is scalability: a global 1,km precipitation field contains approximately $5.1\times {10}^8$ grid cells, and training data must span all climate zones, from tropical convective regimes to polar precipitation.

SpateGAN-ERA5 uses ERA5 as the coarse input (approximately 30,km resolution) and targets a 1,km analysis product derived from satellite observations and rain-gauge networks. The architecture is a space-time GAN: the generator and discriminator operate on sequences of precipitation fields to enforce temporal as well as spatial coherence.

Probabilistic Spatiotemporal Coherent Downscaling

A limitation shared by CorrDiff, RainScaleGAN, and SpateGAN-ERA5 is that they operate on individual time steps: the downscaling is applied frame-by-frame without explicit temporal modelling. This can produce temporal flickering in animated output: successive fine-scale frames may be individually plausible but discontinuous in time.

Schmidt et al. (2025) address this with a probabilistic spatiotemporal coherent downscaling framework. The model generates joint samples ${({𝒙}_{\text{fine},t})}_{t=1}^T$ over a sequence of $T$ time steps, with the diffusion model operating in the joint space-time domain. The key insight is that the fine-scale state at time $t$ is conditionally dependent not only on the coarse input at time $t$ but also on the fine-scale states at nearby times, which constrain the temporal evolution through physical continuity.

Schmidt et al. show that temporal coherence significantly improves the quality of derived products such as diurnal precipitation cycles and wind ramp events (rapid changes in wind speed relevant to power systems), where single-frame downscaling loses the temporal structure.

Zero-Shot Diffusion Posterior Sampling for Downscaling

An attractive alternative to training a purpose-built downscaling model is zero-shot diffusion posterior sampling: use a pre-trained unconditional diffusion model over fine-scale atmospheric fields as a prior, then condition on the coarse observation at inference time using posterior sampling techniques.

Let $p_{\theta }({𝒙}_{\text{fine}})$ be the learned prior over fine-scale fields and let $p({𝒙}_{\text{coarse}}|{𝒙}_{\text{fine}})$ be a known degradation model (e.g., spatial averaging with a Gaussian kernel). The downscaling posterior is: (Posterior)$$\begin{array}{ccc}& p({𝒙}_{\text{fine}}|{𝒙}_{\text{coarse}})\propto p({𝒙}_{\text{coarse}}|{𝒙}_{\text{fine}})p_{\theta }({𝒙}_{\text{fine}}).& \text{<span class="tml-eqn"></span>}\end{array}$$ Sampling from this posterior can be achieved with the diffusion posterior sampling (DPS) algorithm: at each denoising step, the score function of the prior ${\nabla }_{𝒙}\log p_{\theta }({𝒙}_t)$ is augmented with a data-consistency gradient ${\nabla }_{{𝒙}_t}\log p({𝒙}_{\text{coarse}}|{\widehat{𝒙}}_{\text{fine}}({𝒙}_t))$, where ${\widehat{𝒙}}_{\text{fine}}({𝒙}_t)$ is the Tweedie estimate of the clean sample from the noisy sample ${𝒙}_t$.

The zero-shot approach requires no paired coarse–fine training data for the specific downscaling problem at hand, which is highly advantageous for rare events or future climate scenarios where such paired data does not exist. The trade-off is that the degradation model $p({𝒙}_{\text{coarse}}|{𝒙}_{\text{fine}})$ must be specified manually and may not capture the full complexity of the actual coarse model's degradation.

Generative High-Resolution Ensemble for Renewable Energy

A particularly high-value application of generative downscaling is renewable energy forecasting. Wind turbines operate at hub heights of 80–150,m and rotor diameters of 100–260,m, scales far below the resolution of any global weather model. The power output of a wind farm is a nonlinear function of wind speed (proportional to $v^3$ in the rated range), making the tails of the wind speed distribution critically important: an ensemble that correctly captures the probability of extreme winds (whether too low for generation or too high for safe operation) is more valuable than one with a lower mean error.

Wang et al. (2026) present the Generative High-Resolution Ensemble (GHRE) for renewable energy forecasting, a CorrDiff-like architecture that downscales global ensemble forecasts from 25,km to 500,m over wind-farm domains. The model generates an ensemble of fine-scale wind fields from a single coarse ensemble member, effectively multiplying the ensemble size by a factor of $\sim 50$ at negligible additional computational cost.

This result is a compelling practical argument for probabilistic downscaling: beyond capturing uncertainty, it produces better expected value estimates for power generation.

Summary and Comparison of Methods

Exercises

Precipitation Nowcasting: The Gray Zone

On 20 July 2021, in the record hour of 16:00 to 17:00 local time, the city of Zhengzhou, capital of Henan Province, China, received 201.9,mm of rain in a single hour - more than its average rainfall for the whole of July. Over three days the total reached 617.1,mm, close to Zhengzhou's entire average annual rainfall of 640,mm: nearly a year of rain in three days. Metro tunnels flooded, trapping hundreds of commuters. At least 302 people died. The event was not a surprise in hindsight: the ingredients were present, the synoptic setup was identifiable. But short-range operational forecasting failed to anticipate the extreme intensity.

This is the nowcasting problem.

Nowcasting refers to forecasting on timescales of 0–6 hours at spatial resolutions of 1,km or finer. It sits between the domain of numerical weather prediction (which becomes accurate only beyond $\sim$6 hours) and the domain of extrapolation methods (which are accurate only for the first $\sim$30 minutes). This intermediate window - the gray zone - is where precipitation forecasting is hardest, and where the cost of failure is highest.

The Nowcasting Problem

Why Nowcasting Is Hard

Convective Initiation

The hardest part of nowcasting is not predicting where existing rain will move - that is a well-posed advection problem - but predicting where new convection will form. Convective initiation is controlled by sub-kilometre boundary-layer processes: surface heating creating convergence lines, orographic triggers forcing air masses over ridges, land-sea breezes generating coastal convergence. These processes are not resolved by weather radar (which observes precipitation but not the pre-precipitation boundary-layer state) and are highly sensitive to small perturbations. A parcel of air that is 0.1,K warmer than its environment in a nearly neutral boundary layer may trigger deep convection; one that is 0.1,K cooler will not. This extreme sensitivity implies that convective initiation is inherently probabilistic at lead times beyond $\sim$1 hour.

Rapid Error Growth

Atmospheric predictability theory distinguishes the predictability horizon: the timescale beyond which the atmosphere's chaotic dynamics have amplified initial-condition errors to the magnitude of the climatological variability. For synoptic-scale systems, this horizon is $\sim$10 days. For convective-scale features (individual thunderstorms, convective cells), it is $\sim$1 hour. Error doubling times for convective perturbations are on the order of 10–30 minutes - far faster than for large-scale flow.

This means that any deterministic convective-scale forecast at lead times beyond $\sim$2 hours is no longer tracking a specific storm but predicting the distribution of possible storm evolutions. Probabilistic methods are not merely a preference but a theoretical necessity.

Radar Artifacts

Weather radar provides the primary observational input to nowcasting systems. Radar reflectivity is related to precipitation intensity by the Marshall–Palmer relation: (Zpowerlaw)$$\begin{array}{ccc}& Z=a\cdot R^b,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $Z$ is reflectivity in mm${}^6$,m${}^{-3}$, $R$ is rain rate in mm,h${}^{-1}$, and $(a,b)$ are empirical coefficients (typically $a=200,b=1.6$ for stratiform rain). This relation is approximate and varies with drop-size distribution.

Radar observations are corrupted by systematic artifacts: ground clutter (bright returns from buildings and terrain), beam blockage (attenuation behind mountains), beam overshooting (missing near-surface rain at long ranges), bright-band contamination (strong returns from partially melted snowflakes near the 0,${}^{\circ }$C isotherm), and range folding. A nowcasting model trained on raw radar must learn to handle all these artifacts, or pre-processing pipelines must remove them at the cost of information loss.

Classical Approach: Optical Flow Extrapolation

The dominant operational approach to short-range precipitation nowcasting is radar extrapolation: estimate the motion field of the precipitation pattern and advect the current radar image forward in time. Modern operational systems (e.g., PySTEPS - the Probabilistic STEPS) use multi-scale optical flow to estimate the motion field and stochastic perturbations to generate ensemble nowcasts.

Optical-flow nowcasting is accurate for $\tau \lesssim 30$–$60$ minutes: existing precipitation features are advected faithfully as long as their internal structure does not change. Beyond this timescale, the forecast degrades for two reasons: (1) precipitation cells grow and decay independently of advection, invalidating the brightness-constancy assumption; (2) new cells initiate that were absent in the initial radar frame and cannot be extrapolated.

DGMR: Deep Generative Model of Radar

DGMR (Ravuri et al., 2021) was published in Nature and represents the first generative deep learning system for precipitation nowcasting to achieve skill superior to both optical-flow extrapolation and convective-scale NWP at the same computational cost. The paper attracted wide attention both for its meteorological results and for its methodological demonstration that generative models can outperform physics-based systems on high-impact Earth science prediction tasks.

Architecture: Spatial Generator and Temporal Discriminator

DGMR is a conditional GAN with a spatially and temporally structured architecture.

The dual discriminator design is one of DGMR's key architectural contributions. Without the temporal discriminator, the generator can produce individually realistic frames that are discontinuous in time. Without the spatial discriminator, the generator may produce temporally smooth but spatially blurred sequences. Together they enforce both modes of realism simultaneously.

Evaluation and Results

DGMR was trained on the UK Met Office's radar composites at 1,km resolution over England and Wales, covering 2016–2018 as training data. It was evaluated against the operational STEPS system (the PySTEPS predecessor) and the MetNet deterministic deep learning nowcasting system.

The principal skill metric is the Continuous Ranked Probability Score (CRPS), which rewards probabilistic forecasts that assign high probability to the observed outcome. DGMR achieves lower CRPS than STEPS across all lead times from 5 to 90 minutes, with the largest improvement in the 30–90 minute range where extrapolation begins to break down. The paper also reports the Critical Success Index (CSI) and the Fractions Skill Score (FSS), which we define below.

NowcastNet: Physics-Embedded Generative Nowcasting

NowcastNet (Zhang et al., 2023, Nature) advances beyond purely data-driven approaches by explicitly embedding the physical advection equation into the neural network architecture. The guiding insight is that precipitation motion in the 0–2 hour range is dominated by advection - the spatial transport of rain by the wind field - while growth and decay processes are secondary. By encoding advection as an explicit constraint rather than asking the network to learn it from scratch, NowcastNet achieves a more data-efficient and physically interpretable architecture.

Architecture: Explicit Advection Operator

The explicit advection constraint provides a structural inductive bias that significantly reduces the data required to achieve skill on advection-dominated nowcasting (0–2 hours). The stochastic generation term $S$ provides the probabilistic component necessary for the 2–6 hour gray zone.

Performance on Extreme Precipitation

The most striking result in the NowcastNet paper is the comparison on extreme precipitation events. For precipitation thresholds corresponding to the 99th percentile (approximately 16,mm,h${}^{-1}$ in the China Meteorological Administration radar composites used for training), NowcastNet achieves: (CSI)$$\begin{array}{ccc}& {\mathrm{CSI}}_{\text{NowcastNet}}(q_{99})\approx 0.30,\text{vs.}{\mathrm{CSI}}_{\text{HRRR}}(q_{99})\approx 0.04,& \text{<span class="tml-eqn"></span>}\end{array}$$ where HRRR (High-Resolution Rapid Refresh) is the operational 3,km NWP model of NOAA. The contrast is remarkable: NowcastNet achieves a CSI of 0.30 at 4-hour lead time for extreme precipitation, while the state-of-the-art convective-scale NWP achieves only 0.04. The NWP model fails because it cannot maintain convective cells in the correct location at these lead times; NowcastNet uses the advection operator to track existing cells and the generative component to model their growth.

StormDiT: Diffusion Transformer for the Gray Zone

Both DGMR and NowcastNet were designed primarily for the 0–2 hour range and show declining skill in the gray zone (2–6 hours). StormDiT (Sun et al., 2026) targets precisely this window with a Diffusion Transformer architecture that models the full 6-hour nowcasting problem without decomposing it into advection and growth components.

Architecture

StormDiT replaces the convolutional U-Net backbone of standard meteorological diffusion models with a Diffusion Transformer (DiT): a transformer architecture in which the denoising network operates directly on sequences of space-time patch tokens.

The key architectural advantage of StormDiT over GAN-based and U-Net-based diffusion methods is that the self-attention mechanism allows every token to attend to every other token, coupling spatially and temporally distant parts of the storm system without the limited receptive field constraints of convolutions. This is essential for the 2–6 hour gray zone, where the evolution of a storm cell in one part of the domain can be triggered by a convergence line hundreds of kilometres away.

Training and Validation

StormDiT was trained on an archive of 7700+ significant precipitation events extracted from the China Meteorological Administration radar composites at 1,km, 6-minute resolution. Events were selected based on the criterion that the domain-maximum precipitation exceeds 20,mm,h${}^{-1}$ for at least 30 consecutive minutes, ensuring that the training data is dominated by the high-impact events that the model is designed to forecast.

StormDiT achieves FSS$\ge$0.5 (a commonly used skill threshold) at lead times up to 5.5 hours for precipitation thresholds corresponding to 8,mm,h${}^{-1}$. The corresponding lead time for NowcastNet is $\sim$3.5 hours and for DGMR is $\sim$2.5 hours, validating the claim that the DiT architecture provides skill improvement specifically in the gray zone.

CasCast: Cascaded Diffusion for High-Resolution Nowcasting

CasCast (Gong et al., 2024) applies the cascaded diffusion architecture - a sequence of diffusion models operating at progressively finer spatial resolutions - to precipitation nowcasting. The cascade consists of three stages:

1. A coarse-resolution diffusion model generates a 16,km nowcast, capturing synoptic-scale precipitation organisation.

2. A medium-resolution diffusion model (4,km) upscales the 16,km output, adding mesoscale structure.

3. A fine-resolution diffusion model (1,km) adds sub-kilometre convective detail.

At each stage, the conditioning information includes both the radar context and the output of the previous (coarser) stage. This coarse-to-fine structure mirrors the multi-scale organisation of precipitation in the atmosphere, exploiting the physical fact that fine-scale convective cells are embedded within and modulated by larger-scale precipitation systems.

Combating Blur: AANet and Adversarial Autoregressive Models

A systematic failure mode of probabilistic nowcasting models - whether GAN-based or diffusion-based - is blur: when averaged over ensemble members, the probabilistic forecast produces diffuse, low-intensity precipitation fields that lack sharp boundaries. At the ensemble mean level, this is expected (individual members may be sharp but placed at different locations), but when the evaluation focuses on the ensemble mean rather than individual members, blurry forecasts score poorly on standard skill metrics.

AANet (Adversarial Autoregressive Network) addresses this through a deterministic autoregressive architecture augmented with an adversarial loss. The model generates the nowcast sequence autoregressively: at each step, the network receives the current predicted frame and the original context, predicts the next frame, and adds it to the context. The adversarial discriminator operates on the spatial gradients of the predicted frames, penalising blurry outputs at the pixel-edge level.

TrajGRU and the HKO-7 Benchmark

Before DGMR and NowcastNet, the standard benchmark for data-driven precipitation nowcasting was established by Shi et al. (2017) with the TrajGRU architecture and the HKO-7 dataset.

HKO-7 is the Hong Kong Observatory radar dataset, comprising 7 years (2009–2015) of radar composites at 1,km, 6-minute resolution over a $480\times 480$,km domain centred on Hong Kong. The dataset provides a controlled benchmark with diverse precipitation regimes (tropical cyclones, monsoon precipitation, frontal systems) and is publicly available, making it the standard evaluation platform for the community.

TrajGRU extends the convolutional LSTM (ConvLSTM) of Shi et al. (2015) by replacing the fixed convolution in the recurrent update with a trajectory-based warping: the state is advected along learned flow trajectories before the standard GRU update. This provides a structural inductive bias for advection without the explicit PDE constraint of NowcastNet.

TrajGRU remains a competitive deterministic baseline and is the standard against which generative methods are compared in the nowcasting literature.

Nowcasting Pipeline: TikZ Overview

The Value of Probabilistic Nowcasting

Why probabilistic rather than deterministic? The answer is clearest at extreme thresholds, where deterministic forecasts have near-zero skill but probabilistic forecasts can still provide decision-relevant information.

The implication for nowcasting is direct: even a probabilistic forecast that assigns only 10% probability to a flooding event provides actionable decision support for a civil protection agency with a cost-loss ratio of 0.05 (protective measures cost 5% of the damage from an unmitigated flood).

Comparison: Where Generative Models Excel

The key distinction between nowcasting applications and most other generative modelling tasks is that physical correctness is not optional. A generative model for images can hallucinate: the consequence is an aesthetically incorrect photograph. A generative nowcasting model that hallucinates a convective cell that does not exist may trigger an unnecessary evacuation; one that fails to generate a cell that does exist may fail to warn a community at risk. Both errors have real costs. This asymmetry motivates ongoing work on value-aware training objectives that weight the loss according to the decision-theoretic cost of different error types, rather than treating all forecast errors equally.

Exercises

End-to-End Forecasting and Earth System Foundation Models

“The bottleneck is not the forecast model. It never has been. The bottleneck is getting observations into a coherent state of the atmosphere in the first place.” - Peter Bauer, ECMWF Deputy Director of Research

Every AI weather model described in preceding sections-Pangu, GraphCast, FourCastNet, Aardvark-shares a common dependency that is rarely discussed with the excitement reserved for benchmark scores: they all ingest analyses produced by classical data assimilation (DA) systems. These analyses are extraordinarily expensive. ECMWF's 4D-Var assimilation runs on the same supercomputer cluster that produces the ensemble forecasts, consuming roughly as much compute as the forecast itself. Achieving operational-grade analysis requires solving a vast non-convex optimisation problem over hundreds of millions of variables, constrained by the primitive equations and informed by millions of heterogeneous observations arriving from satellites, weather stations, radiosondes, ships, buoys, and commercial aircraft.

The dream of fully machine-learning-driven weather forecasting-where raw sensor data flows directly to predictions without ever touching a physics-based assimilation system-remained elusive until 2025. This section documents the models that broke that barrier and the Earth system foundation models that have begun unifying atmospheric, oceanic, and chemical forecasting under a single pre-training paradigm.

The NWP–DA Bottleneck

To appreciate why end-to-end ML weather prediction is difficult, we must first quantify the DA problem precisely.

The computational cost of solving is dominated by repeated evaluations of the tangent-linear and adjoint of the full nonlinear model $ℳ$. For ECMWF's IFS model, a single 4D-Var analysis requires approximately 60–80 iterations of an inner-loop minimiser, each involving a forward and adjoint model integration. The total wall-clock time is $\sim$30 minutes on a cluster with ${10}^4$ cores-purely to produce the initial condition for a forecast that itself takes only a few minutes to run on AI hardware.

This asymmetry is the NWP–DA bottleneck. The forecast step has been accelerated by two to three orders of magnitude by AI models; the analysis step has not moved at all.

Aardvark: End-to-End Data-Driven Weather Forecasting

The landmark that demonstrated end-to-end ML weather forecasting was Aardvark (Vaughan et al., 2025), published in Nature. Aardvark is the first system to ingest raw, heterogeneous meteorological observations-without any prior data assimilation step-and produce competitive medium-range forecasts of global atmospheric fields and point-specific station forecasts simultaneously.

Architecture overview

The Aardvark pipeline consists of four stages: a SetConv encoder that aggregates heterogeneous observations into a latent grid, a CNN-based latent state processor that propagates the state forward in time, a gridded global decoder for standard NWP output, and a point-specific decoder for station-level predictions.

The SetConv encoder

The core of Aardvark's ability to handle irregular observations is the Set Convolution (SetConv) encoder. Given a set of observation locations ${\{({\ell }_i,{𝐨}_i)\}}_{i=1}^M$, where ${\ell }_i\in {ℝ}^2$ is latitude–longitude position and ${𝐨}_i\in {ℝ}^{d_{\mathrm{obs}}}$ is the observation value vector, SetConv produces a feature map on a regular lat–lon grid by kernel-weighted aggregation.

The SetConv operation is permutation-invariant in the observation index, translation-equivariant under shifts of the grid, and handles $M=0$ (missing observations at a grid point) gracefully by outputting the background climatology embedding. These properties make it robust to the irregular and time-varying nature of real observing networks.

CNN latent state and autoregressive rollout

Once the SetConv encoder has projected observations onto the regular grid, Aardvark uses a convolutional latent state processor inspired by the U-Net architecture. The latent state ${𝐳}_t\in {ℝ}^{H\times W\times D}$ encodes a compressed representation of the atmospheric state, and a learned transition function $f_{\theta }:{𝐳}_t\mapsto {𝐳}_{t+\mathrm{\Delta }t}$ rolls it forward autoregressively.

Dual-decoder design

Aardvark's ability to produce both gridded global forecasts and point-specific station forecasts from a single latent state is enabled by its dual-decoder architecture. The gridded decoder applies a convolutional upsampler to ${𝐳}_{t+k\mathrm{\Delta }t}$ to produce standard NWP output fields (temperature, wind, geopotential, humidity) at each lead time. The station decoder uses cross-attention between station coordinates and the latent grid to produce calibrated predictions at irregular station locations without interpolation artefacts.

Computational efficiency

One of Aardvark's headline results is its computational efficiency relative to classical NWP. The authors demonstrate two orders of magnitude improvement in data assimilation efficiency:

The skill reduction is approximately 5–10% in RMSE for 500 hPa geopotential at 5-day lead, relative to the operational ECMWF system. However, for surface station variables, Aardvark's point-specific decoder actually outperforms interpolated gridded fields from the operational system, because it learns to account for station micro-climatology directly.

Aurora: An Earth System Foundation Model

While Aardvark solved the data assimilation bottleneck, Aurora (Bodnar et al., 2025, Microsoft Research) attacked a different limitation: the proliferation of specialised models. As of 2024, high-quality AI models existed separately for atmospheric weather (Pangu, GraphCast), ocean state (HYCOM-derived neural operators), air quality (separate CNN regressors), and ocean wave height (still classical). Each model required its own training pipeline, its own data curation, and its own validation infrastructure.

Aurora's thesis is that all of these Earth system prediction tasks can be unified under a single pre-training objective, with task-specific fine-tuning achieving state-of-the-art performance on each domain.

Pre-training on 1 million hours of Earth system data

Aurora is pre-trained on over one million hours of Earth system observations and model analyses:

Data source	Hours	Variables	Description
ERA5 reanalysis	720,000	82	Atmospheric fields 1940–2023
MERRA-2 reanalysis	70,000	34	NASA multi-decadal atmosphere
GFS analyses	120,000	25	NOAA operational analyses
HYCOM ocean analyses	60,000	18	Ocean temperature, salinity, currents
CAMS air quality	30,000	14	Aerosols, O${}_3$, NO${}_2$, PM2.5
ERA5-ocean wave	40,000	8	Significant wave height, period
Total	1,040,000	181

Perceiver-based architecture

Aurora's core architectural innovation is the use of a Perceiver IO backbone to create a universal 3D latent representation of the Earth system.

The Perceiver design elegantly decouples the input resolution from the latent computation: inputs can grow (new satellite instruments, finer resolution) without changing the core model. This is the architectural property that makes Aurora truly foundation-model-like: new modalities are encoded by new cross-attention mappings without retraining the latent transformer.

Fine-tuning for diverse Earth system tasks

Aurora is pre-trained on the joint multi-variable prediction task (next-6-hour forecast across all modalities simultaneously). It is then fine-tuned on five downstream tasks:

1. Medium-range weather (1–10 days): global T, U, V, Z, Q at 0.1${}^{\circ }$ resolution. Aurora matches or exceeds GraphCast on all ERA5 benchmarks.

2. Air quality forecasting (1–5 days): PM2.5, NO${}_2$, O${}_3$, SO${}_2$ globally. Aurora delivers 5-day PM2.5 predictions at CAMS operational quality, trained on only 3,000 hours of fine-tuning data.

3. Ocean wave height: significant wave height globally at 12 h intervals. Aurora achieves RMSE comparable to NOAA operational WAVEWATCH III.

4. Tropical cyclone intensity: maximum wind speed and track predictions. Fine-tuned Aurora outperforms operational intensity models at 72 h lead by 12%.

5. Stratospheric ozone: column ozone and vertical profiles. This task benefits most from the joint pre-training, suggesting that learned tropospheric dynamics improve stratospheric representations.

FuXi Weather: ML Data Assimilation and Forecasting

Concurrent with Aardvark, the Chinese Meteorological Administration and Fudan University developed FuXi Weather, an extension of the FuXi cascade model (described in Section GraphCast and Graph Neural Networks) to include a learned data assimilation front-end. FuXi Weather uses a U-Net-style analysis network trained to minimise the error between its output analyses and ERA5 reanalysis, conditioned on observations from AMSU-A, IASI, and ATMS satellite sounders. The approach achieves analysis quality within 5% of the ERA5 RMSE benchmark for key upper-tropospheric variables, while running 800$\times$ faster than the ECMWF 4D-Var.

A critical design choice in FuXi Weather is the use of observation simulation during training: because the true atmospheric state is unknown, the model is trained on synthetic observations generated by applying radiative transfer models to ERA5 fields. This avoids overfitting to specific instrument calibration biases while learning the general physics of atmospheric sounding.

MoTiF: Multi-Source Tropical Cyclone Forecasting

MoTiF (Multi-source Observation Tropical cyclone Integrated Forecasting) addresses the specialised problem of tropical cyclone (TC) prediction from heterogeneous, irregular observation sources. Unlike global models that treat TCs as part of a larger atmospheric state, MoTiF uses a self-supervised multi-source encoder that learns a unified TC vortex representation from:

• Microwave satellite imagery (SSMI/S, AMSR2)

• Infrared geostationary images (GOES-16, Himawari)

• Aircraft reconnaissance observations (Vortex Data Messages)

• Scatterometer surface winds

• Best-track historical records

The self-supervised pre-training objective is masked observation modelling: given a subset of the above sources with some randomly masked, predict the masked observations. This forces the model to learn cross-modal TC structure. Fine-tuned MoTiF achieves TC intensity forecasting errors (in maximum sustained winds) 15% lower than the operational HWRF-B ensemble at 48 h lead time.

Cloud4D: High-Resolution Cloud Property Estimation

Cloud4D extends the end-to-end paradigm to the estimation of three-dimensional cloud microphysical properties from passive satellite imagery. This is an inverse problem: given 2D top-of-atmosphere radiances at multiple wavelengths, infer the full 3D distribution of cloud droplet effective radius, liquid water content, and cloud-top pressure. Classical retrieval algorithms are based on look-up tables derived from radiative transfer simulations and cannot propagate uncertainty.

Cloud4D uses a conditional diffusion model to sample from the posterior $p(\text{3D cloud structure}|\text{satellite radiances})$, trained on high-resolution large-eddy simulations paired with synthetic satellite observations. The model captures the full uncertainty in cloud retrievals and produces an ensemble of 3D cloud fields consistent with the observations.

System	Pub.	Architecture	DA cost	Tasks	Key innovation
Aardvark	2025	SetConv+CNN	$\sim$${10}^{-4}$	Weather+stations	First full e2e; runs on desktop
FuXi Weather	2025	U-Net analysis	$\sim$${10}^{-3}$	Weather	Synthetic obs training
Aurora	2025	Perceiver IO	N/A (pre-trained)	Atm+ocean+chem	Earth system foundation model
MoTiF	2025	Set encoder	$\sim$${10}^{-3}$	Tropical cyclones	Multi-source TC fusion
Cloud4D	2025	Diffusion	N/A	Cloud properties	3D posterior sampling

Exercises

Operational Deployment and Real-World Impact

“We've been saying for years that AI would transform weather forecasting. In 2025, we stopped saying it and started doing it.” - Rick Spinrad, NOAA Administrator, December 2025

The history of numerical weather prediction can be divided into three eras. The first era, from Richardson's 1922 hand-computed forecast to the first machine-computed prediction in 1950, established the mathematical framework. The second era, from 1950 to roughly 2021, was characterised by steady incremental improvement: better models, finer grids, more observations, improved parameterisations. The third era began in 2022 with the first competitive AI forecast models and became an operational revolution by 2024–2026.

This section documents that revolution: the operational deployment of AI weather systems at national meteorological agencies, the commercial sector's rapid adoption, and the measurable impact on human welfare-from hurricane evacuation decisions to agricultural productivity in the developing world.

NOAA Project EAGLE: The US Operational AI Suite

The United States National Oceanic and Atmospheric Administration (NOAA) launched Project EAGLE (Ensemble AI Guidance for Linked Environments) in December 2025, representing the most comprehensive national operational AI weather deployment to date. Project EAGLE comprises three distinct AI-based forecast products:

AIGFS: AI-based Global Forecast System

The AIGFS is a fine-tuned version of GraphCast, trained on NOAA's GDAS (Global Data Assimilation System) analyses from 1990–2024. Key operational specifications:

• Lead time: 0–16 days deterministic forecasts.

• Resolution: ${0.25}^{\circ }$ horizontal, 13 pressure levels.

• Run time: 40 minutes for a complete 16-day forecast, versus $\sim$6 hours for the operational GFS on the same hardware budget.

• Skill: Matches or exceeds operational GFS RMSE for 500 hPa geopotential at Days 1–7; small deficit at Days 8–16.

• Update frequency: 4 times daily (00, 06, 12, 18Z).

The 40-minute wall-clock time enables NOAA to run AIGFS on a single GPU cluster that would previously have been used for one 6-hour model run, leaving the primary supercomputer for ensemble tasks.

AIGEFS: AI-based Global Ensemble Forecast System

The AIGEFS generates a 50-member probabilistic ensemble using a lagged-perturbation technique: initial condition perturbations from the GDAS ensemble Kalman filter are fed into a fine-tuned Aurora model.

• Lead time: 0–24 days probabilistic forecasts.

• Skill extension: AIGEFS shows measurable probabilistic skill at Days 18–24, extending the current operational GEFS skill by 6–9 days. This is equivalent to the improvement achieved by the NWP community over the preceding 15 years.

• Compute cost: AIGEFS 50-member ensemble costs 9% of the equivalent GEFS run on the same hardware.

• Key application: Subseasonal prediction (2–4 weeks), where the probabilistic skill extension is most valuable for agriculture, water management, and energy demand planning.

HGEFS: Hybrid Grand Ensemble Forecast System

The HGEFS is the most ambitious component of Project EAGLE: a hybrid ensemble that combines 31 AI-based members with 30 physics-based members from the operational GEFS, for a total 61-member grand ensemble.

Member type	Count	Model	Cost fraction	Strength
AI deterministic	1	AIGFS	1%	Low-cost high-resolution det.
AI ensemble	30	AIGEFS	8%	Extended probabilistic range
Physics ensemble	30	GEFS	91%	Physical consistency, extremes
Total	61	HGEFS	100%

The HGEFS design reflects a key operational philosophy: AI members are cheap and skillful for most conditions but unreliable for extremes; physics-based members are expensive but physically consistent. The hybrid ensemble covers both regimes.

ECMWF AIFS and the Anemoi Framework

The European Centre for Medium-Range Weather Forecasts has been the world's leading NWP centre for four decades. Its operational adoption of AI marks the most significant institutional validation of the technology.

AIFS (Artificial Intelligence Integrated Forecasting System) became the first AI model in ECMWF's operational forecast suite in early 2024. AIFS is based on a graph neural network architecture trained on 40 years of ERA5 data. Its key operational role is as an ensemble member in the ECMWF ENS rather than as a replacement for the physics-based IFS: it contributes ensemble diversity at negligible cost, improving tail probability estimates without replacing the physically consistent IFS ensemble.

Anemoi: Open-source ML framework

Recognising that AI weather development should not be confined to a handful of large technology companies, ECMWF released Anemoi (Greek: “winds of the north”), an open-source framework for training AI weather models on arbitrary ERA5-derived datasets. Anemoi provides:

• anemoi-datasets: Standardised ERA5 data pipelines with configurable variable sets, resolutions, and temporal splits.

• anemoi-models: A library of graph-based and transformer-based weather model architectures.

• anemoi-training: Distributed training recipes for models from 100M to 10B parameters.

• anemoi-inference: Inference engine for operational deployment with checksum verification and anomaly detection.

Anemoi has enabled a pan-European collaboration among national meteorological services. By mid-2025, the following regional models had been trained using the Anemoi framework:

• Bris (MET Norway): Limited-area AI model for the North Sea and Scandinavia at 2.5 km resolution.

• MeteoSwiss AI: Alpine convection-resolving AI model leveraging high-resolution Swiss radar and radiosonde data.

• FMI-Aurora (Finnish Meteorological Institute): Polar-optimised AI model for Arctic weather prediction.

• DWD-AIFS (Deutscher Wetterdienst): Central European mesoscale model at 1 km resolution for fog and frost prediction.

These regional models demonstrate a key advantage of AI weather: the marginal cost of training a high-resolution limited-area model is dramatically lower than developing an equivalent NWP limited-area system, enabling smaller national services to compete at global accuracy levels.

China: The FuXi–FengWu–Fenghe Ecosystem

China's meteorological AI programme, led by the China Meteorological Administration (CMA) in collaboration with Shanghai AI Laboratory and multiple universities, represents the most integrated national AI weather ecosystem.

FuXi (described in Section GraphCast and Graph Neural Networks) and FengWu are both operational at CMA for medium-range forecasting. NowcastNet provides deep-learning nowcasting (0–6 hour precipitation forecasting) for severe weather warnings.

The most novel component of the Chinese ecosystem is Fenghe (meaning “wind harmony”), the world's first operational meteorological large language model. Fenghe is not a forecast model; it is a risk communication and decision support system. Its key capabilities include:

• Localised risk narratives: Given a set of AI forecast fields, Fenghe generates location-specific natural language risk summaries (e.g., “Flooding risk in Guangzhou will peak between 14:00 and 18:00 local time, with highest risk in Haizhu district”).

• Multi-agency coordination: Fenghe ingests forecast data from multiple agencies (CMA, ECMWF, GFS) and generates consensus language describing areas of forecast uncertainty.

• Public alert generation: Fenghe directly drafts emergency weather alerts in plain language, which are reviewed by human meteorologists before dissemination.

• Agricultural advisory: Daily crop-specific weather advisories for 31 provinces, integrating phenological models with AI weather forecasts.

Fenghe represents the full stack of AI weather deployment: from observation to forecast to risk-based communication to sector-specific decision support.

The Commercial Sector: Horizon AI, Climavision, and Earth-2

The commercial weather sector has adopted AI forecasting with a speed that has surprised even optimistic observers.

Horizon AI and the 2024 Hurricane Season

Horizon AI, a start-up founded by former NOAA and NCAR meteorologists, deployed an ensemble of fine-tuned Pangu and GraphCast models for the 2024 Atlantic hurricane season. The 2024 season was the most active on record (27 named storms, 7 major hurricanes), and provided an unprecedented operational test.

For Hurricane Milton (October 2024), which underwent rapid intensification from Category 1 to Category 5 in under 24 hours before striking Florida's Gulf coast, Horizon AI's models provided the following verified performance:

AXA, NVIDIA Earth-2, and Insurance Risk Assessment

The NVIDIA Earth-2 climate digital twin platform, announced in 2024, is a cloud-based system that allows enterprises to run high-resolution AI weather and climate simulations without owning supercomputing infrastructure. AXA, one of the world's largest insurance groups, became an early adopter for extreme weather risk assessment.

The AXA–Earth-2 partnership uses the following workflow:

1. Earth-2 runs thousands of AI-driven climate simulations to estimate the return periods and severities of extreme events (floods, windstorms, heat waves) at $1\mathrm{km}$ resolution globally.

2. AXA ingests the resulting return-period maps into its catastrophe modelling framework, updating expected loss distributions for insured assets.

3. Operational daily weather forecasts from Earth-2 are used to trigger parametric insurance products (automatic payouts when a forecast wind speed threshold is exceeded at a specified location).

This use case illustrates a key driver of commercial AI weather adoption: the ability to run thousands of simulations is economically transformative for risk industries where pricing depends on tail probability estimates that classical NWP could never provide at scale.

Global Map of Operational Deployments

Case Study: AI Weather for 38 Million Indian Farmers

Nowhere is the real-world impact of improved weather forecasting more tangible than in subsistence agriculture. India has 140 million farming households, most of them rain-dependent smallholders whose crop choices, planting dates, and pest management decisions are intimately tied to weather forecasts. A 24-hour forecast that is 10% more accurate can mean the difference between a protected harvest and a washed-out crop.

The India Meteorological Department (IMD), in partnership with Microsoft Research and ICRISAT (the International Crops Research Institute for the Semi-Arid Tropics), deployed an AI forecast chain in 2025 targeting 38 million farmers in the semi-arid tropics of central and southern India. The system chain is:

1. Global AI forecast: Aurora fine-tuned on IMD analyses provides 1–10 day forecasts at ${0.1}^{\circ }$ resolution.

2. Statistical downscaling: A learned spatial disaggregation model maps ${0.1}^{\circ }$ fields to 5 km resolution, trained on historical IMD station data and AWS (automatic weather stations) installed by the ICRISAT network.

3. Crop advisory generation: A fine-tuned language model (analogous to Fenghe) generates crop-specific advisories in 8 Indian regional languages, tailored to the farmer's crop type, growth stage, and local soil conditions (from a national soil database).

4. Delivery: Advisories are delivered via SMS to feature phones (because 60% of target farmers have no smartphone) and via a voice bot accessible through a free toll number.

An independent evaluation in the 2025 kharif (summer) season showed:

• Uptake: 23 million farmers (60% of target) used the advisory service at least once during the season.

• Crop loss prevention: Participating farmers reported 12% lower crop losses due to unseasonal rainfall compared to non-participating farmers in the same villages.

• Input efficiency: Fertiliser use was 8% lower (through better timing of application relative to rain events), reducing both costs and nitrogen leaching.

• Economic benefit: Estimated INR 28,000 crore (approximately $3.3 billion USD) in prevented crop losses across the season.

This case study illustrates that the economic and humanitarian value of AI weather forecasting is not confined to the advanced industrialised world. The marginal cost of deploying an additional AI-based national weather service is now small enough that developing countries can access forecast quality that was previously available only to the wealthiest nations.

The Economics of AI Weather Forecasting

The economic case for AI weather forecasting rests on three pillars: compute cost reduction, forecast quality improvement, and new capability enablement.

Compute cost reduction

The NOAA Project EAGLE documentation provides the most detailed public cost breakdown. The AIGEFS 50-member ensemble costs 9% of the equivalent GEFS ensemble on the same hardware. Extrapolating:

• NOAA's annual supercomputing budget for weather modelling was approximately $290 million in 2024.

• If 50% of ensemble compute is replaced by AI equivalents at 9% cost, the savings are $\approx \$0.5\times \$290M\times (1-0.09)\approx \$132M$ per year.

• The global NWP community spends approximately $2–3 billion annually on supercomputing. Equivalent AI transition could save $900M–$1.4B per year in compute alone.

Forecast quality improvement value

Independent economic analyses estimate the value of a 1-day extension in useful forecast lead time at approximately $1B–$5B globally, accounting for improvements in aviation routing, energy demand management, agricultural decisions, and emergency management. The 6–9 day skill extension documented for AIGEFS thus represents an estimated economic value of $6B–$45B annually, though these estimates carry large uncertainties.

New capability enablement

The most economically significant new capability is the ability to run large AI ensembles for rare event risk estimation. Running 10,000 AI ensemble members to estimate the probability of a 1-in-100-year flood event costs $\sim$1% of the equivalent physics-based computation. This opens catastrophe modelling, insurance pricing, and long-range risk assessment to AI methods that were previously computationally infeasible.

Exercises

Extremes, Predictability, and Open Problems

The remarkable skill curves celebrated by the machine learning weather prediction community-FourCastNet's five-day 500,hPa geopotential accuracy surpassing the ECMWF medium-range forecast, GraphCast's ten-day headline metrics, GenCast's ensemble spread on the WeatherBench 2 leaderboard-all rest on a foundational assumption: that the metrics being optimized and reported faithfully capture operational value. They do not. Not completely. This section confronts the gap between the leaderboard and the reality faced by forecasters, emergency managers, and the public they serve. That gap is widest precisely where the stakes are highest: in the prediction of extreme events.

The structure of this section is as follows. We first diagnose the performance paradox, the systematic divergence between benchmark metrics and operational skill for high-impact events. We then examine the statistical and physical mechanisms that drive the intensity underprediction pathology common to every current learned weather model. We survey the 2024–2025 hurricane seasons as a natural experiment, examining where track prediction succeeded and where intensity forecasts failed catastrophically. Aurora's comprehensive evaluation across multiple extreme regimes provides a structured case study. We then state eight open problems that, in our assessment, define the research frontier for the next decade. The section closes with a conjecture on the ultimate limits of monolithic data-driven weather modeling and an argument for the hybrid future.

The Performance Paradox: Leaderboard Metrics versus Operational Reality

To make this precise, let ${\widehat{𝐱}}_{t+h}$ denote the forecast issued at time $t$ for lead time $h$, and let ${𝐱}_{t+h}^{\ast }$ denote the verifying analysis. The standard aggregate RMSE is (RMSE Aggregate)$$\begin{array}{ccc}& \mathrm{RMSE}(h)=\sqrt{\frac{1}{|𝒱|}\sum _{(t,\ell )\in 𝒱}\Vert {\widehat{𝐱}}_{t+h}^{(\ell )}-{𝐱}_{t+h}^{(\ell )\ast }{\Vert }^2},& \text{<span class="tml-eqn"></span>}\end{array}$$ where $𝒱$ is the validation set of (time, location) pairs and $\ell$ indexes pressure levels and variables. Now partition the validation set into the “routine” set ${𝒱}_{\mathrm{r}}$ and the “extreme” set ${𝒱}_{\mathrm{e}}$, where ${𝒱}_{\mathrm{e}}$ consists of cases in the top percentile of some impact metric (wind speed, precipitation accumulation, temperature anomaly). Since $|{𝒱}_{\mathrm{e}}|\ll |{𝒱}_{\mathrm{r}}|$, we have (RMSE Decomposed)$$\begin{array}{ccc}& \mathrm{RMSE}(h)\approx {\mathrm{RMSE}}_{\mathrm{r}}(h)+\frac{|{𝒱}_{\mathrm{e}}|}{|𝒱|}({\mathrm{RMSE}}_{\mathrm{e}}(h)-{\mathrm{RMSE}}_{\mathrm{r}}(h)).& \text{<span class="tml-eqn"></span>}\end{array}$$ If $|{𝒱}_{\mathrm{e}}|/|𝒱|\approx 0.01$ (top 1%), then even a model with ${\mathrm{RMSE}}_{\mathrm{e}}\gg {\mathrm{RMSE}}_{\mathrm{r}}$ will show only a 1% degradation in the headline aggregate metric. A model optimized to minimize the aggregate loss has essentially no gradient signal from the extreme events that dominate societal impact.

The aggregate metric bias manifests in a second, more subtle way: model selection. When models are trained and selected using aggregate validation loss, the optimization pressure on rare extremes is quadratically negligible. A model that slightly smooths its predictions (reducing high-frequency variance) achieves lower aggregate RMSE at the cost of substantially degraded extreme event representation. This is the ERA5 smoothing inheritance problem, which we examine next.

ERA5 Smoothing Inheritance: Learning What You Train On

Every major MLWP model trained through 2025 uses ERA5 as its primary training target. ERA5 is a gridded reanalysis product produced by the ECMWF at approximately 31,km horizontal resolution (${0.28125}^{\circ }$ latitude–longitude), generated by assimilating billions of observations into a 4D-Var cycling system. As a reanalysis, ERA5 represents the atmosphere's best estimate given available observations, but it is not equivalent to the true atmosphere.

The spectral structure of ERA5 differs from the true atmosphere in a systematic way that has severe consequences for extremes. Let ${\hat{X}}_k$ denote the spherical harmonic coefficient of ERA5 at total wavenumber $k$, and let ${\hat{X}}_k^{\ast }$ denote the true atmospheric coefficient. The effective spectral transfer function (Spectral Transfer)$$\begin{array}{ccc}& H(k)=\frac{𝔼[|{\hat{X}}_k{|}^2]}{𝔼[|{\hat{X}}_k^{\ast }{|}^2]}& \text{<span class="tml-eqn"></span>}\end{array}$$ satisfies $H(k)\approx 1$ for $k\lesssim 100$ but $H(k)<1$ for $k\gtrsim 200$, with $H(k)\to 0$ as $k\to k_{\max }$ (the Nyquist wavenumber of the ERA5 grid). In physical space, this spectral damping implies that ERA5 systematically underrepresents spatial gradients, peak wind speeds, and precipitation intensities at sub-gridscale and near-gridscale spatial frequencies.

A model trained to minimize mean-squared error against ERA5 therefore learns the ERA5 spectral distribution, including its damping at high wavenumbers. At inference time, even if the model produces physically consistent atmospheric states at coarse scales, it inherits the spectral bias of its training target. Kim et al. (2026) formalize this spectral inheritance with a rigorous test:

Kim et al. (2026) also introduce a butterfly effect test for GenCast, measuring whether ensemble spread grows at rates consistent with the Lorenz attractor's Lyapunov spectrum. They find that GenCast's spread growth rate is qualitatively correct for the first seven days but overestimates predictability (underestimates spread) at lead times beyond ten days, consistent with the spectral inheritance hypothesis: the model's learned dynamics have fewer short-scale instabilities than the true atmosphere, resulting in artificially low ensemble divergence at extended range. Physical consistency tests for mass, momentum, and energy budgets reveal compensating biases that partially cancel in zonal-mean diagnostics but accumulate regionally and seasonally.

The Probabilistic Benchmark: Brenowitz et al. (2025)

The evaluation of probabilistic weather forecasts has a long history in the operational forecasting community, rooted in the proper scoring rule literature (Gneiting and Raftery, 2007) and the ensemble verification tradition of the World Meteorological Organization. Brenowitz et al. (2025) translate this tradition into a rigorous benchmark for MLWP systems, addressing a systematic gap: most published MLWP results report either deterministic skill metrics (RMSE, ACC) or aggregate probabilistic metrics (mean CRPS) without separating the climatological reliability, sharpness, and extreme tail performance.

The Brenowitz benchmark decomposes forecast quality along four orthogonal axes:

1. Reliability: Is the stated probability calibrated? For an ensemble forecast issuing probability $p$ for event $A$, the reliability diagram plots $\hat{p}$ (forecast frequency) against $\bar{p}$ (observed relative frequency). Perfect reliability requires $\hat{p}=\bar{p}$ for all $p$.

2. Resolution: Does the forecast discriminate between cases where $A$ occurs and cases where it does not? Measured by the area under the ROC curve (AUROC) or the Brier skill score.

3. Sharpness: How concentrated is the predictive distribution? A sharp but unreliable forecast is overconfident; sharpness is desirable only jointly with reliability.

4. Extreme skill: For threshold $u$ at the $q$-th climatological quantile, what is the threshold-weighted CRPS (Equation (Twcrps)) as a function of $q\in \{0.9,0.95,0.99,0.999\}$?

A model cannot game all four axes simultaneously. Ensemble mean smoothing improves reliability at the cost of sharpness; overfit ensembles achieve sharpness at the cost of reliability; models that underrepresent variance (due to spectral inheritance) fail on both sharpness and extreme skill. The Brenowitz benchmark reveals that every published MLWP model as of 2025 exhibits measurable reliability degradation beyond day 5, sharpness collapse beyond day 8, and extreme skill that is significantly below the ECMWF ensemble at the 99th percentile threshold for every lead time.

Hurricane Case Studies: Helene, Milton, Beryl, and Lee

The 2024 Atlantic hurricane season provided an exceptionally rich natural experiment for evaluating MLWP systems under operational conditions. Four storms-Helene, Milton, Beryl, and Lee-exhibited contrasting behaviors that stress-test different aspects of learned models.

The Track–Intensity Divergence

The most consistent finding across all four storms is a track–intensity divergence: MLWP models predict track with skill comparable to or exceeding operational NWP, but systematically underpredict peak intensity. This is not merely a resolution effect (though resolution contributes); it reflects the spectral inheritance problem at its most acute.

Track prediction depends primarily on the large-scale steering flow, which is well-represented in ERA5 and well-captured by data-driven dynamics. Intensity prediction depends on inner-core structure at spatial scales of 10–50,km (the eye wall, the radius of maximum wind, the warm core anomaly) that are partially resolved in ERA5 but strongly smoothed. A model inheriting ERA5 spectral statistics learns a smoothed inner core, and therefore cannot reproduce rapid intensification events.

Hurricane Helene (2024): Inland Flooding and Track Success

Helene made landfall in the Florida Big Bend region on 26 September 2024 as a Category 4 hurricane. MLWP models provided 5–7 day track forecasts with mean track errors below 150,km, competitive with operational NWP. The storm's subsequent inland penetration and extreme rainfall event over western North Carolina (200–500,mm accumulations triggering catastrophic flooding) was associated with the interaction of Helene's remnant circulation with complex Appalachian orography. MLWP models underperformed NWP on the rainfall distribution, consistent with their inability to resolve orographic lifting at scales below their effective resolution. The societal impact, over 200 fatalities, was dominated by the inland flood event rather than the coastal wind, a regime where current MLWP skill is lowest.

Hurricane Milton (2024): Rapid Intensification Failure

Milton's rapid intensification from Category 1 to Category 5 in approximately 24 hours (5–7 October 2024) remains one of the most dramatic RI events in Atlantic basin history. MLWP models available for operational guidance at 72-hour lead time forecast a Category 2 or 3 intensity at landfall; the observed peak intensity of 180,kt ($92.5{\mathrm{m}\mathrm{s}}^{-1}$) at sea exceeded all MLWP forecasts by 50–90,kt. This represents a mean absolute error in intensity that is three to five times larger than the track error scaled equivalently. Analysis of the spectral power in Milton's inner-core fields in ERA5 confirms severe smoothing: the ERA5 representation of Milton's eye wall had a warm core anomaly approximately 30% weaker than aircraft reconnaissance observations. A model learning from ERA5-level data was presented with a structurally weakened target; it cannot learn what it was not shown.

Hurricane Beryl (2024): Extended Range Track Success

Beryl formed unusually early in the season (1 July 2024) and struck the Windward Islands, Mexico, and Texas in succession. Its track was predicted with skill at 7–10 day lead times by multiple MLWP systems, providing a proof of concept for the value of machine learning at extended range for the large-scale steering flow. However, peak intensity at each landfall was underpredicted by 15–25,kt, and the rapid weakening over the Yucatan Peninsula followed by reintensification over the Gulf was captured only qualitatively. The track success of Beryl attracted favorable public and policy attention to MLWP, while the intensity errors (which matter for warning decisions) received less media coverage.

Hurricane Lee (2023-24 Analysis): The Non-Landfall Benchmark

Lee (2023) provided a clean benchmark because it underwent explosive RI to Category 5 over the open Atlantic and then recurved without making a direct landfall on a major populated area. Post-season analysis using MLWP hindcasts reveals that the Category 5 peak intensity was forecast by no MLWP system, with all models predicting Category 3 at peak; the recurvature track, however, was accurately forecast at 7-day lead time. Lee thereby exemplifies the core asymmetry: machine learning has solved the where of tropical cyclone motion far better than the how intense.

Aurora (2026) Extreme Event Evaluation

Aurora (Microsoft, 2026) represents the most comprehensive extreme event evaluation published by any MLWP team, covering five regimes: tropical cyclones, atmospheric rivers, continental heatwaves, deep freeze events, and compound flooding. The evaluation methodology follows Brenowitz et al. (2025) in using threshold-weighted scoring, but extends it with conditional skill scores that separate environment-driven predictability from inner-core dynamics.

Atmospheric Rivers

Atmospheric rivers (ARs) are narrow corridors of intense meridional water vapor transport ($\mathrm{IVT}\ge 250\mathrm{k}\mathrm{g}{\mathrm{m}}^{-1}{\mathrm{s}}^{-1}$) responsible for the majority of extreme precipitation on the US West Coast, western Europe, and other mid-latitude continental margins. The integrated vapor transport (IVT) field (IVT)$$\begin{array}{ccc}& \mathrm{IVT}=\frac{1}{g}{\int }_{p_{\mathrm{sfc}}}^{300\mathrm{hPa}}q|𝐮|\mathrm{d}p,& \text{<span class="tml-eqn"></span>}\end{array}$$ where $q$ is specific humidity and $𝐮$ is the horizontal wind vector, depends on column-integrated dynamics at scales of 200–500,km. Aurora predicts AR location (defined by the axis of maximum IVT) with track errors below 50,km at 5-day lead time and below 100,km at 10-day lead time, exceeding operational ensemble skill for the former and matching it for the latter. Intensity (peak IVT magnitude) is again underestimated: Aurora's median IVT magnitude is 15–25% below observations at peak.

Continental Heatwaves

The 2021 Pacific Northwest heatwave (peak ${\mathrm{T}}_{2\mathrm{m}}>{49}^{\circ }\mathrm{C}$) and the 2023 western European heat events were used as retrospective test cases. Aurora (2026) hindcasts initialized 5, 7, and 10 days before peak heat correctly predicted the sign of the anomaly in all cases, but underestimated peak temperatures by 3–6${}^{\circ }$C at 10-day lead time. The spatial footprint of the heat dome was well-predicted. This 3–6${}^{\circ }$C bias corresponds to the difference between a record-breaking heat event and a catastrophic one: it matters enormously for health and infrastructure decisions.

Deep Freeze Events

Sudden stratospheric warming (SSW) events, in which the polar vortex weakens and cold air descends to mid-latitudes, present a different challenge: they are initiated at stratospheric levels (10–50,hPa) that are only partially represented in surface-focused MLWP architectures. Aurora's vertical extent to 0.1,hPa gives it an advantage over models truncated at 50,hPa. On the February 2021 Texas freeze (an SSW-driven event), Aurora (2026) hindcasts predicted the broad temperature anomaly pattern at 7-day lead, with peak cold underestimation of 4–8${}^{\circ }$C, again a systematic warm bias in the extreme tail.

Open Problems in Machine Learning Weather Prediction

We now state eight open problems that define the research frontier. These are not merely engineering challenges; they represent deep scientific questions about the limits of data-driven atmospheric modeling.

The Dominance Conjecture

This conjecture is consistent with the empirical record: NowcastNet specializes on radar nowcasting and underperforms on synoptic forecasting; GraphCast excels at 1–10 day prediction but has no sub-seasonal capability; no MLWP model has been demonstrated to produce physically meaningful climate projections. The conjecture predicts that future MLWP development will organize around model families or foundation models with regime-specific fine-tuning, rather than a single universal weather model.

The Open Problem Frontier Map

Research Challenges

[ enhanced, breakable, colback=exrcoralbg!40, colframe=exrcoral!70, boxrule=0.6pt, arc=3pt, title= Research Challenge A: The Extreme Event Benchmark, coltitle=exrcoral!20!black, fonttitle=, top=10pt, bottom=8pt, left=10pt, right=10pt, before skip=12pt, after skip=12pt ] Specification. Construct a community benchmark for MLWP extreme event skill, addressing the following requirements:

1. Coverage: 10+ extreme event types (TC, AR, heatwave, deep freeze, flash flood, blizzard, derecho, tornado outbreak, extratropical cyclone bomb, dust storm) across all climate zones and seasons.

2. Evaluation: For each event type, define threshold-weighted CRPS at the 90th, 95th, and 99th climatological percentiles, separately for location, timing, and intensity.

3. Baseline: Compare against (a) ERA5 climatology, (b) operational ECMWF ENS, (c) NOAA GEFS, and (d) the best published MLWP model at each horizon.

4. Reproducibility: All events, climatologies, and verification code publicly released with DOI; all submitted systems evaluated on a held-out test set never released publicly.

Expected impact. A standardized extreme event benchmark would redirect research investment toward the highest-impact prediction problems and enable meaningful comparison across model families.

[ enhanced, breakable, colback=defbluebg, colframe=defblue!70, boxrule=0.6pt, arc=3pt, title= Research Challenge B: Enforcing Physical Invariants, coltitle=defblue!20!black, fonttitle=, top=10pt, bottom=8pt, left=10pt, right=10pt, before skip=12pt, after skip=12pt ] Specification. Develop an architectural framework that enforces the following physical invariants in any learned weather model:

1. Mass conservation: ${\partial }_t\rho +\nabla \cdot (\rho 𝐮)=0$ at every grid point and time step.

2. Energy conservation (for adiabatic flow): $\frac{\mathrm{d}}{\mathrm{d}t}(KE+PE+IE)=0$ where $KE=\frac{1}{2}|𝐮{|}^2$, $PE=gz$, $IE=c_{v}T$.

3. Potential vorticity conservation (Ertel's PV): $\frac{\mathrm{D}q}{\mathrm{D}t}=0$ for $q=\frac{1}{\rho }(\omega +2\mathrm{\Omega })\cdot \nabla \theta$.

Benchmark the architecture on: (a) 10-day global drift in total energy, (b) 30-day total water vapor budget closure, (c) zonal-mean PV redistribution over a 100-day simulation. Expected impact. Conservation-enforcing architectures will extend the range of stable autoregressive integration and reduce the systematic climate drift that currently limits MLWP to operational weather timescales.

[ enhanced, breakable, colback=thmgreenbg, colframe=thmgreen!70, boxrule=0.6pt, arc=3pt, title= Research Challenge C: The Hybrid Operational System, coltitle=thmgreen!20!black, fonttitle=, top=10pt, bottom=8pt, left=10pt, right=10pt, before skip=12pt, after skip=12pt ] Specification. Design and prototype a hybrid MLWP plus physics-constrained NWP system with the following architecture:

1. Initialization: Neural data assimilation system producing initial state ${𝐱}_0^{\ast }$ with observation-constrained uncertainty.

2. Short range (0–3 days): High-resolution neural model (25,km) generating $N_s=100$ ensemble members at 6-hourly intervals.

3. Medium range (3–10 days): Reduced-resolution neural model (50,km) with physics-based correction at each 24-hour step to enforce mass and energy closure.

4. Extended range (10–30 days): Statistical post-processing model that blends neural ensemble climatological anomalies with boundary-forced S2S signals.

The system should be evaluated on one full year of quasi-operational retrospective forecasts with full verification suite. Expected impact. The hybrid architecture directly addresses the scale-dependent failure modes identified in the conjecture (Conjecture 1), leveraging neural speed for short range and physics for stability at extended range.

The Hybrid Future: Generative AI and Physics-Constrained NWP

The operational weather prediction systems of 2030 will almost certainly be neither pure NWP nor pure MLWP. The evidence from every operational evaluation conducted between 2022 and 2025 points toward a hybrid architecture that exploits the complementary strengths of the two paradigms.

The hybrid architecture operates as follows. The data assimilation cycle (currently 3D-Var or 4D-Var) is augmented with a neural background error covariance model, replacing the static $𝐁$ matrix with a flow-dependent neural estimate $\widehat{𝐁}({𝐱}_b)$ that adapts to the current synoptic pattern. The resulting analysis ${𝐱}_0^{\ast }$ is passed to both a conventional high-resolution NWP deterministic run and a neural ensemble generator. The neural ensemble provides rapid probabilistic characterization at lower cost than running 50 NWP members, while the deterministic NWP run provides physically consistent verification for ensemble quality control. At medium range, a generative downscaler (CorrDiff or similar) applies to the NWP output, producing high-resolution probabilistic fields without the computational cost of running the full NWP model at fine resolution. Throughout, post-processing neural networks (MOS successors) apply systematic bias corrections learned from the full verification history.

WeatherNext 2 (Google DeepMind, 2025–2026) represents the most advanced publicly described system in this hybrid direction, combining GraphCast dynamics with a diffusion-based ensemble generator and a learned post-processing layer, achieving operational certification in multiple national meteorological services. The system is end-to-end differentiable from observation to verification metric, enabling joint optimization of all components on operational skill scores.

Comprehensive Exercises and Research Challenges

This section provides twelve exercises of varying depth and five research challenges spanning the full mathematical and computational scope of this chapter. The exercises range from fundamental derivations (primitive equations, cost functions) through algorithmic computation (ConvLSTM forward pass, GraphCast message passing) to system-level analysis (operational deployment cost-benefit). The research challenges are open-ended projects suitable for semester-long independent study or early-stage research.

Exercises

Research Challenges

Chapter Summary

This chapter has traced the complete arc of machine learning for weather prediction, from the mathematical foundations of atmospheric dynamics through the current state of the art and to the open problems that define the next decade.