Parameterization of Submesoscale Symmetric Instability in Dense Flows Along Topography

We develop a parameterization for representing the effects of submesoscale symmetric instability (SI) in the ocean interior. SI may contribute to water mass modification and mesoscale energy dissipation in flow systems throughout the World Ocean. Dense gravity currents forced by surface buoyancy loss over shallow shelves are a particularly compelling test case, as they are characterized by density fronts and shears susceptible to a wide range of submesoscale instabilities. We present idealized experiments of Arctic shelf overflows employing the GFDL‐MOM6 in z* and isopycnal coordinates. At the highest resolutions, the dense flow undergoes geostrophic adjustment and forms bottom‐ and surface‐intensified jets. The density front along the topography combined with geostrophic shear initiates SI, leading to onset of secondary shear instability, dissipation of geostrophic energy, and turbulent mixing. We explore the impact of vertical coordinate, resolution, and parameterization of shear‐driven mixing on the representation of water mass transformation. We find that in isopycnal and low‐resolution z* simulations, limited vertical resolution leads to inadequate representation of diapycnal mixing. This motivates our development of a parameterization for SI‐driven turbulence. The parameterization is based on identifying unstable regions through a balanced Richardson number criterion and slumping isopycnals toward a balanced state. The potential energy extracted from the large‐scale flow is assumed to correspond to the kinetic energy of SI which is dissipated through shear mixing. Parameterizing submesoscale instabilities by combining isopycnal slumping with diapycnal mixing becomes crucial as ocean models move toward resolving mesoscale eddies and fronts but not the submesoscale phenomena they host.

