A Data-Driven Approach for Parameterizing Submesoscale Vertical Buoyancy Fluxes in the Ocean Mixed Layer

Abstract

The parameterizations of submesoscale ( $<10$ km) ocean surface flows are critical in capturing the subgrid effects of vertical fluxes in the ocean mixed layer, yet they struggle to infer the full-complexity of these fluxes in relation to the large scale variables that help set them. In this work, we present a data-driven approach for the submesoscale parameterization, utilizing information from the high-resolution submesoscale-permitting MITgcm-LLC4320 simulation (LLC4320). The new parameterization is given by a Convolutional Neural Network (CNN) trained to infer the subgrid mixed layer vertical buoyancy fluxes as a function of relevant large scale variables. In contrast to previous physics-based approaches, such as the Mixed Layer Eddy (MLE) parameterization, here the CNN infers vertical fluxes that are directly computed from the LLC4320 data, where the submesoscales are resolved down to a resolution of approximately 2km. The CNN has significantly high skill compared with the MLE parameterization, which we demonstrate over a wide range of dynamical regimes and resolutions. We find that the improved skill can be attributed to learned physical relationships between submesoscale fluxes and the large scale strain field, currently missing from submesoscale parameterizations in General Circulation Models.

\stackMath\draftfalse\journalname

Journal of Advances in Modeling Earth Systems (JAMES)

Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA Lamont-Doherty Earth Observatory, Columbia University, New York, NY, USA Courant Institute of Mathematical Sciences, New York University, New York, NY, USA

\correspondingauthor

A. [email protected]

{keypoints}

We improve the parameterization of submesoscale vertical buoyancy fluxes using a Convolutional Neural Network (CNN).

The CNN demonstrates high offline skill compared with the physics-based submesoscale parameterization over a wide range of dynamical regimes, seasonality, and resolutions.

We identify strong dependency on the large scale strain field, which is currently missing from submesoscale parameterizations implemented in General Circulation Models.

Plain Language Summary

This works provides a data-driven method for inferring small scale buoyancy fluxes in the upper ocean given from high resolution simulation output. Our method provides high skill compared with the estimates currently used in climate models, suggesting that information from the large scale strain field is key to reducing model biases.

1 Introduction

General Circulation Models (GCMs) and future climate change projections are notoriously sensitive to parameterizations of unresolved phenomena at the ocean-atmosphere interface [IPCC (\APACyear2019), IPCC (\APACyear2021)]. Of particular importance is the ocean mixed layer, where multiscale turbulence modulates the transfer of properties – such as heat, momentum, and carbon – between the atmosphere and ocean interior [<]e.g.,¿[]frankignoul1977stochastic,bopp2015pathways,su2020high,fox2022ocean. Ocean surface submesoscales flows typically appear on the time scale of hours to days and $O(1)$ km in spatial scale [<]e.g.,¿[]mcwilliams2016submesoscale. Submesoscale flows are sandwiched between mesoscale $O(100)$ km flows and boundary layer turbulence $O(1)$ m, and play an important role in vertical transport in the ocean mixed layer. As opposed to mixing and homogenization of mixed layer buoyancy dominated by boundary layer turbulence, submesoscale flows tend to exhibit a positive vertical buoyancy flux and thus have a net restratification effect [Boccaletti \BOthers. (\APACyear2007), Thomas \BOthers. (\APACyear2008), Mahadevan (\APACyear2016), Johnson \BOthers. (\APACyear2016), McWilliams (\APACyear2016), Balwada \BOthers. (\APACyear2021), Taylor \BBA Thompson (\APACyear2023)]. Such positive buoyancy fluxes are at leading order a result of mixed layer instabilities formed along submesoscale fronts, with dynamics attributed to multiscale interactions between the mesoscale flow field and boundary layer turbulence [Fox-Kemper \BOthers. (\APACyear2008), McWilliams (\APACyear2016), Bachman \BBA Klocker (\APACyear2020), Gula \BOthers. (\APACyear2022), Bodner \BOthers. (\APACyear2023)]. It is thus a sensitive interplay between mesoscale flows, submesoscale flows, boundary layer turbulence, and surface forcing that set the stratification in the ocean mixed layer, for which GCM parameterizations struggle to capture in its entirety.

The effects of submesoscale vertical buoyancy fluxes are inferred in many GCMs by the Mixed Layer Eddy (MLE) parameterization [Fox-Kemper \BOthers. (\APACyear2011), Calvert \BOthers. (\APACyear2020), Bodner \BOthers. (\APACyear2023)], which represents the positive vertical buoyancy flux produced by mixed layer instabilities formed along submesoscale fronts [Fox-Kemper \BOthers. (\APACyear2008)]. The MLE parameterization is cast in the form of a streamfunction [<]e.g.,¿[]gent1990isopycnal,Griffies98, ferrari2008parameterization, which represents the net restratification effect prescribed by submesoscale eddy fluxes. Frontal sharpening, known as frontogenesis, helps set the buoyancy gradient from which the mixed layer instabilities are formed. Frontogenesis also contributes vertical buoyancy flux at scales typically smaller than MLE [Hoskins \BBA Bretherton (\APACyear1972), Lapeyre \BOthers. (\APACyear2006), Shakespeare \BBA Taylor (\APACyear2013), Bodner \BOthers. (\APACyear2020), Gula \BOthers. (\APACyear2022), Srinivasan \BOthers. (\APACyear2023)]. The implementation of the MLE parameterization in GCMs [Fox-Kemper \BOthers. (\APACyear2011)] involves a rescaling factor to enhance the subgrid buoyancy gradient, which tends to be suppressed by the smoother GCM [<]e.g.,¿[]stanley2020parameterizing. The frontal width rescaling factor is taken to be proportional to the mixed layer deformation radius [Fox-Kemper \BOthers. (\APACyear2011), Calvert \BOthers. (\APACyear2020)], yet this scaling is found to not satisfactorily hold when frontogenesis or strong surface forcing are present [<]e.g.,¿[]lapeyre2006oceanic,mahadevan2010rapid,callies2018baroclinic. \citeAbodner2023modifying had recently modified the MLE parameterization with a new scaling for frontal width set based on turbulence-induced frontogenesis, linking the MLE parameterization with the Energetics-based Planetary Boundary Layer scheme [<]ePBL;¿reichl2018simplified.

Recent developments in high resolution numerical modeling have presented a new class of submesoscale permitting simulations with kilometer scale resolution [<]e.g.,¿[]rocha2016mesoscale, schubert2019submesoscale,ajayi2021diagnosing. In a comprehensive comparison across eight different submesoscale permitting simulations, \citeAuchida2022cloud find that the MLE parameterization struggles to capture the full range of complexity and magnitude given by submesoscale vertical buoyancy fluxes, even when submesoscale fronts are resolved and the rescaling factor is not necessary. Studies have suggested that such missing physics can be attributed to the relevance of the large scale flow in setting submesoscale vertical buoyancy fluxes, namely, properties of strain, divergence, and vorticity [<]e.g.,¿[]lapeyre2006oceanic,capet2008mesoscale,mcwilliams2017submesoscale, barkan2019role, balwada2021vertical. \citeAzhang2023parameterizing propose a parameterization for submesoscale vertical buoyancy flux enhanced by the effects of the large scale strain field, yet does not include a relationship with boundary layer turbulence. Furthermore, the interaction of submesoscale fronts with surface forcing, and in particular in regards to surface cooling or wind orientation, can induce an Ekman Buoyancy Flux, which emerges as a central component to setting mixed layer stratification, serving as another source for frontogensis, and setting Potential Vorticity conditions for Symmetric Instability [Thomas (\APACyear2005), Thomas \BOthers. (\APACyear2013), Bachman \BOthers. (\APACyear2017), Wenegrat \BBA Thomas (\APACyear2020)]. Thus, even the most recent advancements in submesoscale parameterizations present incomplete physical relationships between submesoscale vertical buoyancy fluxes and the large scale variables that help set it.

Data-driven methods are emerging as powerful tools, with the ability to capture highly complex relationships between variables in turbulent flows. Advances in machine learning parameterizations have yielded promising results in improving physics for subgrid closures. The seminal work by \citeAbolton2019applications demonstrated that Neural Networks (NNs) are able to learn subgrid ocean mesoscale momentum fluxes, while remaining physically-constrained. \citeAzanna2020data and \citeAguillaumin2021stochastic respectively trained a NN to learn the underlying equation and statistics of subgrid mesoscale fluxes as given from data. An example for ocean boundary layer turbulence closures can be found in \citeAsouza2020uncertainty, whom apply a Bayesian approach to estimate empirical parameters in the K-Profile Parameterization [Large \BOthers. (\APACyear1994)] from an ensemble of large eddy simulations. A complimentary approach has been taken in \citeAsane2023parameterizing, where empirical parameters from ePBL are replaced with a NN constrained by observations. In the atmosphere, several comparable approaches have trained NNs to infer vertical buoyancy fluxes in the atmospheric boundary layer [<]e.g.,¿[]yuval2020stable, wang2022non,shamekh2023implicit. Numerous examples for applications exist both in the atmosphere and the ocean for use of machine learning for inference of flow patterns and structures from data [<]e.g.,¿[]chattopadhyay2020predicting, dagon2022machine, xiao2023reconstruction, zhu2023deep.

