Parameterization of Frontal Symmetric Instabilities. I: Theory for Resolved Fronts

A parameterization is proposed for the effects of symmetric instability (SI) on a resolved front. The parameterization is dependent on external forcing by surface buoyancy loss and/or down-front winds, which reduce potential vorticity (PV) and lead to conditions favorable for SI. The parameterization consists of three parts. The first part is a specification for the vertical eddy viscosity, which is derived from a specified ageostrophic circulation resulting from the balance of the Coriolis force and a Reynolds momentum flux (a turbulent Ekman balance), with a previously proposed vertical structure function for the geostrophic shear production. The vertical structure of the eddy viscosity is constructed to extract the mean kinetic energy of the front at a rate consistent with resolved SI. The second part of the parameterization represents a near-surface convective layer whose depth is determined by a previously proposed polynomial equation. The third part of the parameterization represents diffusive tracer mixing through small-scale shear instabilities and SI. The diabatic, vertical component of this diffusivity is set to be proportional to the eddy viscosity using a turbulent Prandtl number, and the along-isopycnal tracer mixing is represented by an anisotropic diffusivity tensor. Preliminary testing of the parameterization using a set of idealized models shows that the extraction of total energy of the front is consistent with that from SI-resolving LES, while yielding mixed layer stratification, momentum, and potential vorticity profiles that compare favorably to those from an extant boundary layer parameterization (Large et al., 1994). The new parameterization is also shown to improve the vertical mixing of a passive tracer in the LES. © 2016 Elsevier Ltd. All rights reserved.


