Predicting permeability from Nuclear Magnetic Resonance (NMR) data is a fundamental yet challenging task in reservoir characterization, primarily due to the uncertainty associated with surface relaxivity ( ρ ) parameters. In this work, we investigate the feasibility of using Machine Learning (ML) to estimate permeability from T_2 distributions and quantify how ρ uncertainty affects predictive accuracy. To address this, we generated a dataset of 15,000 synthetic 3D porous media representing granular sedimentary rock samples. We employed efficient in-house implementations of a Random Walk algorithm (governed by Bloch-Torrey physics) to simulate magnetization decay and obtain T_2 distributions, alongside a Finite Element Method (FEM) solver for the Stokes equations to compute absolute permeability, assuming 100 ρ . In the first experiment, we simulated T_2 distributions by assigning a constant ρ value for all synthetic porous media. In the second experiment, we applied a different ρ value to each synthetic porous medium to emulate real-world uncertainty, representing the scenario where ρ is unknown. The third experiment extends the second by converting the T_2 distributions into surface-to-volume ratio distributions using the specific ρ value assigned in the second experiment to each medium. We systematically compared the Multilayer Perceptron (MLP) performance against the industry-standard Schlumberger-Doll-Research (SDR) model. Overall, the MLP yielded strong predictive performance. The second experiment presented a performance drop for both models, confirming the impact of ρ uncertainty. The main contribution of this work is the systematic quantification of the sensitivity of predictive permeability models to ρ , establishing a controlled benchmark that addresses and reduces existing uncertainties. Additionally, these findings demonstrate that the MLP provides a robust and competitive alternative for permeability estimation in scenarios where ρ is unknown.

We describe CATKE, a parameterization for fluxes associated with small-scale or "microscale" ocean turbulent mixing on scales between 1 and 100 meters. CATKE uses a downgradient formulation that depends on a prognostic turbulent kinetic energy (TKE) variable and a diagnostic mixing length scale that includes a dynamic convective adjustment (CA) component. With its dynamic convective mixing length, CATKE predicts not just the depth spanned by convective plumes but also the characteristic convective mixing timescale, an important aspect of turbulent convection not captured by simpler static convective adjustment schemes. As a result, CATKE can describe the competition between convection and other processes such as shear-driven mixing and baroclinic restratification. To calibrate CATKE, we use Ensemble Kalman Inversion to minimize the error between 21 large eddy simulations (LES) and predictions of the LES data by CATKE-parameterized single column simulations at three different vertical resolutions. We find that CATKE makes accurate predictions of both idealized and realistic LES compared to microscale turbulence parameterizations commonly used in climate models.

NORi is a machine-learned (ML) parameterization of ocean boundary layer turbulence that is physics-based and augmented with neural networks. NORi stands for neural ordinary differential equations (NODEs) Richardson number (Ri) closure. The physical parameterization is controlled by a Richardson number-dependent diffusivity and viscosity. The NODEs are trained to capture the entrainment through the base of the boundary layer, which cannot be represented with a local diffusive closure. The parameterization is trained using large-eddy simulations in an "a posteriori" fashion, where parameters are calibrated with a loss function that explicitly depends on the actual time-integrated variables of interest rather than the instantaneous subgrid fluxes, which are inherently noisy. NORi is designed for the realistic nonlinear equation of state of seawater and demonstrates excellent prediction and generalization capabilities in capturing entrainment dynamics under different convective strengths, oceanic background stratifications, rotation strengths, and surface wind forcings. NORi is numerically stable for at least 100 years of integration time in large-scale simulations, despite only being trained on 2-day horizons, and can be run with time steps as long as one hour. The highly expressive neural networks, combined with a physically-rigorous base closure, prove to be a robust paradigm for designing parameterizations for climate models where data requirements are drastically reduced, inference performance can be directly targeted and optimized, and numerical stability is implicitly encouraged during training.

We describe the user interface, governing equations, and numerical methods underpinning the community ocean modeling software called "Oceananigans". Oceananigans development has been lead by the Climate Modeling Alliance to build a trainable climate model with quantifiable uncertainty. Oceananigans is written in the Julia programming language, which, like similar recent efforts based on modern programming languages, distinguishes it from usual software based on Fortran. Oceananigans can efficiently simulate all scales of ocean motion, ranging from millimeter-scale turbulence in a small box to planetary-scale ocean circulation. Oceananigans design combines (i) a basic structured finite volume algorithm (ii) optimized for high-resolution simulations on GPUs which is (iii) exposed behind a high-level, programmable user interface. This design negotiates a dual mandate for highest-possible performance (to support state-of-the-art applications) and enhanced accessibility (to facilitate adoption and development). The dual mandate aims ultimately to accelerate the progress of Earth system science. Achieving this aim, however, requires a substantial and sustained increase in the collective effort of Oceananigans development.