Plain Language Summary When developing numerical ocean models, processes occurring on scales smaller than the grid size must be approximated in terms of the resolved flow. The term "parameterization" refers to this approximation of small-scale features and is essential for representing turbulent mixing. We consider the effect of a particularly ubiquitous small-scale turbulent process known as symmetric instability (SI). SI occurs throughout the World Ocean and is important in setting oceanic properties through mixing and maintaining energy balance. SI is common in fronts, such as those arising from dense currents known as overflows. Overflows often originate in polar continental margins through cooling and secretion of dense brines as sea ice grows. As the dense waters flow offshore along the seafloor, they become susceptible to small-scale instabilities such as SI. Although crucial for maintaining the density structure of the ocean, SI is presently unresolved in global ocean models. We develop a parameterization for SI using the test case of an Arctic shelf overflow. We test the scheme in various overflow simulations and find it to successfully capture the effects of SI. The need for such a parameterization emerges as models move toward resolving increasingly finer-scale flows but not the small-scale turbulent mixing within them.
YANKOVSKY ET AL. resentations, or parameterizations, for the unresolved turbulent flows. Parameterizing submesoscale turbulence is particularly challenging and urgent in polar oceans-partly because submesoscale phenomena occur on smaller scales at higher latitudes (due to larger f ) and partly because these are the most rapidly changing, climatically significant, and vulnerable regions of our planet (Barnes & Tarling, 2017).
The ocean is dominated by horizontal large-scale current systems and mesoscale flow features, following the paradigm of two-dimensional (2D) turbulence which exhibits an inverse energy cascade to larger scales. In order to maintain an energy equilibrium, mesoscale kinetic energy must be extracted by submesoscale motions, transferring energy downscale to molecular dissipation (Gula et al., 2016;McWilliams et al., 1998). In particular, oceanic fronts-owing to their significant horizontal density gradients, vertical velocity shears, and Ro, Ri ∼ 1-are hotspots for a wide suite of submesoscale processes proposed as conduits for mesoscale energy dissipation (D'Asaro et al., 2011;Molemaker et al., 2010). Numerous theoretical and modeling studies have examined submesoscale turbulence in oceanic fronts stemming from phenomena such as inertial and symmetric instability (SI; Grisouard, 2018;Taylor & Ferrari, 2009), internal wave interactions (Grisouard & Thomas, 2015;Thomas, 2017), mixed layer eddies (Boccaletti et al., 2007;Fox-Kemper et al., 2008), and bottom boundary layer baroclinic instability (Wenegrat et al., 2018). Observations indicate that SI is particularly ubiquitous, occurring in bottom boundary layers (Wenegrat & Thomas, 2020), boundary currents such as the Gulf Stream (Thomas et al., 2013), abyssal flows in the Southern Ocean (Garabato et al., 2019), the Antarctic Circumpolar Current (Ruan et al., 2017;Viglione et al., 2018), and in outflows from the rapidly melting Antarctic ice shelves (Garabato et al., 2017). SI is a globally significant contributor to water mass properties and the energy budget.
Although submesoscale dynamics are crucial components of the ocean circulation, they are unresolved by modern ocean GCMs. Significant work aims to independently develop mesoscale eddy parameterizations as well as subgrid-scale diabatic mixing schemes. However, there have been relatively few attempts to link these processes-that is, represent mesoscale energy loss as a source for irreversible diabatic mixing (an energetic pathway that submesoscale motions may initiate). Mesoscale eddy parameterizations are generally based on the stream function developed by Gent and McWilliams (1990) and Gent et al. (1995), hereinafter referred to as "GM." The premise of GM is to parameterize adiabatic eddy-induced stirring processes by slumping isopycnals according to an eddy diffusivity . The potential energy (PE) released by the isopycnal slumping is not reintroduced into the flow and is assumed to be viscously dissipated without diapycnal mixing-an inaccurate assumption for the real ocean (Tandon & Garrett, 1996). Some studies have sought energetic consistency by (1) reinjecting kinetic energy into the resolved system via a backscatter approach (Bachman, 2019;Jansen & Held, 2014) and (2) parameterizing energy cascade to mixing via Lee waves and internal wave interactions (Eden et al., 2014;Melet et al., 2015;Saenko et al., 2011).
For subgrid-scale diabatic mixing, schemes such as the K-Profile Parameterization (KPP) of Large et al. (1994) are utilized (Roekel et al., 2018). The interior part of KPP represents shear-driven mixing outside the surface mixed layer, similar to the scheme of Pacanowski and Philander (1981); however, both rely on dimensional constants which must be calibrated. Jackson et al. (2008) propose an implicit scheme based on a critical Ri criterion and turbulence decay scale which successfully represents shear-driven, stratified turbulent mixing for various flow scenarios. Similar to early shear mixing schemes, existing SI and submesoscale baroclinic eddy schemes are specialized to certain regions-for example, the forcing-dependent mixed layer SI scheme of Bachman et al. (2017)-and often rely on dimensional parameters. Our aim is to develop a general, implicit, and easily implementable parameterization linking mesoscale energy loss by submesoscale isopycnal slumping with diabatic mixing, capturing the effects of submesoscale SI-driven turbulence. The need for such a scheme emerges as regional and global ocean models approach resolving mesoscale fronts, but not the submesoscale phenomena they host.
The Modular Ocean Model version 6 (MOM6) developed within the Geophysical Fluid Dynamics Laboratory (GFDL) is used in this study. Presently, MOM6 includes parameterizations for the surface and bottom boundary layer, shear mixing according to Jackson et al. (2008), submesoscale mixed layer instabilities according to Fox-Kemper et al. (2011), and transient mesoscale eddies; see Adcroft et al. (2019) for details. However, there is no scheme for representing submesoscale turbulence that may be implemented implicitly for the entire water column-such a parameterization is the objective of this work. We aim to parameterize the effects of pure SI modes, although the resulting scheme may extend to other forms of submesoscale turbulence. We develop the parameterization based on a test case of a 2D symmetrically unstable front arising from a rotating gravity current characteristic of the Arctic Ocean, analogous to the case studied by Yankovsky and , hereinafter referred to as YL2019. Dense gravity currents, also known as overflows, forced by surface buoyancy loss over shallow shelf regions are important contributors to subsurface and abyssal ventilation throughout the World Ocean, yet remain challenging to represent accurately in models (Legg et al., 2009;Snow et al., 2015). Given their characteristic frontal dynamics, complex submesoscale nature, and poor representation in GCMs, dense overflows are a particularly compelling test case for the development of this scheme.
We begin by examining idealized numerical simulations of an overflow that reveal the need for an SI parameterization in a model that resolves a mesoscale front but not the submesoscale dynamics evolving from it. We employ the existing parameterizations in MOM6 and consider two coordinate systems (z* and isopycnal) at various resolutions. In both coordinate systems, when SI is unresolved the water mass modification processes and overflow dynamics are inaccurately represented. We then present the theoretical basis and implementation of the proposed parameterization. Finally, we test and discuss the scheme's performance in z* and isopycnal coordinates. Overall, we find the parameterization to perform remarkably well in representing the effects of submesoscale SI and the resulting turbulence at resolutions that do not explicitly resolve these processes.

Motivation
The motivation for this study stems from a prior work (YL2019) where we identified submesoscale SI as the dominant mechanism leading to turbulent mixing and dissipation of geostrophic energy for a rotating dense overflow. In YL2019, the nonhydrostatic z-coordinate MITgcm  was applied to 2D and three-dimensional (3D) simulations to examine the dynamics of a gravity current representative of shelf overflows originating in the Barents and Kara Seas of the Arctic Ocean. The simulations consisted of an idealized domain with a continental shelf region experiencing negative buoyancy forcing in the form of a heat flux out of the water and a salt flux into the water, representing the effects of cooling and ice formation leading to brine rejection (see Figure 1 of YL2019).
In both 2D and 3D cases, the dense water flows offshore and down the shelfbreak, undergoes geostrophic adjustment, and leads to development of bottom-and surface-intensified jets. The jets descend along the slope through Ekman drainage (Manucharyan et al., 2014), creating a combination of a density front along the topography and geostrophic velocity shear in the vertical (Figure 4 of YL2019 and Figure 1 of this work). SI is initiated, manifesting as small-scale diagonal motions along the front, and leading to secondary Kelvin-Helmholtz shear instability which ultimately creates irreversible mixing and geostrophic energy dissipation. In 3D cases, the jets are baroclinically unstable, but nonetheless SI is prevalent in the bottom boundary and along eddy edges (Figure 12, YL2019). Here, we explore an analogous setup within the hydrostatic MOM6 to test whether the coordinate system and parameterization choices impact the observed dynamics. The results of these simulations demonstrate the need for an SI parameterization.

Model Description
The numerical ocean code used in this study is the GFDL-MOM6. The dynamical core of MOM6 solves the hydrostatic primitive equations formulated in a generalized vertical coordinate form (Adcroft et al., 2019); a variant of the Arbitrary Lagrangian Eulerian method is employed, allowing the use of isopycnal, z*, or hybrid coordinates. Here, we present simulations in z* and isopycnal layer coordinates based on the YL2019 overflow test case. We assume an f-plane with f = 1.43 × 10 −4 s −1 and a nonlinear equation of state (Wright, 1997). Laplacian and biharmonic viscosities, with background values of 1 × 10 −4 m 2 /s and 1 × 10 −4 m 4 /s (respectively) and velocity scales of 1 × 10 −3 m/s, and a Smagorinsky viscosity  with a nondimensional constant of 0.15 are applied. The horizontal isopycnal-height diffusivity and epipycnal tracer diffusivity are set to 1 × 10 −4 m 2 /s and the vertical background diapycnal diffusivity is 1 × 10 −5 m 2 /s. The background values of the horizontal and vertical diffusivities are relatively small and found to have negligible impacts on the flow. The Jackson shear mixing parameterization (Jackson et al., 2008) is used with its default values to represent diapycnal mixing.  Simulations are performed to 80 days, although low-resolution cases are extended to 120 days (steady state is achieved more slowly at lower resolutions). The size of the domain is 80 km in the across-shore x-direction, 2,500 m in depth (z), and for 3D simulations, 100 km in the alongshore y-direction. The 2D nominal resolution case (similar to YL2019) has dx = 125 m and the 3D nominal resolution case has dx = dy = 200 m. In z* coordinates, all cases have 120 vertical layers, with dz = 20.8 m. In isopycnal layer coordinates, there are also 120 layers which are defined linearly in density space. The final potential density distribution (referenced to 0 dbar) of the z* case at 80 days is first computed, then 120 linearly spaced values spanning this range are used to define the isopycnal coordinates (assuming the final density range is independent of coordinate choice). Infinitesimally thin layers represent the densities not present in the initial conditions, accounting for the new density classes created by negative buoyancy forcing in the shelf region. The dense overflow will not be properly resolved if the higher-density classes are unaccounted for in the initial coordinates.
There is a free-slip bottom boundary condition, with linear bottom drag and a dimensionless drag coefficient of 0.003. Other choices for the bottom boundary were explored (including the no-slip condition used in YL2019) but found to produce negligible difference in the results. Boundary conditions are periodic in the y-direction and a sponge is applied in the 10 km offshore edge in x, damping velocities to zero and tracers to their initial values. The model begins from rest and is forced identically to YL2019. A heat flux of 500 W/m 2 out of the water (corresponding to buoyancy forcing of roughly −5 × 10 −6 kg m −2 s −1 ) and a salinity forcing of −3 × 10 −5 kg m −2 s −1 prescribed in terms of an evaporative flux are applied over a 15 km shelf region. As in YL2019, the initial temperature and salinity stratification are based on observations off the Kara and Barents shelves (Rudels et al., 2000). A passive tracer, analogous to a dye, is introduced to track dense fluid as it moves offshore-its values are set to 1.0 at the surface of the forcing region at every time step and damped to zero in the offshore sponge. For a diagram of the simulation domain and initial conditions, see Figure 1 of YL2019. Figure 1 shows the 2D results at 80 days for the z* (left column) and isopycnal layer (right column) coordinate configurations, with vertical coordinate surfaces in black. In the z* case, results are consistent with the nonhydrostatic MITgcm results of YL2019. In the alongshore velocity, we see the bottom-and surface-intensified geostrophic jets formed by the dense outflow being deflected by rotation near the bottom and return flow near the surface. By 80 days, the jets have descended to the bottom of the domain through Ekman drainage and established a velocity shear in the vertical. The tracer concentration and potential density show that the dense water contained within the lower jet has created a dense front adjacent to the slope. The offshore velocity shows the characteristic signature of SI-diagonal velocity beams oriented parallel to the density front. SI sets up small-scale velocity gradients which lead to turbulent dissipation and irreversible mixing; consistent with YL2019.

Results
The primary challenge in the isopycnal layer system is selecting density coordinates to capture both the broad, temporally evolving density structure of the overflow near the surface as well as in the poorly stratified abyssal regions. Due to the surface buoyancy forcing the final density range is much larger than the initial; in linearly spaced density coordinates only 10 layers are initially filled while the remaining 110 are infinitesimally thin and only grow as dense water forms on the shelf. As a result, there is low vertical resolution in regions of low stratification and disproportionately high resolution on the shelf. As is seen in Figure 1, the abyssal ocean has layer thicknesses of nearly 1 km, while many of the high density layers onshore remain infinitesimally thin due to the relatively small volume of dense water and its partitioning into 110 layers. Several other density coordinate schemes were attempted to maximize resolution in various density classes (not feasible in a GCM, where coordinates must be chosen with the entire ocean in mind rather than a local density profile), but all shared the same problem of underresolving either the overflow or the abyss.
Hybrid isopycnal-coordinate models, like the MOM6-based OM4 global ocean model, can avoid the issue of excessively thick layers in weakly stratified water by using a density-like coordinate with an additional compressibility (Adcroft et al., 2019), but we have chosen to use a pure isopycnal coordinate here to illustrate the challenges of representing SI in their most extreme form. Another challenge in the isopycnal coordinates is that certain layers near the surface are filled more rapidly than others, leading to very steep or vertical isopycnals in the shallow shelf region ( Figure 1). As there is no implemented frontal mixing scheme operating in the interior of the water column (the shear mixing scheme only operates on vertical gradients), these horizontal density fronts continue to grow, leading to extreme velocities and numerical divergence. In the abyss, the overly thick layers do not approach resolving submesoscale SI. As a result, the density structure and velocities are erroneous compared to the z* and MITgcm results.
The 3D results shown in Figure 2 further elucidate the problem. Generally, isopycnal-coordinate systems are considered superior for representing overflows (Legg et al., 2006;Winton et al., 1998), as advection in isopycnal coordinates lacks the spurious diapycnal mixing present in z*  and the overflow is able to preserve its density structure as it propagates away from its origin. Comparing the density and passive tracer fields in Figure 2, we see that indeed the overflow is significantly more diffuse in the z* than in the isopycnal layer case. In z*, there are relatively high values of parameterized shear diffusivity adjacent to the slope while in the isopycnal case the values are very low or zero below the near surface. The shear mixing parameterization relies on a Ri criterion to determine where mixing takes place-since vertical gradients are not well captured within the thick isopycnal layers the parameterization is not activated. Thus, although the isopycnal model preserves the density structure of the overflow, there is a lack of representation of water mass modification. The observed lack of frontal mixing motivates the need for parameterization of submesoscale processes, such as SI and its secondary shear instability, that dissipate mesoscale energy and lead to irreversible mixing when resolutions are insufficient to adequately resolve them.

Parameterization for Symmetric Instability
Here, we discuss the relevant theoretical properties of SI and its effects on a geostrophic front, the parameter choices for our scheme, derivation of the stream function, and implementation in the GFDL-MOM6. Our parameterization is aimed at representing the effects of SI in a way that may be implicitly implemented for both surface and deep/interior ocean regions. The scheme is composed of four steps, detailed below.
(1) Identifying unstable regions based on a Richardson number criterion; slumping isopycnals toward a symmetrically stable state. PE released by the isopycnal slumping is calculated.
(2) Assuming conversion of the PE into turbulent kinetic energy (TKE) of the ageostrophic SI perturbations, which grow to finite amplitude, initiate secondary Kelvin-Helmholtz instability, and lead to energy dissipation and diapycnal mixing. Note that in 1, PE is removed from the geostrophic flow, whereas pure SI theory states that KE should be transferred from geostrophic to ageostrophic motions. However, for a geostrophic flow, if PE is reduced, KE will also be reduced and vice versa (based on thermal wind). We aim to represent the end state, which is similar between the two pathways. (3) Calculating diffusivity from the TKE production rate using the Osborn relation (Osborn, 1980). The diffusivity represents the influence of the secondary shear instabilities that SI initiates (Taylor & Ferrari, 2009;Yankovsky & Legg, 2019). (4) Diffusing temperature, salinity, and tracers according to the computed vertical diffusivity. To represent the effect of the secondary instabilities on mixing of momentum, a viscosity component is added to the background viscosity. Here, we employ a Prandtl number of 1 so that the SI viscosity and diffusivity values are equal.

Theory
Pure SI occurs in a flow that is both in hydrostatic and in geostrophic equilibrium (gravitationally and inertially stable), or equivalently, in thermal wind balance. Then, the SI criterion is that Ertel potential vorticity (PV) q, defined as takes an opposite sign to the Coriolis parameter f, so that fq < 0 (Hoskins, 1974). Here, k is the unit vector in the vertical, u is the 3D velocity vector (u, v, w), buoyancy is b = −gρ/ρ 0 , g is gravitational acceleration, ρ is potential density referenced to 0 dbar, and ρ 0 is a reference potential density. For a flow in thermal wind balance Figure 2. Comparison of 60-day fields for the 3D z* (left column) and 3D isopycnal layer (right column) coordinate configurations. First row: alongshoreaveraged potential density with every second vertical layer outlined in black for the z* case and every layer for the isopycnal case. Second row: passive tracer isosurfaces ranging from 1.0 to 0.2 with increments of 0.1 and becoming more transparent as the value decreases. Third row: parameterized shear diffusivity according to the Jackson et al. (2008) shear mixing parameterization.
where u g is geostrophic velocity and ∇ h b is the horizontal buoyancy gradient. Taking ζ a as the vertical component of absolute vorticity, we then rewrite the SI criterion as in Bachman et al. (2017): There are three pure modes of instability that correspond to fq < 0 being satisfied. The first two occur when fζ a N 2 is negative and larger in magnitude than  2 h b . Pure convective instability is the case of N 2 < 0 with fζ a > 0 and pure inertial instability (InI) has N 2 > 0 and fζ a < 0. The third case, pure SI, involves an inertially and convectively stable state (fζ a N 2 > 0) with the second (baroclinic) term  2 h b having a larger magnitude than the first. We may formulate the instability criterion in terms of the balanced Richardson number, equivalent to the Richardson number Ri for a flow in thermal wind balance. The criterion becomes Assuming that planetary vorticity f dominates over the relative vorticity allows us the simplified criterion of Ri B = Ri < 1. Although the criterion Ri B < f/ζ a is more general and applicable to flows with significant relative vorticity we employ Ri B < 1 for simplicity, similar to prior studies on near-geostrophic flows (such as Stone, 1966;Taylor & Ferrari, 2010). Further, using a constant value of 1.0 allows us to neglect nonlinearities within a model time step potentially arising from a criterion value dependent on the flow state. Thus, the criterion for SI we utilize here (further justified in the next section) is that Ri B < 1. An important caveat to keep in mind in addition to the negligible relative vorticity assumption is that fq < 0 on its own is a necessary but not sufficient condition for SI. Some regimes in the bottom boundary layer are characterized by buoyancy fluxes remaining large enough relative to shear production that SI does not develop, even when fq < 0 (Wenegrat & Thomas, 2020).
Real oceanic fronts are often characterized by hybrids of InI and SI, with the pure modes being hard to distinguish as they have similar effects on the flow and their precise definitions vary between studies (Grisouard, 2018). A traditional energetic view defines SI as along-isopycnal motions that grow through extraction of TKE from vertical shear, with a rate given by the geostrophic shear production (GSP) term (Thomas et al., 2013;Wenegrat & Thomas, 2020): An overline denotes a spatial average over the SI scale and primes are deviations from the average. As SI extracts energy from the flow, geostrophic adjustment leads to isopycnal slumping and weakening of the front (Bachman et al., 2017;Salmon, 1998). Examining this process in the surface mixed layer, Haine and Marshall (1998) find that SI is able to restratify on time scales faster than traditional baroclinic instability.
There is also increasing evidence that direct extraction of PE from geostrophic currents is a significant energy source for the growth of InI-SI (Grisouard, 2018;Grisouard & Zemskova, 2020). Bachman and Taylor (2014) consider the linearized primitive equations to solve for growth rates of SI modes. In the hydrostatic limit, the fastest-growing mode is indeed aligned along isopycnals; not the case for the nonhydrostatic limit, where it is shallower than isopycnal slope. Symmetrically unstable slopes form a wedge centered about the isopycnal slope, with SI gaining energy differently depending on the part of the wedge. Figures 1 and 2 in Bachman and Taylor (2014) illustrate the three energetic zones where SI gains energy from (1) geostrophic shear, (2) PE and geostrophic shear, and (3) PE.
Although the precise energetic transfers involved in SI (and its hybrid instabilities) are still an area of active research, here we consider SI to lead to isopycnal slumping and restratification toward a state where Ri B = 1 either by GSP combined with geostrophic adjustment or directly through PE extraction. We choose to represent only the final state (similar between various energetic pathways)-noting that for a geostrophic flow, reduction in isopycnal slope (loss of PE) is tied to a reduction in geostrophic KE. The ageostrophic velocity perturbations of SI also initiate secondary Kelvin-Helmholtz shear instability once they reach finite amplitudes, leading to energy dissipation and small-scale turbulent mixing (Taylor & Ferrari, 2009;Yankovsky & Legg, 2019). In the present parameterization, we consider (1) the initially unstable state defined by Ri B < 1 and (2) the final state by which SI has fully developed, extracted energy from the geostrophically balanced flow leading to isopycnal slumping toward an Ri B = 1 state (directly draining PE, or indirectly removing TKE and leading to geostrophic adjustment), and initiated secondary shear instability with resultant diapycnal mixing.

Parameter Choice
In the first step of the parameterization, we identify regions that are unstable to SI. The two equivalent criteria for instability are As shown in the Motivation (Section 2.2), one of the challenges in isopycnal layer coordinates is the lack of vertical resolution in regions that are poorly stratified. In both 2D and 3D isopycnal cases, the shear mixing parameterization fails to turn on below the well-stratified surface layers, leading to a lack of parameterized water mass modification. The shear mixing parameterization is based on critical Ri values for shear instability, relying solely on vertical density and velocity gradients. However, by using Ri B , this issue is ameliorated as the horizontal density gradients (which are better resolved) are utilized. The Ri B criterion may be formulated using the horizontal buoyancy gradient and isopycnal slope (Equation 9) which are quantities already defined in the model. We therefore propose Ri B as the parameter of choice in identifying unstable regions.
We test this criterion by examining existing 2D z* and isopycnal-coordinate system results to see how Ri and Ri B compare in identifying SI regions. Figure 3 shows a comparison between the two coordinate system results. In the top panel, regions of negative Ertel PV are shown-in z* coordinates the SI is well resolved, with the characteristic negative PV beams first noted in YL2019. The second panel shows regions of the resultant secondary shear instability (Ri < 0.25) which are again well represented in z*. The lower panels show the two Richardson number criteria. In z* coordinates, these give nearly identical results-as expected, since the vertical and horizontal gradients are both well resolved. In the isopycnal layer case, the SI and resultant shear instability are unresolved. Ri B is superior to Ri in identifying regions where the SI should be evolving (along the topography and front, as in z*).

Proposed Stream Function
Here, we present the derivation of the stream function for the proposed SI parameterization. The first step is to slump initially unstable isopycnals toward a state in which Ri B = 1. The isopycnal slope, S, is given by Ri B may be rewritten in terms of S as The criterion for instability in which isopycnal slumping will be implemented is the case of Ri B < 1. If |S| > |f/N| then Ri B < 1 and the system is considered unstable, while if |S| < |f/N| the system is stable. For unstable slopes, the isopycnal will be slumped from S toward the value of f/N. The time scale τ over which the slumping will be applied is chosen to be the ratio of buoyancy frequency to horizontal buoyancy gradient: The time scale sets how quickly we restore the isopycnal to a stable state of Ri B = 1 (Ertel PV to zero). The given form of τ leads to time scales longer than that of the SI growth scale (|f/∇ h b|), as |N| is typically larger than |f|. We find our definition of τ to be consistent with results presented in YL2019 and the high-resolution results of this study. Frontal regions have |N| ∼ 10 −3 , ∇ h b ∼ 1-5 × 10 −7 yielding τ ∼ 0.5-3 h (using |f/∇ h b| yields an order of magnitude smaller τ). Additionally, the present form leads to a stream function which is a function of the isopycnal slope and similar to the stream function already implemented for GM (allowing for increased computational efficiency). We can then approximate the time rate of change in slope magnitude as Note that we now have the isopycnal slope magnitude S multiplied by   h b f N as the rate of change of slope magnitude. This quantity is negative definite if the system is unstable to SI, so that the slope magnitude decreases with time. When implementing the parameterization we include a maximum argument so that if the system is stable, then there will be no change in slope: Note that for stable cases where goes to zero. The rate of change of slope should be positive for negative slopes (magnitude of the negative slope decreases, becoming increasingly positive), and negative for positive slopes. So, the final equation for the rate of change of isopycnal slope is given by The Gent-McWilliams (GM) stream function formulation Gent & McWilliams, 1990;Gent et al., 1995) is expressed as Here, ẑ is the unit vector in z, and κ GM is the isopycnal-height diffusivity parameterizing the effects of mesoscale baroclinic eddies and scales as (Visbeck et al., 1997) Here, α is a scaling factor, l is the lengthscale of the instability, and T is the time scale, which may be taken as the Eady growth rate for baroclinic instability  / (0.3 ) T Ri f . Rewriting the expression for diffusivity we obtain YANKOVSKY ET AL.

