We propose a hybrid method for accurately estimating the score function, i.e., the gradient of the log steadystate density, using a Gaussian Mixture Model (GMM) in conjunction with a bisecting K-means clustering step. Our approach, which we call KGMM, offers a systematic way to combine statistical density estimation with a neural-network-based interpolation of the score, leveraging the strengths of both. We demonstrate its ability to accurately reconstruct the long-time statistical properties of several paradigmatic systems, including potential systems, and chaotic Lorenz-type models, and the Kuramoto-Sivashinsky equation. Numerical experiments show that KGMM yields robust estimates of the score function, even for small values of the covariance amplitude in the GMM, where the standard GMM methods tend to fail because of noise amplification. We compare the performance of KGMM against the conventional Denoising Score Matching (DSM) approach, demonstrating that KGMM achieves more faithful reconstruction of the steady-state distribution for low-dimensional systems at a fraction of the computational cost. These accurate estimates allow us to build effective stochastic reduced-order models that reproduce the invariant measures of the target dynamics.

Many natural systems exhibit cyclo-stationary behavior characterized by periodic forcing such as annual and diurnal cycles. We present a data-driven method leveraging recent advances in score-based generative modeling to construct reduced-order models for such cyclo-stationary time series. Our approach accurately reproduces the statistical properties and temporal correlations of the original data, enabling efficient generation of synthetic trajectories. We demonstrate the performance of the method through application to the Planet Simulator (PlaSim) climate model, constructing a reduced-order model for the 20 leading principal components of surface temperature driven by the annual cycle. The resulting surrogate model accurately reproduces the marginal and joint probability distributions, autocorrelation functions, and spatial coherence of the original climate system across multiple validation metrics. The approach offers substantial computational advantages, enabling generation of centuries of synthetic climate data in minutes compared to weeks required for equivalent full model simulations. This work opens new possibilities for efficient modeling of periodically forced systems across diverse scientific domains, providing a principled framework for balancing computational efficiency with physical fidelity in reduced-order modeling applications.

Understanding subsurface ocean dynamics is essential for quantifying oceanic heat and mass transport, but direct observations at depth remain sparse due to logistical and technological constraints. In contrast, satellite missions provide rich surface datasets-such as sea surface height, temperature, and salinity-that offer indirect but potentially powerful constraints on the ocean interior. Here, we present a probabilistic framework based on score-based diffusion models to reconstruct three-dimensional subsurface velocity and buoyancy fields, including the energetic ocean eddy field, from surface observations. Using a 15-level primitive equation simulation of an idealized double-gyre system, we evaluate the skill of the model in inferring the mean circulation and the mesoscale variability at depth under varying levels of surface information. We find that the generative model successfully recovers key dynamical structures and provides physically meaningful uncertainty estimates, with predictive skill diminishing systematically as the surface resolution decreases or the inference depth increases. These results demonstrate the potential of generative approaches for ocean state estimation and uncertainty quantification, particularly in regimes where traditional deterministic methods are underconstrained or ill-posed.

Recent developments in generative modeling have utilized score-based methods coupled with stochastic differential equations to sample from complex probability distributions. However, these and other performant sampling methods generally require gradients of the target probability distribution, which can be unavailable or computationally prohibitive in many scientific and engineering applications. Here, we introduce ensembles within score-based sampling methods to develop gradient-free approximate sampling techniques that leverage the collective dynamics of particle ensembles to compute approximate reverse diffusion drifts. We introduce the underlying methodology, emphasizing its relationship with generative diffusion models and the previously introduced Föllmer sampler. We demonstrate the efficacy of the ensemble strategies through various examples, ranging from low- to medium-dimensionality sampling problems, including multi-modal and highly non-Gaussian probability distributions, and provide comparisons to traditional methods like the No-U-Turn Sampler. Additionally, we showcase these strategies in the context of a high-dimensional Bayesian inversion problem within the geophysical sciences. Our findings highlight the potential of ensemble strategies for modeling complex probability distributions in situations where gradients are unavailable.

