Insights into the Mixing Efﬁciency of Submesoscale Centrifugal–Symmetric Instabilities

: Submesoscale processes provide a pathway for energy to transfer from the balanced circulation to turbulent dissipation. One class of submesoscale phenomena that has been shown to be particularly effective at removing energy from the balanced ﬂ ow is centrifugal – symmetric instabilities (CSIs), which grow via geostrophic shear production. CSIs have been observed to generate signi ﬁ cant mixing in both the surface boundary layer and bottom boundary layer ﬂ ows along bathymetry, where they have been implicated in the mixing and water mass transformation of Antarctic Bottom Water. However, the mixing ef ﬁ ciency (i.e., the fraction of the energy extracted from the ﬂ ow used to irreversibly mix the ﬂ uid) of these instabilities remains uncertain, making estimates of mixing and energy dissipation due to CSI dif ﬁ -cult. In this work we use large-eddy simulations to investigate the mixing ef ﬁ ciency of CSIs in the submesoscale range. We ﬁ nd that centrifugally dominated CSIs (i.e., CSI mostly driven by horizontal shear production) tend to have a higher mixing ef ﬁ ciency than symmetrically dominated ones (i.e., driven by vertical shear production). The mixing ef ﬁ ciency associated with CSIs can therefore alternately be signi ﬁ cantly higher or signi ﬁ cantly lower than the canonical value used by most studies. These results can be understood in light of recent work on strati ﬁ ed turbulence, whereby CSIs control the background state of the ﬂ ow in which smaller-scale secondary overturning instabilities develop, thus actively modifying the characteristics of mixing by Kelvin – Helmholtz instabilities. Our results also suggest that it may be possible to predict the mixing ef ﬁ ciency with more readily measurable parameters (viz., the Richardson and Rossby numbers), which would allow for parameterization of this effect.


Introduction
Submesoscale currents (roughly defined as having horizontal scales between 0.1 and 10 km) are common in oceanic flows, with significant impacts on global ocean dynamics (McWilliams 2016;Lévy et al. 2018;Buckingham et al. 2019;Wenegrat et al. 2018b). In particular, they are understood to be one of the major pathways for energy in the large scales of the ocean to cascade down to the smallest scales of the flow, which is a necessary condition for that energy to be dissipated (McWilliams 2016).
Recent work has highlighted centrifugal-symmetric instabilities (CSIs) as particularly effective at generating this forward cascade (D'Asaro et al. 2011;Gula et al. 2016b). These CSIs are active both at the surface (Taylor and Ferrari 2010;D'Asaro et al. 2011;Thomas et al. 2013;Gula et al. 2016a;Savelyev et al. 2018) and in the bottom boundary layer and topographic wakes (Allen and Newberger 1998;Dewar et al. 2015;Molemaker et al. 2015;Gula et al. 2016b;Wenegrat et al. 2018a;Wenegrat and Thomas 2020) and, as such, they may be important both for the energetics of global circulation, as well as for the mixing of buoyancy and other tracers. As an example, it has been suggested that mixing by CSIs may lead to significant water mass transformation of Antarctic Bottom Water, possibly affecting the closure of the abyssal overturning circulation Spingys et al. 2021). However, despite evidence of their potential impacts, the characteristics and dynamics of mixing by submesoscale CSIs remain uncertain. Thus, we dedicate this study to the investigation of mixing by submesoscale CSIs.
A common measure of a flow's mixing is given by its mixing efficiency g, which measures the fraction of the energy extracted from the flow that was used to mix the fluid's buoyancy. The value of g is bounded between 0 and 1 and is generally assumed to be g ≈ 0.17 for ocean turbulence (Osborn 1980;Moum 1996;Wunsch and Ferrari 2004;Bluteau et al. 2013;De Lavergne et al. 2016;Mashayek et al. 2017b) which, historically, has fit most measured ocean data reasonably well (Gregg et al. 2018). However, recent investigations have hinted at g varying widely due to submesoscale phenomena. Notably, a recent field study conducted in the Orkney Deep (Spingys et al. 2021) inferred significantly higher mixing efficiencies (with an average value of 0.48) in a location that is likely unstable to CSIs (Garabato et al. 2019). Numerical simulations by Jiao and Dewar (2015) likewise indicated values of g . 0.3, with speculations that the value could be larger if the simulation was run for longer. These results seem to contrast with those by Taylor and Ferrari (2010), which found some forms of CSI are associated with small time-and spatially integrated vertical buoyancy production rates, suggesting small rates of irreversible mixing of buoyancy (Peltier and Caulfield 2003;Caulfield 2021).
With these ideas in mind, we investigate the mixing efficiency of submesoscale CSIs (Haine and Marshall 1998) using large-eddy simulations (LES) of finite-width geophysical setups. This configuration aims to reproduce the natural constraints of oceanic flows (due to rotation, geometry, etc.) and to obtain usefully realistic flow evolution and mixing dynamics. This is in contrast to previous numerical studies with similar lines of investigation, which used more highly idealized setups (i.e., using simplified boundary conditions that do not straightforwardly correspond to oceanic conditions; Shih et al. 2005;Salehipour and Peltier 2015;Maffioli et al. 2016;Garanaik and Venayagamoorthy 2019;Howland et al. 2020) or employed assumptions which can potentially affect mixing patterns [e.g., assuming an infinite-width front (Thomas et al. 2013;Taylor and Ferrari 2010) or a twodimensional flow (Jiao and Dewar 2015)].
We show evidence that, in the submesoscale range of the parameter space, CSIs equilibrate via secondary Kelvin-Helmholtz instabilities. This fact allows us to make direct connections with the literature on the mixing efficiency of turbulence in stratified flows, which provides a framework for explaining the range of mixing efficiencies generated by CSIs. In short, CSIs control the flow's mixing efficiency by modulating the background state for the secondary Kelvin-Helmholtz instabilities, which overturn and create the smaller-scale 3D turbulent motions that ultimately dissipate kinetic energy and mix buoyancy. The result of this cascade is that mixing is more efficient for CSIs dominated by centrifugal modes (i.e., mostly driven by horizontal shear production) and less efficient for symmetrically dominated ones (i.e., driven by vertical shear production).
The paper is organized as follows. In section 2 we briefly review the theory of CSIs and mixing efficiency quantification in oceanic flows. Section 3 defines the problem setup and describes the numerical simulations used in this paper. We describe our results in section 4, making connections with past literature on ocean mixing and turbulence in stratified flows. A discussion and concluding remarks are presented in section 5.