Introduction
The oceanic surface mixed layer is the contact point between the atmosphere and the ocean, and is responsible for communicating atmospheric fluxes into the ocean interior, including fluxes associated with heat, momentum, carbon, oxygen, and other tracers. In addition to influencing transport into the ocean interior, atmospheric fluxes strongly influence the dynamics within the mixed layer itself and drive vertical mixing. At the same time, lateral density gradients tend to increase the vertical stratification within the mixed layer Garrett, 1994, 1995;Hosegood et al., 2006 ). Submesoscale turbulence is the byproduct of the dynamical interplay between lateral density gradients and atmospheric forcing, and has been the focus of many recent studies (e.g. Boccaletti et al., 2007;Capet et al., 2008b-d;; * Corresponding author. E-mail address: sb965@cam.ac.uk (S.D. Bachman). Thomas and Ferrari, 2008;Mahadevan et al., 2010;McWilliams, 2010;Hamlington et al., 2014;Callies et al., 2016 ).
The oceanic submesoscale is typically defined as occupying a range of horizontal scales between 100 m and 10 km. Models with an O(1) km horizontal grid permit many of the important dynamical processes at these scales (e.g., Oschlies, 2002 ), especially the largest fronts and baroclinic instabilities. A more precise definition of submesoscales are those which are marginally constrained by rotation and stratification (the Rossby and Richardson numbers are both O(1) ,  and thus feature both geostrophic and ageostrophic components. Submesoscale dynamics are typically considered to be hydrostatic, which places a lower bound on the range of horizontal scales near the mixed layer depth, or O(100) m, as the aspect ratio approaches unity. Because the Rossby number is near one, the characteristic timescale is 1 / f O(1) day, which is much faster than the mesoscale and results in submesoscales playing a leading role in the evolution of the mixed layer in response to atmospheric forcing. Finally, submesoscale eddies frequently arise as a result of straining by mesoscale eddies, and their own strain fields in turn beget submesoscale frontogenesis Taylor, 2013, 2014 ) and a wide variety of smallscale instabilities that cascade energy down to dissipative scales ( Capet et al., 2008b;Molemaker and McWilliams, 2010 ).
It is the fast evolution of the submesoscale in conjunction with atmospheric forcing that is of principal dynamical interest here. The submesoscale features lateral density gradients that arise via eddy straining, preconditioning the flow to convective instabilities driven by both the wind and surface buoyancy forcing. In particular, when the wind blows in a down-front direction the Ekman transport of dense fluid over light can destabilize the mixed layer. Buoyancy forcing can either exacerbate this destabilization in the case of a positive flux (where the positive z direction is taken to be out of the ocean) or mitigate it for a negative flux. Which among the "zoo" of possible submesoscale instabilities that may arise depends heavily on the local relative vorticity ( Haine and Marshall, 1998;Thomas et al., 2013;Shcherbina et al., 2013 ).
Of the many types of submesoscale instability that arise in the mixed layer, baroclinic instability has thus far received the most attention (e.g. Boccaletti et al., 2007;. From a modeling perspective, the dynamics of baroclinic instability are fairly well-understood and have led to an extensive body of research about how it should be parameterized in eddy-free general circulation models (GCMs). The success of the most well-known of these parameterizations, the Gent and McWilliams (1990) parameterization for mesoscale baroclinic instability, has led to several subsequent papers regarding the nature of tracer transport by subgridscale motions in a variety of dynamical regimes ( Gent et al., 1995;Tandon and Garrett, 1996;Dukowicz and Smith, 1997;Killworth, 1997;Treguier et al., 1997;Griffies, 1998;Griffies et al., 1998;Greatbatch, 1998;Smith and Gent, 2004 , and more). A similar conceptualization to that of Gent and McWilliams (1990) underlies the parameterization of mixed-layer, ageostrophic baroclinic instability at submesoscales Fox-Kemper et al., 2011 ).
A new modeling challenge arises when the resolution of an ocean model permits partial resolution of the eddy field. At this point parameterizations must be carefully constructed and previous parameterizations carefully recast so that they do not either outcompete the resolved eddies for energy (e.g. Henning and Vallis, 2004 ) or double-count the effects of the eddies ( Delworth et al., 2012 ). Unlike Gent and McWilliams (1990) , the parameterization of Fox-Kemper et al. (2011) explicitly depends on the model grid scale and thus exemplifies how physical scalings may allow parameterizations to be recast to adapt to increasing model resolution. Fox-Kemper et al. (2011) can be used even when large mixed layer eddies (hereafter MLE) are resolved. The parameterization for oceanic symmetric instability (hereafter SI) put forth here is intended to be used in models where some fronts and large-scale frontal instabilities (e.g., baroclinic instabilities) are resolved or handled with scale-adaptive parameterizations, but the smaller frontal instabilities, in particular symmetric instabilities, are not resolved. Many high-resolution nested and regional models fall into this category (e.g., Capet et al., 20 08a, 20 08b;Lévy et al., 2010;Sutherland et al., 2011;Mensa et al., 2013;Zhong and Bracco, 2013;Molemaker et al., 2015;Rosso et al., 2015;Siedlecki et al., 2015;Gula et al., 2016 ).
In addition to parameterizations for mesoscale and submesoscale baroclinic instability, many parameterizations exist for small-scale turbulence and mixing. Shear and convective instability are commonly parameterized (for example by Kraus and Turner, 1967;Mellor and Yamada, 1982;Large et al., 1994 ) and, more recently, Langmuir turbulence (e.g., McWilliams and Sullivan, 20 0 0; Smyth et al., 2002;Grant and Belcher, 2009;Van Roekel et al., 2012;Harcourt, 2013;Li et al., 2016 ). Recent work has shown that the mixed layer stratification, energy dissipation, and resolved eddy behavior can be highly sensitive to the details of such parameterizations (e.g. Mukherjee et al., 2016 ). To date, however, no parameterization exists for SI despite observations of its effects on the shear, stratification, and dissipation of kinetic energy in the surface mixed layer ( D'Asaro et al., 2011;Thomas et al., 2016 ).
Numerical simulations with an O(1) km grid that resolve submesoscale fronts and larger MLE typically do not resolve smaller SI. Thus, it is desirable to have a parameterization for the O(100) m SI so that submesoscale simulations of an active field of fronts do not require the added cost of resolving the SI (e.g., , just as KPP or other turbulence parameterizations are used to avoid the cost of resolving O(1) m features (e.g., Hamlington et al., 2014 ). The goal of this paper is to propose a framework for a parameterization that approximates the restratification and tracer mixing by SI in the case where SI modes are unresolved, but the front that undergoes SI is resolved. The performance of this parameterization under different forcing scenarios is evaluated against the results of SI-resolving Large Eddy Simulations (hereafter LES, e.g. Taylor and Ferrari, 2010;Thomas et al., 2013;Hamlington et al., 2014 ) in some representative cases.

Basics of symmetric instability
The principal focus here, SI, arises in baroclinic flows featuring a lateral density gradient and an associated vertically sheared geostrophic flow. SI is typified by overturning circulations about an axis aligned with the geostrophic flow, typically with the flow along density surfaces ( Eliassen, 1949 ). Assuming a geostrophically balanced flow with buoyancy gradients N 2 = ∂ b/∂ z and ∇ h b = ( ∂ b/∂ x, ∂ b/∂ y ) , the SI growth timescale T and horizontal lengthscale L may be estimated for constant shear and stratification ( Stone, 1966 ) as where Ri b is the balanced Richardson number, For typical ocean mixed layer parameters in conditions favorable for SI, 0.01 ≤ U ≤ 0.1 m s −1 , 25 ≤ H ≤ 100 m, and 0.25 ≤ Ri b ≤ 0.95, which imply SI timescales ranging from one minute to an hour and lengthscales of 50-2500 m. These are much smaller than the timescale of 14-18 h and lengthscale of 600 m to 8 km for MLE (e.g. Fox-Kemper et al., 2008 , Eq. (2) and (3)), highlighting the potential importance of SI to the evolution of the mixed layer on fast time scales. The corresponding SI mixing rates can be estimated using in situ measurements of density, wind, and surface heat flux, and are compared against those of MLE in Appendix A .
SI occurs when the sign of the Ertel potential vorticity is anticylonic -i.e., where q is opposite in sign to the Coriolis parameter f , so that fq < 0 ( Hoskins, 1974 ). Here u is the velocity, ˆ k is a unit vector in the vertical, and b = −g ( ρ − ρ 0 ) /ρ 0 is the buoyancy, which is defined in terms of the gravitational acceleration g , density ρ, and constant background density ρ 0 . The condition fq < 0 is most straightforwardly satisfied in convective conditions ( N 2 < 0) which can give rise to gravitational instability, or inertial instability when N 2 > 0 and ζ a < 0, where ζ a = ω a ·ˆ k = f − ∂ u/∂ y + ∂ v /∂ x is the vertical component of the absolute vorticity ω a . These types of instability are not the focus of this study and will not be mentioned further.
In the presence of a lateral density gradient, it is possible that fq < 0 even when N 2 > 0 and ζ a > 0. This is most readily seen by assuming that to leading order the flow is in thermal wind balance so that by which This expression makes clear how lateral buoyancy gradients can "precondition" a flow to instability. When surface forcing reduces the vertical stratification, the destabilizing baroclinic term, −| ∇ h b | 2 , becomes relatively more important and can eventually cause fq to fall below zero. A parallel condition to fq < 0 can be written in terms of the stratification as well as the vorticity, when the flow and vorticity are relatively uniform and Ekman shear and surface wave effects are negligible ( McWilliams and Fox-Kemper, 2013;Hamlington et al., 2014;Haney, 2015;Haney et al., 2015 ). In terms of the Richardson number, it can be shown that fq < 0 corresponds to ( Haine and Marshall, 1998 ). If the geostrophic flow features negligible vertical relative vorticity then f q = f 2 N 2 − | ∇ h b | 2 and Ri b < 1 becomes the sufficient condition for SI to develop Ferrari, 2009, 2010 ). It should be emphasized, however, that fq < 0 is a more general and stringent criterion. SI-resolving simulations when shear and stratification vary substantially show that the potential vorticity criterion holds even when no relevant Richardson number criteria exist ( Haney, 2015;Haney et al., 2015 ). The test cases examined here fall into the simple category where the Richardson number and potential vorticity criterion are identical. Extension to other scenarios where this is not the case (e.g., Hamlington et al., 2014;Haney, 2015;Haney et al., 2015 ) will be tested in future work.
In the surface mixed layer, the anticyclonic PV criterion can be created by surface forcing of PV that is destabilizing ( Thomas, 2005 ). When fq < 0 SI is the most unstable mode for 0.25 < Ri b < 0.95 ( Stone, 1966( Stone, , 1970, and will act to restore the fluid to a marginally stable state ( Thorpe and Rotunno, 1989 ) by mixing in fluid of higher PV from either the thermocline or the surface ( Taylor and Ferrari, 2009 ). The mixed layer is restratified in density as part of this process, achieving a buoyancy frequency of N 2 = | ∇ h b | 2 / f ζ a upon becoming marginally stable to SI when f q = 0 . The SI-neutral state may still be unstable to other, generally slower, forms of instability such as mixed layer baroclinic instability (MLI).
As noted in  , restratification by SI typically exceeds the frontal restratification rate of MLI (see also Appendix A ). MLI is unlike SI in that it is nearly in thermal wind balance in all directions and occurs at larger, more easily resolved, scales. In a freely decaying front, SI first restratifies the front until the front is no longer unstable to SI, which usually coincides with a Richardson number based on the thermal wind shear Ri b = For fq > 0, SI are stabilized, but MLI continues to extract potential energy from the front and restratify. However, when the ocean is forced by convection or downfront winds, fq may be kept below zero for some time, producing forced SI ( Taylor and Ferrari, 2010 ). Due to their scale separation, SI and MLI may form simultaneously with SI developing along sharp fronts and filaments generated by straining associated with mixed layer eddies. Smaller, nonhydrostatic convective or Langmuir instabilities may also co-exist with SI ( Hamlington et al., 2014 ). SI predominately grows by extracting mean kinetic energy from the geostrophic shear. The rate of this growth in the turbulent Down-front winds (thick arrow) drive an Ekman transport perpendicular to the wind vector and across the front (horizontal arrows), carrying dense water over light water. The associated Ekman buoyancy flux ( EBF > 0), in conjunction with atmospheric cooling ( B 0 > 0), reduces the stratification near the surface, destroys PV, and drives the fluid toward an SI -unstable state. The instability extracts energy from the geostrophic shear, and leads to enhanced turbulence and vigorous mixing along isopycnals (colored planes). The resulting circulation consists of overturning cells of alternating orientation (ellipses, with arrows indicating flow direction) and acts to bring light water over dense water in opposition to the EBF , restratifying the mixed layer.
kinetic energy (TKE) budget is given by the geostrophic shear production where the overbar denotes an average over the SI scale. Although for a flow with arbitrary rotation, stratification, and forcing, mixed modes (a combination of baroclinic, symmetric, gravitational, and inertial instabilities, and Langmuir circulations) may develop which derive energy from other sources in the turbulent kinetic energy budget as well ( Li et al., 2012;Thomas et al., 2013 ), the geostrophic shear production is identified as the principle energy source for SI ( Bennetts and Hoskins, 1979;Thomas and Taylor, 2010;Thomas et al., 2013 ). SI also acts as a downscale energy pathway for geostrophic flows. As energy is extracted from the geostrophic shear and the front is weakened, geostrophic adjustment can be expected to lead to a subsequent slumping of isopycnals and a release of mean potential energy. Taylor and Ferrari (2009) show that SI motions undergo a secondary Kelvin-Helmholtz shear instability, which eventually leads to enhanced small-scale dissipation and mixing. SI thus facilitates a conversion of available potential energy associated with the sloping isopycnals and geostrophic flow into TKE. SI can be sustained in the presence of surface forcing which continually contributes to the destruction of PV. The PV evolves according to ∂q ∂t where, using ω a = f ˆ k + ∇ × u , is a PV flux consisting of 1 advective, 2 frictional, and 3 diabatic components. Here F is a frictional or non-conservative body force and Db / Dt is the Lagrangian rate of change of buoyancy. At the surface boundary the frictional and diabatic forcing can supply negative anticyclonic potential vorticity and sustain SI if This combination of frictional and diabatic forcing drives PV toward anticyclonic values ( Fig. 1 ). Thomas (2005) showed that (9) can be related to the surface forcing via an effective buoyancy flux, namely that PV is destroyed if Here (τ w ×ˆ k ) /ρ 0 f is the Ekman transport and B 0 is the surface buoyancy flux driven by atmospheric cooling or evaporation (positive values of either term decrease the buoyancy). F SI is used throughout the rest of this paper as a "forcing" term for SI; although this is not a forcing in the usual sense that it results directly in the formation of SI, it does act as part of the forcing of the PV Eq. (7) and thus tends to drive the flow toward a SIunstable state. The term involving the wind stress is referred to as the Ekman buoyancy flux ( Thomas, 2005;Thomas and Taylor, 2010;D'Asaro et al., 2011 , hereafter EBF), which characterizes the horizontal transport of fluid across the buoyancy gradient in response to wind forcing. A positive EBF occurs during down-front winds, so that the Ekman transport advects fluid down the surface buoyancy gradient, or from the dense side of a front to the less dense.
The relation (10) is fundamental to the parameterization proposed here. The subtitle of the work, "Theory for Resolved Fronts" derives from the need to accurately estimate the covariance of wind components and horizontal buoyancy gradient in the first term of (10) . The parameterization proposed here does not require sufficient model resolution to permit SI, but it does require sufficient resolution to accurately capture the surface PV flux. Fox-Kemper et al. (2011) statistically correct for coarse-resolution models' tendency to underestimate buoyancy gradient magnitudes, but the SI parameterization requires both direction and magnitude estimates of the buoyancy gradient. In this work, it will be assumed that the resolved estimate of (10) is adequate. The case of partiallyresolved SI is discussed by Bachman and Taylor (2014) .
As previously noted, if the PV is low prior to a forcing event, then F SI > 0 will quickly lead to fq < 0. If the initial PV is large and positive, then F SI > 0 will tend to trigger small-scale convective and shear instabilities that will rapidly mix the near surface until nearly neutral ( Ri b ∼ 0.25) stratification and shear results along with geostrophic adjustment ( Tandon and Garrett, 1994 ). Thus, regardless of the initial PV, F SI > 0 will soon result in conditions ripe for SI that may be sustained under continued forcing, or act to accelerate restratification so long as fq < 0. Alongside SI, mixed layer eddies will also act to restratify, but are generally slower to do so than SI whilst fq < 0, but they will persist even after fq > 0 or Ri b > f / ζ a . Under sustained F SI > 0 , SI may remain the dominant restratifying instability, continually combining anticyclonic nearsurface PV with subsurface cyclonic PV.
Based on results from LES of forced SI, Taylor and Ferrari (2010) subdivided the low PV layer into two regions. Near the surface, in a 'convective' layer, turbulent mixing driven by surface cooling and down-front winds are able to maintain weak or even unstable stratification. The interior of the convective layer is associated with a positive vertical buoyancy flux, indicating a transfer of potential energy to kinetic energy. The magnitude of this buoyancy flux is predominantly set by the surface buoyancy flux even in the presence of a down-front wind stress (e.g. Thomas et al., 2013 ). Beneath the convective layer but still within the region of low PV, SI motions dominate the energy budget. In this 'SI layer' the buoyancy flux is generally small but negative, while SI extracts energy from the front through the GSP. The SI layer is also characterized by a stable stratification with Ri b = O (1) and vigorous isoneutral mixing, i.e., mixing tracers along the neutral surfaces of density .
The proposed parameterization in the following section will be developed by considering the basic form of the subgridscale terms in the primitive equations. In particular the variable Reynolds stresses will be parameterized so that the total reduction of resolved energy matches that from LES (e.g. Taylor and Ferrari, 2010;Thomas and Taylor, 2010;Thomas et al., 2013 ). The parameterization must also capture the sensitivity of SI to surface forcing, and be able to distinguish the boundary between the convectivelydominated and SI-dominated sublayers. Finally, the parameterization will also specify a variable isoneutral diffusivity, which mixes resolved tracer gradients in the SI layer and produces the desired profile of the vertical buoyancy flux w b in the convective layer, ensuring proper restratification and mixing of active and passive tracers.

The basic constituents
The goals of the parameterization are to represent the following processes: (1) Appropriate mixing of momentum, buoyancy, and tracers during destabilization by F SI > 0 .
(2) Extraction of energy from the resolved flow by SI.
A successful parameterization should also meet the following conditions: (4) No effect when F SI ≤ 0 or ∇ h b = 0 or fq > 0.
(5) Act only in the SI-unstable part of the surface boundary layer. (6) Maintain energetically consistent boundary conditions on momentum, buoyancy, and PV.
The discussion at the beginning of Section 5 will return to each of these items in turn to summarize how the parameterization meets each of these goals. The performance of the proposed parameterization is tested in Section 4 by comparing against SIresolving large-eddy simulations.
Consider first the Reynolds-averaged hydrostatic Boussinesq equations of motion for a fluid containing a tracer ξ , which may be written Here D/Dt = ∂ /∂ t + ū · ∇ is the material derivative operator, φ is the Boussinesq pressure potential, b is the buoyancy, and F ( u, v ) and D ( b,ξ ) are generic, spatially-variable mixing terms appropriate for each scalar. The overlines denote a spatial coarse-graining over SI scales, but it is assumed that the submesoscale or mesoscale fronts and frontogenetic strain are resolved. The primed quantities refer to deviations from the resolved scales, typified by the SI scales. Taylor and Ferrari (2010) quantify the effects of SI in numerical simulations for forced fronts underneath convective buoyancy fluxes. On subinertial timescales they show that the momentum balance that arises in the SI-unstable layer is where u a = u − u g is the horizontal ageostrophic velocity and the viscous terms ν ∂ 2 ū a ∂z 2 are important only in a thin sublayer near z = 0 .
An expression for the Reynolds stresses in this balance will be used as the foundation for the SI parameterization. In general, one could use this parameterization independently and choose to turn off all other boundary layer parameterizations (e.g. Mellor and Yamada, 1982;Kantha and Clayson, 1994;Large et al., 1994 ) and physics packages when the SI parameterization is on, or to use them to modify an existing parameterization. There are numerical and computational advantages to the latter approach, as the SI parameterization could be designed to use existing procedures (such as the calculation of boundary layer depth).
However, because many boundary layer parameterizations are constructed in a self-consistent manner and may involve specialized algorithms to ensure numerical stability, the details of how to incorporate the SI parameterization into each different scheme would have to be considered on a case-by-case basis. Therefore, the SI parameterization will be presented in a way that will work in a standalone manner. The SI parameterization will be compared against the popular K-profile parameterization ( Large et al., 1994 ) to show how the SI parameterization improves the parameterized mixing and dissipation in SI-unstable conditions.

Parameterization of the SI fluxes and numerical implementation
Surface forcing, through a combination of down-front winds and surface buoyancy loss, reduces PV ( Eq. (10) ) and is the principal driver of mixed-layer SI. Stratification in the mixed layer is then set by a competition between eddy-induced restratification, either by SI or mixed-layer baroclinic instability (e.g. Boccaletti et al., 2007 ), and vertical mixing resulting from the atmospheric forcing. As described in Section 2 , Taylor and Ferrari (2010) found that the competition between restratifying SI and surface forcing results in the formation of two dynamically distinct sublayers within the surface boundary layer -a convective layer of depth h where the surface forcing is sufficient to keep the density unstratified and where the convective buoyancy flux w b > 0 , and a deeper SI-dominated layer of total depth H where the restratifying effect of SI overcomes that of the surface forcing but where PV is still homogenized.
Because of this partitioning of the surface boundary layer, the vertical profiles of the GSP and w b are governed by convective forcing due to surface fluxes and EBF, and also by the difference between h and H . LES simulation results from Taylor and Ferrari (2010) and Thomas et al. (2013) show that the convective layer (0 to −h ) features a positive vertical buoyancy flux, w b , which is approximately linear in z (see also Large et al., 1994 ) and whose magnitude is predominantly set by the surface buoyancy flux, sug-gesting that it can be parameterized as 1 Taylor and Ferrari (2010) further argue that when lateral density gradients are present, the sum of the GSP and w b can be approximated as a linear function of z through the entire surface boundary layer (0 to −H). This sum is related to the surface forcing by where F SI is the sum of the Ekman and surface buoyancy fluxes as written in (10) . By substituting w b con v for w b in (20) and changing the approximation to an equality, Thomas et al. (2013) formed a parameterization for the GSP The GSP ( −u h w · ∂ ū g ∂z ) involves Reynolds stresses that tend to reduce the thermal wind shear by fluxing momentum downgradient. These momentum fluxes are needed in the SI parameterization and will be developed consistently following (21) . Assuming the stresses in the cross-front direction are not meaningfully correlated with the SI and manipulating the trigonometry ( Appendix B ), the Reynolds stresses consistent with SI extracting energy from GSP SI in (21) are The expressions above can be used to obtain two of the three core components of the SI parameterization: the parameterization of the convective layer vertical buoyancy flux is given by (19) , and the SI vertical viscosity can be derived using (21) . An eddy viscosity is obtained by assuming a flux-gradient relationship between the Reynolds stresses and the resolved geostrophic shear, and after substituting (22) into (23) , The above expressions for w b con v and ν SI require suitable definitions of h and H , which should be specified in a way that is consistent with the phenomenology of symmetric instability. Here the SI layer depth is defined according to the phenomenology that the flow is SI-unstable when q is of the opposite sign as f ; therefore, H will be the shallowest depth where a bulk measure of the potential vorticity satisfies the criterion that Here refers to the change in the quantity from the surface to z = −H, and the angle brackets indicate a vertical average over the same depth range. Requiring the criterion on potential vorticity in (25) to use bulk quantities is more numerically stable than using local values of the shear and buoyancy gradient, which would be more prone to noise and the degree to which the model conserves PV.
Once H is known, a quartic equation (e.g. Thomas et al., 2013 ) can be used to solve for the convective layer depth h , where w * = ( B 0 H ) 1 / 3 is the convective velocity, u * = | τ w | /ρ 0 is the friction velocity, θ w is the angle between the wind vector and the geostrophic shear, and c = 14 is an empirical constant. When h / H 1 much of the surface boundary layer is dominated by SI restratification, whereas for h / H 1 the layer is dominated by convective mixing and SI is not expected to be important. These two cases are distinguished by the square-bracketed term, which is large for strong surface forcing ( h / H 1) and small for strong lateral gradients and fronts ( h / H 1). Eq. (26) always has two complex and two real solutions for h / H , and only one real solution falls between 0 and 1 (see Appendix C ). Furthermore, h is positive and no greater than H , so there is no risk of unexpected consequences in the numerical solution of (26) .
The respective forms for w b con v and GSP SI shown above are informed by LES simulations which studied the evolution of forced SI in the presence of destabilizing surface forcing (e.g. Taylor and Ferrari, 2010;Thomas and Taylor, 2010;Thomas et al., 2013 ). A notable limitation of these LES is that, in each case, both EBF ≥ 0 and B 0 ≥ 0, so that it is unclear how (19) and (21) generalize to cases where the Ekman buoyancy flux and surface buoyancy flux are of opposite sign, or when the wind is predominantly in the crossfront direction. Care must be taken when considering cases such as these; for example, if a strong wind blows in the cross-front direction the resulting EBF would be small, and GSP SI as predicted by (21) would likely underestimate the true energy dissipation by the resulting turbulence. Therefore, for now it may be preferable to use the SI parameterization only when the wind is dominantly downfront and when B 0 ≥ 0. This caveat is included in the equation for the GSP parameterization, (21) . Further research is needed to determine appropriate forms for w b con v and GSP SI that span a broader range of surface forcing scenarios, at which point the potential uses for the SI parameterization can be expanded accordingly.
Many ocean models integrate equations for potential temperature, θ , and salinity, s , rather than integrating the buoyancy equation directly, requiring an adaptation of the convective buoyancy flux parameterization in (19) . The LES simulations informing the approximation for w b con v use buoyancy as a state variable rather than θ and s , so details of how the linear structure function in for a thermal expansion coefficient α θ , saline contraction coefficient β s , and reference temperature and salinity values θ 0 and s 0 .
By the linearity of (D.10) , it follows that Defining the surface potential temperature flux as w θ 0 and the surface salinity flux as w s 0 , if it is assumed that the same linear vertical structure in (19) applies to θ and s , so that one recovers the parameterization for w b con v by substitution of (29) into (28) . A similar approach can be made for the convective layer fluxes of other passive tracers which have a nonzero surface flux. Implementation of the SI parameterization then amounts to substituting the appropriate expressions for the convective layer fluxes, as well as those for ν SI and the vertical and along-isopycnal diffusivities ( Section 3.3 ), into the momentum and tracer equations. The full equation set with the SI parameterization included is summarized in Section 3.4 .

Isoneutral and vertical tracer mixing
Taylor and Ferrari (2011) compared large-eddy simulations that were unstable to upright convection and symmetric instability, and found that the vertical mixing rate of passive tracers was greatly reduced in symmetric instability compared to upright convection under the same forcing conditions. Their simulations also included a simple model for light-limited phytoplankton growth, and they found that the reduction in vertical mixing in conditions favorable to symmetric instability resulted in significantly higher phytoplankton concentrations compared with upright convection. These findings suggest that capturing the influence of symmetric instability on the scalar mixing rate is important to accurately model biogeochemical processes at density fronts.
The structures that arise from SI are generally strongly horizontally anisotropic, aligning themselves with isopycnals, and therefore one must be careful in utilizing an isotropic viscosity or diffusivity to parameterize them. We are particularly concerned with transport across density surfaces and transport along density surfaces but across the front. When SI is present, it will greatly enhance the latter. To most simply approximate this cross-front, along-isopycnal transport, it is appropriate to diagnose a scalar vertical diffusivity for the diabatic mixing, and to separately obtain the larger associated cross-front, along-isopycnal mixing coefficients via a tensor rotation.
Previous LES results (e.g. Taylor and Ferrari, 2010;Thomas and Taylor, 2010;Thomas et al., 2013 ) indicate a negative vertical buoyancy flux is present in the SI layer, consistent with mixing by small-scale turbulence driven by shear associated with alongisopycnal SI cells (e.g. Taylor and Ferrari, 2009 ). The vertical diffusivity associated with this mixing, κ SI, v , can be related to ν SI by defining a turbulent Prandtl number such that Note that ν SI , which was defined in Eq. (23) , is a vertical viscosity.
Results from atmospheric boundary layer studies (e.g. Grachev et al., 2007;Anderson, 2009;Kitamura et al., 2013 ) suggest that in a weakly stratified flow Pr T is an increasing function of the local gradient Richardson number, Ri g = N 2 / | ∂ ū /∂z | 2 , and tends to be O(1) for the range of Ri g Ri b < 1 in a SI-unstable flow. Although properties of turbulence in the atmospheric boundary layer are likely to be very different from SI, after evaluating a number of possibilities, a functional form and coefficients informed by the regression analysis of Anderson (2009) is adequate, which is dependent on the local value of Ri b . Fig. 2 shows that diagnostics from previous LES (the details of which are described in Section 4 ) support this scaling for Pr T . In this figure, values of Pr T are saved from the LES by calculating where the overbar indicates a horizontal average over the entire model domain. These values of Pr T are saved along with the corresponding value of Ri b at every depth and at regular intervals in time. The full set of values is then filtered to only include pairs of ( Ri b , Pr T ) after one day of model time has elapsed (to allow SI time to spin up) and in the depth range −1 . 2 h > z > −H m, where an extra twenty percent is added to the convective layer depth to avoid very weakly stratified regions where entrainment is occurring immediately below the convective layer. Values of Ri b are then binned at intervals of 0.05 and the median value of Pr T is calculated for each bin (blue dashed line), and are compared against the scaling from (31) (red line). The median Pr T values generally follow the trend of the scaling, though there is a wide scatter in the range of ( Ri b , Pr T ) pairs. It is unclear whether this relationship between Pr T and Ri b holds in the convective layer, where turbulent transport is largely nonlocal and the destabilizing surface forcing can lead to negative values of Ri b . For now we choose to employ (31) in both the convective and SI layers, which is modified slightly to exclude negative values of Ri b so that Further study is required to clarify whether this scaling for Pr T holds in the convective layer, at which time the expression (33) can be modified accordingly. 2 A principal feature of SI is the formation of individual overturning cells, oriented along tilted isopycnals. To complete the SI parameterization the along-isopycnal mixing associated with these 2 The simulations in Section 4 were also run using constant values of Pr T = 1 and Pr T = 5 to test the sensitivity of changing the value of Pr T . Switching between these choices led to small changes in diagnostics such as the boundary layer depth, PV, and stratification profiles (not shown), but overall the parameterization and its key features remained robust to the specific choice of Pr T . cells should also be represented. Here it is emphasized that this part of the SI parameterization, and the diffusion tensor K SI which is to be defined, is separate from the mixing represented by κ SI, v .
For simplicity, it is assumed that SI motions lead to equal diffusivities in the along-front and cross-front directions.
One may parameterize the eddy flux of a passive tracer ξ in the form of a flux-gradient relationship, where K SI is an eddy diffusion tensor that ensures the proper anisotropy in the horizontal versus vertical directions (e.g. Redi, 1982 ). If it is assumed that the flux is directed along isopycnal surfaces with no diabatic component, in z -coordinates such a tensor can be written where (b x , b y , b z ) = ∇b represent the directional buoyancy gradients, and κ SI, I is an isotropic, along-isopycnal scalar diffusivity.
The eddy diffusion tensor can be simplified if it is assumed that the isopycnal slopes, S = −∇ h b/N 2 , are small (e.g. Cox, 1987;Griffies, 1998;Griffies et al., 1998 ), which occurs when | ∇ h b | < < b z . If restratification by SI is able to maintain Ri b 1 then the isopycnal slopes scale as | S | f 2 /| ∇ h b |. The isopycnal slopes are therefore small unless the front is weak, in which case it is likely that SI plays only a minor role in boundary-layer turbulence. However, it is unclear whether mixing is predominantly alongisopycnal in the convective layer where stratification is weak and isopycnal slopes can be large even for weak fronts, and provisionally this part of parameterization will be developed assuming that K SI can be used in both the convective layer and SI layer. Then because the convective layer isopycnal slopes may be large, and because in small-slope regions the extra terms only add higher-order corrections in any case, it is recommended that the full diffusion tensor be used through the whole boundary layer.
The tensorial structure of K SI specifies that the principal axes of diffusion should be in the along-isopycnal plane, but the magnitude of κ SI, I must be scaled appropriately for a symmetrically unstable mixed layer flow. To this end, suppose such a flow has overturning cells that are nominally tilted along isopycnals. For this purpose, let x denote the along-front direction and y the crossfront direction, and define a horizontal diffusivity associated with lateral fluid parcel displacements within these cells, κ SI,h ∝ v l , for a cross-front velocity perturbation v and displacement l ( Taylor, 1921 ). By (22) , one may scale for the along-front Reynolds stress associated with the SI cells, An appropriate scaling for the along-front velocity perturbation is Since the SI overturning cells are slantwise along isopycnals, a scaling for the cross-front velocity perturbation therefore is where | S | = b y /N 2 is the isopycnal slope. Assuming the displacement is along isopycnals, l ∝ N 2 H/ b y , and the scalings for l and v may be combined to obtain an estimate for the horizontal diffusivity, Here it is understood that this scaling is unique up to multiplication by a dimensionless constant C , which by the Buckingham Pi Theorem ( Buckingham, 1915 ) can be written as a function of only two nondimensional parameters, C = C Ri b , b y / f 2 . For the simulations examined here ( Section 4 ), a simple constant coefficient, C = 1 . 0 , provides a good empirical fit to the LES. This will be the form used throughout the rest of this paper. 3 The preceding scaling relations can now be combined with the structure of K SI to obtain the full SI along-isopycnal diffusion tensor. The horizontal diffusivity can be incorporated by noting that in the SI layer where it is expected that b z > > b x , b y , the dominant horizontal terms lie on the diagonal of the tensor, so that to leading order κ SI,I κ SI,h . The full diffusion tensor then becomes where the minimum operator has been added to ensure that the diffusivity remains bounded even where the local value of Ri b is large, such as might be found at the boundary between the SI layer and pycnocline. As K SI takes the same form as the Redi (1982) tensor for subgridscale baroclinic instability (see Griffies, 1998 ), adding the SI parameterization to an existing tracer transport package should be straightforward. Since the parameterization here presumes that the fronts are resolved, it is unlikely that this parameterization would be used simultaneously with the mesoscale forms of Gent and McWilliams (1990) and Redi (1982) , but the numerical methods should be similar. It might be used simultaneously with a submesoscale baroclinic instability parameterization (e.g., Fox-Kemper et al., 2011 ). In doing so the reader is reminded that the criteria for the SI parameterization to be switched on ( fq bulk < 0) is independent from that for the baroclinic instability parameterization, and that these criteria should be evaluated separately even if their respective transport tensors are combined. Lastly, because of weak stratification and small vertical grid spacing, the modeled surface boundary layer is very sensitive to spurious mixing induced by the discretization of K SI , and thus it is important to ensure numerically that K SI · ∇ b = 0 as closely as possible. Suggested discretization techniques for this tensor are presented in Griffies et al. (1998)

Summary and implementation of the SI parameterization
The SI parameterization has been constructed to simulate the restratification and mixing by SI in a model where it is unresolved, in a way that satisfies the energetics of fully-resolved SI and satisfies the dominant force balance (17) and (18) for a convectively forced, SI-unstable mixed layer. As many different concepts and elements of the parameterization have been introduced piecemeal in the preceding sections, it is useful to now summarize the parameterization and its algorithm.
The SI parameterization represents the downscale transfer of kinetic energy via geostrophic shear production, and the alongisopycnal mixing by SI. With all parts of the parameterization added in and assuming all other SGS parameterizations are included in the forcing and mixing terms, F = (F u , F v ) and D, the momentum and tracer equations now take the form The vertical momentum and continuity equations remain unchanged.
The numerical implementation of the SI parameterization proceeds as follows. The depth of the surface boundary layer, H , is found in each vertical column by finding the shallowest point where fq bulk > 0. The depth of the convective sublayer, h , is then found by solving the quartic equation whose solution is a simple algebraic expression given in Eq. (C.12) . Given the surface wind stress, tracer flux, and local buoyancy gradients, ν SI , κ SI, v , w χ con v , and K SI can now be calculated at each depth, which are then substituted into the momentum and tracer equations. At depths below z = −H, all fluxes and mixing coefficients from the SI parameterization are set to zero, so that any interior mixing is handled by other parameterizations of the modeler's choice. Because the SI parameterization represents a form of boundary layer turbulence which is present only under certain forcing conditions ( EBF, B 0 ≥ 0) and only when fq bulk < 0, it is straightforward to implement it alongside another boundary layer package such as KPP. For example, a standard implementation of KPP would calculate momentum and tracer fluxes in each vertical column of the  Coding the SI parameterization in this way ensures that each vertical column features some sort of boundary layer mixing in the presence of destabilizing forcing. It also allows KPP (or another boundary layer parameterization of choice) to act as the default scheme when the SI parameterization is inactive. Because it is not clear how boundary layer turbulence behaves when the wind has a large cross-front component or the forcing is destabilizing but of opposite sign ( EBF × B 0 < 0 and EBF + B 0 > 0 ), we provisionally suggest turning the SI parameterization off when the wind is not dominantly downfront or B 0 < 0. Further research is needed to understand the dynamics and energetics of boundary layer turbulence in these scenarios, at which point the SI parameterization can be expanded accordingly.
Finally, based on the definition of q bulk in (25) , it is possible that the SI parameterization can become active when N 2 0 and the lateral buoyancy gradient is extremely weak. In this case the isopycnal slope becomes nearly vertical and the turbulence would become more akin to classical upright convection, in which case KPP is the more appropriate choice. These weak-front cases are detected through the calculation for h in (50) , in which case h / H 1. Therefore, for both numerical stability and physical reasons, it is recommended to turn the SI parameterization off when h / H exceeds some threshold value (here chosen to be 0.9). The specific choice of this threshold value required to ensure numerical stability has not been exhaustively tested.

Simulations
A new set of routines have been written for the MIT General Circulation Model ( Marshall et al., 1997 ) which adds the SI parameterization given by (42) -(49) . To test the effects of the parameterization, a series of 2D channel models has been created using the MITgcm and are compared against matching turbulenceresolving LES. The LES are run using a fully nonhydrostatic, spectral flow solver, the details of which can be found in Taylor (2008) and Bewley (2010) .
The three LES test cases used here have been chosen from prior literature (e.g. Taylor and Ferrari, 2010;Thomas and Taylor, 2010;Thomas et al., 2013 ) to evaluate the skill of the SI parameterization at reproducing the vertical profiles of the GSP, buoyancy flux, and state variables from eddy-resolving LES under a range of surface forcing conditions. They are chosen specifically to test the parameterization in cases of surface heat loss and no wind stress ( Taylor and Ferrari, 2010 , simulation 3 D 2 ), wind stress with no heat loss ( Thomas and Taylor, 2010 ), and a case where both types of forcing are large ( Thomas et al., 2013 ). The fixed parameters from each simulation are listed in Table 1 . In these LES the horizontal viscosity and resolution are small enough that all SI modes are resolved, and so that the flow would fully restratify to a state where fq > 0 if the surface forcing were to cease. The MITgcm channel models are initialized with lateral and vertical density gradients to match each LES test case, whose stratification parameters are described in the original papers. The linear equation of state in (27) is used with β s = 0 for simplicity, so that the density is only a function of the potential temperature. The along-front velocity is initialized to be locally in thermal wind balance with the lateral density gradient. The presence of the lateral gradient preconditions the system to SI, and along with a destratifying surface forcing is sufficient to spur the growth of SI in the LES, or activate the SI parameterization in the models which do not resolve SI.
Three versions of each LES are run using the MITgcm, all using coarsened grids in comparison to the LES. One version is run with the SI parameterization active, yet that reverts to KPP in regions that are SI stable. The other two are run with KPP only, where the KPP critical Richardson number in all simulations is set to be 0.3. The two KPP runs are distinguished by whether the optional KPP shear instability parameterization is enabled (see D.1 ), which can have a significant effect on the shear and stratification over time (e.g. Fig. 8 ).
The domain in the MITgcm simulations contains (N X , N Y , N Z ) = (3 , 400 , 100) gridpoints with rigid lids at the vertical boundaries, where the vertical resolution is set according to the depths in Table 1 and the horizontal resolution is set to x = y = 50 0 0 m. The simulations are meant to be approximately 2D, where a small number of gridpoints are added in the along-front dimension to avoid computational issues related to running the MITgcm in true 2D mode. The domain is set in a "double-front" configuration ( Fig. 3 ), where d b dy = M 2 is constant between points 1-150 and reversed in sign ( = −M 2 ) between points 200-350, and is smoothly transitioned over fifty points separating these regions. This setup allows the domain to be periodic in y so that there is no influence from lateral boundaries, and there is sufficient separation between the two fronts that they do not interact over the course of each simulation. The density gradient over the first 150 points matches the values from the LES ( Table 1 ), and represents the area of interest for the comparisons in this study. When referring to the model output mathematically in all figures and equations in this section, an overbar represents a horizontal average over this region of the domain and over one inertial period in time, and primes indicate deviations from this average. Each model contains a mixed layer overlying a "thermocline" where the stratification is high enough to be stable to SI throughout the duration of the simulation. The depth of the mixed layer varies in each model. In the Thomas and Taylor (2010) and Thomas et al. (2013) models the initial mixed layer PV is low enough that the SI parameterization activates immediately. In the Taylor and Ferrari (2010) model, the mixed layer is initially stable to SI, but becomes destabilized over time due to the surface forcing. The horizontal viscosity in all MITgcm simulations is set to ν h = 10 m 2 s −1 , and the background vertical viscosity is ν v = 10 −4 m 2 s −1 . Together these parameter settings do not permit SI modes to be resolved in the domain ( Bachman and Taylor, 2014 ), so that the performance of KPP and the SI parameterization can be measured via comparison with the LES. Thomas et al. (2016) show that even when the instantaneous momentum budget has a nonzero acceleration due to inertial motions, the turbulent Ekman balance in (17) and (18) may still hold as long as the averaging period is longer than the inertial timescale. Because these equations form the basis for the SI parameterization, we opt to compare the models using output averaged over two inertial periods. In this comparison it is assumed that super-inertial, high-frequency variability does not significantly contribute to either the momentum or buoyancy budgets.

Boundary layer depths
Because the wind is downfront in the averaging region of the channel models, the EBF and/or surface heat loss jointly begin destratifying the flow here, starting at the surface and progressing into the interior. The convective and boundary layer depths progressively deepen as the dense, low-PV surface water is mixed downwards. Fig. 4 shows time series of h and H diagnosed in both the LES and MITgcm. In the LES and SI parameterization simulations H is determined using the criterion on q bulk in (25) . The convective layer depth, h , is found in the LES by finding the deepest point where w b > 0 , and in the SI parameterization simulations by solving the quartic Eq. (26) using the value of H at the corresponding time. The KPP boundary layer depths, h KPP , which are diagnosed in the "KPP" simulations are plotted for reference as well. These plots are meant to compare the depth to which the parameterizations mix versus an approximate depth to which the resolved SI mix in the LES, indicating the skill of the parameterizations at predicting the appropriate boundary layer depth.
In the SI parameterization simulations, the criteria for determining H given in (25) and h given in (26) show skill at matching the depth from the LES. The diagnosed KPP boundary layer depth is usually significantly deeper than the LES convective layer depth and tends to be more similar to H . As a result, the tendency is for KPP to parameterize the convective turbulence too deeply. However, h KPP is not well-matched to H either, and overall KPP tends to mix too shallowly in the Taylor and Ferrari (2010) and Thomas and Taylor (2010) simulations and too deeply in the Thomas et al. (2013) simulation. 4 In the Taylor and Ferrari (2010) case h KPP is 4 Cases such as the Thomas and Taylor (2010) simulation where h KPP remains similar to the LES boundary layer depth may tempt the modeler to use the KPP algorithm for calculating H in the SI parameterization rather than the criterion on q bulk . However, because SI dynamics are governed largely by the mixed layer PV, the criterion on q bulk remains more physically relevant than one based on a critical Richardson number, and is more appropriate than simply using KPP to determine H . The KPP boundary layer depths, which are calculated using a bulk Richardson number criterion (see Large et al. (1994) ), are also plotted for reference (red circles and crosses). The KPP boundary layer depth in the Taylor and Ferrari (2010) simulation is strongly affected by the presence of the shear instability component, but this sensitivity is not seen in the other two simulations. This may imply that the extra mixing by the shear instability component plays a strong role in the evolution of the flow below z = −10 m at longer integration times. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) strongly controlled by whether the shear instability component is enabled. When it is enabled, the KPP boundary layer depth remains at a constant value of h KPP = 10 m through most of the simulation, implying that mixing below this depth is controlled largely by enhanced diffusivity (see Appendix D.1 ) rather than the KPP boundary layer scheme. When it is disabled, the KPP boundary layer depth deepens with time throughout the course of the simulation, although it is shallower than the LES boundary layer depth at all times. This behavior of the shear instability component is not seen in the other two simulations, possibly due either to the shorter integration time or to the lack of wind forcing in the Taylor and Ferrari (2010) simulations.

Energetics
Of primary concern in boundary layer parameterizations is whether the parameterization faithfully replicates the energetics of the subgridscale turbulence. A main focus of the SI parameterization is to represent the extraction of mean kinetic energy via the GSP, and the reduction of available potential energy by the convective layer buoyancy flux. To facilitate comparison with KPP, it is necessary to define the effective GSP and w b that is parameterized by KPP. A more detailed description of the KPP parameterization can be found in Appendix D.1 , but for now a brief summary will suffice. In KPP, the vertical flux of momentum by KPP is parameterized as where ν KPP is the KPP vertical viscosity. Assuming that the resolved velocity is approximately in geostrophic balance so that ū ū g , the effective GSP introduced by the KPP momentum flux can be expressed as Likewise, assuming the linear equation of state shown in (27) , the KPP buoyancy flux will be defined as where κ and γ are the vertical diffusivity and nonlocal transport terms specific to each tracer. The KPP GSP and buoyancy flux can contribute significantly to the reduction of mean kinetic and available potential energy in the surface boundary layer, and may reduce the energy at rates inconsistent with fully-developed SI. It is anticipated that the SI parameterization, which is designed specifically to reproduce the energy extraction by resolved SI, will outperform KPP in this regard. The figures in this section show comparisons between both parameterizations and the LES; in all plots the vertical profiles from the LES will be indicated by black lines, the SI parameterization by blue lines, and the KPP simulations by red lines. Fig. 5 shows a comparison of the vertical structures of w b and GSP from each simulation, where the profile of w b from KPP is calculated using (53) and the profile from the SI parameterization is calculated as In each simulation the SI parameterization represents the structure of w b well, inducing a positive buoyancy flux in the convective layer via w b con v and a negative buoyancy flux beneath due to κ SI, v . The convective buoyancy flux is absent in the Thomas and Taylor (2010) simulations, where B 0 = 0 and destabilization only occurs through the downfront wind stress. Nonetheless, outside of the near-surface layer above z = −20 m where there is a positive buoyancy flux in the LES, the structure and magnitude of w b are well-approximated by the SI parameterization. The parameterization is also able to approximate the structure of the GSP reasonably well in the Taylor and Ferrari (2010) and Thomas and Taylor (2010) simulations. It does not match the GSP as well in the Thomas et al. (2013) simulation, largely because the profiles of w b and GSP from the LES deviate slightly from being linear near the bottom of the convective layer (for reference, one may also inspect Fig. 9 from Thomas et al. (2013) ).
The performance of KPP generally compares unfavorably to that of the SI parameterization, particularly in terms of its mixing coefficients, ν KPP and κ θ . The large GSP and large negative values of w b KPP are indicative that the KPP diffusivity and viscosity are too strong. The profiles also can be sensitive to whether the shear instability component is enabled. There is a significant difference in the mixing parameters and fluxes in the Taylor and Ferrari (2010) simulation, where the KPP boundary layer mixing coefficients tend to be much weaker when the component is enabled. This implies that the shear instability component can be effective at mixing fluid from the interior up into the KPP boundary layer. The result is a shallower KPP boundary layer depth, and since the mixing coefficients are proportional to the boundary layer depth, these become smaller as well. The shear instability component has less of an effect in the Thomas and Taylor (2010) and Thomas et al. (2013) simulations. This may be because of the significantly shorter run time (2.5 and 2 days, respectively, versus the 15 days for the Taylor and Ferrari, 2010 simulation), or the strength of the surface forcing, which is over an order of magnitude larger than in Taylor and Ferrari (2010) . The correspondingly larger boundary layer mixing may exert more influence in these simulations relative to the shear instability component. Fig. 6 shows comparisons of the time-integrated GSP and w b from the Taylor and Ferrari (2010) LES against the same parameterized quantities from the MITgcm simulations. In the LES the cumulative dissipation of TKE, SI , (green line, right panel) nearly balances the sum of the GSP and w b , reflecting the role these terms play as a bridge leading to the removal of energy from the mean flow. The parameterizations mimic this process by taking the energy from resolved to unresolved scales; that is, because the GSP and w b associated with SI are not directly resolved in the models, the mean energy is removed directly from the resolved flow instead of being converted into TKE first. The total energy removed by the SI parameterization agrees well with the LES in comparison to KPP, which greatly overestimates the energy lost through the GSP.
The GSP acts to reduce the thermal wind shear, which after geostrophic adjustment is expected to reduce the lateral buoyancy gradient, M 2 , and the isopycnal slope, M 2 / N 2 . 5 Therefore, a potential consequence of removing too much energy via the GSP is excessive shallowing of the isopycnal slope; for water being subducted adiabatically along isopycnals, this may reduce net water mass and tracer exchange by the resolved flow between the mixed layer and ocean interior. Furthermore, excessive reduction of M 2 may lead to erroneous transport by other mixed layer parameterizations whose induced fluxes depend on the lateral buoyancy gradient (e.g. Fox-Kemper et al., 2011 ). The depth-integrated GSP and w b are shown in Fig. 7 . At each time the KPP diffusivity and viscosity are too large, resulting in buoyancy fluxes (top panel) and GSP (middle panel) which are larger in magnitude than the LES and the SI parameterization. In the KPP simulation with shear instability (dashed red line), the sum of both terms offset each other so that the total energy loss is fairly consistent with the LES (bottom panel), even though each term individually is a poor match. The KPP simulation without shear instability is dominated by the GSP, and the total energy loss is significantly too large. This simulation also exhibits a large, spurious oscillation at subinertial timescales which we cannot ex-5 In the models presented here a reduction in the average M 2 does not occur because of the frontal zone configuration of the LES, and because of the lateral homogeneity within the averaging region in the MITgcm models. A reduction in the average M 2 does occur in the MITgcm models if the averaging region is widened to include the areas of the domain where M 2 varies laterally. However, inclusion of these regions is less meaningful in the context of comparing against the LES, and so they are neglected in this analysis.  Taylor and Ferrari (2010 , top row), Thomas and Taylor (2010 , middle row), and Thomas et al. (2013 , bottom row). The plots compare profiles extracted from the LES (black lines) against those from the MITgcm when run with the SI parameterization (blue lines) and KPP, which is run with the shear instability component enabled (dashed red lines) and disabled (solid red lines). Results in the plots are time-averaged over two inertial periods and are shown after the same amount of simulation time has elapsed in each model, which is approximately 15 days for the Taylor and Ferrari (2010) simulations, 2.5 days for the Thomas and Taylor (2010) simulations, and 2 days for the Thomas et al. (2013)  plain, but may be linked to the calculation of the boundary later depth ( Fig. 4 , crosses). It is also clear in the plots that the LES experiences interial oscillations, whose effect on the time-integrated plots ( Figs. 5,6 ,(8)(9)(10) is minimized by averaging over two inertial periods.

Vertical structure of the mean fields
To evaluate the ability of the SI parameterization to reproduce b , N 2 , and q in each layer, it is useful to compare their vertical profiles for each test simulation. Figs. 8 through 10 show vertical profiles of the cross-front (meridional) velocity, along-front (zonal) velocity, buoyancy, Ri b , PV, and GSP + w b from the LES and MITgcm simulations. The profiles are taken after the same amount of simulation time has elapsed in each model, which is approximately 15 days for the Taylor and Ferrari (2010) simulations, 2.5 days for the Thomas and Taylor (2010) simulations, and 2 days for the Thomas et al. (2013) simulations.
Both velocity components (top row) show improvement in the SI parameterization relative to KPP. In particular, the magnitude of the cross-front velocity in the SI parameterization is a better match to the LES than KPP is able to achieve, although neither is able to reproduce the vertical structure in the Thomas and Taylor (2010) simulation. This occurs because the cross-front velocity in these simulations is largely Ekman driven, and therefore sensitive to the vertical variations in the viscosity, which are more skillfully reproduced by the SI parameterization. The along-front velocity, which is dominated by the geostrophic shear, is reasonably wellapproximated by both parameterizations.
Neither parameterization is able to capture local variations in the buoyancy and stratification profiles at all depths. The potential vorticity profiles between the two parameterizations tend to be similar through most of the boundary layer. In all cases the PV from the LES has a negative spike at the surface due to negative N 2 . KPP is able to reproduce this feature in the Thomas and Taylor (2010) and Thomas et al. (2013) simulations. The SI parameterization misses this surface PV signature because κ SI, v reaches its full value in the layer immediately below the surface, nearly completely mixing away regions of negative N 2 . This mixing also is associated with N 2 being too weak in comparison to the LES in the Thomas and Taylor (2010) and Thomas et al. (2013) simulations.
Finally, the bottom right panel shows the sum of w b and GSP, representing the total energy extraction from the mean flow. In this figure a critical point is that the sum of w b and the GSP in the LES has an approximately linear vertical profile, which is nearly matched in the SI parameterization models by design. In comparison, the profiles in the KPP simulations do not match the LES as well, reducing the energy of the front slightly too much in the Thomas et al. (2013) simulation and more severely in the Taylor and Ferrari (2010) and Thomas and Taylor (2010) simulations. This is largely due to the magnitude of the KPP viscosity, which induces a GSP that is significantly larger than is observed in the LES ( Fig. 5 ).

Passive tracers
A proper evaluation of the diffusive component of the SI parameterization, K SI , and in particular the lateral dispersion of tracers, is limited by the horizontal homogeneity of the frontal zone models. A full assessment of this component should be completed using more realistic test cases, and is left to future work. For now, a preliminary evaluation of the diffusive component is tested using three additional simulations seeded from the end state of each of the LES. In these tests, the vertical profiles of buoyancy and velocity from the LES are used as the initial conditions for each simulation in the new set, so that each simulation starts from an identical flow that is SI-unstable. A passive tracer ξ is initialized with where L Z is the domain depth (in meters). The tracer concentration is uniform in both horizontal directions, so that only the vertical tracer gradient is nonzero. Therefore, the tracer mixing by the SI parameterization will be affected by κ SI, v and the vertical part of K SI (note that there is no surface flux of ξ , so that there is no nonlocal parameterization of convective turbulence and w ξ con v = 0 ).
It is expected that turbulent mixing in the LES will be intensified near the surface due to convective instabilities, homogenizing the tracer concentration more quickly than at depth. The dependence of κ SI, I on the GSP, which is larger near to the surface and decays with depth, implies that the mixing by the SI parameterization will also be surface intensified. KPP diagnoses the vertical Fig. 8. Vertical structure comparison from simulations with surface buoyancy loss and wind stress matching the LES in Taylor and Ferrari (2010) . Shown here are the along-front velocity, cross-front velocity, buoyancy, Ri b , potential vorticity, and the EKE production terms ( GSP + w b ). All results are collocated in simulation time, and are taken after 15 days. Results shown in the plots are time-averaged over two inertial periods. Black lines: LES, blue lines: MITgcm with SI parameterization, red lines: MITgcm with KPP (with shear instability -dashed line, without shear instability -solid line). The upper dotted line in each plot is the convective layer depth, z = −h, when these diagnostics are measured; the lower dotted line is the SI layer depth, z = −H. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) structure of its mixing coefficients using a third-order polynomial, which also tends to be skewed stronger toward the surface. The vertical structure and magnitude of the mixing are not matched to SI, however, and are expected to mix the tracer locally at rates that are inconsistent with the LES. In these tests all simulations using KPP have the shear instability component enabled.
Each simulation is run for 48 h of simulated time after initialization. The vertical profiles of the tracer concentration after one, two, and four hours are shown in Fig. 11 , and the profiles after 12, 24, and 48 h are shown in Fig. 12 . In these figures, the black lines represent the tracer concentration in the LES, the blue lines represents the full SI parameterization, the red lines represent KPP, and the green lines represent the SI parameterization without the along-isopycnal diffusion active ( K SI = 0 ).
The rate of tracer mixing is matched to the LES more closely by the SI parameterization than by KPP at nearly every depth. In the Taylor and Ferrari (2010) simulation KPP mixes too quickly down to 20 m over the first four hours, until the tracer concentration is nearly completely homogenized by 12 h. The KPP mix-ing depth also becomes too shallow by this time, and this shallow bias persists through the end of the simulation. The full SI parameterization performs reasonably well in matching both the vertical structure and mixing depth at all times. In the Thomas and Taylor (2010) simulation both parameterizations match the LES well out to 24 h, after which both mix too shallowly. The final profiles at 48 h are very similar, though the SI parameterization profile is a slightly better match to the LES. In the Thomas et al. (2013) simulation the SI parameterization performs markedly better than KPP out to four hours, after which both diverge away from the LES. The LES mixing deepens rapidly in this simulation, which neither parameterization represents well. In all plots the SI parameterization without K SI (green lines) mixes far too slowly, suggesting the presence of additional mixing processes on top of what is parameterized by κ SI, v .
Overall the SI parameterization seems to capture the vertical structure of the tracer profiles, but has difficulty reproducing the correct concentrations due to a shallow mixing bias. Future work is needed to understand an appropriate way to rectify this Fig. 9. Vertical structure comparison from simulations with surface buoyancy loss and wind stress matching the LES in Thomas and Taylor (2010) . Shown here are the along-front velocity, cross-front velocity, buoyancy, Ri b , potential vorticity, and the EKE production terms ( GSP + w b ). All results are collocated in simulation time, and are taken after 2.5 days. Results shown in the plots are time-averaged over two inertial periods. Black lines: LES, blue lines: MITgcm with SI parameterization, red lines: MITgcm with KPP (with shear instability -dashed line, without shear instability -solid line). The upper dotted line in each plot is the convective layer depth, z = −h ; the lower dotted line is the SI layer depth, z = −H. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) bias. It is worth mentioning again that the implementation of an along-isopycnal mixing tensor in this parameterization framework is based on the phenomenology of SI, which results in the formation of overturning cells which are nominally tilted along isopycnals. Here the mixing is envisioned to be adiabatic, though more work is needed clarify the precise direction and magnitude of mixing and the contribution of adiabatic and diabatic components.
This simple test is shown here to motivate further evaluation of the along-isopycnal diffusivity component. Further evaluation of all components is necessary, and could include more complicated forcing scenarios such as variable wind direction and amplitude, unbalanced fronts, or realistic domain boundaries. These are beyond the scope of this initial endeavor.

Conclusions
In this paper a framework for a parameterization for unresolved symmetric instability (SI) has been proposed that is based on the known phenomenology of SI in the presence of destabilizing surface forcing. The parameterization is designed to be sensitive to this surface forcing and adjusts the transport of momentum, buoyancy, and tracer accordingly. The vertical structure of this transport is based on the relative strength of the surface and Ekman buoyancy fluxes, and is set to be consistent with the results of SIresolving LES.
The SI parameterization is designed for use primarily in regional studies and models such as ROMS ( Shchepetkin and McWilliams, 2005 ) where submesoscale fronts are resolved. Only a handful of prototype global models attempt to resolve the submesoscale frontal structures and strain fields pertinent for the growth of SI, and only for very limited duration simulations, although this situation is likely to change as computing power and modeling techniques improve. One principle advantage of the way the SI parameterization is constructed is that it is likely to avoid the problem of competing against resolved SI (e.g. Henning and Vallis, 2004 ) -since both the parameterization and resolved SI "shut off" when the PV becomes positive or fq bulk > 0, the combined effect of the parameterization and partially resolved SI will both drive the flow toward an SI-neutral state. Implementing the SI parameterization is therefore likely to help models avoid the issue of  Thomas et al. (2013) . Shown here are the along-front velocity, cross-front velocity, buoyancy, Ri b , potential vorticity, and the EKE production terms ( GSP + w b ). All results are collocated in simulation time, and are taken after 2 days. Results shown in the plots are time-averaged over two inertial periods. Black lines: LES, blue lines: MITgcm with SI parameterization, red lines: MITgcm with KPP (with shear instability -dashed line, without shear instability -solid line). The upper dotted line in each plot is the convective layer depth, z = −h ; the lower dotted line is the SI layer depth, z = −H. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) arrested restratification when not all SI modes are fully resolved (e.g. Bachman and Taylor, 2014 ).
The SI parameterization consists of three separate components. The first is a specification of the vertical mixing coefficients based on a subinertial turbulent Ekman balance. This part of the parameterization extracts energy from the geostrophic flow at a rate that is consistent with the geostrophic shear production, and thus approximates the net effect of resolved SI. The second component represents a convective vertical buoyancy flux, which is isolated to a near-surface convective layer whose depth is determined by a previously proposed polynomial equation. Finally, the third component is a turbulent diffusivity tensor that represents alongisopycnal mixing by SI. A scaling is derived for the along-isopycnal scalar diffusivity, which is then rotated back into z -coordinates by a tensor transformation ( Redi, 1982 ). Scalings for the mixing coefficients in each component are developed assuming steady forcing and subinertial flow.
At the start of Section 3 , several design goals were described. We revisit these below and comment on the performance of the parameterization with respect to each objective.
Appropriate mixing of momentum, buoyancy, and tracers during destabilization by F SI > 0 .
The vertical mixing component of the parameterization has been developed using the subinertial momentum balance (17) and (18) which arises due to destabilizing surface forcing. The dissipation of the mean kinetic energy by the vertical viscosity, ν SI , is expected to reduce the lateral buoyancy gradient via geostrophic adjustment and thereby increase Ri b . A corresponding vertical diffusivity, κ SI, v , is introduced, representing the diabatic mixing effects of small-scale shear instabilities. The magnitude and vertical structure of the mixing depend on the relative contributions of the Ekman and surface buoyancy fluxes to F SI . Finally, along-isopycnal transport by individual SI overturning cells is parameterized using a diffusivity tensor, K SI .

Extraction of energy from the resolved flow by SI.
The parameterization extracts mean kinetic energy via the vertical viscosity, at a rate which matches the GSP. In the convective layer a positive buoyancy flux, w b con v , is induced which lowers the center of mass of the fluid, reducing the available potential energy.

Along-isopycnal dispersion of tracers.
Along-isopycnal tracer mixing is parameterized by introducing the symmetric tensor K SI . The magnitude of the mixing scales with the GSP, reflecting the ability of stronger surface forcing to generate more efficient mixing.

No effect when
The parameterization is sensitive to the surface forcing and is inactive if F SI = EBF + B 0 ≤ 0 . The activation condition for the parameterization is dependent on a bulk PV criterion which determines H . When ∇ h b approaches zero it is likely that h / H → 1, so the parameterization becomes inactive.
Act only in the SI-unstable part of the surface boundary layer. The parameterization defines the surface boundary layer depth H as the shallowest depth where fq bulk > 0 and sets all mixing coefficients and w b con v to zero below z = −H. If the entire fluid column is SI-stable and fq bulk > 0 at all depths then the parameterization will be completely inactive in that column. The dynamics of the parameterization are designed to match those from SIresolving LES ( Thomas et al., 2013;Hamlington et al., 2014;Haney, 2015;Haney et al., 2015 ) and thus are different in the convective and SI-dominated sublayers. The depth of the convective sub-layer h is determined via the solution of (26) , which is shown in Appendix C .

Maintain consistent boundary conditions on momentum, buoyancy, and PV.
The functional form of GSP SI and w χ con v shown in (47) and (48) , respectively, yield zero flux of momentum and tracer through the surface boundary ( u w = 0 and w χ = 0 ). It is also shown in Appendix D that the parameterization avoids double-counting resolved and parameterized variables.
A set of idealized, frontal zone models has been employed to demonstrate that the parameterized vertical fluxes are consistent with eddy-resolving LES. It has been shown that in the idealized test cases presented here, the SI parameterization compares favorably to KPP at matching the LES vertical profiles of momentum, buoyancy, stratification, PV, GSP, and w b . The profiles obtained using KPP are also shown to be sensitive to whether the KPP shear instability component is active, highlighting the complex interactions and feedbacks that can occur even within the same subgridscale turbulence parameterization.
There are some uncertainties in the formulation of the SI parameterization which cannot be conclusively answered with the current LES test cases, and thus are left to future work. Implementing the along-isopycnal mixing tensor K SI appears to capture the behavior of the LES, but it is presently unclear whether a better, yet still parsimonious, representation in the convective layer exists. Furthermore, the passive tracer mixing tests ( Figs. 11 and 12 ) indicate strong entrainment at the base of the SI layer which is not captured by the parameterization. It is also unclear how the parameterization should behave when the surface forcing is destabilizing overall (EBF + B 0 > 0) , but one of its individual components is stabilizing (i.e. downfront winds and surface heating, or upfront winds and surface cooling). Currently it is recommended that the parameterization only be active when both EBF and B 0 > 0, and further research is needed to understand the dynamical behavior of SI in these cases. Finally, further work is needed to clarify the vertical mixing rate in the convective layer. Currently, the vertical diffusivity κ SI, v is determined as a function of the local value of Ri b using (33) , which is applied in both the convective and SI layers. Previous LES results support this scaling in the SI layer, but whether the same scaling applies in the convective layer remains obscure. It may be possible to incorporate other parameterizations into the SI parameterization which are more skillful at simulating the physics of the turbulence in the convective layer.
In addition to the use of along-isopycnal anisotropic diffusivity ( Section 3.3 ) to represent the transport of tracers by symmetric instability, an anisotropic viscosity may be used in place of the purely vertical viscosity in (24) to correctly account for the geostrophic shear production. Rotating the SI eddy viscosity effect into this direction would more closely resemble the parcel switching analysis of SI (e.g. Haine and Marshall, 1998 ). However, in the simulations considered here, where there is little horizontal shear across the front, such a change would have almost no impact on the results. Furthermore, unlike tracers, vertical momentum redistribution can be carried out non-advectively through pressure (e.g. Harcourt, 2015 ). This effect is of primary importance in Langmuir turbulence, which is dynamically similar to symmetric instability in many ways . Therefore, at this time we choose to retain only a vertical viscosity which can effect the energy budget desired, and leave anisotropic viscosity concerns to await a comparison against LES of SI in a more complex flow.
It should also be noted that the turbulent Ekman balance, (17) and (18) , may be modified by the presence of other dynamics within the surface mixed layer. Other sources of turbulent kinetic energy, such as ageostrophic shear production, can play a significant role in the energetics of the boundary layer (e.g. Thomas and Taylor, 2010 ), but are not accounted for by the SI parameterization. Other physical effects, such as the interaction of SI with Langmuir circulations ( Li et al., 2012;Hamlington et al., 2014;Haney et al., 2015 ) or frontogenetic strain (e.g. Thomas, 2012 ), could be taken into account in future improvements of the parameterization presented here. Strong curvature of the front could also modify the dominant balance, although this effect has been neglected here.
Further experiments are ongoing to evaluate the performance of the parameterization against LES in more realistic and diverse settings. A key challenge in performing these model comparisons is the computational expense of running LES -the multiscale nature of SI dynamics requires that the model domain must be large enough to simultaneously resolve the SI modes while resolving down to the O(1) m scale of convective overturning. Of particular interest in these experiments is exploration of the weak-front limit when ∇ h b is small, and in extending the parameterization to include cases where either the EBF or B 0 are negative.
More comprehensive evaluation and testing of the SI parameterization should include simulations where the surface buoyancy flux, wind direction and magnitude, time-dependence of the wind forcing (in both direction and magnitude), and degree of balance associated with the initial front are varied. As in this paper, the effects of the parameterization on the mixed layer stratification and tracer mixing in each of these scenarios should be compared against a series of SI-resolving LES run with the same domain geometry and initial conditions. It would be interesting to see the effects of the parameterization in simulations such as those used in Rosso et al. (2014) and Hamlington et al. (2014) . Future plans include studying the interaction and joint effect of the SI parameterization with other mixed layer parameterizations (e.g. Large et al., 1994;Fox-Kemper et al., 2011 ) in coarse-resolution reproductions of scenarios such as these, as well as in realistic, submesoscale eddy-permitting simulations of the Subantarctic Front.
one may use the parameterization for GSP SI in (21) to write where the lateral buoyancy gradient is assumed to be measured as near to the surface as possible. The buoyancy flux can be expressed in terms of readily observable quantities ( Steele et al., 2009 ), where g is gravity, ρ is the in situ density, α θ = ρ −1 ∂ρ ∂θ is the effective thermal expansion coefficient, β s = ρ −1 ∂ρ ∂s is the effective haline contraction coefficient, s 0 is the surface salinity, c p is the specific heat of seawater, Q 0 is the net surface heat flux (in W m −2 ), and E and P are the rates of evaporation and precipitation (in m s −1 ).
One may also use (21) to estimate the SI horizontal diffusivity.
Recalling the expression for κ SI, h κ SI, I in (40) , one may again substitute for the GSP to obtain where provisionally C = 0 . 02 ( Section 3.3 ). These scales can be compared against the same quantities as estimated for MLE. Assuming wintertime, non-stormy conditions ( E − P 0) , the MLE and SI scales may be compared for different values of Q 0 . The other parameters needed to calculate B 0 in (A.4) will be taken as α θ = 2 × 10 −4 °C −1 , ρ 0 = 1035 kg m −3 , and c p = 4 × 10 3 J kg −1 °C −1 . The approximate scales are shown in Table A

Appendix B. Transformation of coordinates by counter-clockwise rotation
Consider the counterclockwise rotation of an orthogonal (horizontal) coordinate system by an angle θ . Let the original coordinate system be ( x, y ), and the new coordinate system be denoted by a tilde: ( ˜ x , ˜ y ) . The forward and reverse transformation of a vec-

B1. Transforming the geostrophic shear production
The geostrophic shear production (GSP) is given by Assume that the rotation is applied based on ∇ b , so that the rotation matrix is uniform across the window where Reynolds averaging is performed. Then the rotation commutes with the Reynolds averaging, and substitution of (B.3) and (B.4) and (B.7) and (B.8) into (B.10) yields the system u w from which one can solve for ˜ u h w : If the original coordinate frame is aligned with the front so that (B.20) and the Reynolds stresses in the rotated frame are given by (22) .
Since the SI parameterization only uses the two real roots of (C.1) , the sign of ξ = −4 S 2 + 2 p ± q S is sufficient to determine which pair of roots to take. Consider the term q / S , whose sign is the only difference between (C.4) and (C.5) . Because the real roots will result when ξ ≥ 0, if q / S is positive then ξ is increased in (C.4) and decreased in (C.5) , from which one can conclude (C.4) gives the real roots. Conversely, if q / S is negative then ξ is increased in (C.5) and decreased in (C.4) , and (C.5) must thereby give the real roots.
Begin by considering q and S separately. The expression for q is guaranteed to be positive for α > 0. S is real or imaginary depending on the expression Q − 12 α/Q. Since Q > 0, to show Q − 12 α/Q > 0 is equivalent to showing Q 2 > 12 α. To do so, note plotted as a function of α, which ranges here from 10 −10 to 10 10 .
Then Q 2 > 1728 1 / 6 α 1 / 2 2 = 12 α, (C.11) which implies that S ≥ 0 as well. Then because both q and S are positive the real roots of (C.1) are given by (C.4) . Furthermore, noting that the first two terms in all four roots (C.4) and (C.5) are negative, the only possible positive root (i.e. the one that must be used in the parameterization) is Rigorously establishing bounds on this root is difficult and beyond the scope of this paper. For the purposes of the parameterization, however, it suffices to show that 0 ≤ x 1 ≤ 1 for all choices of α, which confirms that h diagnosed with this formula is guaranteed to be positive and to not exceed H . This is important for numerical reasons, both so that the parameterization activates the correct mixing scheme within each sublayer and so that calculating x 1 does not require extra "fail-safes" against unexpected behavior. Fig. C.13 shows x 1 evaluated as a function of α, where α ranges from 10 −10 to 10 10 . The evidence from the numerical solution ( Fig. C.13 ) suggests that x 1 is bounded between 0 and 1 for all possible values of α.
In practice, the value of x 1 can be calculated simply by using the algebraic formula in (C.12) . The authors have noted that numerical solvers may have trouble with this procedure when α is either very large or very small, and may be unable to find a solution. Therefore, it may suffice to simply set x 1 = 1 − α −1 / 3 when α > 10 6 and x 1 = α 1 / 4 when α < 10 −6 , where a numerical examination has shown the asymptotic limits to be excellent approximations (less than 1% relative error).

Appendix D. Prevention of double-counting
The transport of momentum and buoyancy across density surfaces can play a major role in setting the shear and stratification in the mixed layer, which in turn can affect air-sea exchange and the mixing and subduction of tracers. Subgridscale parameterizations run the risk of interfering with a model's resolved transport characteristics, and must therefore be carefully constructed so as not to incur problems related to "double-counting" (e.g. Henning and Vallis, 2004;Delworth et al., 2012 ). In the case of the SI parameterization, it is important to ensure that the parameterization does not act as a spurious source of momentum or buoyancy.
SI acts to redistribute momentum in the vertical direction, but does not accelerate or decelerate the depth-integrated flow. Appropriate boundary conditions on the SI parameterization must be enforced to reflect this. To see this, consider the Reynolds-averaged momentum equations, (11) and (12) , with the SI parameterization included. These will be written here as where for notational simplicity all other frictional and subgridscale forcing terms in the momentum equation are wrapped into a new variable F . Removing the balanced part of the flow and assuming that the turbulent stresses at the bottom of the SI-unstable layer go to zero, integrating in the vertical yields (D.2) In order for the SI parameterization to not act as a spurious source of momentum requires that the parameterized Reynolds stress term vanish at the surface. This is achieved automatically in (21) since GSP SI = 0 at z = 0 , and is consistent with the kinematic boundary condition that w = 0 at the surface. Another point to consider in the construction of the SI parameterization is its effect on the vertical budget of an arbitrary tracer, For the SI parameterization to not act as a spurious source of tracer requires w ξ SI to vanish at the surface. This condition is satisfied for κ SI, v and K SI since both are proportional to GSP SI , which is zero at the surface. The convective layer flux w ξ con v vanishes according to (48) , so that all three terms comprising w ξ SI vanish at the surface. Note that if one sets ξ = b, this derivation also shows that the SI parameterization does not affect the mixed layer buoyancy budget. It has thus been shown that the SI parameterization satisfies the requirement that it should not act as a net source of momentum or buoyancy, and avoids the problem of "double-counting" resolved and parameterized variables.

D1. The K-profile parameterization (KPP)
In hydrostatic models parameterizations such as the popular KPP scheme ( Large et al., 1994 ) are needed to represent vertical fluxes driven by convectively unstable and wind shear conditions. Here a brief overview of KPP is warranted in order to motivate comparison with the SI parameterization. In models where KPP is enabled, the vertical fluxes for a generic tracer χ are parameterized as w χ KPP = −κ χ ∂ χ ∂z − γ χ (D.6) for z > −h KPP , where h KPP is the surface boundary layer depth diagnosed by KPP, κ χ is a depth-dependent vertical diffusivity, and γ χ is a nonlocal transport term which is nonzero only in convectively unstable conditions. 6 The generic form of the KPP diffusivity coefficients is where the vertical dependences of the velocity scales, w * , and nondimensional structure functions, G , are designed to smoothly taper the mixing coefficients to zero at the surface and to match the interior mixing coefficients at the bottom of the planetary boundary layer. Note that the velocity scale is unique to each scalar, while the structure function is not. The nonlocal terms have the generic form where C s is a nondimensional scale factor and w χ 0 indicates the flux of χ through the ocean surface. Note that the respective forms of κ χ and γ χ allow (D.6) to be rewritten in the form By the linearity of (D.10) , it follows that w b = g α θ w θ − β s w s , (D.11) so that, given the general form for each vertical flux in (D.6) , the KPP buoyancy flux can be written where w b 0 = B 0 as defined in Section 2 . Likewise, the vertical flux of momentum by KPP is parameterized as (D.13) with no contribution from a nonlocal term. In KPP the viscous coefficient ν KPP has a vertical structure that can be distinct from those of κ s and κ θ , and is tapered to match the value and first derivative of the model viscosity for the ocean interior. Assuming that the resolved velocity is approximately in geostrophic balance so that ū ū g , the effective GSP introduced by the KPP momentum flux can be expressed as (D.14) The KPP GSP and buoyancy flux can contribute significantly to the reduction of mean kinetic and available potential energy in the surface boundary layer, and may reduce the energy at rates inconsistent with fully-developed SI.
Finally, KPP includes an optional parameterization for unresolved mixing due to shear instabilities, which becomes active when the local gradient Richardson number where Ri 0 is a critical value typically defined to be 0.7. When this condition is met, the vertical diffusivity for all scalars (including momentum) is increased by an amount dependent on the value of where the value κ 0 = 5 × 10 −3 m 2 s −1 is chosen to fall within the range of maximum observed diffusivities reported for the seasonal thermocline ( Peters et al., 1988 ). This additional diffusivity is calculated where necessary before the KPP boundary layer scheme is called, and is overridden by the boundary layer diffusivities for z > −h KPP . The inclusion or exclusion of this shear instability component can have a significant effect on the mixed layer momentum and tracer profiles over time. Results from simulations using KPP with and without the shear instability component are shown in Section 4 to demonstrate these effects.