10.1029/2020MS002264
The expression for the GM stream function then takes the form We formulate the expression for the proposed SI stream function Ψ SI in an analogous way to GM to facilitate its implementation into the model: A lengthscale, R, is chosen to equal the horizontal grid spacing (dx or dy depending on the component). We assume that the submesoscale SI we aim to parameterize occurs at and/or below the grid scale, initiates secondary shear instability, and results in energy dissipation and mixing at the grid scale. After applying the stream function based on Equation 19, we compute the change in PE due to the isopycnal slumping. We assume that this PE is converted to TKE of the finite amplitude SI motions that initiate a forward energy cascade leading to local dissipation and diapycnal mixing by secondary shear instability: We then derive a diffusivity κ from ΔTKE by assuming a balance between the rate of TKE production and the loss of TKE to viscous turbulent dissipation (ϵ) and mixing (κN 2 ): According to the Osborn model (Osborn, 1980), the turbulent mixing rate in a stratified fluid may be estimated as Here, ϵ is TKE dissipation and Γ is a mixing efficiency based on a flux Richardson number R f ; common practice for mixing in a stratified fluid involves assuming a constant R f = 0.17, yielding Γ = 0.2 (Bluteau et al., 2013). We follow the convention of assuming Γ = 0.2 and employ the Osborn relation to rewrite the dissipation term: Solving for κ we obtain The scaling factor β ranges from 0 to 1, where 0 assumes purely viscous energy dissipation with no associated diapycnal mixing (as in GM), while 1 assumes that all of the energy is converted to TKE of the SI and leads to local diapycnal mixing through the resulting secondary shear instability. In the case where β = 0, our parameterization is therefore similar to GM, with the difference that it slumps isopycnals on smaller scales determined by the Ri B criterion. In the case where β = 1, the TKE of the submesoscale SI is transferred entirely to local mixing. By convention, using          Γ 1 Γ with Γ = 0.2 yields β = 0.167-this value will be used throughout the present study. Note that this simplification assumes the mixing efficiency, which is not yet entirely understood for SI, corresponds directly to that of secondary shear instabilities (for which Γ = 0.2 is also a simplification). As in Melet et al. (2012), we scale β by N 2 /(N 2 + Ω 2 ), where Ω is the angular velocity of the Earth, to ensure κ remains bounded when stratification is small. In the final step of the scheme, temperature, salinity, and passive tracers are diffused vertically based on the calculated diffusivity. We impose a Prandtl number (Pr) criterion, setting Pr = 1 so that the momentum diffusivity (viscosity) is equal to the computed tracer diffusivity, accounting for mixing of momentum in shear-driven instability. In our scheme, Pr is a user-defined parameter; a more physically justified estimate may be obtained by using the relation between the dissipation and diffusivity: In the simulations presented here, the traditional Richardson number Ri is poorly constrained and we avoid employing it as a parameter within the SI scheme. However, the value of Ri B and Ri in the high-resolution z* case is in the 0.2-1.0 range within symmetrically unstable regions, so that ν is indeed comparable to κ.
Thus, the parameterization has two distinct components-isopycnal slumping and vertical mixing-representing unique processes initiated by SI (removal of energy from a thermally wind balanced front and the cascade of that energy into dissipation/mixing by secondary shear instability). The proper representation of the net effect of SI relies on the two components implemented together.

