A Gradient Based Subgrid‐Scale Parameterization for Ocean Mesoscale Eddies

Mesoscale eddies play an important role in transport of heat and biogeochemical tracers in the global ocean circulation. Resolving these energetic eddies, however, is challenging in ocean general circulation models (OGCM) because it requires a horizontal grid spacing of ≲1/8° that is computationally expensive. As a result, we are required to parameterize mesoscale eddy effects on large‐scale ocean flows. In this work, we introduce a new subgrid‐scale (SGS) model that is developed based on a Taylor series expansion of resolved variables to parameterize subgrid mesoscale eddy transports and momentum fluxes in OGCM. We have performed an a priori study to evaluate the performance of our new gradient model using high‐resolution ocean simulations. Our results show that the gradient model well predicts the actual SGS thickness fluxes in the zonal and meridional directions in coarse‐resolution simulations with the grid spacing ≳1/4°. The unresolved kinetic energy at the ocean surface is also skillfully estimated. More importantly, unlike current mesoscale eddy parameterizations, which are mainly developed based on an assumption of flat bottom ocean, our new SGS model can capture the structure of unresolved standing meanders at the ocean surface. We have also developed a dynamic procedure for setting in non‐dimensional parameters in our new parameterization through a non‐ad hoc and tuning‐free method. Overall, this work suggests that implementing the gradient model in OGCM can improve the model accuracy with an affordable computational cost in eddy‐permitting and non‐eddying simulations.

