NeuralGCM: Hybrid Differentiable GCM

• NeuralGCM is a hybrid general circulation model that integrates a differentiable dynamical core with neural network closures to represent unresolved physics.
• It leverages online training with ERA5 data and stochastic parameterizations to enhance forecast skill, climate simulation, and ensemble calibration.
• The model bridges conventional dynamics and machine learning, offering improved physical consistency while highlighting challenges in land-surface coupling and stratospheric processes.

NeuralGCM, short for Neural General Circulation Model, is a hybrid atmospheric general circulation model that combines a differentiable dynamical core with neural-network parameterizations for unresolved physics. It was introduced as the first fully differentiable hybrid GCM and has since been evaluated across deterministic weather forecasting, probabilistic ensemble forecasting, multi-decade atmosphere-only climate simulation, precipitation modeling, seasonal prediction, climate-change stress tests, and differentiable storyline construction for extremes (Kochkov et al., 2023). In the subsequent literature, NeuralGCM is treated both as a concrete model family and as a reference architecture within the broader class of deep-learning Earth system models, with particular emphasis on the interaction between explicit large-scale dynamics, learned closures, stability, and out-of-distribution generalization (Pallotta et al., 8 Apr 2026, Baxter et al., 6 Oct 2025, Zhang et al., 2 Oct 2025).

1. Model class and scientific position

A recurrent misconception is to regard NeuralGCM as a purely data-driven emulator. Across the literature, it is instead described as a hybrid or “grey-box” GCM: a traditional dynamical core advances the resolved atmospheric state, while neural components replace or augment subgrid-scale parameterizations such as convection, cloud microphysics, radiation, boundary-layer turbulence, and surface fluxes (Kochkov et al., 2023, Rucker et al., 31 Oct 2025). This architectural split is central to the model’s identity, because it preserves explicit dynamical balances while allowing learned closures to absorb systematic deficiencies in conventional parameterization schemes.

The model was proposed in the context of two longstanding gaps in ML-based atmospheric modeling. First, purely learned weather models had demonstrated high deterministic forecast skill but had not shown comparable ensemble calibration or stable multi-year climate integrations. Second, offline ML parameterizations had not established robust behavior when embedded in long-horizon atmospheric simulation. NeuralGCM addresses both issues by training learned physics online in the context of a differentiable solver, rather than as an external surrogate detached from prognostic integration (Kochkov et al., 2023).

The published record positions NeuralGCM between conventional atmosphere models and fully autoregressive neural emulators such as ACE2. In comparative studies, this hybrid construction is repeatedly associated with physically consistent circulation responses, improved dynamical structure, and stable long integrations, while also exposing characteristic limitations: missing land-surface feedbacks in atmosphere-only configurations, incomplete slow-mode fidelity in the stratosphere, and uncertain extrapolation to perturbed climates outside the training distribution (Zhang et al., 2 Oct 2025, Rucker et al., 31 Oct 2025, Baxter et al., 6 Oct 2025).

2. Dynamical core, state representation, and coupling

Across studies, NeuralGCM is consistently formulated as a prognostic time-stepping system in which resolved dynamics and learned subgrid tendencies are added within each update. Representative reported forms include

$$x_{t+\Delta t}=M_{\rm dyn}\!\left(x_t+\Delta t\,\phi(x_t;\theta)\right),$$

and

$$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$

where the neural component predicts subgrid increments or tendencies and the core advances the atmospheric state (Tian et al., 2024, Rucker et al., 31 Oct 2025).

The original formulation uses a fully differentiable spectral-dynamical core on pseudo-spectral Gaussian grids at TL63, TL127, or TL255 resolution, with 32 terrain-following $\sigma$ levels and semi-implicit IMEX Runge–Kutta time integration (Kochkov et al., 2023). Other evaluations describe closely related configurations at 2.8° horizontal resolution with 37 standard pressure or sigma levels, fixed global-mean surface pressure, or approximately 30 levels extending to about 1 hPa, together with short core time steps of roughly 5–15 minutes and outer 6-hour rollout intervals for training and forecasting (Baxter et al., 6 Oct 2025, Chen et al., 11 Feb 2026, Zhang et al., 2 Oct 2025). A tangent-linear and adjoint study expresses the dynamical core on a centroidal-Voronoi mesh with C-grid staggering and split-explicit time stepping, emphasizing the differentiable composition of physical-core and neural-physics Jacobians rather than long-climate evaluation (Tian et al., 2024). These differences indicate that the literature reports several task-specific views or versions of the model rather than a single immutable discretization.