Theoretical background
A brief review of CSIs and mixing efficiency follows, and the reader is directed to other works for further details (Haine and Marshall 1998;Bluteau et al. 2013;Gregg et al. 2018;Caulfield 2021).

a. CSI theory
Centrifugal-symmetric instabilities (sometimes referred to simply as symmetric instabilities) are defined here as those that emerge when qf , 0 (Haine and Marshall 1998). Here f is the Coriolis frequency and q is the Ertel potential vorticity (PV): where = is the gradient operator, b is the buoyancy, u is the velocity vector, andk is the unit vector in the vertical (z) direction. When the flow is in thermal wind balance Eq. (1) can be rewritten asq 2 is the balanced Richardson number, and u b is the (horizontal) velocity component in thermal wind balance. (The subscript b is used to indicate an assumption of thermal wind balance throughout this paper.) This essentially reduces the instability criterion toq b , 0. Given a dynamical definition of the submesoscale range as Ro b ∼ Ri b ≈ 1 (McWilliams 2016), it can be seen from Eq. (2) that submesoscale flows are particularly likely to be unstable to CSIs. It is useful to characterize CSIs based on their primary source of kinetic energy (Thomas et al. 2013). For the purposes of this paper we focus on the horizontal shear production rates ( SP h , associated with the centrifugal modes; · denotes any volume averaging procedure) and vertical shear production rates ( SP y , associated with symmetric modes). Hence, a straightforward way to characterize CSIs is by estimating their ratio, which (assuming a background flow with uniform gradients, and that CSI-unstable parcels move in paths whose angle with the horizontal direction is small) can be approximated as [Thomas et al. 2013, their Eq. (42)] where R SP is the ratio of horizontal to vertical shear production rates. Thus the larger R SP , the more centrifugally dominated a CSI (and opposite for symmetrically dominated CSIs). In all cases in this paper the second term in parenthesis is small and we approximate Eq. (3) by which can be understood as the ratio of the two nonunitary terms in Eq. (2). Hence, wheneverq b is negative primarily due to Ro b being sufficiently negative (i.e., due to the horizontal shear), we call the ensuing instability a centrifugally dominated CSI (R SP . 1). Similarly, whenq b , 0 due to Ri b being small (i.e., due to the vertical shear), we say the ensuing instability is symmetrically dominated (R SP , 1). We note that, while CSIs are generally understood to grow using the kinetic energy of the balanced flow through the shear production rate terms (Haine andMarshall 1998), Wienkers et al. (2021) showed that in the limiting case of Ro b = 0 and for fronts with relatively shallow isopycnal slopes, symmetrically dominated CSIs can grow primarily at the expense of the potential energy of the balanced flow, making the vertical buoyancy flux term dominant. Although these limiting cases may be relevant for some broad fronts, we do not focus on this part of the parameter space in this study.
In the general case (a CSI where both or either Ro b and Ri b contribute toq b being negative), and again assuming a balanced background flow, the linear inviscid growth rate v of the instability is (Haine and Marshall 1998) v which reveals that, at a fixed latitude, a given instance of CSI will grow faster the more negativeq b is.

b. Mixing efficiency theory
We focus on the mixing efficiency g, which we define as (Peltier and Caulfield 2003;Venayagamoorthy and Stretch 2010;Salehipour and Peltier 2015) where « k = 2n S ij S ij is kinetic energy dissipation rate [S ij (u i /x j 1 u j /x i )/2 is the strain rate tensor, and n is the total flow viscosity, which may include an eddy viscosity] and « p is the irreversible mixing rate of buoyancy (see appendix A for details about its calculation). Note that there are many definitions of g in the literature [see Gregg et al. (2018) for a review], but we choose Eq. (6) because it specifically considers only irreversible processes, making it more accurate. We also consider the cumulative mixing efficiency (Gregg et al. 2018;Peltier and Caulfield 2003;Mashayek and Peltier 2012a,b;Caulfield 2021): Since both « k and « p eventually go to zero after a sufficiently long time, C(t) approaches an asymptotic value as t → '. This makes C a better approach to quantify the cumulative mixing of a given instability over its lifetime. Dimensional analysis indicates that the mixing efficiency (either g or C) of a given flow depends on several parameters, including the flow's Froude number, Richardson number, turbulence and molecular Prandtl numbers, and buoyancy Reynolds number (Salehipour and Peltier 2015;Maffioli et al. 2016;Khani 2018;Caulfield 2021), although it remains unclear which ones are the most important or what is the functional shape of these dependencies (Caulfield 2021). Here we focus on the buoyancy Reynolds number (Shih et al. 2005), which, based on recent literature, seems to organize results from idealized numerical simulations reasonably well (Shih et al. 2005;Bouffard and Boegman 2013;Salehipour and Peltier 2015). It was proposed in part because it is easier to estimate in field campaigns than a more traditional Reynolds number, serving as a proxy for the intensity of turbulence. It can be written as where N 2 0 is a constant background stratification and n mol is the molecular viscosity of the fluid.

Problem setup
We use a numerical setup that approximates geophysical flows while allowing the Rossby and Richardson numbers of the flow to be easily varied. In this section we describe that setup in detail, including the numerical tools used for the simulations.