Implementation Within the GFDL-MOM6
The four steps of the proposed SI parameterization are summarized in Figure 4. Isopycnal slumping according to Equation 19 defines the SI stream function in the same form as the GM stream function implemented into the mesoscale eddy closure module in the MOM6 source code. Ψ SI is added to the module by the same methodology as Ψ GM . The zonal (x-direction) and meridional (y-direction) components of the stream function are first computed independently. As derived, Ψ SI goes to zero in the symmetrically stable limit, where In the unstable case, to prevent division by zero as N → 0 we modify  h b N by adding an extra term in the denominator (here written for the x-direction, analogous in y): We justify this correction term by noting that isopycnal slopes are assumed to be much smaller than 1 according to the hydrostatic assumption employed in MOM6. Generally,     / h b b z, and the correction term is insignificant.
The zonal and meridional transports are computed for each model layer and limiting is applied based on the mass available in the two neighboring grid cells. The SI stream function has the effect of decreasing the slope of isopycnals, thus releasing PE (taken as positive). The PE release is computed at each layer interface and every horizontal grid cell. In localized regions with negative values and columns where the net PE release is negative, such as convectively unstable regions, the PE values are zeroed out and the parameterization is not applied. In the next step, we assume that the PE is converted to TKE of the small-scale slantwise motions associated with the SI and then dissipated/mixed by secondary shear instability. The fraction of TKE that leads to local diapycnal mixing is controlled by a user-defined parameter β, which is set to 0.167 here.
In YL2019, SI and the resultant shear-driven mixing are indeed found to be highly efficient at dissipating geostrophic energy and initiating water mass modification in the absence of nonlocal turbulent processes such as breaking internal waves. In the final step, an existing MOM6 subroutine is used to convert the TKE into vertical diffusivity at each model grid point where parameterized slumping occurred. The salinity, temperature, and tracers are diffused according to these values. Note that in MOM6 the vertical diffusivity is split into an along-isopycnal and a diapycnal component; the diapycnal component typically dominates (thus we use "diapycnal" and "vertical" interchangeably here). We assume Pr = 1 (Pr may be modified as an input parameter), adding a viscosity equal to the derived diffusivity to the background viscosity value at each point where the parameterization was applied.
Our scheme's stream function Ψ SI is implemented (and proportionally limited) alongside Ψ GM , allowing both the GM and SI schemes to operate simultaneously and in a scale-aware manner. In the limit where eddies are unresolved (where in principle the SI scheme should be off), the value of κ GM (a mesoscale diffusivity) will be significantly larger that κ SI (a submeoscale diffusivity). Additionally, κ SI may be modified by the user to a larger value to test the effects of representing energy loss to mixing even when mesoscale eddies are unresolved (not explored in this work).