Reported prognostic states include combinations of temperature, winds, humidity, surface pressure, geopotential, divergence, vorticity, cloud liquid water, and cloud ice (Kochkov et al., 2023, Whittaker et al., 12 Jun 2025, Yuval et al., 2024). In the original weather-and-climate paper, the core state per grid cell and vertical level comprises divergence $\delta$, vorticity $\zeta$, temperature perturbation $T'$, log-surface-pressure $\ln p_s$, and specific humidities $q$, $q_{c_l}$, and $q_{c_i}$ (Kochkov et al., 2023). Other task-specific formulations use states such as $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$0 or $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$1, reflecting the same hybrid principle under different diagnostic or implementation conventions (Tian et al., 2024, Whittaker et al., 12 Jun 2025).

The learned physics modules are likewise described at multiple levels of detail. The original model uses a single-column Encode–Process–Decode ResNet with a 384-dimensional latent state, five residual blocks of 3-layer MLPs, standardized inputs, location embeddings, and encoder/decoder networks that map between ERA5 pressure-level data and the model’s sigma-coordinate state (Kochkov et al., 2023). Later descriptions refer to a feedforward convolutional-block network or a feed-forward deep neural network with 4–6 layers, 256–512 neurons per layer, residual connections, and layer normalization (Tian et al., 2024, Baxter et al., 6 Oct 2025). A precipitation-optimized variant adds an explicit precipitation network with latent size 64 and three residual blocks, predicting hourly precipitation while diagnosing evaporation from a column water-budget constraint (Yuval et al., 2024).

A further defining feature is stochasticity. The stochastic NeuralGCM injects random perturbations into learned parameterizations to generate ensemble spread (Kochkov et al., 2023, Asch et al., 17 Jun 2026). In the original model, 20 independent Gaussian random fields are injected into the encoder and physics module, with learned spatial and temporal scales, enabling an unbiased two-sample CRPS loss during training (Kochkov et al., 2023).

3. Training data, objectives, and probabilistic calibration

Training regimes vary by application, but the dominant supervision source is ERA5 reanalysis. The original model was trained on ERA5 from 1979–2019, hourly on 37 pressure levels at 0.25°, conservatively regridded to the model grid, with a deterministic curriculum from 6-hour to 5-day rollouts and an analogous stochastic curriculum (Kochkov et al., 2023). Other studies report ERA5-based training over 1979–2017, 1979–2018, 1979–2019, 1981–2020, or 2001–2018, depending on whether the target task is weather prediction, climate simulation, seasonal forecasting, or forced-response evaluation (Rucker et al., 31 Oct 2025, Zhang et al., 30 Apr 2025, Whittaker et al., 12 Jun 2025, Chen et al., 11 Feb 2026).

Observed precipitation has become an additional training target in later work. In the precipitation-optimized configuration, IMERG V07 precipitation is conservatively regridded to 2.8° and incorporated jointly with ERA5 state variables; the model is then selected for both stability in decade-long simulations and low validation precipitation CRPS (Yuval et al., 2024). A related uncertainty-quantification study states that NeuralGCM uses ERA5 temperature and wind fields from 1979–2018 and is fine-tuned against IMERG precipitation estimates from 2000–2018, with stochastic closures tuned to match observed precipitation statistics (Asch et al., 17 Jun 2026).

The training losses also differ by configuration. The original deterministic NeuralGCM minimizes a weighted sum of filtered rollout MSE, spectral power mismatch, and batch-mean bias,

$$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$2

with decoder fine-tuning via plain grid-space MSE on 1-day to 24-hour rollouts (Kochkov et al., 2023). The stochastic version uses CRPS in nodal and spectral space (Kochkov et al., 2023). Other descriptions emphasize mean-squared error on one-step-ahead tendencies or full-state forecasts, sometimes augmented by conservation penalties, weight decay, or weak physics regularization (Zhang et al., 2 Oct 2025, Tian et al., 2024, Rucker et al., 31 Oct 2025, Chen et al., 11 Feb 2026). One climate-change evaluation explicitly notes soft penalties for mass and water-vapor conservation blocks and $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$3 regularization on network weights, optimized with Adam at learning rates of order $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$4 (Zhang et al., 2 Oct 2025).