a. Initial conditions
We start our simulations with a thermal-wind-balanced front configuration given by the following equations: where u, y, and w are the eastward (x), northward (y), and upward (z) velocity directions; u 0 is a velocity constant; N 0 denotes a constant Brunt-Väisälä frequency; f 0 is the Coriolis frequency; f y , F y , and f z are nondimensional functions of y and z given by where s y and s z are decay length scales; and y 0 serves as a length offset. For the purposes of our paper s z = 80 m, and y 0 is always set to be half the length of our domain in the y direction (4 km; see section 3b). The equations above define a Gaussian-shaped front centered at y 0 with a vertically constant vertical shear of u 0 f y /s z a width s y , and a superimposed spatially uniform background stratification N 2 0 . A vertical cross section of the front showingq b can be seen in Figs. 1a and 1b for two different sets of parameters (details are given in section 3b). Recall that CSIs emerge in the regions whereq b (shown in the color map) is negative.
If we exclude molecular viscosity and diffusivity (on the basis of a high-Reynolds-number assumption), our setup can be fully defined with the parameters u 0 , s y , s z , N 2 0 , f, the eddy viscosity n e , and the eddy diffusivity of buoyancy k e . Application of dimensional analysis produces five nondimensional parameters: in addition to Rossby, Richardson, and Reynolds numbers.
Here d is the aspect ratio and Pr t is the turbulent Prandtl number. For simplicity we only consider the case Pr t = 1 (which in our case is a turbulent Prandtl number since we use C H O R E T A L . eddy diffusivity closures) and, given the uncertainty of d in real oceanic conditions, we assume that the aspect ratio does not affect results as strongly as the Rossby or Richardson numbers. We thus report d, but do not make efforts to explore its range. To use representative values to characterize our simulations, we use Ro r and Ri r , which we refer to as reference Rossby and Richardson numbers, to characterize the parameter space. They are defined as the Rossby and Richardson numbers at the point of the domain whereq b is initially (i.e., at t = 0) the lowest. Recall that this corresponds to the point with the fastest linear growth rate for CSIs according to Eq. (5), making Ro r and Ri r relevant quantities of the flow evolution. For our setup, this point always lies at z = 0 but the y location is found numerically given the challenge of obtaining a closed-form expression for it from Eqs. (9)-(14). The reference point is shown as white circles in Figs. 1a and 1b. A parameter space of Ro r 2 1/Ri r is shown in Fig. 2, where the color map shows values ofq b at the reference point.
Finally, following previous literature (Shih et al. 2005;Salehipour and Peltier 2015), we use the buoyancy Reynolds number [properly defined for LES cases in Eq. (19)] to diagnose the turbulence intensity related to the stabilizing effect of stratification. We focus our exploration of parameter space on the Rossby and Richardson numbers and we use the buoyancy Reynolds number as a diagnostic quantity.

b. Simulations
We use the Julia package Oceananigans (Ramadhan et al. 2020) to run a series of numerical simulations with Eqs. (9)-(14) as initial conditions. Oceananigans uses a finite volume discretization based on that of MITgcm (Marshall et al. 1997), and we run it with a fifth-order weighted essentially non-oscillatory advection scheme and a third-order Runge-Kutta time-stepping method. The bulk of our simulations are three-dimensional (3D) LES (whose parameters can be found in Table 1) where we solve the filtered nonhydrostatic incompressible Boussinesq equations where k is the unit vector in the vertical (z) direction, r 0 is a reference density, p is the pressure, n mol = 10 26 m 2 s 21 is the molecular viscosity, t is the subgrid-scale (SGS) stress tensor, k mol = 1.5 3 10 27 m 2 s 21 is the molecular diffusivity, and k is the SGS buoyancy flux. We also run three auxiliary two-dimensional (2D) simulations that solve Eqs. (17) and (18), but dropping the SGS closure and molecular diffusivities in favor of a constant eddy diffusivity (the purposes of these simulations is made clear in section 4 and their parameters can be found in Table 2). The two-dimensional domains retain all three velocity components despite only formally including the y and z directions, in what is sometimes called 2.5D setup (Kämpf 2010).
All simulations are bounded in the y and z directions, and the 3D simulations are periodic in the x (alongfront) direction. In all cases a buoyancy gradient of db/dz N 2 0 was imposed at the top and bottom boundaries (in order to minimize initial dissipation of buoyancy before the onset of turbulence) and all other nonperiodic boundary conditions imposed zero fluxes for the momentum components and the buoyancy scalar. No-flux boundary conditions at the top and bottom of the domain were also tested and found to not affect our findings. A constant background rotation rate f was imposed on the domain for each simulation and sponge layers were included on both ends of the y direction with a width of 1/16 of the domain length each to absorb internal waves and simulate open boundaries.
For the 2D setups, a constant isotropic eddy diffusivity was used with its value set to be as low as possible while still producing well-resolved simulations. Resolvedness was verified both by inspecting the small scales of the flow visually and by ensuring that an effective Kolmogorov scale [(n 3 e /« k ) 1/4 , estimating the smallest scales of motion in our 2D simulations] was always at least ≈30% larger than the grid spacing (further refinement produced no significant change in the results). In the 3D simulations we used a constant-coefficient Smagorinsky model closure (Smagorinsky 1963) with a modification that reduces the eddy viscosity in stably stratified regions (Lilly 1962). In appendix B we follow Khani (2018) and analyze results from extra simulations at different resolutions and show that our results are well converged at the resolution at which we present them.
The simulation parameters for the main (3D LES) runs are given in Table 1 and their location in the Ro r 2 1/Ri r parameter space can be seen in Fig. 2, where each symbol corresponds to a simulation. Relevant simulation parameters for the auxiliary 2D runs are given in Table 2. Their values for Ro r and Ri r are exactly the same as those of their 3D counterparts, so they do not expand the exploration of the parameter space. Finally, the code used to generate the results used in this paper is publicly available in Chor et al. (2022).