Results
Here, we present 2D and 3D results testing the effects of the proposed SI parameterization. If the vertical and/or horizontal grid spacing is insufficient to resolve the velocity shears arising through SI, then the existing shear mixing scheme fails to capture the relevant frontal mixing. In order to test our scheme, we first consider high-resolution 2D z* coordinate cases in which the SI and water mass modification processes are well represented. We degrade the resolution and test whether our scheme yields improvements as the flow evolves. We then move to isopycnal layer coordinates where the vertical resolution in the abyssal ocean is even lower than the degraded z* cases. Finally, we examine the effects of the parameterization in full 3D simulations in z* and isopycnal coordinates.

2D Parameterization Results
We begin by performing 2D simulations in z* coordinates, identical to the setup described in the Motivation (Section 2.2). The nominal resolution (dx = 125 m and dz = 20.8 m) z* case presented in Figure 1 performed well in resolving the SI and was consistent with the dynamics identified in YL2019. To further ensure numerical convergence and accuracy of these results, we perform an "ultra-high" resolution case with dx = 31.25 m and dz = 5.2 m (4 times more points in x and z) without SI parameterization. This case is shown in the leftmost column of Figure 5. The potential density and tracer distributions are nearly identical to the nominal resolution results, the jets have the same magnitudes and locations, and the SI is resolved along the topography. This case is taken to be the benchmark for comparing with lower-resolution cases. We next examine a "low"-resolution simulation with dx = 250 m and dz = 41.7 m (2 times fewer points in x and z relative to the nominal case) with the parameterization off and on (middle and rightmost columns, respectively). This case takes longer to reach steady state as the offshore flow in the shallow shelf region is poorly resolved, so the simulations are extended to 120 days. In the case without SI parameterization, we see that the potential density and tracer signatures are less pronounced along the topography compared to the high-resolution case. The geostrophic jets are weaker, and the SI is underdeveloped due to the smaller shear and density gradients arising from poor resolution. Significant improvement is apparent in the case with SI parameterization. Due to isopycnal slumping and parameterized frontal mixing, the dense shelf water is able to propagate off the shelf more readily, allowing the dynamics to evolve. The jet structure, density, and tracer concentration are nearly identical to the high-resolution case, though the overflow is slightly diffuse and oriented more parallel to the slope.
To further test the SI parameterization, we examine an extreme "ultra-low"-resolution case where there are 4 times fewer points (relative to nominal) in x and z, so that dx = 500 m and dz = 83.3 m-SI is fully unresolved. Figure 6 shows the resulting fields at 120 days. In the case without SI parameterization, the tracer and density fields show a buildup of dense water on the shelf. Due to the shear mixing scheme relying on vertical gradients, there is no mechanism by which to mix across the front and the dense water is unable to propagate offshore. The geostrophic jets are entirely absent as the relevant dynamics do not evolve. In the case with SI parameterization, there is cross-front mixing that allows the dense water to propagate offshore, undergo geostrophic adjustment, and evolve according to YL2019. Though the SI is not resolved, the jets and density gradients are-the SI scheme leads to diapycnal mixing and maintains evolution of the dynamics. Although there is spurious diapycnal mixing and the final density and tracer values are slightly lower in the overflow compared to the higher-resolution case, the SI parameterization overall performs remarkably well.
We now examine the diffusivity values that the SI parameterization imparts to the flow, considering the nominal, low-, and ultra-low-resolution cases. We consider diffusivity from two sources: the shear mixing parameterization and the newly implemented SI parameterization. When the SI parameterization is turned off, the shear mixing diffusivity is the sole component (aside from spurious mixing). In Figure 7, we show the total (SI plus shear) diffusivity values with the parameterization on and off. Without the parameterization, we see that as the resolution is degraded, the diffusivity value decreases. When the SI scheme is turned on, the total diffusivity is consistent with the case without SI parameterization at the nominal resolution. Interestingly, when resolution is degraded with the SI scheme on, the diffusivity maintains consistent values in the correct locations. Thus, the SI parameterization is able to represent the effects of SI-driven turbulent mixing even at low resolutions where the SI and vertical gradients are poorly resolved and the shear mixing parameterization fails.
In Figure 8, we examine the relative magnitudes of the shear (left column) and SI (right column) diffusivity components for four different resolutions with the SI parameterization turned on. For the ultra-high-resolution case, diffusivity is strongly dominated by the shear component, as SI and partly its secondary shear instabilities are resolved. At the nominal resolution, the SI diffusivity is slightly higher as some of the smallest scale velocity shears captured in the ultra-high-resolution case are no longer resolved. However, the shear diffusivity component is still dominant. As we move to low and ultra-low resolutions, the diffusivity becomes primarily set by the parameterized SI component. SI is no longer resolved and the frontal mixing is accounted for by representing the effects of SI and its secondary shear instabilities. Note that with the SI scheme off in the ultra-low-resolution case there was no parameterized shear diffusivity (Figure 7), while with the SI scheme on there are small values of shear diffusivity along the topography coexisting with the SI parameterization. Even though SI is dominant, the scheme allows the correct dynamics to evolve and leads to some shear regions along the topography that are dissipated by the shear mixing parameterization. In other words, the SI scheme also helps improve the performance of the shear mixing scheme at the lowest YANKOVSKY ET AL.
10.1029/2020MS002264 16 of 28 Figure 5. 2D z* simulation xz-slices of potential density, offshore velocity, alongshore velocity, and passive tracer concentration for the ultra-high-resolution case, low-resolution case with the SI parameterization off, and low-resolution case with the SI parameterization on (left, middle, and right column, respectively). Results are shown at 80 days for the ultra-high-resolution case and 120 days for the low-resolution cases. The black lines indicate where coordinate surfaces are defined; in the ultra-high-resolution case, every eighth vertical level is shown. Figure 6. 2D z* simulation xz-slices of potential density, offshore velocity, alongshore velocity, and passive tracer concentration with the SI parameterization off (left) and on (right) at the ultra-low resolution (dx = 500 m and dz = 83.3 m). Results are shown at 120 days and black lines indicate where coordinate surfaces are defined. resolutions. Overall, the SI diffusivity at ultra-low resolution has a very similar structure to the shear diffusivity at the ultra-high-resolution case.
To further test the performance of the SI scheme, we consider its effects on the isopycnal simulations where SI is unresolved-an extreme limit of the degraded resolution z* simulations as the thickness of the isopycnal layers in the abyssal ocean is on the order of 1,000 m. Figure 9 shows potential density and tracer concentration for cases with the parameterization off and on (left and right columns, respectively). In the case with SI parameterization, at 80 days, there is a slightly higher concentration of dense water present in the abyssal portion of the domain and the shelfbreak contains layers that are more uniformly filled and undergoing mixing. In the case without SI parameterization, there is less mixing across density fronts and the initially infinitesimally thin layers are not uniformly filled. By 120 days, the difference is more pronounced, with the case with SI parameterization having smoother higher-density isopycnal layers present along the bottom. The tracer concentration also shows a smoother structure that is closer to the high-resolution z* cases.
In Figure 10, the two diffusivity components for the isopycnal layer case with SI parameterization and the shear (or total) diffusivity for the case without SI parameterization are shown. At 80 days, the deeper portion of the domain without SI parameterization has diffusivities near zero; with the SI parameterization, both SI and shear diffusivities are elevated. By 120 days, the case with SI parameterization has higher diffusivities along the shelf but lower diffusivities downslope. The source of the diffusivity is primarily the SI component. The main difference in diffusivity values is in the shelf region, where the SI scheme allows for more uniform mixing, particularly within the infinitesimally thin layers. The increased diffusivity within frontal regions allows for a more uniform distribution of dense water among the linearly defined isopycnal layers and thus more realistic propagation of dense water offshore. Without the SI scheme, a few isopycnal layers become filled preferentially on the shelf, nearly vertical isopycnals are established with an unmixed density front in the horizontal direction, and dense water fails to propagate offshore (similar to the ultra-low-resolution z* case without SI parameterization). Additionally, note that when imposing Pr = 1, the viscosity values stemming from the SI scheme will be on the order of 1 × 10 −4 m 2 /s, similar to the background value YANKOVSKY ET AL.    . xz-Slices of potential density (upper three rows) and passive tracer concentration (bottom row) for the isopycnal layer 2D cases; the case with SI parameterization off is shown on the left and the case with SI parameterization on is shown on the right. Upper two rows are at 80 days; lower two rows are at 120 days. prescribed in the model. The scheme's viscosity contribution primarily serves to dissipate momentum to maintain numerical stability and has little influence on the dynamics.
Results for all the 2D cases are summarized in Figure 11. Here, the potential density anomaly (final minus initial density) is weighted by the tracer concentration and plotted as a function of x for all the presented 2D cases. Solid lines indicate simulations where the SI parameterization is off, dashed lines parameterization on, and scaled topography is shown in gray shading. At ultra-high-and nominal resolutions, the parameterization does very little, as the SI is already well resolved. The nominal resolution case exhibits a slightly lower density anomaly due to heightened spurious mixing. The nominal case lies at the precise resolution where the shear mixing and SI scheme criteria are both well satisfied and thus, combined with spurious mixing, slightly overmix relative to the ultra-high case. As we move to lower resolutions, we see a large improvement with the SI scheme on. At the low resolution, the case with SI parameterization has a significantly higher weighted density anomaly than the case without and is much closer to the ultra-high-resolution results (in fact, better than the nominal resolution case with SI parameterization). Likewise for the ultra-low case-the density structure of the overflow is significantly improved with the SI scheme. In isopycnal coordinates, the parameterization also adds a physically based diffusivity component, bringing the tracer-weighted density anomaly closer to the higher-resolution z* results. Overall, the SI scheme helps preserve the relevant dynamics and overflow structure when the submesoscale range is unresolved.
10.1029/2020MS002264 21 of 28 Figure 10. Diffusivities (SI and shear components) for the 2D isopycnal-coordinate cases with the SI parameterization on (top two rows) and off (bottom row) at 80 and 120 days. The first two columns show the region along the slope and the third is a closeup of the shelfbreak.