Probabilistic post-processing has become a separate line of work. For NeuralGCM ensembles, online conformal prediction wraps raw ensemble quantiles in an adaptive interval

$$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$5

with the scalar padding updated by $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$6 (Asch et al., 17 Jun 2026). At 5-day lead time and target 90% coverage, the raw NeuralGCM 95th–5th percentile interval achieved area-weighted empirical coverage of approximately 85.6%, increasing to approximately 90.0% after conformal adjustment; extreme-event coverage improved by up to 15–20 percentage points in many regions, while global CRPS changed by less than 1% (Asch et al., 17 Jun 2026).

4. Forecast, climate, and process-level performance

In deterministic weather forecasting, NeuralGCM-0.7° is reported to be comparable to GraphCast at 1–3 days and to outperform ECMWF-HRES in RMSE for $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$7, $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$8, and $$x(t+\Delta t)=x(t)+\Delta x_{\rm dyn}+\Delta x_{\rm sub},$$9 on 2020 forecasts (Kochkov et al., 2023). In probabilistic forecasting, NeuralGCM-ENS at 1.4° shows lower ensemble-mean RMSE and lower CRPS than ECMWF-ENS across 1–15 days, with spread–skill ratio near 1 (Kochkov et al., 2023). The same study emphasizes stable power spectra, improved geostrophic balance relative to pure ML models, realistic derived fields, and water-budget consistency through diagnosed $\sigma$0 (Kochkov et al., 2023).

For climate simulation, the original atmosphere-only experiments report 35 two-year runs at 1.4° forced by observed SST and sea ice, with global $\sigma$1 RMSE of 0.16 K versus 0.45 K for climatology, and realistic seasonal cycles in precipitable water, total kinetic energy, Hadley circulation, monsoon winds, and eddy kinetic energy (Kochkov et al., 2023). At 2.8°, 22 stable 40-year AMIP-like runs reproduce the global warming trend, and the spatial 850-hPa bias is smaller than all 22 CMIP6 AMIP models considered in that comparison (Kochkov et al., 2023). Emergent tropical cyclone statistics are also reported: a 1.4° run yields 83 tropical cyclones versus 86 in ERA5, whereas an X-SHiELD reference regridded to 1.4° yields 40 (Kochkov et al., 2023).

A PMP-based evaluation broadens this picture. In AMIP mode, with outputs interpolated to a common 2.5°×2.5° grid, NeuralGCM shows portrait-plot nRMSE values of approximately 0.90 for ta200, 0.85 for ta850, 0.95 for ua200, 0.92 for ua850, 0.96 for va200, 0.94 for va850, and 0.88 for psl (Pallotta et al., 8 Apr 2026). ENSO skill is mechanically high because prescribed SSTs drive the Niño 3.4 index, with correlation greater than 0.95. Extratropical mode amplitude ratios are near unity: approximately 1.02 for SAM, 0.98 for NAM, 1.00 for NAO, and 0.97 for PNA (Pallotta et al., 8 Apr 2026). MJO propagation is weaker than observed, with NeuralGCM EWR = 2.05 versus observed $\sigma$2 and EOR around 0.90, indicating slight underestimation of eastward MJO power (Pallotta et al., 8 Apr 2026).

Circulation-variability diagnostics reach similar conclusions. NeuralGCM reproduces MJO, Kelvin, and Rossby peaks in tropical wave spectra, underestimates high-frequency Kelvin and inertio-gravity power by about 20–30%, and matches the MJO band within 10% of ERA5 (Baxter et al., 6 Oct 2025). Its 250-hPa eddy-momentum-flux co-spectrum closely follows ERA5, with peak flux near $\sigma$3 m s$\sigma$4 and latitude about $\sigma$5S, although amplitudes are slightly overestimated by 10–15% (Baxter et al., 6 Oct 2025).