a. Time evolution of the mixing efficiencies
All our simulations go through qualitatively similar evolutions: 2D primary instabilities (CSIs) develop quickly in the initially unstable (q b , 0) region, followed by the sudden onset of the secondary instabilities creating 3D turbulence and releasing internal waves, followed by a longer decay of the turbulence. We focus for now on simulations CIfront1 and FIG. 2. Simulations (circles and triangles) on top of the Ro r 2 1/Ri r parameter space, where Ro r and Ri r are the reference Rossby and Richardson numbers, respectively, which are defined shortly after Eq. (16). The dark gray areas denote regions of negative Ri r (thus impossible to achieve in a stably stratified environment such as ours) and light gray areas denote regions stable to CSIs (q b . 0). The dot-dashed line corresponds to Ro r Ri r = 21, which theoretically separates centrifugally dominated (to the left of the line) from symmetrically dominated CSIs (to the right of the line) for a thermal-wind-balanced environment. Note that simulation SIfront5 does not appear in the plot but is located at the same point as simulation SIfront4.  SIfront4 (representative of centrifugally and symmetrically dominated CSIs, respectively) to illustrate that process in this section, and encourage readers to refer to the animations that are available in the supplemental material to gain more intuition. Results for these simulations are shown in Figs. 3a and 3b, which show the kinetic energy dissipation rate and the mixing rate of buoyancy. After a quiescent start (indicated by low values of « k and « p ; refer to appendix A for details about the calculation of « p ), primary CSIs (which are mostly 2D in the y-z plane) develop within one inertial period. Between one and three inertial periods the shear from the primary instabilities becomes sufficiently strong to generate secondary instabilities (see section 4b) that mediate the transition to full 3D turbulence; this roughly coincides with the first peak in « k . The ensuing turbulent flow can be seen in Figs. 1c-f. Note that simulation SIfront4 reaches the onset of turbulence earlier than simulation CIfront1 because it has lower initial values ofq b (see Fig. 2), which translates into a faster growth rate for the CSIs per Eq. (5) (Haine and Marshall 1998).
Internal waves are generated in all our simulations during the emergence of the secondary instabilities (which is explosive in nature). However, the total amount of energy radiated via internal waves (as quantified by the energy dissipated in the sponge layers) is never larger than around 1/1000 of the kinetic energy dissipated by the instabilities, which qualitatively matches the findings of Kloosterziel et al. (2007) for centrifugal instabilities. Interestingly, more waves are visible on the lighter side of the front compared to the heavier side. This can be seen in Figs. 1c and 1d (the portion of the domain shown does not include the sponge layers).
In Figs. 3c and 3d we show two measures of mixing efficiency: the instantaneous mixing efficiency g (dashed lines) and the cumulative mixing efficiency C (solid lines). The pattern of the instantaneous measure is significantly noisier than the cumulative one, with abrupt changes in g over short times (this is especially true for simulation CIfront1). This variability suggests caution in extrapolating instantaneous mixing efficiencies from observations as a means of characterizing the integrated mixing of a given flow throughout its lifespan. For the purposes of our analysis, we overcome this limitation by using the cumulative mixing efficiency C [Eq. (7)]. As t → ', C(t) converges to a value C ' , which we take to be representative of the total mixing of the flow. In practice a good approximation for C ' can be obtained by taking C at around 12 inertial periods (after its value has approximately converged in all our simulations), which we adopt as our approach.
We plot results for C ' in Fig. 4 as a function of the ratio between the average horizontal and vertical shear production rates (results for C ' are the larger, bolder symbols).   . 1) to have higher mixing efficiencies than symmetrically dominated ones. We also plot the maximum value of the instantaneous mixing efficiency g max for each simulation as smaller symbols with slight transparency. The same pattern is evident, with the mixing efficiency increasing as the modes become more centrifugally dominated. Figure 4 shows a clear pattern in which values of C ' for symmetrically dominated CSIs (R prim SP , 1) are lower than the canonical value of 0.17 (shown as a dashed gray line for reference), while values for centrifugally dominated CSIs (R prim SP . 1) are higher. Additionally, values of g max for centrifugally dominated CSIs can reach even higher values. This large range of values in Fig. 4 is in qualitative agreement with previous indications that mixing efficiencies of submesoscale CSIs can significantly deviate from the commonly used value (Taylor and Ferrari 2010;Spingys et al. 2021).
We note that the mixing efficiencies found in these simulations are somewhat more moderate than those reported from 2D simulations using constant eddy viscosities, where it has been argued that centrifugal instabilities may generate g ≈ 1 (Jiao and Dewar 2015). We are able to reproduce similar results for our basic frontal configuration when using a similar 2D constant-viscosity setup (i.e., low Reynolds number direct numerical simulation, matching the simulations used in the study), but not when using LES closures and in 3D. Furthermore, even in 2D constant-viscosity simulations, using C as a metric (instead of g) also indicates more moderate mixing efficiencies since the largest values of g happen after most of the turbulence has ceased (see, e.g., Jiao and Dewar 2015, their Fig. 14) and thus contributes little to the total mixing performed by the flow. These numerical results (and the connections to Kelvin-Helmholtz mixing discussed below) thus are taken to indicate that, while centrifugal instabilities can generate significantly enhanced mixing efficiencies, some of the prior results indicating CSIs generating near perfectly efficient mixing may have been reflective of numerical methods (mainly associated with the fact that 2D simulations cannot properly represent the small-scale 3D eddies that ultimately dissipate energy and mix buoyancy), and not entirely representative of high-Reynolds number oceanic flows.