In this work, we introduce a data-driven approach for parameterizing submesoscale- vertical buoyancy fluxes. The main goal is to train a Convolutional Neural Network (CNN) on realistic simulated data to learn an improved functional relationships between submesoscale vertical buoyancy fluxes and the large scale physics that help set it. The data used to train and test the CNN is sampled from the global ocean MITgcm-LLC4320 [<]hereafter LLC4320,¿menemenlis2021pre simulated at a high resolution of $1/48^{o}$ , developed ahead of the Surface Water and Ocean Topography (SWOT) altimeter mission. The LLC4320 output has been widely studied for submesoscale applications, which cumulatively have demonstrated that submsoscale energetics and dynamics are captured relatively well down to its effective resolution [<]e.g.,¿[]rocha2016mesoscale, rocha2016seasonality, su2018ocean,torres2018partitioning, dong2020seasonality,uchida2022cloud, gallmeier2023evaluation. We describe the LLC4320 data, preprocessing, and choices of CNN inputs and output in section 2. Results of the CNN prediction on unseen data are presented and compared with the MLE parameterization in section 3, together with sensitivity and generalization tests for robustness. In section 4, we apply two complimentary methods to explain the relationship learned by the CNN between the physical variables used as inputs and the submesoscale vertical buoyancy flux. Discussion and concluding remarks are given in section 5.

2 Data and methods

2.1 Designing the learning problem

We choose to include input features that are correlated, or have known analytical relationships, with submesoscale vertical buoyancy fluxes. We thus leverage the physical relevance demonstrated by the MLE parameterization, and extract variables which appear in the \citeAfox2011parameterization and \citeAbodner2023modifying formulations. We also include correlated large scale velocity derivatives as demonstrated in \citeAbarkan2019role, balwada2021vertical,zhang2023parameterizing. The input features (table 1) consist of the depth-averaged horizontal buoyancy gradient magnitude, $|\overline{\nabla b}^{z}|$ , where overbar represents the filtering operator described in section 2.2, and superscript $z$ represents a depth averaging operator over the mixed layer depth; Coriolis parameter, $\overline{f}$ ; mixed layer depth, $\overline{H_{ML}}$ ; surface heat flux , $\overline{Q^{*}}$ ; surface wind stress magnitude, ${|\tau|}=\sqrt{\overline{\tau_{x}}^{2}+\overline{\tau_{y}}^{2}}$ ; boundary layer depth, $\overline{H_{BL}}$ ; depth-averaged strain magnitude, $\overline{\sigma}^{z}=\sqrt{(\overline{u_{x}}^{z}-\overline{v_{y}}^{z})^{2}+(\overline{v_{x}}^{z}+\overline{u_{y}}^{z})^{2}}$ ; depth-averaged vertical vorticity, $\overline{\zeta}^{z}=\overline{v_{x}}^{z}-\overline{u_{y}}^{z}$ ; depth-averaged horizontal divergence, $\overline{\delta}^{z}=\overline{u_{x}}^{z}+\overline{v_{y}}^{z}$ . Note that we exclude the Brunt-Väisälä frequency, $\overline{N^{2}}$ , [Fox-Kemper \BOthers. (\APACyear2008)] as its average quantity in the mixed layer contains identical information as $H_{ML}$ , and thus would not contribute to the CNN training.

Formally, we define our 9 input features as,

\mathbf{X}=(|\overline{\nabla b}^{z}|,\overline{f},\overline{H_{ML}},\overline{Q^{*}},{|\tau|},\overline{H_{BL}},\overline{\sigma}^{z},\overline{\zeta}^{z},\overline{\delta}^{z}),

(1)

and a single output as,

\mathbf{Y}=\overline{w}^{z}\overline{b}^{z}-\overline{wb}^{z}

(2)

For compactness, we will frequently refer to the output as $\overline{w^{\prime}b^{\prime}}^{z}:=\overline{w}^{z}\overline{b}^{z}-\overline{wb}^{z}$ , where $()^{\prime}$ represents the submesoscale. The CNN provides the subgrid fluxes as a function, $S$ , of the large scale variables, such that $S(\mathbf{X})\rightarrow\mathbf{Y}$ . Figure 1 illustrates a schematic of the CNN with 9 input features and 1 output. The processing involved in computing all inputs and output from the LLC4320 is described in more detail below.

Refer to caption — Figure 1: Schematic of CNN method with 9 input features and one output feature (table 1) at a $1/4^{o}$ resolution.

CNN inputs, $\mathbf{X}$
Depth-averaged buoyancy gradient magnitude	$\|\overline{\nabla b}^{z}\|$
Coriolis parameter	$\overline{f}$
Mixed layer depth	$\overline{H_{ML}}$
Surface heat flux	$\overline{Q^{*}}$
Surface wind stress magnitude	$\overline{\|\tau\|}=\sqrt{\overline{\tau_{x}}^{2}+\overline{\tau_{y}}^{2}}$
Boundary layer depth	$\overline{H_{BL}}$
Depth-averaged strain magnitude	$\overline{\sigma}^{z}=\sqrt{(\overline{u_{x}}^{z}-\overline{v_{y}}^{z})^{2}+(\overline{v_{x}}^{z}+\overline{u_{y}}^{z})^{2}}$
Depth-averaged vertical vorticity	$\overline{\zeta}^{z}=\overline{v_{x}}^{z}-\overline{u_{y}}^{z}$
Depth-averaged horizontal divergence	$\overline{\delta}^{z}=\overline{u_{x}}^{z}+\overline{v_{y}}^{z}$
CNN Output, $\mathbf{Y}$
Depth-averaged subgrid vertical buoyancy flux	$\overline{w^{\prime}b^{\prime}}^{z}:=\overline{w}^{z}\overline{b}^{z}-\overline{wb}^{z}$

Table 1: Input and output features used in the CNN method. Overbar represents the filtering operator described in 2.2, and superscript

z

represents a depth averaging operator over the mixed layer depth.

2.2 Input and output features

The LLC4320 is a $1/48^{o}$ Massachusetts Institute of Technology general circulation model (MITgcm), named after its Latitude‐Longitude polar Cap (LLC) grid with 4320 points on each of the 13 tiles. The LLC4320 is initialized from the Estimating the Circulation and Climate of the Ocean (ECCO), Phase II project, and is forced at the surface by atmospheric reanalysis, at 6 hourly temporal resolution. The LLC4320 simulation also includes a synthetic surface pressure field to mimic tidal forcing present in the real ocean. A total of 14 months are available at hourly frequency from September 2011 to November 2012 [Forget \BOthers. (\APACyear2015), Menemenlis \BOthers. (\APACyear2008), Menemenlis \BOthers. (\APACyear2021)].

Figure 2 illustrates a snapshot of the global surface horizontal buoyancy gradient from LLC4320 coarsened to $1/4^{o}$ , with boxes highlighting the sampled regions used in this study. Table 2 lists the coordinates of each of the 12 sampled regions of approximately $15^{o}\times 15^{o}$ . Buoyancy gradients are a key variable in the MLE parameterization, which helps motivate the choices of regions selected for training. We include a mix of regions with strong variability and others more quiescent [<]e.g.,¿[]torres2018partitioning in effort to train the CNN on a range of dynamical regimes influenced by seasonality and local flow properties.

Region	Latitudinal Range	Longitudinal Range
1. Gulf Stream	(30,45)	(-60,-45)
2. South Atlantic	(-30,-15)	(-25,-10)
3. Equator Atlantic	(-8,8)	(-30,-15)
4. Malvinas Current	(-55,-40)	(-60,-45)
5. California Current	(30,45)	(-140,-125)
6. South Pacific	(-45,-30)	(-140,-125)
7. Kuroshio Extension	(25,40)	(145,160)
8. North Pacific	(10,25)	(-180,-165)
9. Southern Ocean, New Zealand	(-60,-45)	(-175,-160)
10. Agulhas Current	(-55,-40)	(20,35)
11. Indian Ocean	(-25,-10)	(70,85)
12. Arabian Sea	(0,15)	(55,70)

Table 2: Coordinate range of sampled regions from the LLC4320 corresponding to the blue boxes in figure 2.

From each of the regions listed in table 2, the LLC4320 output is preprocessed in preparation of the CNN training, in which we compute the CNN input and output features listed in table 1.

In order to filter the synthetic fast-varying wave field from the LLC4320 output, the temporal resolution of each variable is converted from hourly into an average over 12 hour window periods. The potential density anomaly $\sigma_{0}$ is computed from the potential temperature and salinity fields, with reference pressure of 0 dbar and $\rho_{0}=1000$ kg m^-3. The mixed layer depth, H_ML, is defined as the depth at which $\sigma_{0}$ increased by $0.03$ kg m^-3 from its value at 10m depth [de Boyer Montégut \BOthers. (\APACyear2004)]. The buoyancy field is then defined as $b=-g\sigma_{0}/\rho_{0}$ , where $g=9.81$ m/s² is the gravity acceleration. Velocities and wind stresses are interpolated to commute with the tracer grid. The boundary layer depth $H_{BL}$ is already an LLC4320 output computed from the K‐profile Parameterization criteria for which the Richardson number exceeds the critical value of 0.3 [<]KPP,¿[]large1994oceanic.