3D Parameterization Results
Here, we present 3D simulations, identical to the 2D setup but with a 100 km extent in the y-direction. We examine two cases: nominal resolution with dx = dy = 200 m (taken as the benchmark) and coarser resolution with dx = dy = 1 km. In both cases, there are 120 layers in the vertical, defined in z* and in layer coordinates identically as for the 2D simulations (linearly in z or density space). As in YL2019, the dynamics of the 3D case are generally similar to the 2D case, with the formation of bottom-and surface-intensified geostrophic jets. Unlike the 2D case, the jets become baroclinically unstable and eddies lead to stirring and propagation of dense water offshore and down the slope. However, the dominant mechanism driving irreversible mixing and water mass modification is still SI, found along the topography and at eddy edges. Thus, in simulations where eddies and fronts are resolved, there is a need to parameterize turbulent SI-driven mixing. Figure 12 shows the alongshore-averaged potential density, passive tracer, total diffusivity, and SI diffusivity contribution for the 1 km horizontal resolution results with the SI parameterization on/off in z* and isopycnal layer coordinates. As shown in Figure 2, z* coordinates capture shear mixing processes well in 3D cases but are also characterized by spurious diapycnal mixing. Thus, they are not the ideal coordinate choice for representing overflows. The SI parameterization does little to change the structure of the overflow in z* coordinates, as mixing is already well accounted for. The main difference between the cases with/without the YANKOVSKY ET AL.
10.1029/2020MS002264 22 of 28 Figure 11. Tracer-weighted potential density anomaly at steady state for the 2D cases in z* and isopycnal layer coordinates, at various resolutions indicated by the legend (4 times is ultra-high resolution and 0.25 times is ultra-low resolution). Results are shown at 80 days for the nominal and high-resolution cases, and at 120 days for lower resolutions. The shaded region shows topography, scaled so that the surface is at 1.0. Solid lines have the SI parameterization off, and dashed lines on. The lower panel shows a magnified view of the region enclosed by the black outline in the upper panel.
SI parameterization is that the diffusivity is SI dominated rather than shear dominated when the scheme is turned on. There is a slightly improved tracer distribution due to the combination of isopycnal slumping and frontal mixing, particularly in the uppermost 1,000 m where there is more dense shelf water present and it is mixed more uniformly across the topography.
Isopycnal coordinates have significantly smaller parameterized mixing along the topography and in the abyssal portions of the domain (Figure 2). The shear mixing parameterization fails to capture the vertical gradients when the isopycnal layers become too thick while the region near the surface exhibits significant mixing due to the well-resolved vertical gradients. When the SI scheme is turned on, the tracer distribution becomes more evenly mixed across the topography. The isopycnal slumping component of the scheme helps preserve a thicker overflow along the slope (particularly, in the uppermost 1,000 m) by propagating dense water offshore, and the diffusivity component introduces heightened mixing within density fronts and ed-YANKOVSKY ET AL.
10.1029/2020MS002264 23 of 28 Figure 12. One-kilometer resolution results at 60 days for the 3D z* and isopycnal layer coordinate configurations, with the SI parameterization on/off. Alongshore-averaged potential density, passive tracer concentration, total diffusivity, and SI diffusivity (for the cases with SI parameterization) are shown (top to bottom). dies. As in z* cases, the diffusivity in the layered isopycnal case also becomes SI dominated when the scheme is on. Overall, the SI scheme is successful in helping preserve a realistic density structure of the overflow throughout the water column and provides a physically justified and energetically consistent mechanism for water mass modification where the shear mixing parameterization is poorly activated. Figure 13 shows the tracer-weighted density anomaly for the 3D cases with the SI scheme on/off in z* and isopycnal layer coordinates (scaled topography is shown in gray shading). In the nominal resolution z* case, we see that the parameterization leads to a lower tracer-weighted density anomaly consistent with increased parameterized mixing. The lower-resolution z* case has a significantly lower density anomaly due to the spurious diapycnal mixing, and the parameterization does little to change the water mass transformation. The isopycnal cases without the parameterization have a similar tracer-weighted density anomaly signature to the nominal resolution z* results (see Figure 12 for examining the differences between the two coordinate systems). With the parameterization turned on, parameterized mixing increases and the density anomaly decreases. We consider the isopycnal layer simulation at nominal resolution with the parameterization on to be the most realistic case, as it lacks spurious diapycnal mixing and incorporates an SI-driven diffusivity. The SI parameterization adds a physically justified diffusivity throughout the symmetrically unstable region adjacent to the topography and allows for more accurate dynamics by properly representing mixing in frontal regions. Thus, the SI scheme is successful in improving the representation of density structure and mixing experienced by overflows, particularly in isopycnal cases.
10.1029/2020MS002264 24 of 28 Figure 13. Tracer-weighted potential density anomaly at steady state for the 3D z* and isopycnal layer cases with the SI parameterization on/off at the nominal resolution and at 1 km horizontal resolution. The shaded region shows topography, scaled so that the surface is at 1.0. Solid lines have the SI parameterization off and dashed lines on. The lower panel shows a magnified view of the region enclosed by the black outline in the upper panel. The benchmark which we believe captures the overflow most realistically is the nominal resolution "Layer Par. ON" case.