b. The nature of the secondary instabilities
The variations in mixing efficiency across CSI simulations indicate changes in the mixing generated by secondary instabilities during the equilibration process. In this section we therefore identify the secondary instabilities that mediate the transition from CSI modes to full 3D turbulence. Previous work by Taylor and Ferrari (2009) has shown that pure symmetric instabilities (CSIs in the absence of any centrifugal modes) equilibrate via Kelvin-Helmholtz instabilities (KHIs). Griffiths (2003) likewise inferred KHIs as the equilibration mechanism for centrifugal instability, although this was based on simulations that did not directly resolve overturning motions.
A significant difference between centrifugal and symmetric instabilities}which might be hypothesized as the source of the enhanced mixing efficiencies}is that the fastest growing linear centrifugal modes cross isopycnals, whereas the symmetric modes do not (Thomas et al. 2013). This suggests the possibility that buoyancy advection by centrifugal instabilities adds a gravitational instability component to the equilibration process, which is known to have higher mixing efficiencies than KHIs (Gayen et al. 2013;Wykes and Dalziel 2014). Given the uncertainty in the equilibration mechanism of centrifugal instabilities, we therefore focus in this section on the onset of the secondary instabilities, where we find that they indeed are of the Kelvin-Helmholtz type, as shown below.
Thus, for the purposes of this section, we run three 2D simulations with a constant eddy diffusivity (2D_CIfront1, 2D_CIfront3, and 2D_SIfront4; see Table 2) that are otherwise identical to simulations CIfront1, CIfront3, and SIfront4. The use of a constant eddy diffusivity avoids possible artificial changes in the energetics and dynamics due to the SGS model [e.g., see Piomelli et al. (1990), which showed that SGS models can underpredict the growth rate of perturbations and delay the transition to turbulence], and the two-dimensionality is designed to save computational resources since we anticipate both the primary (CSI) and secondary (KHI) instabilities to be 2D in the y-z plane (Peltier and Caulfield 2003;Rahmani et al. 2014). Note that, while KHIs have been recently shown to produce 3D instabilities after the formation of the billows but before their full overturning (Mashayek and Peltier 2012a,b;Thorpe 2012), the initial billows are always 2D as shown by Squire's theorem (Squire 1933) and we stop our instability analysis before the 3D effects become important.
We focus the analysis on a small portion of the domain (to avoid interference by the edges of CSI modes and other features of the flow) and quantify the horizontal and vertical shear production rates separately, as well as the buoyancy production rate and the Richardson number. The subdomains used, however, were verified to be representative of the turbulence transition of the CSI modes as a whole.
Results are shown in Figs. 5a-f for the centrifugally dominated Simulation 2D_CIfront1. Figures 5a-e show the evolution of Ri in snapshots as time progresses, with Ri values between 0 and 1/4 shaded gray in order to indicate areas that are susceptible to KHIs (Miles 1961;Howard 1961). It is clear in the first panels that a large horizontal portion of the subdomain is susceptible to KHIs as indicated by the gray-shaded areas. The light white-bluish areas indicate that a portion of the domain also has slightly negative Ri (their magnitudes are mostly smaller than 0.05), which is a consequence of the primary centrifugal modes crossing isopycnals, as expected. In Fig. 5d, we see undulations qualitatively characteristic of KHIs before a decay into turbulence in Fig. 5e. Figure 5f shows subdomain means (denoted by · s ; the subdomain being the rectangular domain portion shown in Figs. 5a-e) of the vertical buoyancy flux and shear production rate components for Simulation 2D_CIfront1 (details about this calculation can be found in appendix C, section b). At early times, all the averages are approximately zero, but the secondary instability growth is dominated by vertical shear production. Note that the buoyancy production rate is actually negative (implying energy moving from kinetic form to potential form) despite portions of the subdomain having slightly unstable stratification. Thus, despite buoyancy advection generating regions potentially susceptible to gravitational instability, the primary energy source for the secondary instabilities remains vertical shear production. The same analysis applied to Simulation 2D_CIfront3 produces similar results (Fig. 5g).
These characteristics are expected of KHIs (Peltier and Caulfield 2003), which strongly suggests that centrifugally dominated CSIs equilibrate through secondary KHIs, as argued by Griffiths (2003). The same analysis for symmetrically dominated CSIs (using Simulation 2D_SIfront4 and shown in Fig. 5g) produces very similar results (albeit with smaller regions of unstable stratification due to the alignment of the symmetric modes with isopycnals) and identical conclusions}consistent with earlier analysis by Taylor and Ferrari (2009). To make sure that this feature is not specific to our frontal setup, we ran the same analysis for a centrifugally unstable interior jet similar to the one considered by Jiao and Dewar (2015), but with a higher resolution than used in that paper. The results (not shown) are again extremely similar to the ones just described, indicating that CSIs in the submesoscale portion of the parameter space (regardless of being symmetrically or centrifugally dominated) equilibrate via KHIs that emerge from the vertical shear of the primary modes before the onset of gravitational instability.

