Mesoscale eddies modulate the stratification, mixing and dissipation pathways, and tracer transport of oceanic flows over a wide range of spatiotemporal scales. The parameterization of buoyancy and momentum fluxes associated with mesoscale eddies thus presents an evolving challenge for ocean modelers, particularly as modern climate models approach eddy-permitting resolutions. Here we present a parameterization targeting such resolutions through the use of a subgrid mesoscale eddy kinetic energy budget (MEKE) framework. Our study presents two novel insights: (1) both the potential and kinetic energy effects of eddies may be parameterized via a kinetic energy backscatter, with no Gent-McWilliams along-isopycnal transport; (2) a dominant factor in ensuring a physically-accurate backscatter is the vertical structure of the parameterized momentum fluxes. We present simulations of 1/2$^\circ$ and 1/4$^\circ$ resolution idealized models with backscatter applied to the equivalent barotropic mode. Remarkably, the global kinetic and potential energies, isopycnal structure, and vertical energy partitioning show significantly improved agreement with a 1/32$^\circ$ reference solution. Our work provides guidance on how to parameterize mesoscale eddy effects in the challenging eddy-permitting regime.

Regional patterns of sea level rise are affected by a range of factors including glacial melting, which has occurred in recent decades and is projected to increase in the future, perhaps dramatically. Previous modeling studies have typically included fluxes from melting glacial ice only as a surface forcing of the ocean or as an offline addition to the sea surface height fields produced by climate models. However, observational estimates suggest that the majority of the meltwater from the Antarctic Ice Sheet actually enters the ocean at depth through ice shelf basal melt. Here we use simulations with an ocean general circulation model in an idealized configuration. The results show that the simulated global sea level rise pattern is sensitive to the depth at which Antarctic meltwater enters the ocean. Further analysis suggests that the response is dictated primarily by the steric response to the depth of the meltwater flux.

There is growing interest in discovering interpretable, closed-form equations for subgrid-scale (SGS) closures/parameterizations of complex processes in Earth system. Here, we apply a common equation-discovery technique with expansive libraries to learn closures from filtered direct numerical simulations of 2D forced turbulence and Rayleigh-Benard convection (RBC). Across common filters, we robustly discover closures of the same form for momentum and heat fluxes. These closures depend on nonlinear combinations of gradients of filtered variables (velocity, temperature), with constants that are independent of the fluid/flow properties and only depend on filter type/size. We show that these closures are the nonlinear gradient model (NGM), which is derivable analytically using Taylor-series expansions. In fact, we suggest that with common (physics-free) equation-discovery algorithms, regardless of the system/physics, discovered closures are always consistent with the Taylor-series. Like previous studies, we find that large-eddy simulations with NGM closures are unstable, despite significant similarities between the true and NGM-predicted fluxes (pattern correlations > 0.95). We identify two shortcomings as reasons for these instabilities: in 2D, NGM produces zero kinetic energy transfer between resolved and subgrid scales, lacking both diffusion and backscattering. In RBC, backscattering of potential energy is poorly predicted. Moreover, we show that SGS fluxes diagnosed from data, presumed the ‘truth’ for discovery, depend on filtering procedures and are not unique. Accordingly, to learn accurate, stable closures from high-fidelity data in future work, we propose several ideas around using physics-informed libraries, loss functions, and metrics. These findings are relevant beyond turbulence to closure modeling of any multi-scale system.

The understanding and representation of energetic transfers associated with ocean mesoscale eddies is fundamental to the development of parameterizations for climate models. We investigate the influence of eddies on flow vertical structure as a function of underlying dynamical regime and grid resolution. We employ the GFDL-MOM6 in an idealized configuration and systematically consider four horizontal resolutions: 1/4, 1/8, 1/16, and 1/32 degree. We analyze the distributions of potential and kinetic energy, decomposed into barotropic and baroclinic, and eddy and mean parts. Kinetic energy increases and potential energy decreases as resolution increases and captures more baroclinically-unstable modes. The dominant trend in vertical structure is an increasing fraction of kinetic energy going into the barotropic mode, particularly its eddy component, as eddies are increasingly resolved. We attribute the increased baroclinicity at low resolutions to inaccurate representation of vertical energy fluxes, leading to suppressed barotropization and energy trapping in high vertical modes. We also explore how the underlying dynamical regime influences energetic pathways. In cases where large-scale flow is dominantly barotropic, resolving the deformation radius is less crucial to accurately capturing the flow’s vertical structure. We find the barotropic kinetic energy fraction to be a useful metric in assessing vertical structure. In the highest-resolution case, the barotropic kinetic energy fraction correlates with the scale separation between the deformation scale and the energy-containing scale, i.e. the extent of the eddy-driven inverse cascade. This work suggests that mesoscale eddy parameterizations should incorporate the energetic effects of eddies on vertical structure in a scale-aware, physically-informed manner.