We restrict our method to learning the depth-averaged submesoscale vertical buoyancy fluxes. This approach remains close to the MLE parameterization, which was developed as a bulk formula for the depth-averaged submesoscale vertical buoyancy fluxes [Fox-Kemper \BOthers. (\APACyear2011), Bodner \BOthers. (\APACyear2023)], and a vertical structure function determines the shape of the parameterization at depth (eq. (S4) in the supplementary material). To this end, all 3D variables are averaged over the mixed layer depth, which we denote by superscript $z$ . Horizontal filtering is applied to all input and output variables to bring the high-resolution, fine-scale LLC4320 data, to the low-resolution large-scale, such that it is comparable to a GCM resolution. We have tested our method over resolutions of $1^{o},{1/2}^{o},{1/4}^{o},{1/8}^{o},{1/12}^{o}$ . To achieve this, we include both spatial and temporal filtering. In space, we we apply a top-hat filter. As an example for the ${1/4}^{o}$ resolution case, the top-hat filter is an average over $12$ grid points of the original LLC4320 grid [Loose \BOthers. (\APACyear2022)]. An additional 7-day temporal filter is applied such that all large-scale variables are assumed to vary on timescales larger than a week, e.g., larger than the submesoscale. Note that the 7-day temporal filter is larger than the inertial period with the added intention of reducing noise from tides and internal waves that may have otherwise contaminated the sampled data [Torres \BOthers. (\APACyear2018), Jones \BOthers. (\APACyear2023)]. To review, for any variable $\xi$ , the filtering procedure is defined as both in horizontal space $<\xi>_{h}$ and in time $\widetilde{\ \ \xi\ \ }^{t}$ . For compactness, we denote the combination of both averaging operators as $\overline{\xi}:=\widetilde{<\xi>}^{t}_{h}$ , which is applied to all input and output features. If $\xi$ is a three-dimensional variable, we also include the superscript $\overline{\xi}^{z}$ which corresponds to a vertical average over the mixed layer depth.

Thus, the sub-grid process we are parameterizing is defined by a combined temporal and spatial filter, which is the fast-varying (within the 12 hour to 7 day range) and small-scale (below the spatial filter scale) mixed layer buoyancy flux. This quantity is designed to be the CNN output, formally written as the total turbulent stress tensor, $\overline{w^{\prime}b^{\prime}}^{z}:=\overline{w}^{z}\overline{b}^{z}-\overline{wb}^{z}$ . The CNN output of mixed-layer-averaged vertical buoyancy fluxes is thus akin to an estimate of submesoscale fluxes. We demonstrate this in the variance-preserving co-spectrum of $w$ and $b$ in the mixed layer, shown in figure 3. It can be seen that the maximum covariance is predominantly in the submesoscale range. Special cases indicate that more complex dynamics influence mixed layer buoyancy fluxes, either by strong mesoscale activity, such as in the Arabian Sea, or by turbulent activity, such as in the Equatorial Atlantic. Naturally, the filter scale (illustrated by the grey lines in figure 3b) will impact the properties of the learned subgrid flux, which we explore in section 3. We have included all regions to gain a variety of dynamical regimes in our training data. In sections 3.1 and 3.2 we test the ability of the CNN to generalize over different seasons and regions, respectively.

2.3 CNN architecture and training

All features listed in table 1 are normalized by a global mean and standard deviation computed over all regions. To train the CNN, we randomly select 80% of the, approximately, 10,000 samples given from all regions combined, after temporal and spatial filtering is applied. The remaining 20% is left unseen by the CNN and is used only to test the prediction of the trained CNN and compare with the target LLC4320 data and MLE parameterization in section 3 below.

We use a CNN architecture for regression inspired by applications for mesoscale eddy parameterizations [Bolton \BBA Zanna (\APACyear2019), Guillaumin \BBA Zanna (\APACyear2021), Perezhogin \BOthers. (\APACyear2023)]. A hyper-parameter sweep over the number of hidden layers, kernel size, learning rate, and weight decay, was used to find the best performing CNN. The CNN is trained over 100 epochs while minimizing the Mean Squared Error (MSE) loss (shown in figure S2). Results presented here are based on a CNN with a kernel size of 5X5 in the first layer, 7 hidden convolutional layers with kernel size of 3x3, a learning rate of $2\times 10^{-4}$ , and weight decay of $0.02$ , which gives an MSE of $0.23$ in normalized units. The total number of learnable parameters is approximately $300,000$ .

3 CNN prediction of subgrid submesoscale fluxes

Once training is complete, the relationship between the CNN inputs and output is optimized over the learnable parameters. In other words, the CNN has learned a functional form, $S(\mathbf{X})$ , between the input features, $\mathbf{X}$ , and submesoscale vertical buoyancy fluxes, $\mathbf{Y}$ . In this section, we examine the extent to which the relationships learned by the CNN generalize onto data that was not included in the training process. For this purpose, we evaluate the skill of the CNN compared with the target LLC4320 data as well as the \citeAbodner2023modifying version of the MLE parameterization. In subsections 3.1 and 3.2 we test whether our results are sensitive to the choice of training on a subset of the timeseries or sampled regions, respectively.

Figure 4 displays a single sample of the $1/4^{o}$ submesoscale vertical buoyancy fluxes, $\overline{w^{\prime}b^{\prime}}^{z}$ , as given by the LLC4320 data (4a), CNN (4b), and MLE parameterization (4c). The sample is from the Arabian Sea region during the month of September. For illustration purposes, the LLC4320 with filter scales corresponding to $1/12^{o},1/8^{o},1/2^{o},1^{o}$ , is shown in figure S3. The majority of the fluxes appear to be positive, which is the bulk restratification effect infered by the MLE parameterization. However, the CNN captures much more of the fine-scale structure as well as the sign, including intermittent negative fluxes.

To examine the statistics beyond a single sample, we compute the joint histogram of the $1/4^{o}$ LLC4320 submesoscale fluxes and those given by the CNN. The joint histogram is computed over the entire unseen test dataset, which contains $20\%$ of random samples over all regions, and provides a metric over several orders of magnitude. The LLC4320-CNN joint histogram is compared with that of the LLC4320 and the MLE parameterization in figure 5. In the case of positive fluxes, the CNN prediction remains close to the LLC4320 data, as can be seen by the alignment along the one-to-one grey line, a significant improvement on the MLE parameterization. For the negative fluxes the one-to-one alignment is less pronounced, likely due to the significantly smaller number of negative samples, however the CNN prediction is an improvement on the MLE parameterization which does not infer negative fluxes by construction.

We examine the ability of the CNN in the $1/4^{o}$ resolution experiment to capture seasonality in figure 6, which shows the seasonal cycle of the spatially-weighted-average of $\overline{w^{\prime}b^{\prime}}^{z}$ decomposed by region. In each panel, the LLC4320 target data is compared with the MLE parameterization and the CNN predictions on unseen data. Once again, in virtually all regions the CNN prediction outperforms the MLE parameterization, particularly where fluxes appear to be strongest during the winter and spring months.

The CNN and MLE parameterization skill in terms of $R^{2}$ values decomposed by region and across all resolution experiments are shown in figure 7. The CNN prediction skill remains above that of the MLE in all filter scale experiments and in all regions. The lowest resolution tends to have the highest skill and it decreases as resolution increases. Note that the MLE parameterization is a bulk formula for $\overline{wb}^{z}$ , which is not as sensitive to the different resolutions due to the frontal rescaling factor. As a reminder, figure 3 illustrates how different filter scales define the amount of variability captured in the data. The variability of $\overline{wb}^{z}$ varies both in scale and by location, which impacts the learned output of the CNN in the different regions. The CNN prediction skill is found to be especially sensitive in the high-resolution experiments, where it performs well in some regions (e.g. the New Zealand region), but not in others (e.g. the Indian Ocean region). In the low resolution experiments, the fields tend to be smoother, as much of the variability is averaged out, thus presenting an easier learning problem for the CNN.

To better understand the dependency of our method on the training data, in the following subsections we perform two sensitivity tests by holding out parts of the training data, retraining the CNN, and examining CNN prediction skill on the unseen regional or seasonal data.

3.1 Holding out seasonality from training

Submesoscale fluxes tend to be strongest during months of winter and spring [<]¿[e.g.]callies2015seasonality,johnson2016global. The strong seasonality is also observable in the climatology presented in figure 6 with respect to the associated hemisphere. We perform two experiments in which we hold out winter and summer months from the training data, to examine the ability of the CNN to generalize on unseen seasonal variability, with skill given in terms of the $R^{2}$ value. We thus create two new training and test datasets to better understand the overall sensitivity of our method to submesoscale seasonality:

•

Winter held out refers to training data which excludes from the time series the months of January, February, March from all regions in the Northern Hemisphere, and July, August, September from regions in the Southern Hemisphere. Note that we have removed equatorial regions from the analysis entirely as the submesoscale equatorial seasonality is less trivial. The remainder of the time series– e.g. spring, summer, fall– is used to train the CNN, and predictions are made on the unseen winter data.
•

Summer held out is same as the above, where we now exclude July, August, September from the Northern Hemisphere and January, February, March from the Southern Hemisphere. Equatorial regions are once again excluded.