We introduce an approach for analyzing the responses of dynamical systems to external perturbations that combines score-based generative modeling with the generalized fluctuation-dissipation theorem. The methodology enables accurate estimation of system responses, including those with non-Gaussian statistics. We numerically validate our approach using time-series data from three different stochastic partial differential equations of increasing complexity: an Ornstein-Uhlenbeck process with spatially correlated noise, a modified stochastic Allen-Cahn equation, and the 2D Navier-Stokes equations. We demonstrate the improved accuracy of the methodology over conventional methods and discuss its potential as a versatile tool for predicting the statistical behavior of complex dynamical systems.

We present a method to downscale idealized geophysical fluid simulations using generative models based on diffusion maps. By analyzing the Fourier spectra of fields drawn from different data distributions, we show how a diffusion bridge can be used as a transformation between a low-resolution and a high-resolution dataset, allowing for new sample generation of high-resolution fields given specific low-resolution features. The ability to generate new samples allows for the computation of any statistic of interest, without any additional calibration or training. Our unsupervised setup is also designed to downscale fields without access to paired training data; this flexibility allows for the combination of multiple source and target domains without additional training. We demonstrate that the method enhances resolution and corrects context-dependent biases in geophysical fluid simulations, including in extreme events. We anticipate that the same method can be used to downscale the output of climate simulations, including temperature and precipitation fields, without needing to train a new model for each application and providing a significant computational cost savings.

