NORi: An ML-Augmented Ocean Boundary Layer Parameterization

TL;DR

NORi addresses the challenge of representing shallow-ocean boundary layer turbulence by coupling a simple physics-based eddy-diffusivity closure with neural-network residuals that capture nonlocal entrainment. It uses a posteriori calibration with Neural Ordinary Differential Equations to optimize trajectory-based losses, yielding stable, long-time integrations and strong generalization across forcing regimes. In single-column tests, NORi matches LES and competes with or outperforms classic closures like k-ε and CATKE, while enabling larger time steps and reduced data needs. The large-scale Centennial double-gyre test demonstrates NORi’s potential for climate-scale applications, emphasizing its data-efficient design and the importance of physics-informed structure for robust long-time ocean simulations.

Abstract

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.

Figures (16)

• Figure 1: Schematic of the vertical profile of buoyancy as a function of depth below the ocean surface---buoyancy is negatively proportional to the density anomaly generated by a change in temperature and salinity. At time $t=0$, the buoyancy decreases linearly with depth. After some time, in response to buoyancy loss at the surface, a well-mixed buoyancy layer develops in the upper ocean. In the absence of entrainment through the base of the mixed layer, the buoyancy profile smoothly connects to the stratified interior. Entrainment results in a further deepening of the boundary layer and the development of a sharp buoyancy jump at its base.
• Figure 2: Large-eddy simulations (LES) for free convection and pure wind stress scenarios in a horizontally doubly-periodic domain of size $(L_x, L_y, L_z) = (128, 128, 128) \, m$ with a grid resolution of $0.5m$ and a Coriolis parameter of $f = 8d-5s^{-1}$. The top row shows a convective turbulence LES driven by surface cooling, while the bottom row shows a shear turbulence LES driven by surface wind stress. The first and second columns show snapshots of the buoyancy and vertical velocity fields. The third and fourth columns show the time evolution of the horizontally-averaged temperature and salinity profiles characterized by a deepening mixed layer in response to the surface forcings.
• Figure 3: The top panel shows the calibrated diffusivity and viscosity values as a function of the local gradient Richardson number $Ri$ in the base closure model. In the convective range, the viscosity and diffusivity are constant, while they decrease linearly towards a background value in the shear range. The lower panels show the vertical profiles of momentum, temperature, salinity, and buoyancy from LES simulations forced with four different air-sea fluxes and the corresponding predictions from the base closure (the dashed lines represent the initial profiles). From left to right: a wind-forced example with no rotation and a convection example with rotation used to train the base closure; a wind-forced example with rotation and a wind-forced example with heating and precipitation used for validation.
• Figure 4: Training and validation losses of NORi. Top row: mean and individual losses against epoch over the final integration horizon of $43.3hours$ for the training suite (Tables \ref{['table free convection LES training suite']}, \ref{['table wind + convection LES training suite']}, and \ref{['table misc LES training suite']}) in the left panel and the validation suite (Table \ref{['table LES validation suite']}) in the right panel. All losses are normalized only once at epoch 0—an epoch is defined as one iteration over the entire training suite. The gray lines are the individual losses for each training/validation case, while the colored lines are their means. Bottom row: mean training and validation losses over the three stages of "curriculum learning," where each stage covers a longer integration period: $15hours$, $23.3hours$, and $43.3hours$. The training loss is re-normalized at the start of each stage, as explained in \ref{['section appendix loss scalings']}. The epoch with the lowest training loss at each stage is used to initialize network weights for the next stage. The final model weights are selected from the model with the lowest validation loss in the final stage.
• Figure 5: Temperature, salinity, and buoyancy profiles for selected training (columns 1 through 3) and validation (columns 4 through 6) examples generated with the calibrated NORi model (black lines), the base closure (orange lines) as well as area-averaged large-eddy simulation (LES) solutions (green lines). The profiles are computed $1.75days$ after the initial conditions (dashed lines).