The upper panel in figure 8 shows the results of the two seasonality experiments. $R^{2}$ values of the CNN prediction on the unseen summer are found to be relatively skillful across resolutions and over the different regions included in the seasonality experiment. Contrarily, CNN predictions on winter months do not yield as skillful results, with consistently equal or lower $R^{2}$ value. Although not entirely surprising, these results suggest that the CNN is able to generalize and make skillful predictions when the strongest fluxes, generally exhibited during winter and spring months, are included in the training data.

3.2 Holding out regions from training

We perform another set of experiments to test the ability of the CNN to generalize on regions that are not included in the training data. We thus generate 12 new datasets that correspond to removing one region at a time from the training dataset. We retrain the CNN in 12 different experiments, and make predictions on a different unseen region each time.

The lower panel in figure 8 illustrates that $R^{2}$ values of the CNN on the unseen regional data remain consistent with those found on the full training set (figure 7) across resolutions and over all regions. This suggests that the training data covers a wide enough range of dynamical regimes that enables generalization of the CNN on regions not included in training, an especially important result given that a fairly small number of regions were included in training compared with the full ocean. Thus, the learned relationships between the input features and $\overline{w^{\prime}b^{\prime}}^{z}$ by the CNN are realizable across dynamical regimes, especially in the low resolution experiments.

4 Feature importance

We have shown that the CNN improves on the MLE parameterization, but an important remaining question is why? What relationships are learned between the physical variables used as inputs and $\overline{w^{\prime}b^{\prime}}^{z}$ that lead to better predictions by the CNN? With such complex and nonlinear relationships, it is difficult to decipher which input feature is most important and why. Many methods exist that help explain and interpret the dependency of CNN outputs to its inputs. Here, we have chosen two complimentary methods that help gain insight on the learned local and non-local relationships and the importance of individual inputs to $\overline{w^{\prime}b^{\prime}}^{z}$ .

4.1 Impact of input feature on CNN prediction skill

To test the dependency of the CNN prediction on certain input features, we perform a set of ablation experiments, where we remove one input feature at a time, retrain the CNN, and examine the resulting prediction skill in terms of $1-R^{2}$ value. If removing a certain input feature results in a very low $R^{2}$ value, or high $1-R^{2}$ value, it indicates that the CNN strongly depends on said input feature. The top panel in figure 9 displays the $1-R^{2}$ value given from the ablation experiments for each input features, and compared across all resolutions.

Notably, removing strain, $\overline{\sigma}^{z}$ , as an input feature results in high $1-R^{2}$ values, demonstrating the largest reduction in skill consistently across all resolutions. Interestingly, there appears to be very little sensitivity to the removal of any other input feature, including those used by the MLE parameterization. These results strongly suggest that the primary reason the CNN predictions surpass those of the MLE parameterization, are due to the newly-captured non-local relationship between submesoscale vertical buoyancy fluxes and the large scale strain field. Note that the relevance of strain to submesoscale fluxes is not a new result [<]e.g.,¿[]balwada2021vertical,zhang2023parameterizing, but these findings emphasize the relevance of strain to improving submesoscale $\overline{w^{\prime}b^{\prime}}^{z}$ parameterizations.

4.2 Sensitivity of output to relative to input features

We next apply a complimentary method to the ablation experiment above. We compute the Jacobian, defined as the gradients of the CNN output with respect to the input features, $\nabla_{\mathbf{X}}S(\mathbf{X})$ . Note that unlike the ablation experiment, where we examined the CNN skill on the full output domain, the Jacobian is a local metric, considering only the sensitivity of a single output point to a single point in the input feature map. The Jacobian is especially useful to evaluate the point-wise sensitivity of the output to each input feature by taking gradients along the CNN weights [<]e.g.,¿[]ross2023benchmarking. We compute the Jacobian over the entire unseen test dataset, and examine its average values for each input feature, thus providing a metric for how sensitive, on average, the CNN output is to each input feature. We contrast the Jacobian with the $R^{2}$ values given by the ablation experiments in figure 9, where for the Jacobian, a high score indicates that the CNN output, $\overline{w^{\prime}b^{\prime}}^{z}$ , is sensitive to local changes in a certain input. We find that the highest-ranked input feature, for which $\overline{w^{\prime}b^{\prime}}^{z}$ is most sensitive to, is the mixed layer depth, $H_{ML}$ , which is generally a 1D, local, physical property, determined by surface forcing. The sensitivity to mixed layer depth is followed by sensitivity to boundary layer depth, $H_{BL}$ , the buoyancy gradient and vorticity. Note that $\overline{w^{\prime}b^{\prime}}^{z}$ does not appear to be sensitive to local-changes in surface heat flux, surface wind stress, or Coriolis, which is likely due to these fields being smoother in the LLC4320 at the scales relevant for the Jacobian. Despite strain being the most important feature in the previous section, it is only in the $1^{o}$ resolution experiment that the Jacobian exhibits sensitivity of $\overline{w^{\prime}b^{\prime}}^{z}$ to vorticity, divergence, and strain, indicating that these fields are significantly non-local unless the largest filter is applied.

To further understand the relevance of locality, in figure 10, we examine the Jacobian of the output center point with respect to the full domain of each input feature. Figure 10a shows an example for the buoyancy gradient input feature, where the shaded area illustrates the CNN’s receptive field needed to predict a single output point. Averaging over that halo, we examine the fraction of Jacobian over the number of grid points, which can be thought of as the percentage of sensitivity for each input feature that is being captured by the CNN. Figure 10b is an example for the $1/4^{o}$ resolution case, where on average, we find that 7 grid points away from the center is sufficient for capturing 90% of the Jacobian fraction, e.g. the 90% of sensitivity between the output and input features. This relatively local receptive field is also found to be consistent across resolutions (figure 10c), despite the varying importance of input features found previously.

5 Discussion and Conclusions

The parameterization for submesoscale vertical buoyancy fluxes plays a key role in setting stratification in the ocean mixed layer, and as such contributes to the exchange between the ocean and atmosphere systems. In this work, we introduce an improved parameterization based on a data-driven approach, where a CNN is trained to learn submesoscale vertical buoyancy fluxes, $\overline{w^{\prime}b^{\prime}}^{z}$ , given by large scale variables that help set it. The subgrid flux, $\overline{w^{\prime}b^{\prime}}^{z}$ , is inferred by the CNN as a function of 9 large-scale input features with known relevance to submesoscale vertical buoyancy flux: $\overline{\nabla b}^{z}$ , $f$ , $H_{ML}$ , $N^{2}$ , $Q^{*}$ , $\tau$ , $H_{BL}$ , $\overline{\sigma}^{z}$ , $\overline{\delta}^{z}$ , $\overline{\zeta}^{z}$ (see table 1). The data used for training is given from 12 regions sampled from the global high-resolution LLC4320 simulation output. The CNN is trained over a random selection of $80\%$ of all data, while the remaining $20\%$ is unseen by the CNN and is used for testing. We perform five resolution experiments of $1/12^{o},1/8^{o},1/4^{o},1/2^{o},1^{o}$ and compare with the \citeAbodner2023modifying formulation of the MLE parameterization. We consistently find that the CNN predictions improves on the MLE parameterization, with higher $R^{2}$ values across all regions, seasons, and resolutions tested in this study. We additionally perform several sensitivity experiments, where we test the CNN’s ability to generalize on regional or temporal data held out during training. It is found that the CNN, in particular in the low resolution experiments, is able to make skillful predictions on unseen data as long as it is trained on seasons where submesoscales are most active, which generally corresponds to winter and spring months. We have also shown that the CNN is able to generalize on most regions that were held out during training, with dependency on the dominant scales of variability in $\overline{w^{\prime}b^{\prime}}^{z}$ .

The significant improvement on the MLE parameterization indicates that the CNN has learned meaningful relationships between the input features and $\overline{w^{\prime}b^{\prime}}^{z}$ that are likely not captured by the MLE parameterization, and are able to generalize over widely different dynamical regimes. Thus, we applied two complimentary explainability methods which enable a closer look and the relationships between the CNN output, $\overline{w^{\prime}b^{\prime}}^{z}$ , and the CNN input features. The point-wise dependency is computed by the Jacobian of the output with respect to inputs along the CNN weights. We find that a strong dependency between $\overline{w^{\prime}b^{\prime}}^{z}$ and the mixed layer depth emerges. In contrast, from a set of ablation experiments, we find that including strain as an input feature significantly improves skill in terms of $R^{2}$ values. To summarize, we find that the CNN exhibits strong dependency on the local relationship between $\overline{w^{\prime}b^{\prime}}^{z}$ and the mixed layer depth, a 1D property driven by surface forcing, and non-locally on the large scale strain field, a variable currently missing from the MLE paremeterization in GCMs. An interesting application here would be to test whether this result holds in other submesoscale permitting simulations, such as those compared in \citeAuchida2022cloud. Other parameterizations such as \citeAzhang2023parameterizing have suggested a theoretical formulation that includes a relationship between submesoscale vertical buoyancy flux with the strain field. An equation discovery approach [<]e.g.,¿[]zanna2020data, may enable a closer comparison with \citeAzhang2023parameterizing, and whether a strong relationship between strain and submesoscale fluxes as well as a new shape function emerge in a similar fashion.