Discussion and Conclusions
We developed a parameterization capturing the effects of submesoscale SI-driven turbulence in various oceanic regions where SI may be an important contributor to water mass modification and energetics. The scheme relies entirely on nondimensional parameters and is simple enough to be implemented implicitly into the GFDL-MOM6. Implicit schemes for turbulent mixing are essential in ocean models with large time steps or isopycnal models where layers may be quite thin (Jackson et al., 2008). The premise of the parameterization is to rely on a balanced Richardson number criterion (Ri B < 1) to identify symmetrically unstable regions, taking into account horizontal density gradients rather than solely vertical shear and density gradients as the traditional Ri does. We then slump the isopycnals toward a stable state, assume the released PE corresponds to the TKE of SI motions which is dissipated by secondary shear instability, and apply a vertical diffusivity and equivalent viscosity (Pr = 1) from the TKE based on the Osborn relation (Osborn, 1980). The scheme applies over the entire water column, taking into account SI arising in bottom boundary layers and intermediate depths.
The energetic transfer we represent in this scheme is that of mesoscale energy being removed through submesoscale SI and transferred to dissipation and diapycnal mixing by shear instability. The question of SI energetics is presently an area of active research. Pure SI is sometimes defined in terms of its fastest-growing mode-along-isopycnal motion fed by GSP leading to geostrophic adjustment (indirectly releasing PE)although linear stability analysis shows existence of SI modes which directly extract PE from the flow. Further, real oceanic fronts are often characterized by hybrids of SI and InI, which are challenging to distinguish and have similar effects of directly draining PE (Grisouard, 2018). Our parameterization generalizes to a variety of plausible energetic transfers, as we consider only the final state by which SI has evolved and returned Ri B to 1 by slumping isopycnals. In addition to InI-SI hybrids, there are other submesoscale turbulent phenomena that such a scheme may be extended to. Processes characterized by forward cascade of mesoscale energy into local diapycnal mixing, such as submesoscale baroclinic eddies (also exhibiting isopycnal slumping combined with energy loss to mixing/dissipation), may be similarly parameterized.
Broadly, our scheme sets up communication between a GM-like isopycnal slumping and irreversible diabatic mixing. The criterion for slumping is chosen specifically for SI but may be easily modified to occur on length/time scales of other submesoscale features (or even mesoscale phenomena). Our scheme has diffusivity κ SI and stream function Ψ SI by default operating on submesoscales; when eddies are parameterized κ GM and Ψ GM dominate in magnitude (as they operate on mesoscales). However, the user may explore setting higher values of κ SI to directly simulate removal of mesoscale energy even in cases where eddies are unresolved. This was not tested in the present study but offers the possibility of exploring a GM-like scheme within coarser-resolution regional models and GCMs without the inaccurate assumption of purely viscous energy dissipation made by GM.
We implemented the parameterization within the GFDL-MOM6 and tested it for the case of a rotating gravity current characteristic of the Arctic continental shelf regions, analogous to YL2019. We first considered the effects of the SI scheme in 2D z* coordinate cases. At high resolutions, SI is well resolved and the shear mixing parameterization successfully captures turbulent mixing. As the resolution is degraded, the shear diffusivity becomes smaller as the vertical gradients and SI cease to be resolved. When the SI scheme is applied, the SI diffusivities become comparable to the shear diffusivities of the high-resolution cases. The SI scheme introduces the correct magnitude of mixing and allows the dynamics to evolve more realistically. The isopycnal slumping and frontal mixing near the surface allow the 2D isopycnal case to have an improved density structure and parameterized mixing along the topography, where the shear mixing scheme previously failed. In 3D coordinates, the SI scheme is also successful. It allows for more uniform mixing in the SI regions along the slope, slumps isopycnals and creates mixing in the near-surface density fronts, and allows dense water to move offshore more easily.
Numerous regions throughout the World Ocean are characterized by frontal dynamics conducive to SI development. The gravity current scenario considered here for the Arctic may be extended to other rotating, buoyancy-driven flows. For instance, the Antarctic shelves similarly experience dense water formation leading to overflows that contribute to Antarctic Bottom Water; properly representing mixing within these flows is crucial for constraining ocean circulation. Other systems susceptible to SI-driven turbulence in-clude outflows and fronts arising in the Antarctic marginal ice zone, the Antarctic Circumpolar Current, bottom boundary layers, freshwater frontal systems (both in polar and lower-latitude regions), and western boundary currents. Our focus on the Arctic is a particularly challenging scenario for the development of this scheme, as submesoscale phenomena occur on the smallest scales near the poles.
The parameterization is aimed at models that resolve mesoscale features but not the submesoscale phenomena they host. Presently, the aim is high-resolution regional models as well as low-latitude regions of GCMs, where mesoscales are captured. We recommend turning on the SI scheme for simulations with grid spacings smaller than the Rossby deformation radius, where mesoscale eddies are at least partially resolved. However, due to its computationally inexpensive and scale-aware properties, the scheme may be turned on at coarser resolutions. The scheme does not significantly affect the run-time of simulations (less than 3% change) and is scale aware so that even if turned on globally it will only influence dynamics in the relevant regions. As GCMs approach higher resolutions, our scheme will become increasingly crucial in the dynamically significant polar regions.
Thus, we have successfully developed and implemented a submesoscale mixing parameterization that captures the effects of geostrophic energy dissipation and secondary shear mixing initiated by submesoscale SI. We found that in cases where mesoscale features such as eddies and fronts are resolved, but the submesoscales are not, the scheme improves representation of frontal processes and accurately captures water mass modification through diapycnal mixing. Future work involves testing the parameterization in flow scenarios other than the rotating gravity current case considered here. The instability criterion may be refined to incorporate relative vorticity effects and scenarios in boundary layers where fq < 0 is an insufficient condition for instability. A more sophisticated approach may also be developed to better isolate convectively unstable regions where the PE released by the isopycnal slumping is of the wrong sign. Presently, we simply assume that if the PE release for a column is of the wrong sign, the parameterization is not applied for that entire column (even if the negative column integral only arises from a convectively unstable surface layer). Additionally, a more complex diffusivity calculation may be applied that does not assume complete local dissipation. We used an existing module within MOM6 which converts TKE to diffusivity. Instead, the TKE may be passed to the Jackson shear mixing parameterization to calculate diffusivities over a nonlocal region where shear instability may be important. Overall, even with the simplifications made, the proposed SI scheme performs remarkably well and provides an important contribution in capturing the effects of submesoscale turbulence that has previously not been considered.