The resolution of climate models is limited by computational cost. Therefore, we must rely on parameterizations to represent processes occurring below the scale resolved by the models. Here, we focus on parameterizations of ocean mesoscale eddies and employ machine learning (ML), namely relevance vector machines (RVM) and convolutional neural networks (CNN), to derive computationally efficient parameterizations from data, which are interpretable and/or encapsulate physics. In particular, we demonstrate the usefulness of the RVM algorithm to reveal closed-form equations for eddy parameterizations with embedded conservation laws. When implemented in an idealized ocean model, all parameterizations improve the statistics of the coarse-resolution simulation. The CNN is more stable than the RVM such that its skill in reproducing the high-resolution simulation is higher than the other schemes; however, the RVM scheme is interpretable. This work shows the potential for new physics-constrained interpretable ML turbulence parameterizations for use in ocean climate models.

There is large uncertainty in the future sea level change at regional scales under anthropogenic global warming. This study uses a novel design of ocean-only general circulation model (OGCM) experiments to investigate the ocean’s response to surface buoyancy and momentum flux perturbations, as part of the Flux-Anomaly-Forced Model Intercomparison Project (FAFMIP), and compares with results from coupled, atmosphere-ocean GCM (AOGCM) experiments. Much of the inter-model spread is driven by the response to surface heat flux perturbations. In a multi-model ensemble of OGCMs forced with identical surface heat flux perturbations, regional sea level and ocean heat content changes demonstrate considerable disagreement, especially in the North Atlantic. Spread in both residual mean advection and diapycnal diffusion changes contribute to much of the multi-model disagreement over regional heat content change. Residual mean advection changes are related to the large spread in simulated Atlantic meridional overturning circulation (AMOC) weakening (20-50%). We find approximately 10% more AMOC weakening in response to surface heat flux perturbations in AOGCMs relative to OGCMs with consistent ocean models. This enhanced AMOC weakening is driven by an atmosphere-ocean feedback which amplifies the surface heat flux perturbation. In the North Pacific, there is little agreement amongst the ensemble over which processes lead to ocean warming, with varying contributions from residual mean advection and diapycnal diffusion. For the Pacific basin, the atmosphere-ocean feedback reduces sea surface temperature (SST) warming by 0.5°C. In the Southern Ocean, the atmosphere-ocean feedback is not generally important for buoyancy and momentum flux perturbations.

Recently, a growing number of studies have used machine learning (ML) models to parameterize computationally intensive subgrid-scale processes in ocean models. Such studies typically train ML models with filtered and coarse-grained high-resolution data and evaluate their predictive performance offline, before implementing them in a coarse resolution model and assessing their online performance. In this work, we systematically benchmark the online performance of such models, their generalization to domains not encountered during training, and their sensitivity to dataset design choices. We apply this proposed framework to compare a large number of physical and neural network (NN)-based parameterizations. We find that the choice of filtering and coarse-graining operator is particularly critical and this choice should be guided by the application. We also show that all of our physics-constrained NNs are stable and perform well when implemented online, but generalize poorly to new regimes. To improve generalization and also interpretability, we propose a novel equation-discovery approach combining linear regression and genetic programming with spatial derivatives. We find this approach performs on par with neural networks on the training domain but generalizes better beyond it. We release code and data to reproduce our results and provide the research community with easy-to-use resources to develop and evaluate additional parameterizations.

Coupled climate simulations that span several hundred years cannot be run at a high-enough spatial resolution to resolve mesoscale ocean dynamics. These mesoscale dynamics backscatter to macroscales. Recently, several studies have considered Deep Learning to parameterize subgrid forcing within macroscale ocean equations using data from idealized simulations. In this manuscript, we present a stochastic Deep Learning parameterization that is trained on data generated by CM2.6, a high-resolution state-of-the-art coupled climate model with nominal resolution 1/10° . We train a Convolutional Neural Network for the subgrid momentum forcing using macroscale surface velocities from a few selected subdomains. At each location and each time step of the coarse grid, rather than predicting a single number, we predict the mean and standard deviation of a Gaussian probability distribution. This approach requires training our neural network to minimize a negative log-likelihood loss function rather than the Mean Square Error, which has been the standard in applications of Deep Learning to the problem of parameterizations. Each prediction of the mean subgrid forcing can be associated with an uncertainty estimate and can form the basis for a stochastic subgrid parameterization. Offline tests show that our parameterization generalizes well to the global oceans, and a climate with increased CO2 levels, without further training. We test our stochastic parameterization in an idealized shallow water model. The implementation is stable and improves some statistics of the flow. Our work demonstrates the potential of combining Deep Learning tools with a probabilistic approach in parameterizing unresolved ocean dynamics.

We describe an idealized primitive equation model for studying mesoscale turbulence and leverage a hierarchy of grid resolutions to make eddy-resolving calculations on the finest grids more affordable. The model has intermediate complexity, incorporating basin-scale geometry with idealized Atlantic and Southern oceans, and with non-uniform ocean depth to allow for mesoscale eddy interactions with topography. The model is perfectly adiabatic and spans the equator, and thus fills a gap between quasi-geostrophic models, which cannot span two hemispheres, and idealized general circulation models, which generally have diabatic processes and buoyancy forcing. We show that the model solution is approaching convergence in mean kinetic energy for the ocean mesoscale processes of interest, and has a rich range of dynamics with circulation features that emerge only due to resolving mesoscale turbulence.