We have shown that the CNN improves on the MLE parameterization in an offline setting. A next natural step is to explore the implications of better captured $\overline{w^{\prime}b^{\prime}}^{z}$ in a GCM and compare with the MLE parameterization online. We have designed our method to correspond with the existing implementation of the MLE parameterization in GCMs, where $\overline{w^{\prime}b^{\prime}}^{z}$ in (S3) can simply be replaced with the CNN. A relatively small receptive field of 7 grid points is found to be sufficient at capturing relationships between the input features and $\overline{w^{\prime}b^{\prime}}^{z}$ , which suggests that a smaller network may aid future implementation efforts in GCMs [C. Zhang \BOthers. (\APACyear2023)]. A decomposition may be preferred to distinguish the bulk restratification effect with the intermittent negative fluxes, and will allow a more natural relationship with vertical buoyancy fluxes already estimated in boundary layer turbulence parameterizations [Large \BOthers. (\APACyear1994), Reichl \BBA Hallberg (\APACyear2018)]. The exact formulation, implementation, and evaluation of impact on climate variables is left for future work.

Acknowledgements

AB was supported by a grant from the Simons Foundation: award number 855143. AB, DB and LZ received M²LInES research funding through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures program. We thank members of the M²LInES project for support and constructive feedback during the formulation of ideas, in particular, Pavel Perezhogin, Chris Pedersen, Ryan Abernathey, Carlos Fernandez-Granda, and Fabrizio Falasca. The authors would also like to thank the Pangeo project for providing open-source code which enabled timely analysis for working with the LLC4320 data.

References