c. The role of the secondary instabilities in the mixing efficiency
Given that the transition to 3D turbulence is mediated by KHIs, it is now possible to connect these geophysical flows with results on turbulence in stratified flows. We note that, despite the connection made here between submesoscale CSIs and stratified turbulence results, the main purpose of this work is not to provide a detailed analysis of the parameter space dependence of small-scale turbulence in stratified environments. Rather, our main goal is to use established relations from previous literature to uncover mechanisms by which submesoscale flows mix buoyancy and dissipate kinetic energy.
Many recent investigations focus on the buoyancy Reynolds number Re b (Shih et al. 2005;Bouffard and Boegman 2013;Salehipour and Peltier 2015;Mashayek et al. 2017a) to explain the mixing efficiencies of overturning motions in stratified environments, which we found to be a good predictor of g in our simulations. In experimental settings and in direct numerical simulations, Re b is well defined as presented in Eq. (8), but it needs to be adapted for use with our LES, where the eddy viscosity varies in time and space. We thus define the subgridscale buoyancy Reynolds number as where · q denotes an average over the region whereq b , 0 at t = 0. 1 In principle, any consistent averaging procedure can be used to define Re sgs b , but we choose · q due to the changing distribution of unstable areas in our setup across the parameter space (see, e.g., Figs. 1a and 1b). We chose to use N 2 0 here instead of db/dz q since we want to characterize the background stratification against which the overturning motions need to do work without the influence of the locally unstable stratification generated by the overturning motions themselves. Quantitative aspects of the results change slightly when evaluating Eq. (19) using db/dz q , but the qualitative characteristics remain unchanged.
Note that, apart from the different averaging operators in the numerator, Eq. (19) is mathematically equivalent to Eq. (8) and it produces the same numerical results if all scales of the flow are well sampled and represented down to the Kolmogorov scale (e.g., in direct numerical simulations). However, for the case of LES, the strain rate tensor S ij cannot capture the smallest scales (whose effects are instead parameterized via the subgrid-scale viscosity n e ), and thus Eq. (19) produces smaller numerical values than Eq. (8). As a result, while low values of Re b indicate a low-turbulence regime where molecular viscosity effects dominate (Brethouwer et al. 2007), that need not be the case for low values of Re sgs b , where these regimes can still be turbulent. For this reason the magnitude of Re sgs b is expected to differ from Re b (as a function of the resolved scales of the LES), making direct comparisons of magnitude between Re sgs b and Re b from DNS not possible. We thus opt to maintain separate nomenclature for Re sgs b and Re b , using Re b only when the molecular viscosity is explicitly used. We explore an estimation of Re b for our simulations in appendix D.
Comparisons between these metrics are further confounded as previous studies investigating Re b rely on very well-controlled numerical simulations, where turbulent regions are more easily predicted and averaging procedures can be performed straightforwardly. This is not the case for our simulations, which are significantly more realistic, contributing to a patchy turbulence pattern (e.g., Figs. 1e and 1f). This pattern is more representative of real ocean turbulence, but again excludes quantitative comparisons of buoyancy Reynolds number magnitude with more idealized studies and between different configurations (Mashayek et al. 2017b;Howland et al. 2020;Caulfield 2021)}see also discussion in Mashayek et al. (2021, section 7). Given these facts and our focus on making a connection between submesoscale flows and turbulence in stratified environments (as opposed to investigating the dynamics of stratified turbulence per se), we concentrate analyses on the relative variations of g with Re sgs b , without focusing on absolute values of Re sgs b . Similar to previous studies of KHI (Shih et al. 2005;Salehipour and Peltier 2015), we compare instantaneous mixing efficiencies 1 Note that · q is different from the previously introduced · and · s , which denote an average over the whole domain and an average over the rectangular subdomain shown in Figs. 5a-e, respectively. g with instantaneous values of Re sgs b in Fig. 6a. To ensure that only cases with significant 3D turbulence are taken into account, we only consider times after the peak in the dissipation rate « k (indicating a transition to full 3D turbulence) and discard points where « k q is smaller than 10 210 m 2 s 23 . Figure 6a shows a pattern where, for small values of Re sgs b , g does not depend on Re sgs b , followed by a power-law dependence for larger values which follows an approximate 21/2 slope (as evidenced by comparing it with the solid line). In our simulations, these characteristics persist regardless of the averaging procedure used, as long as most of the averaging volume is turbulent. Furthermore, both the region of approximately constant g for Re sgs b ≈ 1 and the region of power-law dependence match well with previous findings for KHI (Shih et al. 2005;Salehipour and Peltier 2015) despite the significant difference between setups.
The agreement between our data and simulations of idealized KHIs is further evidence that the mixing efficiencies of these submesoscale instabilities are ultimately controlled by the small-scale overturning motions of the flows that emerge as a consequence of CSIs. This suggests that CSIs control the mixing efficiency by adjusting the background state from which KHIs emerge, namely, the stratification, vertical shear, and dissipation rate, which directly modulate the Richardson and buoyancy Reynolds numbers. Along with the Prandtl number, this sets all three nondimensional parameters necessary to characterize overturning motions in stratified flow (Mashayek et al. 2017a, section 2.2). Although a Froude number is also necessary if one includes the kinetic energy as a relevant quantity (Caulfield 2021, section 2.4), we found no clear Froude number dependence for the mixing efficiency in our results. This is in contrast to some other results which found that it is the preferred parameter for organizing values of g (Maffioli et al. 2016;Garanaik and Venayagamoorthy 2019). This discrepancy may be due to specific dynamics of the frontal configuration used here or simply due to our simulations spanning too small a range of the parameter space for any dependence to be evident. Further investigation of this topic is outside the scope of the present paper and we leave it for future work.
We further find an inverse relation between Re sgs b and 2 Ro q Ri q in our simulations, shown in Fig. 6b ( This result can explain the pattern of mixing efficiencies seen in Fig. 4, where symmetrically dominated CSIs (where 2 Ro q Ri q , 1) tend to have higher buoyancy Reynolds number than centrifugally dominated CSIs (where 2 Ro q Ri q . 1). This relation can be used to plot g as a function of Ro q Ri q in Fig. 6c, where we also see that points collapse rather well. This comparison is similar to that in Fig. 4, and we see that the result again indicates that centrifugally dominated CSIs tend toward higher values of g, and the opposite for symmetrically dominated CSIs.
It is worth mentioning that the range of values for N 2 0 is significantly larger than the range of values of other parameters in our simulations (see Table 1). As such, the collapse of points in Figs. 4 and 6b could be partially explained by changes in N 2 0 . While our frontal configuration is such that large variations in N 2 0 are needed to cover the submesoscale range of the Ro r 2 1/Ri r parameter space (without relying on unrealistic values of other parameters), this is not necessarily always the case in the ocean. Even though preliminary investigations with an interior jet geometry produced a similar result as that shown in Fig. 4 (not shown)}indicating that the relation between the mixing efficiency and Ri r Ro r is not a function of the frontal configuration considered here}a more thorough investigation is necessary to better assess its generality in other flow configurations. We leave this investigation for future studies.