The Earth system is often modeled as a dynamical system in what has come to be known as Earth Systems Models. When used to study anthropogenically forced climate change, these models are forced in such a way that they are not in a statistically stationary state. Yet, statistical statements are still made about the Earth system using only a single trajectory by taking temporal averages. At each moment in time, one draws a sample from a different distribution, raising questions about the utility of temporal averages, in stark contrast to the utility of temporal averages in ergodic systems. This work follows in the tradition of using a toy model to examine properties present in the Earth system. We aim to examine how we can make meaningful statistical statements in non-stationary systems when only dealing with a single trajectory. We use the Lorenz equations with a time-varying parameter as a starting point for comparing ensemble averages to temporal averages. We find that, in so far as the control parameter induces a slow and smooth change in the dynamics, the resulting statistics of ensemble averages compare well to temporal averages.

This is the second part of a two-part paper. We apply the methodology of the first paper (Souza, J. Fluid Mech., vol. 997, 2024, A1) to construct a data-driven finite-volume discretization of the Liouville/Fokker-Planck equation of a high-dimensional dynamical system, i.e. the compressible Euler equations with gravity and rotation evolved on a thin spherical shell. We show that the method recovers a subset of the statistical properties of the underlying system, steady-state distributions of observables and autocorrelations of particular observables, as well as revealing the global Koopman modes of the system. We employ two different strategies for the partitioning of a high-dimensional state space, and explore their consequences.

Current eddy‐permitting and eddy‐resolving ocean models require dissipation to prevent a spurious accumulation of enstrophy at the grid scale. We introduce a new numerical scheme for momentum advection in large‐scale ocean models that involves upwinding through a weighted essentially non‐oscillatory (WENO) reconstruction. The new scheme provides implicit dissipation and thereby avoids the need for an additional explicit dissipation that may require calibration of unknown parameters. This approach uses the rotational, “vector invariant” formulation of the momentum advection operator that is widely employed by global general circulation models. A novel formulation of the WENO “smoothness indicators” is key for avoiding excessive numerical dissipation of kinetic energy and enstrophy at grid‐resolved scales. We test the new advection scheme against a standard approach that combines explicit dissipation with a dispersive discretization of the rotational advection operator in two scenarios: (a) two‐dimensional turbulence and (b) three‐dimensional baroclinic equilibration. In both cases, the solutions are stable, free from dispersive artifacts, and achieve increased “effective” resolution compared to other approaches commonly used in ocean models.

Analyzing Partial Discharge (PD) signals is crucial to assessing the health of insulation in high-voltage systems. Nevertheless, noise often distorts these signals, hindering the ability to obtain precise information. This paper proposes a novel deep-learning approach using two denoising autoencoders (DAEs) to learn data representations and eliminate noise during reconstruction. By leveraging DAEs’ capacity to capture essential features within the latent space, this method enhances the analysis of PD signals and yields more accurate results. This paper investigates the effectiveness of two deep-learning architectures for denoising partial discharge signals in high-voltage insulation systems. Experimental results carried out on a PD dataset demonstrated the efficiency of the Linear AE model in removing noise in sets A, B, and C suggesting that DAEs hold great promise in PD signal denoising.

Nuclear magnetic resonance (NMR) techniques are key in the study of porous reservoir rocks. They can provide valuable insight into the pore size distribution of the pore space of a given rock sample due to its dependence on the magnetic fluid/matrix interaction. The pore space is often studied at the μm scale through the use of micro-CT images, which are often composed of hundreds of millions of voxels, posing significant challenges to numerical simulations. In this paper, we present an image-based, fully explicit, and matrix-free finite element implementation for the simulation of NMR relaxation process that is capable of handling such large 3D problems in single GPUs. The chosen explicit time-integration scheme uses a lumped capacitance formulation and stabilization via hyperbolization, and it is capable of handling arbitrary time-step sizes with controllable error levels. The image-based representation of the pore space is used for a memory-efficient, matrix-free formulation of the time integration using massively parallel processes on a single GPU. In addition, we propose the substitution of a global digital roughness correction factor that depends on the porous space’s geometry for a problem-independent local correction factor, based on nodal neighborhoods. We show that the numerical scheme converges with successive refinements as expected and that our local correction coefficient is capable of estimating the correct S/V parameter of several different classical geometries. We tested our formulation against an image-based Random Walk simulation of four digital rock core samples, achieving good agreement between them. We manage to simulate a giga-voxel image-based model on a personal use GPU (less than 10GB of memory use) in 33 min with our FEM implementation.

The Earth system is often modeled as a dynamical system in what has come to be known as Earth System Models. When used to study anthropogenically forced climate change, these models are forced in such a way that they are not in a statistically stationary state. Yet, statistical statements are still made about the Earth climate system using only a single trajectory by taking temporal averages. At each moment in time, one draws a sample from a different distribution, raising questions about the utility of temporal averages, in stark contrast to the utility of temporal averages in ergodic systems. This work follows in the tradition of using a toy model to examine properties present in the Earth climate system. We aim to examine how we can make meaningful statistical statements in non-stationary systems when only dealing with a single trajectory. We use the Lorenz equations with a time-varying parameter as a starting point for comparing ensemble averages to temporal averages. We find that, in so far as the control parameter induces a slow and smooth change in the dynamics, the resulting statistics of ensemble averages compare well to those of temporal averages.

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.