Ajayi \BOthers. (\APACyear2021) \APACinsertmetastarajayi2021diagnosing{APACrefauthors}Ajayi, A., Le Sommer, J., Chassignet, E\BPBIP., Molines, J\BHBIM., Xu, X., Albert, A.\BCBL \BBA Dewar, W. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleDiagnosing cross-scale kinetic energy exchanges from two submesoscale permitting ocean models Diagnosing cross-scale kinetic energy exchanges from two submesoscale permitting ocean models.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems136e2019MS001923. \PrintBackRefs\CurrentBib
Bachman \BOthers. (\APACyear2017) \APACinsertmetastarbachman2017parameterization{APACrefauthors}Bachman, S\BPBID., Fox-Kemper, B., Taylor, J\BPBIR.\BCBL \BBA Thomas, L\BPBIN. \APACrefYearMonthDay2017. \BBOQ\APACrefatitleParameterization of frontal symmetric instabilities. I: Theory for resolved fronts Parameterization of frontal symmetric instabilities. i: Theory for resolved fronts.\BBCQ \APACjournalVolNumPagesOcean Modelling10972–95. \PrintBackRefs\CurrentBib
Bachman \BBA Klocker (\APACyear2020) \APACinsertmetastarbachman2020interaction{APACrefauthors}Bachman, S\BPBID.\BCBT \BBA Klocker, A. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleInteraction of jets and submesoscale dynamics leads to rapid ocean ventilation Interaction of jets and submesoscale dynamics leads to rapid ocean ventilation.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography50102873–2883. \PrintBackRefs\CurrentBib
Balwada \BOthers. (\APACyear2021) \APACinsertmetastarbalwada2021vertical{APACrefauthors}Balwada, D., Xiao, Q., Smith, S., Abernathey, R.\BCBL \BBA Gray, A\BPBIR. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleVertical fluxes conditioned on vorticity and strain reveal submesoscale ventilation Vertical fluxes conditioned on vorticity and strain reveal submesoscale ventilation.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography5192883–2901. \PrintBackRefs\CurrentBib
Barkan \BOthers. (\APACyear2019) \APACinsertmetastarbarkan2019role{APACrefauthors}Barkan, R., Molemaker, M\BPBIJ., Srinivasan, K., McWilliams, J\BPBIC.\BCBL \BBA D’Asaro, E\BPBIA. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleThe role of horizontal divergence in submesoscale frontogenesis The role of horizontal divergence in submesoscale frontogenesis.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography4961593–1618. \PrintBackRefs\CurrentBib
Boccaletti \BOthers. (\APACyear2007) \APACinsertmetastarboccaletti2007mixed{APACrefauthors}Boccaletti, G., Ferrari, R.\BCBL \BBA Fox-Kemper, B. \APACrefYearMonthDay2007. \BBOQ\APACrefatitleMixed layer instabilities and restratification Mixed layer instabilities and restratification.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography3792228–2250. \PrintBackRefs\CurrentBib
Bodner \BOthers. (\APACyear2023) \APACinsertmetastarbodner2023modifying{APACrefauthors}Bodner, A\BPBIS., Fox-Kemper, B., Johnson, L., Van Roekel, L\BPBIP., McWilliams, J\BPBIC., Sullivan, P\BPBIP.\BDBLDong, J. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleModifying the Mixed Layer Eddy Parameterization to Include Frontogenesis Arrest by Boundary Layer Turbulence Modifying the mixed layer eddy parameterization to include frontogenesis arrest by boundary layer turbulence.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography531323–339. \PrintBackRefs\CurrentBib
Bodner \BOthers. (\APACyear2020) \APACinsertmetastarbodner2020perturbation{APACrefauthors}Bodner, A\BPBIS., Fox-Kemper, B., Van Roekel, L\BPBIP., McWilliams, J\BPBIC.\BCBL \BBA Sullivan, P\BPBIP. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleA perturbation approach to understanding the effects of turbulence on frontogenesis A perturbation approach to understanding the effects of turbulence on frontogenesis.\BBCQ \APACjournalVolNumPagesJournal of Fluid Mechanics883. \PrintBackRefs\CurrentBib
Bolton \BBA Zanna (\APACyear2019) \APACinsertmetastarbolton2019applications{APACrefauthors}Bolton, T.\BCBT \BBA Zanna, L. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleApplications of deep learning to ocean data inference and subgrid parameterization Applications of deep learning to ocean data inference and subgrid parameterization.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems111376–399. \PrintBackRefs\CurrentBib
Bopp \BOthers. (\APACyear2015) \APACinsertmetastarbopp2015pathways{APACrefauthors}Bopp, L., Lévy, M., Resplandy, L.\BCBL \BBA Sallée, J\BHBIB. \APACrefYearMonthDay2015. \BBOQ\APACrefatitlePathways of anthropogenic carbon subduction in the global ocean Pathways of anthropogenic carbon subduction in the global ocean.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters42156416–6423. \PrintBackRefs\CurrentBib
Callies \BBA Ferrari (\APACyear2018) \APACinsertmetastarcallies2018baroclinic{APACrefauthors}Callies, J.\BCBT \BBA Ferrari, R. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleBaroclinic instability in the presence of convection Baroclinic instability in the presence of convection.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography48145–60. \PrintBackRefs\CurrentBib
Callies \BOthers. (\APACyear2015) \APACinsertmetastarcallies2015seasonality{APACrefauthors}Callies, J., Ferrari, R., Klymak, J\BPBIM.\BCBL \BBA Gula, J. \APACrefYearMonthDay2015. \BBOQ\APACrefatitleSeasonality in submesoscale turbulence Seasonality in submesoscale turbulence.\BBCQ \APACjournalVolNumPagesNature communications616862. \PrintBackRefs\CurrentBib
Calvert \BOthers. (\APACyear2020) \APACinsertmetastarcalvert2020impact{APACrefauthors}Calvert, D., Nurser, G., Bell, M\BPBIJ.\BCBL \BBA Fox-Kemper, B. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleThe impact of a parameterisation of submesoscale mixed layer eddies on mixed layer depths in the NEMO ocean model The impact of a parameterisation of submesoscale mixed layer eddies on mixed layer depths in the nemo ocean model.\BBCQ \APACjournalVolNumPagesOcean Modelling154101678. \PrintBackRefs\CurrentBib
Capet \BOthers. (\APACyear2008) \APACinsertmetastarcapet2008mesoscale{APACrefauthors}Capet, X., McWilliams, J\BPBIC., Molemaker, M\BPBIJ.\BCBL \BBA Shchepetkin, A. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleMesoscale to submesoscale transition in the California Current System. Part II: Frontal processes Mesoscale to submesoscale transition in the california current system. part ii: Frontal processes.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography38144–64. \PrintBackRefs\CurrentBib
Chattopadhyay \BOthers. (\APACyear2020) \APACinsertmetastarchattopadhyay2020predicting{APACrefauthors}Chattopadhyay, A., Hassanzadeh, P.\BCBL \BBA Pasha, S. \APACrefYearMonthDay2020. \BBOQ\APACrefatitlePredicting clustered weather patterns: A test case for applications of convolutional neural networks to spatio-temporal climate data Predicting clustered weather patterns: A test case for applications of convolutional neural networks to spatio-temporal climate data.\BBCQ \APACjournalVolNumPagesScientific reports1011317. \PrintBackRefs\CurrentBib
Dagon \BOthers. (\APACyear2022) \APACinsertmetastardagon2022machine{APACrefauthors}Dagon, K., Truesdale, J., Biard, J\BPBIC., Kunkel, K\BPBIE., Meehl, G\BPBIA.\BCBL \BBA Molina, M\BPBIJ. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleMachine Learning-Based Detection of Weather Fronts and Associated Extreme Precipitation in Historical and Future Climates Machine learning-based detection of weather fronts and associated extreme precipitation in historical and future climates.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Atmospheres12721e2022JD037038. \PrintBackRefs\CurrentBib
de Boyer Montégut \BOthers. (\APACyear2004) \APACinsertmetastarde2004mixed{APACrefauthors}de Boyer Montégut, C., Madec, G., Fischer, A\BPBIS., Lazar, A.\BCBL \BBA Iudicone, D. \APACrefYearMonthDay2004. \BBOQ\APACrefatitleMixed layer depth over the global ocean: An examination of profile data and a profile-based climatology Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans109C12. \PrintBackRefs\CurrentBib
Dong \BOthers. (\APACyear2020) \APACinsertmetastardong2020seasonality{APACrefauthors}Dong, J., Fox-Kemper, B., Zhang, H.\BCBL \BBA Dong, C. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleThe seasonality of submesoscale energy production, content, and cascade The seasonality of submesoscale energy production, content, and cascade.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters476e2020GL087388. \PrintBackRefs\CurrentBib
Ferrari \BOthers. (\APACyear2008) \APACinsertmetastarferrari2008parameterization{APACrefauthors}Ferrari, R., McWilliams, J\BPBIC., Canuto, V\BPBIM.\BCBL \BBA Dubovikov, M. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleParameterization of eddy fluxes near oceanic boundaries Parameterization of eddy fluxes near oceanic boundaries.\BBCQ \APACjournalVolNumPagesJournal of Climate21122770–2789. \PrintBackRefs\CurrentBib
Forget \BOthers. (\APACyear2015) \APACinsertmetastarforget2015ecco{APACrefauthors}Forget, G., Campin, J\BHBIM., Heimbach, P., Hill, C., Ponte, R.\BCBL \BBA Wunsch, C. \APACrefYearMonthDay2015. \BBOQ\APACrefatitleECCO version 4: An integrated framework for non-linear inverse modeling and global ocean state estimation Ecco version 4: An integrated framework for non-linear inverse modeling and global ocean state estimation.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development8103071–3104. \PrintBackRefs\CurrentBib
Fox-Kemper \BOthers. (\APACyear2011) \APACinsertmetastarfox2011parameterization{APACrefauthors}Fox-Kemper, B., Danabasoglu, G., Ferrari, R., Griffies, S., Hallberg, R., Holland, M.\BDBLSamuels, B. \APACrefYearMonthDay2011. \BBOQ\APACrefatitleParameterization of mixed layer eddies. III: Implementation and impact in global ocean climate simulations Parameterization of mixed layer eddies. iii: Implementation and impact in global ocean climate simulations.\BBCQ \APACjournalVolNumPagesOcean Modelling391-261–78. \PrintBackRefs\CurrentBib
Fox-Kemper \BOthers. (\APACyear2008) \APACinsertmetastarfox2008parameterization{APACrefauthors}Fox-Kemper, B., Ferrari, R.\BCBL \BBA Hallberg, R. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleParameterization of mixed layer eddies. Part I: Theory and diagnosis Parameterization of mixed layer eddies. part i: Theory and diagnosis.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography3861145–1165. \PrintBackRefs\CurrentBib
Fox-Kemper \BOthers. (\APACyear2022) \APACinsertmetastarfox2022ocean{APACrefauthors}Fox-Kemper, B., Johnson, L.\BCBL \BBA Qiao, F. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleOcean near-surface layers Ocean near-surface layers.\BBCQ \BIn \APACrefbtitleOcean Mixing Ocean mixing (\BPGS 65–94). \APACaddressPublisherElsevier. \PrintBackRefs\CurrentBib
Frankignoul \BBA Hasselmann (\APACyear1977) \APACinsertmetastarfrankignoul1977stochastic{APACrefauthors}Frankignoul, C.\BCBT \BBA Hasselmann, K. \APACrefYearMonthDay1977. \BBOQ\APACrefatitleStochastic climate models, Part II Application to sea-surface temperature anomalies and thermocline variability Stochastic climate models, part ii application to sea-surface temperature anomalies and thermocline variability.\BBCQ \APACjournalVolNumPagesTellus294289–305. \PrintBackRefs\CurrentBib
Gallmeier \BOthers. (\APACyear2023) \APACinsertmetastargallmeier2023evaluation{APACrefauthors}Gallmeier, K\BPBIM., Prochaska, J\BPBIX., Cornillon, P., Menemenlis, D.\BCBL \BBA Kelm, M. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleAn evaluation of the LLC4320 global ocean simulation based on the submesoscale structure of modeled sea surface temperature fields An evaluation of the llc4320 global ocean simulation based on the submesoscale structure of modeled sea surface temperature fields.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development Discussions1–42. \PrintBackRefs\CurrentBib
Gent \BBA Mcwilliams (\APACyear1990) \APACinsertmetastargent1990isopycnal{APACrefauthors}Gent, P\BPBIR.\BCBT \BBA Mcwilliams, J\BPBIC. \APACrefYearMonthDay1990. \BBOQ\APACrefatitleIsopycnal mixing in ocean circulation models Isopycnal mixing in ocean circulation models.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography201150–155. \PrintBackRefs\CurrentBib
Griffies (\APACyear1998) \APACinsertmetastarGriffies98{APACrefauthors}Griffies, S\BPBIM. \APACrefYearMonthDay1998. \BBOQ\APACrefatitleThe Gent–McWilliams skew flux The gent–mcwilliams skew flux.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography285831–841. \PrintBackRefs\CurrentBib
Guillaumin \BBA Zanna (\APACyear2021) \APACinsertmetastarguillaumin2021stochastic{APACrefauthors}Guillaumin, A\BPBIP.\BCBT \BBA Zanna, L. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleStochastic-deep learning parameterization of ocean momentum forcing Stochastic-deep learning parameterization of ocean momentum forcing.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems139e2021MS002534. \PrintBackRefs\CurrentBib
Gula \BOthers. (\APACyear2022) \APACinsertmetastargula2022submesoscale{APACrefauthors}Gula, J., Taylor, J., Shcherbina, A.\BCBL \BBA Mahadevan, A. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleSubmesoscale processes and mixing Submesoscale processes and mixing.\BBCQ \BIn \APACrefbtitleOcean Mixing Ocean mixing (\BPGS 181–214). \APACaddressPublisherElsevier. \PrintBackRefs\CurrentBib
Hoskins \BBA Bretherton (\APACyear1972) \APACinsertmetastarhoskins1972atmospheric{APACrefauthors}Hoskins, B\BPBIJ.\BCBT \BBA Bretherton, F\BPBIP. \APACrefYearMonthDay1972. \BBOQ\APACrefatitleAtmospheric frontogenesis models: Mathematical formulation and solution Atmospheric frontogenesis models: Mathematical formulation and solution.\BBCQ \APACjournalVolNumPagesJournal of the atmospheric sciences29111–37. \PrintBackRefs\CurrentBib
IPCC (\APACyear2019) \APACinsertmetastarSROCC{APACrefauthors}IPCC. \APACrefYear2019. \APACrefbtitleSpecial Report on the Ocean and Cryosphere in a Changing Climate Special report on the ocean and cryosphere in a changing climate [Book]. \APACaddressPublisherCambridge, United Kingdom and New York, NY, USACambridge University Press. {APACrefDOI} 10.1017/9781009157964 \PrintBackRefs\CurrentBib
IPCC (\APACyear2021) \APACinsertmetastarRN1{APACrefauthors}IPCC. \APACrefYear2021. \APACrefbtitleClimate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change Climate change 2021: The physical science basis. contribution of working group i to the sixth assessment report of the intergovernmental panel on climate change (\BVOL In Press) [Book]. \APACaddressPublisherCambridge, United Kingdom and New York, NY, USACambridge University Press. {APACrefDOI} 10.1017/9781009157896 \PrintBackRefs\CurrentBib
Johnson \BOthers. (\APACyear2016) \APACinsertmetastarjohnson2016global{APACrefauthors}Johnson, L., Lee, C\BPBIM.\BCBL \BBA D’Asaro, E\BPBIA. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleGlobal estimates of lateral springtime restratification Global estimates of lateral springtime restratification.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography4651555–1573. \PrintBackRefs\CurrentBib
Jones \BOthers. (\APACyear2023) \APACinsertmetastarjones2023using{APACrefauthors}Jones, C\BPBIS., Xiao, Q., Abernathey, R\BPBIP.\BCBL \BBA Smith, K\BPBIS. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleUsing Lagrangian filtering to remove waves from the ocean surface velocity field Using lagrangian filtering to remove waves from the ocean surface velocity field.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems154e2022MS003220. \PrintBackRefs\CurrentBib
Lapeyre \BOthers. (\APACyear2006) \APACinsertmetastarlapeyre2006oceanic{APACrefauthors}Lapeyre, G., Klein, P.\BCBL \BBA Hua, B\BPBIL. \APACrefYearMonthDay2006. \BBOQ\APACrefatitleOceanic restratification forced by surface frontogenesis Oceanic restratification forced by surface frontogenesis.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography3681577–1590. \PrintBackRefs\CurrentBib
Large \BOthers. (\APACyear1994) \APACinsertmetastarlarge1994oceanic{APACrefauthors}Large, W\BPBIG., McWilliams, J\BPBIC.\BCBL \BBA Doney, S\BPBIC. \APACrefYearMonthDay1994. \BBOQ\APACrefatitleOceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization.\BBCQ \APACjournalVolNumPagesReviews of geophysics324363–403. \PrintBackRefs\CurrentBib
Loose \BOthers. (\APACyear2022) \APACinsertmetastarloose2022gcm{APACrefauthors}Loose, N., Abernathey, R., Grooms, I., Busecke, J., Guillaumin, A., Yankovsky, E.\BDBLothers \APACrefYearMonthDay2022. \BBOQ\APACrefatitleGCM-filters: A Python package for diffusion-based spatial filtering of gridded data Gcm-filters: A python package for diffusion-based spatial filtering of gridded data.\BBCQ \APACjournalVolNumPagesJournal of Open Source Software770. \PrintBackRefs\CurrentBib
Mahadevan (\APACyear2016) \APACinsertmetastarmahadevan2016impact{APACrefauthors}Mahadevan, A. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleThe impact of submesoscale physics on primary productivity of plankton The impact of submesoscale physics on primary productivity of plankton.\BBCQ \APACjournalVolNumPagesAnnual review of marine science8161–184. \PrintBackRefs\CurrentBib
Mahadevan \BOthers. (\APACyear2010) \APACinsertmetastarmahadevan2010rapid{APACrefauthors}Mahadevan, A., Tandon, A.\BCBL \BBA Ferrari, R. \APACrefYearMonthDay2010. \BBOQ\APACrefatitleRapid changes in mixed layer stratification driven by submesoscale instabilities and winds Rapid changes in mixed layer stratification driven by submesoscale instabilities and winds.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans115C3. \PrintBackRefs\CurrentBib
McWilliams (\APACyear2016) \APACinsertmetastarmcwilliams2016submesoscale{APACrefauthors}McWilliams, J\BPBIC. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleSubmesoscale currents in the ocean Submesoscale currents in the ocean.\BBCQ \APACjournalVolNumPagesProceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences472218920160117. \PrintBackRefs\CurrentBib
McWilliams (\APACyear2017) \APACinsertmetastarmcwilliams2017submesoscale{APACrefauthors}McWilliams, J\BPBIC. \APACrefYearMonthDay2017. \BBOQ\APACrefatitleSubmesoscale surface fronts and filaments: Secondary circulation, buoyancy flux, and frontogenesis Submesoscale surface fronts and filaments: Secondary circulation, buoyancy flux, and frontogenesis.\BBCQ \APACjournalVolNumPagesJournal of Fluid Mechanics823391–432. \PrintBackRefs\CurrentBib
Menemenlis \BOthers. (\APACyear2008) \APACinsertmetastarmenemenlis2008ecco2{APACrefauthors}Menemenlis, D., Campin, J\BHBIM., Heimbach, P., Hill, C., Lee, T., Nguyen, A.\BDBLZhang, H. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleECCO2: High resolution global ocean and sea ice data synthesis Ecco2: High resolution global ocean and sea ice data synthesis.\BBCQ \APACjournalVolNumPagesMercator Ocean Quarterly Newsletter31October13–21. \PrintBackRefs\CurrentBib
Menemenlis \BOthers. (\APACyear2021) \APACinsertmetastarmenemenlis2021pre{APACrefauthors}Menemenlis, D., Hill, C., Henze, C., Wang, J.\BCBL \BBA Fenty, I. \APACrefYearMonthDay2021. \APACrefbtitlePre-SWOT Level-4 Hourly MITgcm LLC4320 Native 2km Grid Oceanographic Version 1.0. Pre-swot level-4 hourly mitgcm llc4320 native 2km grid oceanographic version 1.0. \APACaddressPublisherVer. \PrintBackRefs\CurrentBib
Perezhogin \BOthers. (\APACyear2023) \APACinsertmetastarperezhogin2023generative{APACrefauthors}Perezhogin, P., Zanna, L.\BCBL \BBA Fernandez-Granda, C. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleGenerative data-driven approaches for stochastic subgrid parameterizations in an idealized ocean model Generative data-driven approaches for stochastic subgrid parameterizations in an idealized ocean model.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2302.07984. \PrintBackRefs\CurrentBib
Reichl \BBA Hallberg (\APACyear2018) \APACinsertmetastarreichl2018simplified{APACrefauthors}Reichl, B\BPBIG.\BCBT \BBA Hallberg, R. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleA simplified energetics based planetary boundary layer (ePBL) approach for ocean climate simulations. A simplified energetics based planetary boundary layer (epbl) approach for ocean climate simulations.\BBCQ \APACjournalVolNumPagesOcean Modelling132112–129. \PrintBackRefs\CurrentBib
Rocha, Chereskin\BCBL \BOthers. (\APACyear2016) \APACinsertmetastarrocha2016mesoscale{APACrefauthors}Rocha, C\BPBIB., Chereskin, T\BPBIK., Gille, S\BPBIT.\BCBL \BBA Menemenlis, D. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleMesoscale to submesoscale wavenumber spectra in Drake Passage Mesoscale to submesoscale wavenumber spectra in drake passage.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography462601–620. \PrintBackRefs\CurrentBib
Rocha, Gille\BCBL \BOthers. (\APACyear2016) \APACinsertmetastarrocha2016seasonality{APACrefauthors}Rocha, C\BPBIB., Gille, S\BPBIT., Chereskin, T\BPBIK.\BCBL \BBA Menemenlis, D. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleSeasonality of submesoscale dynamics in the Kuroshio Extension Seasonality of submesoscale dynamics in the kuroshio extension.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters432111–304. \PrintBackRefs\CurrentBib
Ross \BOthers. (\APACyear2023) \APACinsertmetastarross2023benchmarking{APACrefauthors}Ross, A., Li, Z., Perezhogin, P., Fernandez-Granda, C.\BCBL \BBA Zanna, L. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleBenchmarking of machine learning ocean subgrid parameterizations in an idealized model Benchmarking of machine learning ocean subgrid parameterizations in an idealized model.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems151e2022MS003258. \PrintBackRefs\CurrentBib
Sane \BOthers. (\APACyear2023) \APACinsertmetastarsane2023parameterizing{APACrefauthors}Sane, A., Reichl, B\BPBIG., Adcroft, A.\BCBL \BBA Zanna, L. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleParameterizing Vertical Mixing Coefficients in the Ocean Surface Boundary Layer using Neural Networks Parameterizing vertical mixing coefficients in the ocean surface boundary layer using neural networks.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2306.09045. \PrintBackRefs\CurrentBib
Schubert \BOthers. (\APACyear2019) \APACinsertmetastarschubert2019submesoscale{APACrefauthors}Schubert, R., Schwarzkopf, F\BPBIU., Baschek, B.\BCBL \BBA Biastoch, A. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleSubmesoscale impacts on mesoscale Agulhas dynamics Submesoscale impacts on mesoscale agulhas dynamics.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems1182745–2767. \PrintBackRefs\CurrentBib
Shakespeare \BBA Taylor (\APACyear2013) \APACinsertmetastarshakespeare2013generalized{APACrefauthors}Shakespeare, C\BPBIJ.\BCBT \BBA Taylor, J\BPBIR. \APACrefYearMonthDay2013. \BBOQ\APACrefatitleA generalized mathematical model of geostrophic adjustment and frontogenesis: uniform potential vorticity A generalized mathematical model of geostrophic adjustment and frontogenesis: uniform potential vorticity.\BBCQ \APACjournalVolNumPagesJournal of fluid mechanics736366–413. \PrintBackRefs\CurrentBib
Shamekh \BOthers. (\APACyear2023) \APACinsertmetastarshamekh2023implicit{APACrefauthors}Shamekh, S., Lamb, K\BPBID., Huang, Y.\BCBL \BBA Gentine, P. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleImplicit learning of convective organization explains precipitation stochasticity Implicit learning of convective organization explains precipitation stochasticity.\BBCQ \APACjournalVolNumPagesProceedings of the National Academy of Sciences12020e2216158120. \PrintBackRefs\CurrentBib
Souza \BOthers. (\APACyear2020) \APACinsertmetastarsouza2020uncertainty{APACrefauthors}Souza, A\BPBIN., Wagner, G., Ramadhan, A., Allen, B., Churavy, V., Schloss, J.\BDBLothers \APACrefYearMonthDay2020. \BBOQ\APACrefatitleUncertainty quantification of ocean parameterizations: Application to the K-profile-parameterization for penetrative convection Uncertainty quantification of ocean parameterizations: Application to the k-profile-parameterization for penetrative convection.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems1212e2020MS002108. \PrintBackRefs\CurrentBib
Srinivasan \BOthers. (\APACyear2023) \APACinsertmetastarsrinivasan2023forward{APACrefauthors}Srinivasan, K., Barkan, R.\BCBL \BBA McWilliams, J\BPBIC. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleA forward energy flux at submesoscales driven by frontogenesis A forward energy flux at submesoscales driven by frontogenesis.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography531287–305. \PrintBackRefs\CurrentBib
Stanley \BOthers. (\APACyear2020) \APACinsertmetastarstanley2020parameterizing{APACrefauthors}Stanley, Z., Grooms, I., Kleiber, W., Bachman, S., Castruccio, F.\BCBL \BBA Adcroft, A. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleParameterizing the Impact of Unresolved Temperature Variability on the Large-Scale Density Field: Part 1. Theory. Parameterizing the impact of unresolved temperature variability on the large-scale density field: Part 1. theory.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems1212e2020MS002185. \PrintBackRefs\CurrentBib
Su \BOthers. (\APACyear2020) \APACinsertmetastarsu2020high{APACrefauthors}Su, Z., Torres, H., Klein, P., Thompson, A\BPBIF., Siegelman, L., Wang, J.\BDBLHill, C. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleHigh-frequency submesoscale motions enhance the upward vertical heat transport in the global ocean High-frequency submesoscale motions enhance the upward vertical heat transport in the global ocean.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans1259e2020JC016544. \PrintBackRefs\CurrentBib
Su \BOthers. (\APACyear2018) \APACinsertmetastarsu2018ocean{APACrefauthors}Su, Z., Wang, J., Klein, P., Thompson, A\BPBIF.\BCBL \BBA Menemenlis, D. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleOcean submesoscales as a key component of the global heat budget Ocean submesoscales as a key component of the global heat budget.\BBCQ \APACjournalVolNumPagesNature communications91775. \PrintBackRefs\CurrentBib
Taylor \BBA Thompson (\APACyear2023) \APACinsertmetastartaylor2023submesoscale{APACrefauthors}Taylor, J\BPBIR.\BCBT \BBA Thompson, A\BPBIF. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleSubmesoscale dynamics in the upper ocean Submesoscale dynamics in the upper ocean.\BBCQ \APACjournalVolNumPagesAnnual Review of Fluid Mechanics55103–127. \PrintBackRefs\CurrentBib
Thomas (\APACyear2005) \APACinsertmetastarthomas2005destruction{APACrefauthors}Thomas, L\BPBIN. \APACrefYearMonthDay2005. \BBOQ\APACrefatitleDestruction of potential vorticity by winds Destruction of potential vorticity by winds.\BBCQ \APACjournalVolNumPagesJournal of physical oceanography35122457–2466. \PrintBackRefs\CurrentBib
Thomas \BOthers. (\APACyear2008) \APACinsertmetastarthomas2008submesoscale{APACrefauthors}Thomas, L\BPBIN., Tandon, A.\BCBL \BBA Mahadevan, A. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleSubmesoscale processes and dynamics Submesoscale processes and dynamics.\BBCQ \APACjournalVolNumPagesOcean modeling in an Eddying Regime17717–38. \PrintBackRefs\CurrentBib
Thomas \BOthers. (\APACyear2013) \APACinsertmetastarthomas2013symmetric{APACrefauthors}Thomas, L\BPBIN., Taylor, J\BPBIR., Ferrari, R.\BCBL \BBA Joyce, T\BPBIM. \APACrefYearMonthDay2013. \BBOQ\APACrefatitleSymmetric instability in the Gulf Stream Symmetric instability in the gulf stream.\BBCQ \APACjournalVolNumPagesDeep Sea Research Part II: Topical Studies in Oceanography9196–110. \PrintBackRefs\CurrentBib
Torres \BOthers. (\APACyear2018) \APACinsertmetastartorres2018partitioning{APACrefauthors}Torres, H\BPBIS., Klein, P., Menemenlis, D., Qiu, B., Su, Z., Wang, J.\BDBLFu, L\BHBIL. \APACrefYearMonthDay2018. \BBOQ\APACrefatitlePartitioning ocean motions into balanced motions and internal gravity waves: A modeling study in anticipation of future space missions Partitioning ocean motions into balanced motions and internal gravity waves: A modeling study in anticipation of future space missions.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans123118084–8105. \PrintBackRefs\CurrentBib
Uchida \BOthers. (\APACyear2022) \APACinsertmetastaruchida2022cloud{APACrefauthors}Uchida, T., Le Sommer, J., Stern, C., Abernathey, R., Holdgraf, C., Albert, A.\BDBLothers \APACrefYearMonthDay2022. \BBOQ\APACrefatitleCloud-based framework for inter-comparing submesoscale permitting realistic ocean models Cloud-based framework for inter-comparing submesoscale permitting realistic ocean models.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development Discussions1–32. \PrintBackRefs\CurrentBib
Wang \BOthers. (\APACyear2022) \APACinsertmetastarwang2022non{APACrefauthors}Wang, P., Yuval, J.\BCBL \BBA O’Gorman, P\BPBIA. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleNon-local parameterization of atmospheric subgrid processes with neural networks Non-local parameterization of atmospheric subgrid processes with neural networks.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems1410e2022MS002984. \PrintBackRefs\CurrentBib
Wenegrat \BBA Thomas (\APACyear2020) \APACinsertmetastarwenegrat2020centrifugal{APACrefauthors}Wenegrat, J\BPBIO.\BCBT \BBA Thomas, L\BPBIN. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleCentrifugal and symmetric instability during Ekman adjustment of the bottom boundary layer Centrifugal and symmetric instability during ekman adjustment of the bottom boundary layer.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography5061793–1812. \PrintBackRefs\CurrentBib
Xiao \BOthers. (\APACyear2023) \APACinsertmetastarxiao2023reconstruction{APACrefauthors}Xiao, Q., Balwada, D., Jones, C\BPBIS., Herrero-González, M., Smith, K\BPBIS.\BCBL \BBA Abernathey, R. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleReconstruction of surface kinematics from sea surface height using neural networks Reconstruction of surface kinematics from sea surface height using neural networks.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems1510e2023MS003709. \PrintBackRefs\CurrentBib
Yuval \BBA O’Gorman (\APACyear2020) \APACinsertmetastaryuval2020stable{APACrefauthors}Yuval, J.\BCBT \BBA O’Gorman, P\BPBIA. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleStable machine-learning parameterization of subgrid processes for climate modeling at a range of resolutions Stable machine-learning parameterization of subgrid processes for climate modeling at a range of resolutions.\BBCQ \APACjournalVolNumPagesNature communications1113295. \PrintBackRefs\CurrentBib
Zanna \BBA Bolton (\APACyear2020) \APACinsertmetastarzanna2020data{APACrefauthors}Zanna, L.\BCBT \BBA Bolton, T. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleData-driven equation discovery of ocean mesoscale closures Data-driven equation discovery of ocean mesoscale closures.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters4717e2020GL088376. \PrintBackRefs\CurrentBib
C. Zhang \BOthers. (\APACyear2023) \APACinsertmetastarzhang2023implementation{APACrefauthors}Zhang, C., Perezhogin, P., Gultekin, C., Adcroft, A., Fernandez-Granda, C.\BCBL \BBA Zanna, L. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleImplementation and Evaluation of a Machine Learned Mesoscale Eddy Parameterization into a Numerical Ocean Circulation Model Implementation and evaluation of a machine learned mesoscale eddy parameterization into a numerical ocean circulation model.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2303.00962. \PrintBackRefs\CurrentBib
J. Zhang \BOthers. (\APACyear2023) \APACinsertmetastarzhang2023parameterizing{APACrefauthors}Zhang, J., Zhang, Z.\BCBL \BBA Qiu, B. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleParameterizing Submesoscale Vertical Buoyancy Flux by Simultaneously Considering Baroclinic Instability and Strain-Induced Frontogenesis Parameterizing submesoscale vertical buoyancy flux by simultaneously considering baroclinic instability and strain-induced frontogenesis.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters508e2022GL102292. \PrintBackRefs\CurrentBib
Zhu \BOthers. (\APACyear2023) \APACinsertmetastarzhu2023deep{APACrefauthors}Zhu, R., Li, Y., Chen, Z., Du, T., Zhang, Y., Li, Z.\BDBLWu, L. \APACrefYearMonthDay2023. \BBOQ\APACrefatitleDeep learning improves reconstruction of ocean vertical velocity Deep learning improves reconstruction of ocean vertical velocity.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters5019e2023GL104889. \PrintBackRefs\CurrentBib