Precipitation-specific optimization materially alters the model’s performance envelope. On 24-hour accumulated precipitation over the 2020 test year and latitudes $\sigma$6, a precipitation-trained NeuralGCM shows CRPS about 10–20% lower than the ECMWF ensemble on days 1–15 and ensemble-mean RMSE of about 0.8 mm day$\sigma$7 versus approximately 1.0 mm day$\sigma$8 for ECMWF ENS (Yuval et al., 2024). In 2002–2014 climate statistics, global mean-precipitation MAE is 0.45 mm day$\sigma$9, compared with approximately 0.74 mm day$\delta$0 across 37 CMIP6 AMIP models; over land, Rx1day bias MAE is about 1.8 mm day$\delta$1 versus about 4.0 mm day$\delta$2 for a representative CMIP6 model, and peak-time MAE of the diurnal cycle is about 1.2 h versus approximately 2.1 h for ERA5 (Yuval et al., 2024).

5. Stress tests, seasonal prediction, and extreme-event applications

NeuralGCM has been examined under increasingly demanding use cases beyond standard weather verification. For seasonal prediction, one study prescribes simplified boundary forcings in which SST and sea ice follow their climatological annual cycle while persisting anomalies present at initialization (Zhang et al., 30 Apr 2025). In 5-month hindcasts for 1990–2023, July–November ACCs for 500-hPa geopotential height and 1000-hPa temperature exceed 0.7 in tropical main development regions through July–September and remain statistically significant through October, while NRMSE stays below unity (Zhang et al., 30 Apr 2025). Seasonal tropical-cyclone frequency correlations with IBTrACS are reported as $\delta$3 in the North Atlantic, $\delta$4 in the East Pacific, and $\delta$5 in the West Pacific for deterministic/stochastic configurations; track-density pattern correlation reaches $\delta$6 (Zhang et al., 30 Apr 2025).

Heatwave storyline studies exploit a different property: full differentiability. For the 2021 Pacific Northwest event, a nudged or ensemble-based NeuralGCM reproduces the observed heatwave with daily-max area-mean peak metrics of RMSE $\delta$7 K, Bias $\delta$8 K, and $\delta$9 for NGCM-ENS, compared with RMSE $\zeta$0 K, Bias $\zeta$1 K, and $\zeta$2 for E3SM (Duan et al., 2024). In a pseudo-global-warming experiment for 2050, NeuralGCM projects $\zeta$3 K and $\zeta$4 K over the heatwave core, while a fully spun-up E3SM future-land configuration projects substantially larger warming, illustrating the effect of omitted land feedbacks (Duan et al., 2024).

Differentiability also enables direct optimization of initial conditions for extreme storylines. Using an 11-day forecast initialized on 21 June 2021 and 75 Adam steps, one study constructs a trajectory whose 5-day average peak temperature is $\zeta$5C above the most extreme member of a 75-member stochastic-NeuralGCM ensemble; the optimized trajectory reaches 2 m temperatures up to $\zeta$6C versus approximately $\zeta$7C in the control ensemble and exhibits enhanced 500-hPa wavenumber-3 amplitude and blocking (Whittaker et al., 12 Jun 2025).

Climate-change stress tests show both strengths and unresolved structural issues. Under a +2 K uniform-SST benchmark, NeuralGCM produces global-mean $\zeta$8 K at 1000 hPa, land-enhanced warming up to about +2.5 K, global precipitation increase of about +6% or roughly 3% K$\zeta$9, an upper-tropospheric warming maximum of about +3.0 K near 200 hPa, mid-tropospheric warming of about +2.3 K at 500 hPa, lower-stratospheric cooling of about $T'$0 K, and subtropical jet strengthening up to +4 m s$T'$1 (Zhang et al., 2 Oct 2025). A separate hierarchical evaluation including P2K and P4K finds $T'$2 K K$T'$3 in a 36-year P4K free run and $T'$4 K K$T'$5 in a 20-member one-year +4 K ensemble, with $T'$6 and about $T'$7 % K$T'$8, respectively; the zonal-mean warming structure is captured qualitatively in the troposphere but degrades in the UTLS, where pattern correlation falls from 0.62 to 0.39 when upper levels are included (Chen et al., 11 Feb 2026).