Leveraging recent work on data-driven methods for constructing a finite state space Markov process from dynamical systems, we address two problems for obtaining further reduced statistical representations. The first problem is to extract the most salient reduced-order dynamics for a given timescale by using a modified clustering algorithm from network theory. The second problem is to provide an alternative construction for the infinitesimal generator of a Markov process that respects statistical features over a large range of time scales. We demonstrate the methodology on three low-dimensional dynamical systems with stochastic and chaotic dynamics. We then apply the method to two high-dimensional dynamical systems, the Kuramoto-Sivashinky equations and data sampled from fluid-flow experiments via Particle Image Velocimetry. We show that the methodology presented herein provides a robust reduced-order statistical representation of the underlying system.

This is the first of a two-part paper. We formulate a data-driven method for constructing finite-volume discretizations of an arbitrary dynamical system's underlying Liouville/Fokker-Planck equation. A method is employed that allows for flexibility in partitioning state space, generalizes to function spaces, applies to arbitrarily long sequences of time-series data, is robust to noise and quantifies uncertainty with respect to finite sample effects. After applying the method, one is left with Markov states (cell centres) and a random matrix approximation to the generator. When used in tandem, they emulate the statistics of the underlying system. We illustrate the method on the Lorenz equations (a three-dimensional ordinary differential equation) saving a fluid dynamical application for Part 2 (Souza, J. Fluid Mech., vol. 997, 2024, A2).

The parameterization of fluxes associated with representing unresolved dynamics in turbulent flows, especially in the atmosphere and ocean (which have a vast range of scales), remains a challenging task. This is especially true for Earth system models including complex biogeochemistry and requiring very long simulations. The problem of representing the dependence of the mean flux of a passive tracer in terms of the mean has a very long history; in this study, we take a somewhat different approach. We use a formalism showing that the mean flux will be a functional of the mean gradients, a formalism that can be used to calculate the structure of the functional which is non-local in both space and time. Two-dimensional turbulent simulations are used to explore the weight of nearby (in space or time) gradients. We also use stochastic velocities and iterated maps to show that the results are similar. The functional formalism provides an understanding of when non-locality needs to be considered and when a local eddy diffusivity can be a reasonably good approximation. Furthermore, the formalism provides guidance for the development of data-driven parameterizations.

NMR techniques are key in the study of porous reservoir rock, both experimentally and numerically. The T_2 relaxation process, the most common application of NMR, measures the loss of coherence of transversal magnetization and strongly depends on the fluid/matrix interaction—thus providing useful insights into the pore size distribution of a rock sample. The pore space is often studied at the μm scale through the use of micro-CT images, which are formed by stacks of images with hundreds of thousands of pixels each, posing significant challenges to numerical simulations. In this paper, we present an image-based, fully explicit, and matrix-free finite element implementation for the simulation of the T_2 relaxation process that is capable of handling such large problems. The utilized mathematical model considers relaxation due to bulk effects and surface relaxivity, not taking into account the effects of magnetic field gradients. We propose the usage of stable time marching schemes that use hyperbolization as means of acquiring stability with large time-steps. We compare the numerical performance of different time-marching schemes, showing that the application of the Leap-Frog method in a hyperbolized form of the equation can give the best trade-off between memory use and numerical convergence. Additionally, we show that the use of a lumped mass matrix allows for a fully explicit and simpler implementation while adding negligible amounts of numerical error.

We analyse a class of stochastic advection problems by conditionally averaging the passive tracer equation with respect to a given flow state. In doing so, we obtain expressions for the turbulent diffusivity as a function of the flow statistics spectrum. When flow statistics are given by a continuous-time Markov process with a finite state space, calculations are amenable to analytic treatment. When the flow statistics are more complex, we show how to approximate turbulent fluxes as hierarchies of finite state space continuous-time Markov processes. The ensemble average turbulent flux is expressed as a linear operator that acts on the ensemble average of the tracer. We recover the classical estimate of turbulent flux as a diffusivity tensor, the components of which are the integrated autocorrelation of the velocity field in the limit that the operator becomes local in space and time.

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.

Climate models must simulate hundreds of future scenarios for hundreds of years at coarse resolutions, and a handful of high-resolution decadal simulations to resolve localized extreme events. Using Oceananigans.jl, written from scratch in Julia, we report several achievements: First, a global ocean simulation with breakthrough horizontal resolution -- 488m -- reaching 15 simulated days per day (0.04 simulated years per day; SYPD). Second, Oceananigans simulates the global ocean at 488m with breakthrough memory efficiency on just 768 Nvidia A100 GPUs, a fraction of the resources available on current and upcoming exascale supercomputers. Third, and arguably most significant for climate modeling, Oceananigans achieves breakthrough energy efficiency reaching 0.95 SYPD at 1.7 km on 576 A100s and 9.9 SYPD at 10 km on 68 A100s -- the latter representing the highest horizontal resolutions employed by current IPCC-class ocean models. Routine climate simulations with 10 km ocean components are within reach.