2 of 26 studies have suggested that an inverse energy cascade occurs in the ocean on scales larger than the Rossby deformation radius (Bachman et al., 2017;Kjellsson & Zanna, 2017;Rhines, 1975;Scott & Wang, 2005;Tulloch et al., 2011), and both numerical and observational work have attempted to infer the inter-scale energy interactions from quasi geostrophic (QG) turbulence (Held et al., 1995;Xu & Fu, 2011), submesoscale instability (Molemaker et al., 2010), or internal wave generation, interactions, and breaking (Khani & Waite, 2020;Kunze, 2019;Nagai et al., 2015), but a complete understanding is an active research area. As with other anisotropic turbulent flows in the atmosphere and ocean, such as atmospheric boundary layer turbulence (e.g., Porté-Agel et al., 2000), stratified turbulence (Khani, 2018;Khani & Waite, 2015), and buoyancy-driven turbulence (e.g., Moeng & Sullivan, 1994), SGS parameterization of mesoscale eddies is clearly sensitive to model resolution (e.g., Hallberg, 2013;. The character of mesoscale eddies in ocean general circulation models (OGCM) considerably changes with latitude, stratification and ocean depth (Hallberg, 2013), but in general, we may categorize them into three model resolution ranges: (a) non-eddying (1° and coarser) in which the entire process of baroclinic instability must be parameterized, (b) eddy permitting (1/4°-1/2°), in which meanders and eddies are marginally resolved at the grid scale but these are actively damped by the model dissipation scheme; and (c) eddy-resolving (1/6° and finer, alternatively described as eddy rich given the acknowledgment that resolving all scales of importance is likely impossible), in which eddies are captured along with at least some of their interactions and secondary instabilities. In each range, the dynamics at and below the grid spacing have strikingly different behaviors and likely require a different parameterization approach. Simulations on basin or global scales using mesoscale eddy-resolving models are being run for longer and longer duration, but there are still many orders of magnitude remaining to the scales where turbulence and viscous dissipation occur (Fox-Kemper et al., 2019;Griffies et al., 2009). For this reason, "resolution-aware" SGS parameterizations are currently of great interest (see e.g., Bachman et al., 2017Bachman et al., , 2020Fox-Kemper & Menemenlis, 2008;Grooms et al., 2015;Jansen et al., 2019;Porta Mana & Zanna, 2014). The resolution-aware scheme in this context refers to the large-eddy simulation property that the SGS parameterization directly varies with the grid spacing, unlike the Reynolds-averaged scheme where unresolved fluxes are always averaged over the spatial domain. In this work, we introduce a gradient-based model, which is developed using a Taylor series expansion of resolved variables, to parameterize ocean mesoscale eddies. We do an a priori testing to study the performance of the gradient model using high-resolution ocean simulations with different bathymetry (in an a priori study, the performance of a SGS parameterization is evaluated using explicit filtering). The rest of this paper is in the following order: background on ocean mesoscale parameterizations, along with our new model development are introduced in Section 2. Section 3 includes methodology and simulations setup. Results and discussion are given in Sections 4-7, followed by concluding remarks and future research work in Section 8.

Background and Theoretical Development
The governing equations of motion under the Boussinesq approximation with the hydrostatic balance for the stacked shallow water model is written as follows 3 of 26 Here, u = (u, v), h, and C are the zonal and meridional velocities, thickness and passive tracer, respectively, and = ( ∕ ∕ ) is the lateral gradient operator; f and ′ = ( +1 − ) ∕ 0 are the Coriolis parameter and reduced gravity; ρ denotes density layers and ρ 0 is a reference density. The wind stress is shown by τ x ; i and j change from 1, 2 to show zonal and meridional directions, respectively; n denotes a density layer, η b is the bottom of the ocean and N is the total number of density layers (isopycnals). Unit vectors in the zonal and vertical directions are shown by ̂ and ̂ , respectively.

Current Eddy Parameterizations in Ocean Mesoscale
Possibly the most widely used SGS parameterization for non-eddying ocean models is the one proposed by Gent and McWilliams (1990, hereafter GM), which represents the influence of unresolved baroclinic eddies on resolved motions by dissipation of available potential energy (APE), and subsequently, lateral stirring of tracers is described by an isopycnal advection process.
If we employ a spatial filtering (or averaging) operator ⋅ to the mass continuity Equation 2, we have where ̃ ℎ =̃ ℎ −̃ h is not known and needs to be parameterized using resolved motions. In the GM model, the unresolved volume transport ̃ ℎ is parameterized as a lateral diffusion of isopycnal layer thickness. This allows the total SGS transport below (and including) a particular density layer n to be written as an overturning streamfunction ψ, and to be related to the isopycnal slope ̃ of the layer elevation: where the bar ⋅ sign denotes time average, and κ G is the GM eddy diffusivity coefficient. At the ocean surface, the GM model is constrained by the assumption of no barotropic SGS transport (Appendix A). The "SGS" label in this context means the unresolved fluxes prescribed by the parameterization, but in general SGS means all effects of unresolved motions that might appear in the numerical model's governing equations (e.g., the SGS thickness flux appears in the mass continuity equation and SGS/Reynolds stress in the momentum equation).
The GM overturning streamfunction implementation imposes the constraint that the SGS volume transport integrates to zero over the water column, consistent with the parameterization's motivation in the baroclinic instability process (indicated as the barotropic limit in Appendix A). As noted by Khani et al. (2019), this barotropic constraint at the ocean surface implies that the transport by unresolved topographic features/standing meanders in ocean simulations is not parameterized in the GM model. The GM parameterization has been modified and improved over the years, in particular through the desire to implement a variable and dynamically-motivated coefficient κ G .
For example, Visbeck et al. (1997) have proposed a GM coefficient itself based on isopycnal slope ( ∝  ̃ , where  is the local buoyancy frequency). Alternatively, the GM coefficient can be written as the product of a SGS velocity scale (derived from the eddy KE E sgs ) and a mixing length, that is, ∝ E 1∕2 sgs , where various forms of scaling for ℓ mix have been suggested, including a scale based on the Rossby deformation radius (the injection scale of baroclinic instability) or the Rhines scale (the large-scale end of the inverse cascade) (Canuto et al., 2018(Canuto et al., , 2019. However, these improved GM-based subgrid mesoscale eddy parameterizations still do not allow for the contributions of unresolved barotropic standing meanders at the ocean surface. It is also a near-universal characteristic of GM-based parameterizations that they act to purely dissipate potential energy (PE) as long as κ G > 0, in line with the underlying hypothesis of extraction of APE to KE through the baroclinic process.
Although APE dissipation by the mesoscale field is generally accepted, it has been pointed out that the lateral unresolved momentum fluxes (SGS stresses) may include a negative eddy viscosity term or other types of backscatter (energizing motions at the grid scale) approach (Bachman, 2019;. Attempts have been made to implement large-scale energizing terms dynamically (Berloff, 2016), stochastically (Grooms et al., 2015;Porta Mana & Zanna, 2014) or through principles of energy conservation (Bachman, 2019;Jansen et al., 2019). Progress in this direction has been encouraging, but in all cases 4 of 26 new complications have been added to the eddy parameterizing problem, including: (a) extra coefficients need to be tuned (e.g., Grooms et al., 2015;, (b) model parameters need to be determined using eddy-resolving simulations that add an extra computational cost (e.g., Berloff, 2016;Grooms et al., 2015;Porta Mana & Zanna, 2014), (c) stability issues in numerical schemes that leads to setting a very small Courant-Friedrichs-Lewy (CFL) condition, resulting in an inefficient time-stepping scheme (e.g., Jansen et al., 2019;. More importantly, the energizing terms in these models are practically employed to fix the highly dissipative nature of the SGS momentum models (e.g., regular or biharmonic Smagorinsky dissipation) in ocean simulations that excessively remove KE from scales around grid spacing. As a result, these energizing terms may not be designed to represent the inverse cascade mechanism for quasi-geostrophic flows, they rather alleviate intemperate dissipation around grid spacing (see e.g., Augier & Lindborg, 2013;Khani & Porté-Agel, 2022;Khani & Waite, 2020). Furthermore, a few nonlinear mesoscale eddy parameterizations, for which the mixing length hypothesis includes a non-Newtonian relationship between unresolved stress and rate-of-strain tensors, have recently been introduced to improve the state of mesoscale eddy parameterization in the context of ocean modeling (see e.g., Anstey & Zanna, 2017;Bachman et al., 2018).

Gradient Parameterization for Ocean Mesoscale
The gradient (or nonlinear) model is firstly introduced for parameterizing SGS momentum fluxes in large-eddy simulations (Clark et al., 1979;Leonard, 1975;Meneveau & Katz, 2000). The Gradient model is motivated by the fact that it can approximate actual subgrid momentum stresses by avoiding extra computational expenses of explicit filtering in the similarity model (as elaborated in Meneveau & Katz, 2000, employing the similarity model in addition to the purely dissipative Smagorinsky model can improve the performance of SGS parameterizations of unresolved momentum stresses). In this work, we construct a gradient model for parameterizing unresolved mesoscale eddy transport to benefit coarse-resolution ocean simulations using similarity assumption between resolved and subgrid fluxes. Concretely, we can expand the resolved velocity ̃ and thickness h at a given point x over a distance r, which is on the order of the filter scale Δ f , using a Taylor series where the tilde sign ⋅ shows resolved (coarse-grained) variables by a filter function, G is a filter function that is employed over the lateral (latitude-longitude) domain D. The filter function has a dimension of 1/L 2 where the length scale is set by the filter scale Δ f for coarse-graining from high-resolution grid spacing Δ to coarse-resolution filter scale Δ f (see below); i = {1, 2} denote zonal and meridional directions, and j is a dummy index (according to the Einstein's summation convention); n denotes the density layer. In general, G(x) can be anisotropic and non-homogeneous. Using Equations 8 and 9, and also the Taylor series expansion of total transport ̃ ℎ , we can mathematically approximate the SGS eddy transport n =̃ ℎ −̃ h in the following form (Appendix B) where the non-dimensional parameter c e , in general, depends on the filter scale, domain geometry and flow physics.
Similarly, we can develop a thickness-weighted gradient model for unresolved momentum fluxes using the Taylor series expansion of resolved and thickness-weighted velocities. In fact, the SGS momentum flux =̃ * −̃ ̃ * can be approximated as follows (see e.g., Khani & Porté-Agel, 2017b;Khani & Porté-Agel, 2022;Khani & Waite, 2020) where ( ⋅) * =h(⋅)∕h is a thickness-weighted filtering operator (see Appendix C in Khani et al., 2019, for more information); k is a dummy index. Noteworthy, the diagonal components of the unresolved momentum flux can be used to estimate the SGS KE K , written as 5 of 26 where i is a dummy index in Equation 12. It is worthwhile mentioning that the thickness weighting operator in our work is based on the spatial filtering approach. There are certain mathematical properties associated with this operator: for example, m * * ≠m * , and the residual of thickness weighting has non-zero thickness-weighted quantity (i.e., m′ * ≠ 0 where m ′ = m −m * ). Here, m is a scalar/vector variable. Alternative thickness-weighting operators are defined in other research work with their associated mean and residual properties (see e.g., Greatbatch, 1998;Smith, 1999;Young, 2012;Loose et al., 2022, for more information).
SGS parameterizations based on the mixing length hypothesis (such as the GM model) contain fairly poor local correlation between the parameterized SGS fluxes (and their divergences) with the actual SGS fluxes (and their divergences) (this weak local correlation has been shown in isotropic and stratified turbulence, see e.g., Clark et al., 1979;Khani & Porté-Agel, 2017a). Unlike mixing-length-based models, the gradient model has shown high local correlations between the actual and parameterized SGS fluxes (and their divergences) (see e.g., Khani & Porté-Agel, 2017b;Meneveau & Katz, 2000). So, it is promising to improve the parameterization of unresolved mesoscale eddy fluxes using the gradient model. Furthermore, the gradient model is not a purely dissipative scheme (unlike the GM-based model, the SGS term for the gradient model is not solely diffusive to remove energy at small resolved scales; see Appendix C). As a result, the SGS APE flux based on the gradient model includes local energizing and diffusivity exchanges based on the dynamics of resolved large-scale and parameterized unresolved-scale motions. These unique characteristics provide a framework for subgrid mesoscale eddy parameterization to include an energizing SGS term through the dynamics of the flow. Consistent with our approach, a data-driven ocean mesoscale momentum parameterization has been recently developed using a combination of gradients of resolved variables in the momentum budget (Zanna & Bolton, 2020). This data-driven SGS momentum model, which is very similar to our gradient parameterization, has improved the performance of coarse-resolution simulations by skillfully capturing large-scale ocean dynamics (Zanna & Bolton, 2020).

Simulation Setup and Methodology
Simulations from Modular Ocean Model (MOM6, see Adcroft et al., 2019) are analyzed in this study. We consider a strictly adiabatic configuration with no diapycnal mixing that allows a full equilibrium integration in 80 years of simulation time. Six constant density layers are considered from top to bottom with ρ = 1,025.5, 1,027, 1,027.5, 1,028 and 1,028.1 (kg/m 3 ), respectively. Also, the initial layer thicknesses from top to the bottom are 150, 250, 600, 1,000, 1,000 and 1,000 (m), respectively. A domain in the Southern Hemisphere, spanning on a sphere from latitudes of 65°S to 15°S and longitudes of 0°E to 90°E, is considered. The domain is periodic (reentrant) in the zonal direction. A zonally-symmetric wind stress, which varies with latitude, is distributed across the top 20 (m) of the domain. The wind stress is easterly at low latitudes and wind direction changes at θ = 33.3°S, a latitude below which the wind stress is westerly. Three ocean bathymetry configurations are considered: a flat bottom channel, a channel with bottom topography and a Southern Hemisphere idealized configuration (Neverworld). Figure 1 shows the wind stress distribution, along with the three bathymetry configurations.We have performed eddy-resolving simulations with the horizontal grid spacing of (1∕8 • ) for all configurations; in addition, an extra simulation with the horizontal grid spacing of (1∕16 • ) is performed for the Neverworld configuration. The Coriolis parameter f = 2Ω sin θ, with Ω = 7.2921 × 10 −5 (s −1 ) is considered. The model setup is the same as that is reported in Khani et al. (2019) (simulations are slightly different: here, we run an extra Neverworld simulation with a horizontal grid spacing of 1/16°, and also the time averaging is performed over years 70-80 in all simulations). In this work, we mainly focus on developing an a priori study for testing the performance of the gradient model in mesoscale eddy parameterization using eddy-resolving ocean simulations. This model setup is useful to test our new SGS parameterization for mesoscale eddy and momentum fluxes because it provides a range of simulations on configurations with and without bottom topography and/or continental barriers.
In the current work, we set non-dimensional parameters c e and c m using the trial-and-error method. We have found that c e and c m can be influenced by bottom topography if the scale of topography ℓ t ≲ Δ f . In addition, the wind stress direction that can make different types of instabilities in areas with easterly or westerly winds (see e.g., Khani et al., 2019;Tulloch et al., 2011) could affect non-dimensional parameters. Tables 1 and 2 show the non-dimensional parameters c e and c m where they mainly change with Δ f . We found that results are less biased 6 of 26 when we use slightly different c m for tuning the zonal and meridional SGS momentum fluxes ( Table 2 reports the mean of them for SGS KE). In some cases when Δ f ≳ ℓ t , c e might be different in the easterly versus westerly wind region (Table 1). The GM coefficient κ G in the flat bottom channel is set to vary with Δ f as shown in Table 3. It is worthwhile mentioning that the GM parameterization here is a "baseline" model, with a constant κ G , which most ocean simulations at non-eddying regime with Δ f ≳ 1° employ. GM and gradient models yield very different SGS spatial and temporal patterns while constant non-dimensional parameters can only change the magnitude of eddy/ momentum fluxes and not their structures (see below). The GM model in this work is a reference parameterization

A Priori Testing for SGS Parameterization
Using a spatial filtering approach, we can diagnose subgrid mesoscale eddy fluxes by employing a coarse-resolution filter function G to the results of eddy-resolving simulations (see e.g., Khani et al., 2019, for more information). Concretely, the coarse-filter resolved variable m can be obtained as follows where m is the actual variable at high resolution, A is the area of a grid cell, and G is the top-hat filter Here, A f is the area of the filter stencil and N f = Δ f /Δ is an integer. Using this approach, we can compute the actual and parameterized subgrid mesoscale eddy transport: Hereafter, superscripts "a" and "g" denote the actual and gradient-based-parameterized contributions, respectively. Furthermore, the actual and parameterized eddy-induced overturning streamfunction n ψ can be obtained as follows where superscript GM denotes the GM-based eddy parameterization. For overturning streamfunction and SGS transport vectors n ψ = n ψ 1 e x + n ψ 2 e y and = 1 + 2 , respectively, we mainly consider the meridional transport with 2 =̃ ℎ −̃ h in Equation 17, and meridional derivative of isopycnal layer in the GM model Similarly, we can compute the actual and gradient-based-parameterized thickness-weighted SGS KE, given as Here, u n represents geostrophic velocity that is calculated from the gradient of the Montgomery potential. We have tested a different version of the actual SGS KE: . We employ different filter widths (scales) to coarse-grain the eddy-resolving simulations at 1/16° and 1/8° into eddy-permitting and non-eddying simulations at 0.5°, 0.75°, 1°, 1.5°, 2°, 3°, 4°, and 5°. At each filter scale, we will compare the actual SGS fluxes with those parameterized from the gradient model to evaluate the performance of the new subgrid mesoscale eddy parameterization. We compare the results of and , K and K , and , and 12 and 12 . We evaluate these results over the latitude-longitude plane at different density layers. We also provide comparisons between zonally-integrated actual and parameterized meridional volume transports, SGS KE and subgrid PE transfers.
The non-dimensional parameters vary with resolution in the a priori study. This is due to the effect of the filter width on the SGS fluxes when we do explicit filtering since there is an aligned correlation between filter-resolved and sub-filter motions at a given filter width in high-resolution simulations. In the a posteriori study, non-dimensional parameters have to be set only by information from resolved scales because there is no SGS fluxes from explicit filtering in non-eddying simulations. Therefore, the non-dimensional parameters are more resolution (scale) independent in the a posteriori runs. This behavior also demonstrates that tuning online simulations in the a posteriori study using the non-dimensional parameters that are trained using the a priori study in offline simulations (with explicit filtering) may not yield stable numerical solutions (as shown in e.g., Zanna & Bolton, 2020). In Section 7.1, we introduce a dynamic method to determine non-dimensional parameters in the a posteriori coarse-resolution ocean simulations (see below).

Subgrid Mesoscale Eddy Transport
The overview of eddy-resolving simulations has been presented in detail in our previous work (see Section 2.1 in Khani et al., 2019, for more information). In this work, hence, we mainly focus on studying the performance of the gradient parameterization, which is an alternative model for representing subgrid mesoscale eddy fluxes, using high-resolution ocean simulations. Figure 2 shows the time-averaged and zonally-integrated meridional subgrid eddy transport in the flat bottom channel configuration at different filter scales. Different lines show the SGS eddy transport that is carried by different density layers numbering from 1 for the ocean surface layer to 6 for the ocean bottom layer. Solid lines depict the actual SGS volume transport 2 =̃ ℎ −̃ h , while dashed lines represent the parameterized SGS volume transport 2 = Δ 2 ̃ ∕ h ∕ . Here, subscript 2 denotes the meridional transport with meridional velocity v, and j = 1, 2. At the coarse-filter scale Δ f = 5°, SGS meridional eddy transports in ocean surface and bottom layers are significant (around 10 [Sv] as shown in Figure 2a; note that 1 [Sv] = 1 × 10 6 [m 3 /s]). Interestingly, SGS meridional eddy transport that is parameterized by gradient model (dashed lines) can reproduce the actual SGS meridional eddy transport (solid lines) with high accuracy even at very coarse filter scale of 5°. As we decrease the filter scale (increase the horizontal resolution from Figure 2a toward Figure 2e), the SGS eddy transport decreases because more eddy structures are explicitly resolved at smaller filter scale (higher resolution), and this decrease is at the expense of increased "resolved" eddy transport. For example, at the filter scale Δ f = 1°, magnitudes of SGS eddy transport is less than 4 [Sv] and 2 [Sv] in regions with westerly and easterly wind stresses, respectively (as shown in Figure 2e).

Gradient Model
In the channel with bottom topography configuration, the zonally-integrated SGS meridional volume transport is comparative at the coarse filter scale Δ f = 5°. The SGS transport decreases with increased filter scale from 5° to 1° (Figure 3). The parameterized time-averaged and zonally-integrated SGS meridional eddy transport ⟨ 2 ⟩ can reasonably model the actual SGS transport ⟨ 2 ⟩ in the presence of significant topography even at very coarse resolution (e.g., Δ f = 5°) in ocean surface and bottom layers (compare dashed and solid lines in Figure 3). Here, the angle bracket 〈⋅〉 denotes an integration over the zonal direction. The thick solid black line in Figure 3 shows the vertical sum of zonally-integrated SGS meridional transport over the water column. This sum exhibits the transport carried by unresolved standing meanders, generated from bottom topography, at the upper ocean through a barotropic contribution. Note that this sum is zero in the flat bottom channel configuration (as shown in Figure 2). It is explained in Khani et al. (2019) that the volume transport by unresolved standing meanders can not be parameterized using GM-based parameterizations because the GM 9 of 26 model is only developed assuming a flat bottom ocean. In other words, the vertical sum of the SGS meridional volume transport over the water column has been set to zero since the ocean surface is flat (i.e., ∕ ≈ 0 in the GM model, where denotes the resolved thickness displacement at the ocean surface; see also Appendix A). The gradient model, however, can skillfully parameterize the unresolved standing meanders; there- Time-averaged and zonally-integrated unresolved geostrophic meridional volume transport for flat bottom channel at different filter scales, ranging from 5° to 1°. Actual ⟨ 2 ⟩ = ⟨̃ ℎ −̃ h ⟩ and gradient model parameterized are shown with solid and dashed lines, respectively. Here, n denotes density layers where from top to the bottom are numbered from 1 to 6 with cyan, magenta, blue, gray, red and green colors, respectively; black lines shows the vertical sums of the actual and parameterized geostrophic volume transports over the water column (these are zero in the flat bottom channel). Small vertical lines show the location where a layer vanishes and the layer beneath is exposed to the surface. Time averaging has been performed from years 70-80 (over the last 10 years of simulations). The latitudes of zero wind is indicated. Note that 1 [Sv] = 1 × 10 6 (m 3 /s). 10 of 26 fore, the vertical sum of SGS meridional transport over the water column in the gradient model can capture the structure of unresolved meanders at the ocean surface through the barotropic transport. It is clear that the vertical sum of zonally-integrated parameterized SGS meridional transport is fairly well overlapped with that of actual SGS transport at different resolutions (compare thick black dashed and solid lines in Figure 3). topography at different filter scales, ranging from 5° to 1°. Actual ⟨ 2 ⟩ = ⟨̃ ℎ −̃ h ⟩ and gradient model parameterized ⟩ transports are shown with solid and dashed lines, respectively. Here, n denotes density layers where from top to the bottom are numbered from 1 to 6, respectively; black lines shows the vertical sums of the actual and parameterized geostrophic volume transports over the water column. Small vertical lines show the location where a layer vanishes and the layer beneath is exposed to the surface. Time averaging has been performed from years 70-80 (over the last 10 years of simulations). The latitudes of zero wind is indicated. Note that 1 [Sv]= 1 × 10 6 (m 3 /s).

of 26
In the Neverworld configuration, in which the ocean bathymetry includes continental barriers in addition to rough bottom topography, time-averaged and zonally-integrated parameterized SGS meridional volume transport ⟨ 2 ⟩ can fairly well estimate the actual SGS transport ⟨ 2 ⟩ at different density layers (see dashed and solid lines in Figure 4; they are likely overlapped at each density layer, in particular, away from continental lands). It is noteworthy mentioning that due to the explicit filtering difficulties near the land in isopycnal coordinate models (see e.g., Section 4.1 in Khani et al., 2019), the zonal integration of resolved variables' differentiation can show some wiggles in SGS transport near boundaries. This noisy transport is visible around latitudes 47°S and 33°S where continental shelves exist in the Neverworld configuration (see Figure 1d). Also, these wiggles are more significant near boundaries when Δ f ≳ 2° (not shown). As a result, we need to perform some treatment when isopycnals intersect with shelves at relatively coarse resolution ocean simulations when the scale of topography/ land is not fully resolved (some of these treatments are discussed in Khani et al., 2019; Aluie, 2019, but further ⟩ transports are shown with solid and dashed lines, respectively. Here, n denotes density layers where from top to the bottom are numbered from 1 to 6, respectively. Small vertical lines show the location where a layer vanishes and the layer beneath is exposed to the surface. Time averaging has been performed from years 70-80 (over the last 10 years of simulations). The latitudes of zero wind is indicated. Note that 1 [Sv] = 1 × 10 6 (m 3 /s).

of 26
research work is required to be done for near boundary parameterization in ocean models; this topic is beyond the scope of this work).

GM Model
As it is explained in Section 3.1, the GM model parameterizes SGS overturning streamfunction, which defines the unresolved volume transport below a given density layer in a flat bottom ocean (see Equations 17 and 18). The GM parameterization we employed in this study is the "baseline" model with the constant GM coefficient that is given in Table 3 for different filter widths. Khani et al. (2019) have diagnosed the global average κ G ≈ 1,000 (m 2 /s) as a reasonable estimate for the flat bottom channel configuration when Δ f ≳ 3° using coarse-grained ocean simulations. The time-averaged and zonally-integrated SGS streamfunction (vertical sum of SGS meridional volume transport below a density layer) along with the GM parameterization for SGS streamfunction are shown in Figure 5. In the westerly wind region (latitudes <−33.3°), zonally-integrated SGS meridional transports by top layers (layers 2, 3, and 4) in the GM model are well overlapped with the actual SGS streamfunctions at upper ocean, while the GM-parameterized and actual SGS meridional transports show some inconsistencies in SGS transports carried by layers 5 and 6 at the bottom of the ocean ( Figure 5). In the easterly wind region (latitudes >−33.3°), the parameterized zonally-integrated SGS meridional transports by the GM model at upper ocean carried by layers 2 and 3 are remarkably weaker than those of actual SGS streamfunctions. However, the GM-parameterized zonally-integrated SGS meridional transport is consistent with that of actual transport at the bottom of the ocean in the easterly region ( Figure 5). Moreover, as we increase the filter scale (decrease the horizontal grid spacing) from 5° to 1°, the SGS streamfunction decreases (Figures 5a-5e). Overall, the GM model underestimates the SGS meridional volume transport in easterly and westerly regions. Figure 6 shows the field of actual and parameterized SGS meridional transport at the surface layer in the latitude-longitude plane with Δ f = 4° (top and bottom panels show the flat bottom channel and channel with topography, respectively). Unresolved mesoscale eddies are mainly dominated around θ = −33.3°, a latitude in which the wind stress changes its direction (i.e., wind stress shear exists at 33.3°S; see Figure 1a). Moreover, eddy patterns can be influenced by the layer vanishing that occurs at around θ = −40° in the ocean top layer. The horizontal scale of unresolved (subgrid) mesoscale eddies in the flat bottom channel configuration is relatively small compared with those in channel with bottom topography case due to an efficient transformation of APE to KE in the presence of significant topography (isopycnal slopes are shallower in the presence of bottom topography as shown in Mazloff et al., 2013;Abernathey & Cessi, 2014;Khani et al., 2019). Noteworthy, SGS eddies are prominent when wind stress changes its direction ( Figure 6). The gradient model is able to estimate the scale of SGS eddies in cases with significant topography, an indication that this new model can capture unresolved standing meanders at the ocean surface even at coarse-resolution models (Figures 6c and 6d). Figure 7 shows the time-averaged SGS meridional eddy transport in the latitude-longitude plane at different density layers in the Neverworld configuration with 2° grid spacing (actual and gradient-parameterized SGS transports are shown in left and right panels, respectively). At the ocean surface, the SGS meridional transport is significant in the Western boundary current area (Figure 7a). Also, the subgrid standing eddy structures are visible in the latitude ranges 40°S to 50°S that show meridional eddy transport in Agulhas region (the South tip of the African continent in Figure 7a). Interestingly, the gradient model can successfully capture subgrid meridional eddy transports in Western boundary currents and also at the South tip of Africa (Figure 7b). Actual SGS meridional transport is weaker over the second density layer (below surface layer), while SGS eddies have similar structures. Noteworthy, the gradient model can reproduce similar SGS meridional transport below the top layer (Figures 7c and 7d). Near bottom of the ocean over the density layer 4, the SGS meridional transport is even weaker in comparison with that in top layers, but we can see remarkable SGS fluxes near the Drake Passage area (latitude band 50°S to 60°S in Figures 7e and 7f). Moreover, jet streams across the Southern boundary and along the zonal direction are seen at the bottom layers (Figure 7e). The gradient model can well parameterize the SGS meridional fluxes at the bottom of the ocean (Figures 7e and 7f).

Diagnostics of Eddy Transports in a Latitude-Longitude Map
The actual and parameterized SGS zonal eddy transport over the latitude-longitude plane in the Neverworld configuration at 2° grid spacing is shown in Figure 8. At the surface layer, unresolved eddy structures are mainly important in the Western boundary current, Agulhas region and Antarctic Circumpolar Current (ACC). The 13 of 26 SGS zonal transport in basins is notably influenced by unresolved mesoscale eddy structures at the surface layer (Figures 8a and 8b). Zonal transports of SGS eddies become relatively weaker below the surface layer and also in the bottom of the ocean (in particular at basins). Nevertheless, mesoscale eddy structures are still dominant in Here, i and n denote density layers. Sum of transports below density layers from top to the bottom are numbered from 1 to 6, respectively. SGS GM transports below the surface layer, which is shown by the cyan line, is zero due to the GM zero-barotropic transport constraint (Appendix A). Small vertical lines show the location where a density layer vanishes and the eddy-induced circulation is taken care by the layer beneath. Time averaging has been performed from years 70-80 (over the last 10 years of simulations). The latitudes of zero wind is indicated. Note that 1 [Sv] = 1 × 10 6 (m 3 /s). 14 of 26 zonal transports at the Western boundary, Agulhas and ACC and Drake Passage region in ocean sub-surface and bottom layers (Figures 8c-8f). The gradient model is able to parameterize the unresolved zonal fluxes in surface and bottom of the ocean at coarse-resolution simulations (compare parameterized SGS fluxes in right panels with those from actual SGS eddies in left panels of Figure 8).

Subgrid Kinetic Energy and Momentum Fluxes
The unresolved surface KE field at different filter scales, ranging from 0.5° to 1.5°, in the Neverworld configuration is shown in Figure 9. With Δ f = 0.5° at the eddy-permitting resolution, the SGS KE is mainly significant near the Western boundary region, however, eddy structures are also visible in the south tip of the African continent and along the ACC (Figures 9a and 9b). At non-eddying resolution with Δ f = 1° and 1.5°, where most of the energetic mesoscale eddies are not resolved, the unresolved KE becomes stronger in the Western boundary and Agulhas regions (Figures 9c-9f). Moreover, SGS KE are important in basins at coarse-resolution simulations with Δ f = 1.5°, where eddy structures influence KE and ocean circulation in the South Atlantic and Indian Ocean gyres (Figures 9e and 9f). Figure 10 shows zonally-averaged unresolved KE (SGS KE) for the Neverworld configuration at different filter scales (resolutions). The unresolved KE is significant at the ocean surface; SGS KE decreases in subsurface layers 2 and 3, and the contribution of SGS KE is relatively small in layers 5 and 6 at the bottom of the ocean (Figure 10). This trend is in line with the structure of SGS meridional volume transports in Figure 4. The unresolved KE is large at large filter scale Δ f (i.e., coarse resolution) because eddies with scales smaller than Δ f are not resolved in simulations with coarse resolution. As Δ f decreases (in high resolution simulations), n K sgs becomes smaller because most of the mesoscale eddies are resolved at finer filter scales (see Figures 10a-10c for which Δ f changes from 1.5° to 0.5°).

of 26
The subgrid horizontal momentum flux 1 τ 12 at the surface layer for the Neverworld configuration is shown in Figure 11. The weight of SGS horizontal momentum fluxes are mostly considerable around the Western boundary current, in the South tip of the Africa continent and in the ACC region. At eddy-permitting resolution with Δ f = 0.5°, where the horizontal momentum transports are better resolved, SGS fluxes are weaker than SGS momentum transports in non-eddying resolutions with Δ f = 1° and 1.5°, in which horizontal momentum fluxes are mainly unresolved (Figure 11). Also, the gradient model constructs similar patterns of unresolved momentum transport as those are shown by actual SGS horizontal momentum fluxes at different spatial resolution (compare left vs. right panels in Figure 11). Noteworthy, the non-dimensional parameter c m for the horizontal 16 of 26 SGS transport increases with resolution ( Figure 11). For actual and parameterized n τ 12 terms, we have found that the thickness-weighted resolved momentum fluxes ̃ ̃ * ≈̃ * ̃ and ̃ ̃ * ≈ ̃ * ̃ (not shown); therefore, we define 1 τ 12 fluxes as shown in Figure 11.
The gradient model can approximate unresolved KE structures fairly well at eddy-permitting and non-eddying ocean resolutions with Δ f ≳ 0.5° (compare right and left panels in Figure 9; see also dashed and solid lines in Figure 10). The estimation of unresolved KE, a.k.a. mesoscale eddy kinetic energy (MEKE), has been controversial in eddy-permitting and non-eddying ocean models, where we are required to solve an explicit trans- port equation for subgrid eddy KE in OGCM (see MEKE model in e.g., Jansen et al., 2019;Perezhogin, 2020). The gradient model suggests an alternative approach to estimate MEKE by employing a filter function at resolutions coarser than the grid spacing (i.e., Δ f > Δ) over the lateral domain to compute the thickness-weighted local eddy KE at each time step. For example, computing GM diffusivity coefficient based on the MEKE model can benefit from estimating the SGS KE using the thickness-weighted gradient model.  18 of 26

Subgrid Potential Energy Transfer
We can evaluate the performance of mesoscale eddy parameterizations by diagnosing subgrid PE transfer in the PE budget based on the mass continuity equation in the stacked shallow water model. The actual, gradient and GM PE transfer are given as follows (Appendix C) (24) Figure 12 shows the zonally-integrated actual and parameterized PE transfer at different spatial resolution in the flat bottom channel. Interestingly, the zonally-integrated SGS transfer is mainly negative (subgrid term removes PE). The gradient model well captures structures and magnitudes of integrated SGS PE transfer in line with those from the actual subgrid PE (Figures 12a and 12b). However, the magnitudes of the GM PE transfer are about 2.3 (1.5) times larger than those from the actual subgrid PE in simulations with the grid spacing of 1° (2°) (Figures 12c and 12d). The GM PE transfer structures are also a little different than those from actual PE transfer a least in upper layers (see e.g., layers 2 and 3 in Figures 12c and 12d).The patterns of gradient-parameterized and actual subgrid PE transfer at coarser resolution (Δ f ≳ 1°) are also very similar (not shown).

Discussion and Implementation
For a new SGS parameterization, it is essential to provide a systematic approach for setting in correlation coefficients. In practice though, this step could be overlooked and the approach is only limited to a trial-and-error method (e.g., most backscatter parameterizations for the KE budget employ trial-and-error to set correlation coefficients; see e.g., Leith, 1990;Bachman, 2019;Perezhogin, 2020). Correlation coefficients, in general, could depend on different factors, such as the model configuration, forcing term, physical 19 of 26 scales and grid spacing. In an a priori testing with explicit filtering, this situation becomes more complicated because correlation/non-dimensional coefficients could also be controlled by filter functions and filter widths (see e.g., Bachman, 2019;Chow et al., 2005;Khani & Porté-Agel, 2017a). These outcomes clearly demonstrate that c e and c m depend on the dynamics of the ocean flows and can vary with the model resolution and ocean bathymetry. As a result, we propose a dynamic SGS method (see e.g., Germano et al., 1991;Khani & Waite, 2015;Lilly, 1992) to set the non-dimensional coefficient c e in parameterization of mesoscale eddy fluxes using the gradient model. Note that scale-dependent non-dimensional parameters are not inconsistent with the self-similarity assumption of the gradient model because (a) self-similarity in this context refers to the structure of SGS transports while correlation coefficients (non-dimensional parameters) can only change their magnitudes (see Equations 10 and 11, and Tables 1 and 2); (b) SGS transports, by definition, depend on the grid spacing (Δ f in Equations 10 and 11), and so it is not inconsistent with the self-similarity assumption if c e and c m also vary

A Dynamic Method to Evaluate c e for Mesoscale Eddy Fluxes
If we employ a second coarse spatial filter in the horizontal direction, so called the test filter indicated by the hat sign ⋅ , we can write actual and parameterized fluxes associated with the test filter in the following forms (in analogy to Equations 15 and 16) where Δ is the test filter width, which is coarser than Δ f (i.e., Δ > Δ ), and i, j = {1, 2} stand for zonal and meridional directions, respectively. We have also assumed that c e does not change at coarser filter due to a scale similarity assumption (see e.g., Lilly, 1992). If we also directly employ the second filter Δ to Equation 15 and 16, and subtract Equations 25 and 26 from those, respectively, we can write

of 26
where unlike the first terms in the right-hand-side of Equation 15 and 25, all variables in Equations 27 and 28 are explicitly known at filter scale Δ f . We can directly evaluate c e by minimizing the error between n A i and n B i using a regression method (note that the goal in SGS modeling is to parameterize n A i using n B i ). If we define Q to be the square of the error, that is, we can use a least squares approach by setting ∂Q/∂c e = 0 to estimate c e as follows In an a posteriori study of eddy-permitting and non-eddying ocean simulations, we will use Equation 30 to set the non-dimensional parameter c e in the gradient model. The grid spacing Δ will indicate the resolution shown by ⋅ sign, and we will need to employ an explicit horizontal test filtering function with Δ > Δ at each density layer and time step to compute coarse-grained variables at test filter scale shown by ⋅ sign. As a result, we can dynamically calculate c e at each time step over the latitude-longitude plane using Equation 30. In an a posteriori testing for the MOM6 code, we will implement a subgrid mesoscale parameterization based on the gradient model, in which the time and space dependent coefficient c e is computed dynamically. In the mean APE budget, the new model can either locally dissipate APE motions at the grid scale Δ or energize mean APE, depending the sign of Π shown by Equation 23. The APE diffusive part of the gradient-based SGS model removes APE mainly at grid scale Δ. The locally energizing APE part reinforces large-scale APE (baroclinic mode) that ultimately results in large-scale APE and KE interactions through resolved motions in eddy-permitting (or non-eddying) ocean simulations. Although the net SGS PE transfer is diffusive (see Section 6), we may still face some numerical stability issues in the a posteriori study due to locally energizing SGS PE fluxes. We think the dynamic method will alleviate these numerical instabilities. Nevertheless, we might be required to reduce or clip large values of locally energizing SGS PE fluxes when we set non-dimensional parameters.

Conclusions
In this work, we have proposed a new subgrid mesoscale eddy parameterization (the gradient model) for OGCM. The new model is developed based on the Taylor series expansion of resolved velocity and thickness to parameterize SGS volume transport and momentum fluxes in the zonal and meridional directions. We have also employed the gradient model to estimate subgrid mesoscale KE. An a priori study on the gradient model using eddy-resolving ocean simulations at different ocean bathymetry and filter scales is performed. It has been shown that the gradient model can reproduce actual SGS fluxes with high accuracy even at non-eddying filter scales Δ f ≳ 1°. Zonally-integrated meridional SGS transport, latitude-longitude maps of meridional and zonal SGS fluxes, and maps of SGS KE are skillfully captured by the gradient model in configurations with and without bottom topography and continental barriers. In particular, unlike the Gent and McWilliams (1990) model that is basically developed for a flat bottom ocean (Appendix A), the gradient model can reproduce the unresolved standing meanders at the ocean surface. We have also diagnosed the subgrid PE transfer in the a priori study. It is shown that the net SGS PE transfer is diffusive in the channel and Neverworld bathymetry at different spatial resolution ( Figure 12). Overall, the net estimated subgrid PE transfer in the gradient model is consistent with the actual SGS PE transfer, whereas the baseline GM model overestimates the net subgrid PE fluxes (Figure 12).
We have introduced a dynamic method to estimate non-dimensional parameters in the gradient model based on physics of the ocean flow. This method will help in an a posteriori study of gradient model to set in the model coefficients through a non-ad hoc and tuning-free procedure in eddy-permitting and non-eddying ocean simulations. For future work, we plan to implement the gradient model in GFDL-MOM6 code to parameterize SGS eddy and momentum transports and to estimate SGS KE. We will test our new model with idealized

of 26
where c e is a non-dimensional parameter, and k is a dummy index.
Noteworthy, the last term on the right-hand-side of Equation B3 will survive if the filter function includes non-zero odd moments (e.g., in a diffusion-based spatial filter, see e.g., Grooms et al., 2021). For this case, Equation B4 will be modified as follows where  and  are linear functions of Δ f . Overall, the second term in the right-hand-side of Equation B6 represents anisotropic features of the spatial filter function while the leading order of this term is still ∼ Δ 2 because both  and  are linear functions of Δ f . Similar to c e in Equation 30, we can determine the anisotropic correlation coefficients  and  using the dynamical method over the spatial domain (we leave this topic for future investigation). Note that the accuracy of the gradient model with diffusion-based filtering is  ( Δ 3 ) as shown in (B6).

of 26
we can show that the SGS PE transfer in GM model in Equation C6 is equivalent to a SGS PE transfer based on interface diffusion when the bottom of the ocean is flat. In the presence of bottom topography, however, the GM model based on thickness diffusion in Equation C4 may result in unphysical structures at ocean interfaces (Bleck et al., 1992;Holloway, 1997). As a result, we may employ the SGS PE transfer in the GM model based on the interface diffusion in configurations with a bathymetric map (as used in e.g., Loose et al., 2022).
With a gradient-based parameterization, the PE budget can be written as follows (again, we assume periodic boundary conditions in the lateral direction) where ̃ = (1∕2) ̃ ∕ + ̃ ∕ is the strain-rate tensor. It is clear that, with a positive κ g , the sign of SGS term in the right-hand-side (C7) is still not determined and it depends on the zonal and meridional slopes of thickness and the rate-of-strain tensor (i.e., energy interactions based on the flow dynamics). Assuming κ g > 0, a negative Π shows removing of PE (diffusive) while a positive Π shows energizing resolved motions by reinforcing large-scale PE. Overall, the gradient-model Equation C7 can hold both negative and positive SGS PE transfer, whereas the GM-model Equation C6 only shows negative subgrid PE transfer. The parameterized SGS terms Π and Π should be compared with the actual subgrid PE transfer which also can be negative and positive depending on the sign of SGS eddy diffusivity and the zonal and meridional slopes of the thickness.