Supplementary Material

The formula for $\boldsymbol{\Psi}_{MLE}$ in \citeAfox2008parameterization is provided by a scaling for,

\overline{w^{\prime}b^{\prime}}^{z}\propto\frac{H_{ML}|\overline{\nabla_{H}b|}^{z}}{|f|},

(S3)

where $H_{ML}$ is the mixed layer depth, $f$ is the Coriolis parameter, $w$ is vertical velocity, $b$ is buoyancy, and $\nabla_{H}b$ is the horizontal buoyancy gradient. We follow the notation in \citeAfox2008parameterization, where the horizontal spatial resolution of the GCM is denoted $\overline{(\cdot)}$ and $(\cdot)^{\prime}$ is the unresolved subgrid variable. Superscript $z$ represents a vertical average over the mixed layer depth. . The scaling for submesoscale vertical buoyancy flux represents the bulk extraction of potential energy by MLEs within the mixed layer. A shape function $\mu(z)$ estimates the depth, $z$ , at which the MLE fluxes are activated,

\mu(z)=\max\left(0,\left[1-\left(\frac{2z}{H_{ML}}+1\right)^{2}\right]\left[1+\frac{5}{21}\left(\frac{2z}{H_{ML}}+1\right)^{2}\right]\right)

(S4)

where $\mu(z)$ is set to vanish at the surface and below the mixed layer $H_{ML}$ .