Regional trend benchmarks further delimit the regime of skill. Over 1981–2014, NeuralGCM captures Arctic amplification with Arctic 850-hPa warming of $T'$9 K decade$\ln p_s$0 versus 0.25 in ERA5, and it captures tropical upper-tropospheric warming at 250 hPa with $\ln p_s$1 K decade$\ln p_s$2 versus 0.20 in ERA5 (Rucker et al., 31 Oct 2025). It also captures heat-extreme trends in Western Europe and the Midwest US, but underestimates Southwest US TXx trends at $\ln p_s$3 K decade$\ln p_s$4 versus 0.35 in ERA5 and fails to capture Southwest US drying, with 700-hPa humidity trend $\ln p_s$5 g kg$\ln p_s$6 decade$\ln p_s$7 versus $\ln p_s$8 in ERA5 (Rucker et al., 31 Oct 2025).

6. Failure modes, controversies, and development directions

The most persistent limitations concern slow modes, linearized operators, land-surface coupling, and extrapolation. In circulation-variability benchmarks, NeuralGCM fails to produce a consistent quasi-biennial oscillation: no ensemble member shows a stable $\ln p_s$9 month period, the 50-hPa amplitude remains in the 20–30 m s$q$0 range versus 40 m s$q$1 in ERA5, and vertical Hovmöller diagnostics show no coherent downward migration (Baxter et al., 6 Oct 2025). Southern Annular Mode propagation is partially captured in autocorrelation structure, yet the observed $q$2 day spectral peak is absent; the ensemble-mean spectrum is red-noise-like (Baxter et al., 6 Oct 2025). These studies attribute the problem to short rollout losses, bottom-heavy weighting, and limited stratospheric vertical resolution (Baxter et al., 6 Oct 2025).

Data-assimilation studies expose a different weakness. Tangent-linear and adjoint models derived from NeuralGCM exhibit patchy temperature sensitivities, spurious vertical lobes, and elevated high-wavenumber noise, with vertical-level variance in temperature adjoints 2–5 times larger than that of MPAS-A and spectral power concentrated at vertical wavenumbers $q$3 (Tian et al., 2024). The implication is that naïvely inserting NeuralGCM TL/AD operators into 4D-Var or EnKF systems can contaminate gradients and ensemble covariances with unphysical small-scale structure (Tian et al., 2024).

Out-of-distribution climate forcing remains unresolved rather than absent. The +2 K and +4 K warming studies show that NeuralGCM can reproduce several robust large-scale responses, particularly precipitation scaling and amplified tropical upper-tropospheric warming, but also that it can depart from expected radiative and stratospheric behavior (Zhang et al., 2 Oct 2025, Chen et al., 11 Feb 2026). Reported issues include unavailable TOA radiation diagnostics in some configurations, occasional individual blow-ups, underpredicted poleward jet shifts, UTLS cooling in the deep tropics, excessive polar upper-level warming, and biased eddy heat- and momentum-flux convergence under warming (Zhang et al., 2 Oct 2025, Chen et al., 11 Feb 2026). A plausible implication is that the hybrid architecture is strongest where training data constrain tropospheric balances and weakest where the learned closures must infer poorly sampled stratospheric or radiative processes.

Land-surface omission is another recurrent constraint. Several studies explicitly state that the atmosphere-only model lacks interactive soil moisture, evapotranspiration, vegetation, and land-energy feedbacks (Duan et al., 2024, Pallotta et al., 8 Apr 2026, Rucker et al., 31 Oct 2025). This is linked to underestimation of future heatwave amplitude relative to land-coupled references and to difficulty reproducing regional drying and near-surface heat-extreme trends in arid regions (Duan et al., 2024, Rucker et al., 31 Oct 2025).

The published development directions are correspondingly concrete. They include extending training to perturbed climates, adding land-surface and coupled ocean feedbacks, refining learned precipitation and cloud processes, strengthening conservation and stability constraints, increasing stratospheric resolution or loss weighting, and incorporating multi-scale objectives that penalize monthly-to-annual drift rather than only 6-hour to 5-day errors (Pallotta et al., 8 Apr 2026, Baxter et al., 6 Oct 2025, Zhang et al., 2 Oct 2025). Taken together, these recommendations do not negate NeuralGCM’s established capabilities; they define the present boundary between its demonstrated competence in fast, large-scale atmospheric dynamics and the remaining requirements for robust long-horizon Earth-system prediction.