Dynamical cores used to study the circulation of the atmosphere employ various numerical methods ranging from finite-volume, spectral element, global spectral, and hybrid methods. In this work, we explore the use of Flux-Differencing Discontinuous Galerkin (FDDG) methods to simulate a fully compressible dry atmosphere at various resolutions. We show that the method offers a judicious compromise between high-order accuracy and stability for large-eddy simulations and simulations of the atmospheric general circulation. In particular, filters, divergence damping, diffusion, hyperdiffusion, or sponge-layers are not required to ensure stability; only the numerical dissipation naturally afforded by FDDG is necessary. We apply the method to the simulation of dry convection in an atmospheric boundary layer and in a global atmospheric dynamical core in the standard benchmark of Held and Suarez (1994, https://doi. org/10.1175/1520-0477(1994)075<1825:apftio>2.0.co;2).Plain Language Summary Numerical models cannot explicitly represent all degrees of freedom that characterize atmospheric flows due to limitations in computing power. One must allocate the available computational degrees of freedom to reduce the degradation of the solution. In this work, we explore the use of the discontinuous Galerkin numerical method, a hybrid approach that combines the accuracy of spectral methods with the flexibility of finite volume methods. We apply it to idealized dry atmospheric simulations and show that the method is robust and incorporates physical principles to best account for unresolved processes.

Generated images of score-based models can suffer from errors in their spatial means, an effect, referred to as a color shift, which grows for larger images. This paper investigates a previously-introduced approach to mitigate color shifts in score-based diffusion models. We quantify the performance of a nonlinear bypass connection in the score network, designed to process the spatial mean of the input and to predict the mean of the score function. We show that this network architecture substantially improves the resulting quality of the generated images, and that this improvement is approximately independent of the size of the generated images. As a result, this modified architecture offers a simple solution for the color shift problem across image sizes. We additionally discuss the origin of color shifts in an idealized setting in order to motivate the approach.

The small-scale microphysical processes governing the formation of precipitation particles cannot be resolved explicitly by cloud resolving and climate models. Instead, they are represented by microphysics schemes that are based on a combination of theoretical knowledge, statistical assumptions, and fitting to data ("tuning"). Historically, tuning was done in an ad hoc fashion, leading to parameter choices that are not explainable or repeatable. Recent work has treated it as an inverse problem that can be solved by Bayesian inference. The posterior distribution of the parameters given the data-the solution of Bayesian inference-is found through computationally expensive sampling methods, which require over O 105 evaluations of the forward model; this is prohibitive for many models. We present a proof of concept of Bayesian learning applied to a new bulk microphysics scheme named "Cloudy," using the recently developed Calibrate-Emulate-Sample (CES) algorithm. Cloudy models collision-coalescence and collisional breakup of cloud droplets with an adjustable number of prognostic moments and with easily modifiable assumptions for the cloud droplet mass distribution and the collision kernel. The CES algorithm uses machine learning tools to accelerate Bayesian inference by reducing the number of forward evaluations needed to O 102. It also exhibits a smoothing effect when forward evaluations are polluted by noise. In a suite of perfect-model experiments, we show that CES enables computationally efficient Bayesian inference of parameters in Cloudy from noisy observations of moments of the droplet mass distribution. In an additional imperfect-model experiment, a collision kernel parameter is successfully learned from output generated by a Lagrangian particle-based microphysics model.

Microphysics methods for climate models and numerical weather prediction typically track one, two, or three moments of a droplet size distribution for various categories of liquid, ice, and aerosol. Such methods rely on conversion parameters between these categories, which introduces uncertainty into predictions. While higher‐fidelity options such as bin and Lagrangian schemes exist, they require too many degrees of freedom for climate modeling applications and introduce numerical challenges. Here we introduce a flexible spectral microphysics method based on collocation of basis functions. This method generalizes to a linear bulk scheme when using few basis functions and to a smoothed bin scheme with more degrees of freedom. Tested in an idealized box setting, the method improves spectral accuracy for droplet collision‐coalescence and may eliminate the need for precipitation autoconversion rates required by bulk methods; furthermore, it generalizes well to multimodal distributions with less complexity than a bin method. The potential to extend this collocation representation to multiple hydrometeor classes suggests a path forward to unify liquid, ice, and aerosol microphysics in a single, flexible, computational framework for climate modeling.

Reliable projections of sea‐level rise depend on accurate representations of how fast‐flowing glaciers slip along their beds. The mechanics of slip are often parameterized as a constitutive relation (or “sliding law”) whose proper form remains uncertain. Here, we present a novel deep learning‐based framework for learning the time evolution of drag at glacier beds from time‐dependent ice velocity and elevation observations. We use a feedforward neural network, informed by the governing equations of ice flow, to infer spatially and temporally varying basal drag and associated uncertainties from data. We test the framework on 1D and 2D ice flow simulation outputs and demonstrate the recovery of the underlying basal mechanics under various levels of observational and modeling uncertainties. We apply this framework to time‐dependent velocity data for Rutford Ice Stream, Antarctica, and present evidence that ocean‐tide‐driven changes in subglacial water pressure drive changes in ice flow over the tidal cycle.

Tropical rainfall to first order responds to variations in Earth’s orbit through shifts of the intertropical convergence zone (ITCZ) and changes in zonally averaged rainfall intensity. Here, a conceptual model is developed that represents both processes and their response to orbital insolation variations. The model predicts the seasonal evolution of tropical rainfall between 30°S and 30°N. Insolation variations impact seasonal rainfall in two different ways: thermodynamically, leading to variations in rainfall intensity through modulation of the water vapor content of the atmosphere; and dynamically, leading to shifts of the ITCZ through modulation of the global atmospheric energy budget. Thermodynamic and dynamic effects act together to shape the annual-mean response of tropical rainfall to changes in Earth’s orbit. The model successfully reproduces changes in annual-mean rainfall inferred from paleo-proxies across several glacial–interglacial cycles. It illuminates how orbital precession and variations of Earth’s obliquity affect tropical rainfall in distinct ways near the equator and farther away from it, with spectral signatures of precession and obliquity variations that shift with latitude. It also provides explanations for the observed different phasings of rainfall minima and maxima near the equator and away from it. For example, the model reproduces a phase shift of ~10 ka between rainfall records from caves in northern Borneo (4°N) and from China (approximately 30°N). The model suggests that such phase shifts arise through a different weighting of ITCZ shifts and variations in rainfall intensity, thus providing insight into the mechanisms that drive tropical rainfall changes on orbital time scales.

In the zonal mean, the ITCZ lies at the foot of the ascending branch of the tropical mean meridional circulation, close to where the near-surface meridional mass flux vanishes. The ITCZ also lies near the energy flux equator (EFE), where the column-integrated meridional energy flux vanishes. This latter observation makes it possible to relate the ITCZ position to the energy balance, specifically the atmospheric net energy input near the equator and the cross-equatorial energy flux. Here the validity of the resulting relations between the ITCZ position and energetic quantities is examined with reanalysis data for the years 1979–2014. In the reanalysis data, the EFE and ITCZ position indeed covary on time scales of seasons and longer. Consistent with theory, the ITCZ position is proportional to the cross-equatorial atmospheric energy flux and inversely proportional to atmospheric net energy input at the equator. Variations of the cross-equatorial energy flux dominate seasonal variations of the ITCZ position. By contrast, variations of the equatorial net energy input, driven by ocean energy uptake variations, dominate interannual variations of the ITCZ position (e.g., those associated with ENSO).

We examine how tropical zonal mean precipitation biases in current climate models relate to the atmospheric energy budget. Both hemispherically symmetric and antisymmetric tropical precipitation biases contribute to the well‐known double‐Intertropical Convergence Zone (ITCZ) bias; however, they have distinct signatures in the energy budget. Hemispherically symmetric biases in tropical precipitation are proportional to biases in the equatorial net energy input; hemispherically antisymmetric biases are proportional to the atmospheric energy transport across the equator. Both relations can be understood within the framework of recently developed theories. Atmospheric net energy input biases in the deep tropics shape both the symmetric and antisymmetric components of the double‐ITCZ bias. Potential causes of these energetic biases and their variation across climate models are discussed.

The ocean’s overturning circulation is inherently three-dimensional, yet modern quantitative estimates of the overturning typically represent the subsurface circulation as a two-dimensional, two-cell streamfunction that varies with latitude and depth only. This approach suppresses information about zonal mass and tracer transport. In this article, the authors extend earlier, zonally averaged overturning theory to explore the dynamics of a “figure-eight” circulation that cycles through multiple basins. A three-dimensional residual-mean model of the overturning circulation is derived and then simplified to a multibasin isopycnal box model to explore how stratification and diabatic water mass transformations in each basin depend on the basin widths and on deep and bottom-water formation in both hemispheres. The idealization to multiple, two-dimensional basins permits zonal mass transport along isopycnals in a Southern Ocean–like channel, while retaining the dynamical framework of residual-mean theory. The model qualitatively reproduces the deeper isopycnal surfaces in the Pacific Basin relative to the Atlantic. This supports a transfer of Antarctic Bottom Water from the Atlantic sector to the Pacific sector via the Southern Ocean, which subsequently upwells in the northern Pacific Basin. A solution for the full isopycnal structure in the Southern Ocean reproduces observed stratification differences between Atlantic and Pacific Basins and provides a scaling for the diffusive boundary layer in which the zonal mass transport occurs. These results are consistent with observational indications that North Atlantic Deep Water is preferentially transformed into Antarctic Bottom Water, which undermines the importance of an adiabatic, upper overturning cell in the modern ocean.

We show that the climate phenomena of El Nino and La Nina arise naturally as states of macro-variables when our recent causal feature learning framework (Chalupka 2015, Chalupka 2016) is applied to micro-level measures of zonal wind (ZW) and sea surface temperatures (SST) taken over the equatorial band of the Pacific Ocean. The method identifies these unusual climate states on the basis of the relation between ZW and SST patterns without any input about past occurrences of El Nino or La Nina. The simpler alternatives of (i) clustering the SST fields while disregarding their relationship with ZW patterns, or (ii) clustering the joint ZW-SST patterns, do not discover El Nino. We discuss the degree to which our method supports a causal interpretation and use a low-dimensional toy example to explain its success over other clustering approaches. Finally, we propose a new robust and scalable alternative to our original algorithm (Chalupka 2016), which circumvents the need for high-dimensional density learning.

The ITCZ lies at the ascending branch of the tropical meridional overturning circulation, where near-surface meridional mass fluxes vanish. Near the ITCZ, column-integrated energy fluxes vanish, forming an atmospheric energy flux equator (EFE). This paper extends existing approximations relating the ITCZ position and EFE to the atmospheric energy budget by allowing for zonal variations. The resulting relations are tested using reanalysis data for 1979-2014. The zonally varying EFE is found as the latitude where the meridional component of the divergent atmospheric energy transport (AET) vanishes. A Taylor expansion of the AET around the equator relates the ITCZ position to derivatives of the AET. To a first order, the ITCZ position is proportional to the divergent AET across the equator; it is inversely proportional to the local atmospheric net energy input (NEI) that consists of the net energy fluxes at the surface, at the top of the atmosphere, and zonally across longitudes. The first-order approximation captures the seasonal migrations of the ITCZ in the African, Asian, and Atlantic sectors. In the eastern Pacific, a third-order approximation captures the bifurcation from single-to double-ITCZ states that occurs during boreal spring. In contrast to linear EFE theory, during boreal winter in the eastern Pacific, northward cross-equatorial AET goes along with an ITCZ north of the equator. EFE and ITCZ variations driven by ENSO are characterized by an equatorward (poleward) shift in the Pacific during El Nino (La Nina) episodes, which are associated with variations in equatorial ocean energy uptake.

This work addresses the question of the stability of stratified, spatially periodic shear flows at low Péclet number but high Reynolds number. This little-studied limit is motivated by astrophysical systems, where the Prandtl number is often very small. Furthermore, it can be studied using a reduced set of “low-Péclet-number equations” proposed by Lignières [“The small-Péclet-number approximation in stellar radiative zones,” Astron. Astrophys. 348, 933–939 (1999)]. Through a linear stability analysis, we first determine the conditions for instability to infinitesimal perturbations. We formally extend Squire’s theorem to the low-Péclet-number equations, which shows that the first unstable mode is always two-dimensional. We then perform an energy stability analysis of the low-Péclet-number equations and prove that for a given value of the Reynolds number, above a critical strength of the stratification, any smooth periodic shear flow is stable to perturbations of arbitrary amplitude. In that parameter regime, the flow can only be laminar and turbulent mixing does not take place. Finding that the conditions for linear and energy stability are different, we thus identify a region in parameter space where finite-amplitude instabilities could exist. Using direct numerical simulations, we indeed find that the system is subject to such finite-amplitude instabilities. We determine numerically how far into the linearly stable region of parameter space turbulence can be sustained.

Motivated by experiments on sheared suspensions that show a transition between ordered and disordered phases, we here study the long-time behavior of a sheared and overdamped two-dimensional system of particles interacting by repulsive forces. As a function of interaction strength and shear rate we find transitions between phases with vanishing and large single-particle diffusion. In the phases with vanishing single-particle diffusion, the system evolves towards regular lattices, usually on very slow time scales. Different lattices can be approached, depending on interaction strength and forcing amplitude. The disordered state appears in parameter regions where the regular lattices are unstable. Correlation functions between the particles reveal the formation of shear bands. In contrast to single-particle densities, the spatially resolved two-particle correlation functions vary with time and allow to determine the phase within a period. As in the case of the suspensions, motion in the state with low diffusivity is essentially reversible, whereas in the state with strong diffusion it is not.