Discussion and conclusions
Centrifugal-symmetric instabilities (CSIs) have been suggested to play a potentially important role in generating turbulent mixing in the ocean surface boundary layer and along topography at the bottom. Thus, we have used LES to investigate several geophysical flows that are unstable to submesoscale CSIs, with the goals of investigating the mechanisms by which the submesoscales produce buoyancy mixing and dissipation, and systematically examining the resulting mixing efficiencies. All simulations in this paper follow a similar evolution: primary CSIs quickly develop in the domain, increase the vertical shear, which prompts the emergence of secondary instabilities (which we showed to be Kelvin-Helmholtz instabilities, KHIs) that mediate the transition to small-scale turbulence, which ultimately dissipates kinetic energy and mixes buoyancy.
We showed that CSIs can generate a wide range of mixing efficiencies (0.05 # C ' # 0.3), which can depart significantly from the community-standard value of 0.17 (Gregg et al. 2018;Caulfield 2021), suggesting caution in the use of a single mixing efficiency value for parameterizations where submesoscale turbulence is active. This variation in mixing efficiency is a consequence of the submesoscales, with centrifugally dominated CSIs tending to have higher instantaneous and cumulative mixing efficiencies than symmetrically dominated instabilities (see Fig. 4). This pattern of mixing efficiencies due to CSIs can be well reproduced using only the Richardson and Rossby numbers (Ri r Ro r ; Fig. 4), suggesting a potential strategy for improving parameterized estimates of mixing due to submesoscale instabilities.
In all simulations considered here KHIs mediate the transition to turbulence, allowing us to explain the observed patterns in mixing efficiency by leveraging results from the stratified turbulence literature. Specifically, we show that variations in mixing efficiency can be understood as the result of CSIs setting the background state on which KHIs grow. CSIs modulate the strength of vertical shear, stratification, and turbulence intensity which have been shown to influence the mixing efficiency of KHIs through the Richardson and buoyancy Reynolds numbers (along with the Prandtl number; Mashayek et al. 2017a;Caulfield 2021). Notably, we were able to reproduce the dependency of the instantaneous mixing efficiency g on the buoyancy Reynolds number Re b using the surrogate parameter Re sgs b , shown in Fig. 6a (note that absolute values of Re sgs b and Re b are not directly comparable given our use of a LES; see appendix D ). The satisfactory collapse of points reproducing a result that is well known in the stratified turbulence literature is evidence that these small overturning instabilities are what ultimately sets the mixing efficiency, providing a direct connection between submesoscale dynamical processes and stratified turbulence. We believe this to be one of the primary contributions of this paper, since it is likely that this control mechanism for CSIs extends beyond the portion of the parameter space explored here, providing the community with additional tools to analyze observations and develop parameterizations.
These results provide a potential explanation of recent observational findings of elevated mixing efficiencies in conditions susceptible to centrifugally dominated CSI in the Orkney Deep (Garabato et al. 2019;Spingys et al. 2021), as well as low mixing efficiencies in simulations of forced-symmetric instability in the surface boundary layer where the stratification remains small (Taylor and Ferrari 2010). We note, however, that, despite this qualitative agreement, we do not find mixing efficiencies as large as implied by some previous work on CSIs (Spingys et al. 2021). This may be a result of the uncertainty in observational estimates, or that the mixing efficiency of CSIs can vary over an even wider range as a consequence of other parameters or flow geometries not varied here. For example, in weak fronts with Ro b ≈ 0 and shallow isopycnal slopes, CSIs can grow by extracting potential energy from the balanced flow (Wienkers et al. 2021), potentially introducing a gravitational component to the energetics that may contribute to higher rates of buoyancy mixing. However, to the extent that our finding of KHIs mediating the transition to turbulence is general for CSIs, we expect our results to be robust, as they depend on the local background state felt by the growing KHI modes, and not directly on the geometry or parameters at the submesoscale.
Finally, evidence that CSIs are common in both the surface and bottom boundary layer suggests the variations in mixing efficiency shown here may be an important aspect of largerscale ocean dynamics and circulation (Allen and Newberger 1998;Taylor and Ferrari 2010;D'Asaro et al. 2011;Thomas et al. 2013;Gula et al. 2016a;Savelyev et al. 2018;Dewar et al. 2015;Molemaker et al. 2015;Gula et al. 2016b;Wenegrat et al. 2018a;Wenegrat and Thomas 2020). The case of abyssal flows offers a particularly compelling example, as observations suggest the possibility of CSIs generated by flow along bottom topography (Ruan et al. 2017;Spingys et al. 2021). Centrifugally dominated instabilities}generated preferentially in regions of steep slopes and strong stratification (Wenegrat and Thomas 2020)}in particular provide a route for the efficient mixing of buoyancy, and hence may contribute to abyssal water mass transformation, a key component of the global overturning circulation. Quantification of the integrated effect of CSIs in both the surface and bottom boundary layer, and the variations of mixing efficiency documented here, remains an open target for future study.
Data availability statement. The numerical model simulations upon which this study is based are too large to archive or to transfer. Instead, all the code used to generate the results is available at https://doi.org/10.5281/zenodo.6509962 (Chor et al. 2022).

Calculation of the Rate of Irreversible Mixing of Buoyancy
We calculate the irreversible mixing of buoyancy based on the theory of Winters et al. (1995). An evolution equation for background potential energy for a control volume can be written as (Winters et al. 1995) where E b is the average background potential energy (the portion of the potential energy unavailable for conversion into kinetic form) per unit mass, S adv and S diff are the advective and diffusive fluxes of E b across the volume's boundaries. The sponge layers used in our simulations do not directly modify the buoyancy, so they do not appear in Eq. (A1). The term « p is the average irreversible mixing of buoyancy (due to diapycnal mixing within the control volume), and it appears as a nonnegative rate of change in Eq. (A1) because potential energy lost due to internal mixing is irreversibly stored as background potential energy (Winters et al. 1995;Winters and D'Asaro 1996). The term S adv is identically zero for our simulations due to the boundary conditions. The term S diff , on the other hand, is nonzero for our domain but its effect on E b was found to be negligibly small. Thus we assume S diff ≈ 0. This allows us to simplify Eq. (A1), leading to our equation for « p which is similar to Eq. (18) of Winters et al. (1995). Thus, in order to apply Eq. (A2), we estimate E b by adiabatically sorting the buoyancy field b at every time step to arrive at a reference state that minimizes horizontal buoyancy gradients (Winters et al. 1995). Although this approximation is the main source of error in our calculation of « p , we found that the error is small enough to be neglected.

Grid Resolution Analysis
To verify that our results were not dependent on grid resolution, we ran extra simulations of setups CIfront1 and SIfront4 (representative of centrifugally and symmetrically dominated simulations, respectively) with varying grid resolutions. These simulations were exactly the same as the ones whose results are presented in main text, except for the spacings Dx, Dy, and Dz, which were increased by factors of 2, 4, and 8 (while keeping the ratios Dx/Dz = Dy/ Dz constant).
Khani (2018) compared LES of idealized stratified turbulent flows against 3D direct numerical simulations (resolving all scales of the flow without the need of a turbulence closure) of the same flows and found that LES produced the correct result when their grid spacing was approximately equal or smaller than the Ozmidov length scale. Thus, we compare our grid spacing with the Ozmidov length, defined as (Khani 2018) where « k is average of « k over turbulent regions of the flow (defined as regions where « k . 10 210 m 2 s 23 , although the precise value of the threshold was verified to not significantly impact results), and the factor 2p follows the definition in Khani (2018) and is included since it is often the corresponding wavenumber which appears in applications (Khani and Waite 2014). We plot the quantity Dz/L O as a function of Dz in Fig. B1a for simulations CIfront1 and SIfront4, where the simulations whose results are presented in section 4 correspond to the finest resolutions in the panels. It is clear that in all simulations meet or exceed the threshold identified by Khani (2018), with the finer simulations being well within the necessary threshold. It is worth mentioning that, as seen in Fig. B1a, symmetrically dominated CSIs can be well represented with much coarser grids in our setup due to the smaller stratification, contributing to larger L O . We also show the cumulative mixing efficiency C ' in

a. Shear production terms for the primary instabilities
The general definition of the shear production terms comes from the turbulent kinetic energy prognostic equation and reads where U i = (U, V, W) is a Reynolds-averaged velocity vector about which the turbulent fluctuations u i are calculated, and summation is implied for the i index only (Stull 1988). In this case we want to consider the rate at which shear of the average flow transfers energy to the primary instabilities, namely, the CSIs. Thus, the fluctuations u i should ideally capture the CSIs only. Given the nature of our setup, this is challenging to achieve with directional averages (recall that CSIs are mainly 2D in nature, so averaging in the x direction would not achieve this result). Hence we consider an ensemble average over many realizations of this flow and make the assumption that such an average of the flow velocities is well approximated by the flow velocities at the initial condition. We then approximate U i as the flow velocities at the initial condition [given by Eqs. (9)-(14)], simplify Eq. (C1) accordingly, and define horizontal and vertical shear production rate terms for the primary instabilities as SP prim h 2 u y U y , SP prim y 2 u w U z : According to this definition the shear production rate is zero at t = 0 and starts to evolve as the instabilities start to develop. We quantify the value of the shear production rate terms at a time t = 15/v max , where v max is the maximum growth rate for CSIs [Eq. (5)]. This choice of time captures a well-developed CSI before the onset of full 3D turbulence. Different choices of time were investigated (including some based not on v max but on the evolution of « k ) and the results were found to be robust.

b. Shear production terms for the secondary instabilities
For this section, the purpose of the analysis is to capture the rate of energy input into the secondary instabilities by the CSIs. Ideally, it is then necessary to capture only the secondary instabilities in the fluctuation terms u i , and the background flow (with the CSIs) should be captured in the U j terms. Similarly to the primary instabilities analysis, the best approach we found is to consider an ensemble average that we assume to be well approximated by the state of the flow at a time t = t 1 in which the primary instabilities are well developed, but the secondary instabilities still have not started emerging. This choice is done manually, since a programmatic way to choose t 1 consistently across simulations could not be found. We found, however, that the precise choice of time does not alter the results significantly as long as the two aforementioned criteria are observed and as long as we consider a portion of the domain that isolates the emergence of secondary instabilities.
For these calculations U i = (U, V, W) Þ 0 (since they correspond to CSIs), and for a 2.5D setup (without an x direction) we can define these shear production rate terms for the secondary instabilities as SP second h 2u y U y 2 y 2 V y 2 w y W y , [defined in Eq. (19)] produces smaller values in a LES than it otherwise would if all scales of motion were resolved. However, provided that the subgrid-scale closure accurately represents the effects of eddies smaller than the grid scale, the values and distribution of « k should roughly coincide with real-ocean values and Re b can be used.
For this analysis we apply Eq. (8) with the averaging operator · q and calculate Re b for our simulations as Re b « k q n mol N 2 0 : (D1) Figure D1 shows the instantaneous mixing efficiency g as a function of Re b for all simulations in this work. Despite the larger scatter and the more ambiguous decay slope for Re b տ 10 2 , the qualitative behavior of points is similar to that of Fig. 6a: a plateau of g ≈ 0.3 and a subsequent decay of g with increasing Re b . Moreover, we see that the range of values of Re b and location of the peak in g roughly match that of previous idealized studies (Shih et al. 2005;Salehipour and Peltier 2015). This result, together with the resolution sensitivity test in appendix B, suggests that the mixing efficiency in the LES follows the same trends as in previous direct numerical simulations.