Mixed layer sub-mesoscale parameterization – Part 1: Derivation and assessment

. Several studies have shown that sub-mesoscales (SM ∼ 1km horizontal scale) play an important role in mixed layer dynamics. In particular, high resolution simulations have shown that in the case of strong down-front wind, the re-stratiﬁcation induced by the SM is of the same order of the de-stratiﬁcation induced by small scale turbulence, as well as of that induced by the Ekman velocity. These studies have further concluded that it has become necessary to include SM in ocean global circulation models (OGCMs), especially those used in climate studies. The goal of our work is to derive and assess an analytic parameterization of the vertical tracer ﬂux under baroclinic instabilities and wind of arbitrary directions and strength. To achieve this goal, we have divided the problem into two parts: tracer for ocean codes mesoscales, dynamic equations including the non-linear terms for which we employ developed and assessed in previous work. We present a detailed analysis for down-front and up-front winds with the following results:


Abstract.
Several studies have shown that sub-mesoscales (SM ∼1km horizontal scale) play an important role in mixed layer dynamics. In particular, high resolution simulations have shown that in the case of strong down-front wind, the re-stratification induced by the SM is of the same order of the de-stratification induced by small scale turbulence, as well as of that induced by the Ekman velocity. These studies have further concluded that it has become necessary to include SM in ocean global circulation models (OGCMs), especially those used in climate studies.
The goal of our work is to derive and assess an analytic parameterization of the vertical tracer flux under baroclinic instabilities and wind of arbitrary directions and strength. To achieve this goal, we have divided the problem into two parts: first, in this work we derive and assess a parameterization of the SM vertical flux of an arbitrary tracer for ocean codes that resolve mesoscales, M, but not sub-mesoscales, SM. In Part 2, presented elsewhere, we have used the results of this work to derive a parameterization of SM fluxes for ocean codes that do not resolve either M or SM.
To carry out the first part of our work, we solve the SM dynamic equations including the non-linear terms for which we employ a closure developed and assessed in previous work. We present a detailed analysis for down-front and up-front winds with the following results: (a) down-front wind (blowing in the direction of the surface geostrophic velocity) is the most favorable condition for generating vigorous SM eddies; the de-stratifying effect of the mean flow and re-stratifying effect of SM almost cancel each other out, Correspondence to: V. M. Canuto (vcanuto@giss.nasa.gov) (b) in the up-front wind case (blowing in the direction opposite to the surface geostrophic velocity), strong winds prevents the SM generation while weak winds hinder the process but the eddies amplify the re-stratifying effect of the mean velocity, (c) wind orthogonal to the geostrophic velocity. In this case, which was not considered in numerical simulations, we show that when the wind direction coincides with that of the horizontal buoyancy gradient, SM eddies are generated and their re-stratifying effect partly cancels the de-stratifying effect of the mean velocity. The case when wind direction is opposite to that of the horizontal buoyancy gradient, is analogous to the case of up-front winds.
In conclusion, the new multifaceted implications on the mixed layer stratification caused by the interplay of both strength and directions of the wind in relation to the buoyancy gradient disclosed by high resolution simulations have been reproduced by the present model.
The present results can be used in OGCMs that resolve M but not SM.

Introduction
Recently, there has been a considerable interest in submesoscales (SM) which are oceanic structures with sizes O(1 km) and a life time of the order of days. If one considers that the highest resolution O(1/10 0 ) in stand-alone OGCMs (ocean global circulation models) can represent structures of about 10 km which is 10 times larger than SM sizes and that OGCMs employed in thousand years runs for climate studies have in general a 1 0 resolution (corresponding to structures 100 times larger than SM sizes), it seems clear that a good deal of important physical processes have thus far gone unrepresented in many OGCMs.
The parameterization of SM cannot be constructed by analogy with mesoscales, M. Indeed, while mesoscales are 680 V. M. Canuto and M. S. Dubovikov: Mixed layer sub-mesoscale parameterization -Part 1 characterized by a small Rossby number Ro=ζ /f 1 (where ζ and f are the relative and planetary vorticities respectively), mixed layer SM are characterized by Ro∼1 (Thomas et al., 2008). Thus, the dynamics of M and SM are quite different and so are the final results for the fluxes of interest. For this reason, the parameterizations for M that were suggested for the deep ocean cannot be extrapolated to describe SM in the mixed layer and a new parameterization is required.
Most of the present knowledge about mixed layer SM comes from high resolution numerical simulations (Levy et al., 2001;Thomas and Lee, 2005;Mahadevan, 2006;Mahadevan and Tandon, 2006;Klein et al., 2008;Thomas et al., 2008;Capet et al., 2008, C8;Levy et al., 2010;Mahadevan et al., 2010, MTF). These studies have revealed many interesting features of SM, e.g., their contribution to the vertical mixing of buoyancy and tracers in the upper ocean. Among the most salient effects of SM on global ocean properties is the well documented tendency to re-stratify the mixed layer (Spall, 1995;Nurser and Zhang, 2000;C8;MTF). The effect of SM on deep convection has been recently demonstrated by Levy et al. (2009, Fig. 9) who point out a better agreement with the new mixed layer data by Boyer Montegut et al. (2004) south of the WBC (Western Boundary Current). Even in non-convective regimes, one can expect a significant cancellation between SM and small scale fluxes leading to the mixed layer re-stratification (e.g., C8, Fig. 12; Klein et al., 2008); Hosegood et al. (2008) estimated that SM contribute up to 40% of the re-stratification process. In addition, as noted by Lapeyre et al. (2006) and Klein et al. (2008), the surface layers re-stratification is compensated by de-stratification of the ocean interior pointing to an interesting dynamical connection between surface and interior processes. Another important effect of SM concerns the location of the WBC that is shifted south by 4 0 and whose off shore extension penetrates further to the east, in better agreement with observations (Levy et al., 2009). An earlier study by Treguier et al. (2005), who found a significant increase (from ∼30 Sv to ∼70 Sv) in the barotropic transport in the Gulf Stream when moving from 1 0 to 1/6 0 resolution, was recently confirmed and the inclusion of SM further increased the transport by ∼50 Sv. Finally, the structure of the MOC (meridional overturning circulation) was also significantly affected by SM not so much in its intensity as to its location (Levy et al., 2009, Fig. 12).
This brief summary of some of the results of very high resolution (1/54 0 , ∼2 km) regional studies highlights the importance of SM and thus the question arises as to how much of SM physics is actually accounted for by global OGCMs. Even today's highest resolution global ocean models ∼1/10 0 (Maltrud and McClean, 2005;Sasaki et al., 2008) are not able to capture the SM field and much less are in the position to do so the OGCMs coupled to an atmospheric model used in climate studies where the resolution is O(1 0 ) but generally lower. The significant global processes revealed in going from 1 0 resolution (∼100 km) and 1/54 0 resolution (∼2 km) are presently absent in such global models especially in climate studies. Therefore, a reliable parameterization of SM in terms of the resolved fields has become necessary to ensure the physical completeness of mixed layer mixing processes.
As MTF stressed, "parameterization of the circulation induced by SM eddies in the presence of wind forcing is required in climate models in order to simulate the restratification correctly". To accomplish such a goal, a parameterization: (1) must be valid for an arbitrary tracer since a model for buoyancy only is insufficient because it cannot describe an important ingredient such as CO 2 , (2) must include a wind of arbitrary direction and intensity since MTF have shown that its effect on SM fluxes is large and because forcing in future climates is likely to be quite different from today's, (3) must reproduce simulation data and finally, (4) to be usable in climate codes, it must be expressed only in terms of resolved fields, that is, the parameterization must be averaged over mesoscale fields. No parameterization presently available satisfies these criteria.
We begin by considering the model independent dynamic equation for an arbitrary mean tracer 1 : where the SM horizontal and vertical tracer fluxes are defined as follows: In Eqs. (1a,b), an overbar indicates averages over submesoscales. To derive a parameterization for OGCMs that do not resolve either SM or M, one must further average the present results over the mesoscale fields, a problem we have studied elsewhere . The method we employ to derive the SM parameterization is analytical and thus it can be followed and checked in detail. The final result for the vertical SM flux for an arbitrary tracer under arbitrary buoyancy and wind conditions, is also expressed analytically. The model includes non-linear interactions. To obtain the desired results, one carries out three steps: first, one solves in Fourier space the SM dynamic equations describing the u , w , τ fields; second, one constructs the second-order correlation function w τ in kspace and third, one integrates the results over all wave vectors to obtain the fluxes in physical space in terms of the resolved fields, to be used in Eq. (1a). The procedure was first worked out for the linear case by Killworth (1997) and for the non-linear case by Dubovikov (2005, 2006, CD5, 6). Though the dynamic equations describing the SM velocity and temperature fields are formally the same as those describing mixed layer mesoscales that were discussed in , in the present case they must be solved in the regime appropriate to SM, namely for a Rossby number Ro=ζ /f =O(1) rather than Ro 1, as in the case of mesoscales. In addition, SM are trapped in the ML while mesoscales extend throughout the entire water column and form coherent structures (Provenzale, 1999).
The key difficulty in solving the SM dynamic equations is the presence of non-linear terms whose closure is expressed in Eq. (3a) below. Since the latter is a key ingredient of the present model and since the original derivation (Canuto and Dubovikov, 1997) is somewhat involved, in Appendices A and B we have attempted to find a way to present a more physical approach to Eq. (3a) with the goal of highlighting the physical rather than the technical features of Eq. (3a).
In addition to the derivation of the closure relations (Eq. 3a), there is the issue of the assessment of Eq. (3a) when applied to flows different than the present one so as to justify their use in the present context. Such an assessment was made using data from freely decaying flows, 2-D flows, rotating flows, unstably stratified flows, shear driven flows, DNS data, etc. and the results were in good agreement with the data (see Canuto et al., 1999 and references therein). Even so, we consider such an assessment necessary but not sufficient for the credibility of the parameterization of SM fluxes derived below. The additional requirement consists in assessing the model predictions against results from SM resolving simulations. A first simulation corresponds to a system forced only by baroclinic instabilities and no wind (Fox-Kemper et al., 2008, FFH) while a second one corresponds to a flow under realistic wind and buoyancy forcing (C8, 0.75 km resolution; MTF, 1 km resolution). They will be discussed in Sects. 5-7.
The following two conditions must be further satisfied by an SM parameterization: (a) it must reproduce existing data and (b) it must predict new features to be assessed when such data become available. In this context, it must be mentioned that our work was posted as an OS Discussions (Canuto and Dubovikov, 2009, CD9) before Dr. A. Mahadevan kindly sent us the MTF manuscript. Our model would have been falsified had its predictions turned out to be inconsistent with MTF data. However, the model predictions in CD9 not only did not contradict the simulation data, but called attention to the same qualitative SM effects as MTF did in their paper.
To make the SM parameterization usable in OGCMs, we looked for analytical solutions of the SM dynamic equations and to achieve that goal, we introduced the assumption that the fluxes are mostly contributed by their spectra in the vicinity of their maxima. Though this introduces errors of several tens of a percent, the advantages of obtaining analytic expressions for the vertical tracer flux in terms of resolved fields in the presence of both frontogenesis and Ekman pumping, was worth exploring. Another approximation which has helped us obtain analytical results follows from the assumption that the SM kinetic energy K SM exceeds K= u 2 /2 where u is the baroclinic component of the mean velocity (we call attention to the fact that K is considerably smaller than the mean kinetic energy), that is, the condition of applicability of the present treatment is predicated on the assumption: which we shall check several times in the following. Furthermore, following Killworth (2005), we adopt the approximation that due to the mixed layers strong mixing, one can neglect τ z in the SM equations. Anticipating our main result, the vertical flux that enters Eq. (1a) will be shown to have the following form: where u + S plays the role of a bolus velocity. Since τ z is small, one may make the analogy with the mesoscale bolus velocity more complete by adding to Eq. (1d) the term w + S τ z , where w + S is found from the continuity condition ∂ z w + S +∇ H ·u + S =0 (Killworth, 2005).
The organization of the paper is as follows. In sec.2 we discuss the dynamic equations for the SM fields in the ML and apply the turbulence closure model to the non-linear terms; in Sect. 3 we present the form of ∂ z F V and F V which we derive in Appendix C; in Sect. 4 we derive the explicit form of the SM kinetic energy K SM in terms of the resolved fields that, together with the results of the previous section, completes the problem of expressing ∂ z F V in terms of resolved fields in the presence of both frontogenesis and Ekman pumping. In Sect. 5 we study the case of a strong wind when the Ekman velocity exceeds the geostrophic one. We shall show that when a strong wind blows in the direction of the geostrophic velocity or of ∇ H b, it tends to de-stratify the mixed layer but at the same time it generates SM that tend to re-stratify the mixed layer. On the other hand, when the wind blows in directions opposite to geostrophic velocity or to ∇ H b, it re-stratifies the mixed layer, an effect that is strengthened by the re-stratifying effect of SM. In Sect. 6 we compare the model results with the data from the SM resolving simulations of Capet et al. (2008). In Sect. 7 we compare the model results for the no-wind case. In Sect. 8, we present some conclusions.

Sub-mesoscales dynamic equations near the surface
Consider an arbitrary tracer field τ and separate it into mean and fluctuating partsτ =τ +τ . The dynamical equation for the SM tracer field τ is obtained by subtracting the equation for the mean tracer τ from that of the total field τ . Since this procedure is well known and entails only algebraic steps with no physical assumptions, we cite only the final result (the notation is explained in footnote 1): must be noted that in Eqs. (2a) no closure has been used for the non-linear terms. Without the non-linear terms, Eqs. (2a) formally coincide with those describing mixed layer mesoscales tracer fields studied by Killworth (2005). The difference in representing M and SM lies in the scales over which the averages (represented by an overbar in Eqs. 2a), is taken: in the case of mesoscales, averages are over scales exceeding mesoscales while in the case of submesoscales, averages are meant to be over scales smaller than mesoscales but larger than submesoscales. Furthermore, in describing mesoscales, one has Ro 1 and Ri 1, whereas in the case of sub-mesoscales, both Ro, Ri∼O(1). Following Killworth (2005), we neglect the terms containing τ z and τ z in which case the first of Eqs. (2a) simplifies to: Without the non-linear term, Eq. (2b) coincides with Eq. (2) of Killworth (2005) for the mesoscale buoyancy field. Within the same approximation, the equation for the horizontal SM velocity is given by: where e z is the unit vector along z axis. Next, we Fourier transform Eqs. (2b,c) in horizontal planes and time. Following Killworth (1997Killworth ( , 2005, we keep the same notation u , τ for the submesoscale fields in the k−ω space and assume that when Eqs. (2b,c) are Fourier transformed, the mean fields u and ∇ H τ are constant in time and horizontal coordinates. We thus obtain: where we have added the continuity equation that provides the z-derivative of w . We recall that τ , u and the nonlinear terms are functions of the horizontal wave vector and frequency (k, ω) and z, whileū is a function of z only and ∇ H τ is z independent. Equations (2e) form a closed system whose solution provides the necessary ingredients to construct the vertical flux (Eq. 1b) provided one has a closure for the non-linear terms, a problem discussed in Appendices A and B with the result that, in the vicinity of k = k 0 where the SM energy spectrum E(k) has its maximum, the non-linear terms Q H have the following forms: where the scale =k −1 0 may be interpreted as the SM horizontal length scale. As it was shown in detail in CD5, k 0 is obtained from the solution of the eigenvalue problem which is derived from the eddy dynamical equations (Eq. 2e). In the limit of a strong non-linearity represented by: where K is the kinetic energy of the baroclinic component of the mean velocity (defined in Eq. 4b), the solution of the eigenvalue problem yields the following result: where r S is the Rossby deformation radius of the mixed layer (ML) of depth h and where N is the buoyancy frequency in the ML. Relation (Eq. 3c) is in qualitative agreement with other evaluations of the submesoscale length scale discussed in the literature (Boccaletti et al., 2007;Thomas et al., 2008;Fox-Kemper and Ferrari, 2008). In Fig. 9 of Fox-Kemper et al. (2008), the authors, using simulation data, plot SM length scales defined using different variables and in unit of the Stone length scale which is of the order of the deformation radius. Equation (2e) together with Eq. (3a), represent a stochastic Langevin equation which has played a major role in turbulence modeling studies (Kraichnan, 1971;Leith, 1971;Herring and Kraichnan, 1971;Chasnov, 1991). The advantage of the Langevin equation is that it is linear in the fluctuating fields and thus allows one to compute second-order moments while the original Eqs. (2b,c) are non-linear and do not allow an analytical computation of such correlation functions. The key problem is to find a model for the non-linear terms Q's that leads to a Langevin equation whose correlation functions are sufficiently close to those of the original Eqs. (2b,c). This is the closure problem for the nonlinear terms. In CD5, we used the closure (Eq. 3a) derived by Canuto and Dubovikov (1997) and solved the eigenvalue problem to which the mesoscale dynamic equations were shown to reduce. Closure (Eq. 3a) has a simple interpretation within the mixing length approach. In fact, the first two relations are quite standard with χ −1 being the characteristic time scale while the third relation containing the characteristic length scale (Eq. 3c) and velocity, is the only possible combination that leads to a time scale.
The advantage of the Langevin equation is that it allows us to express all SM fields in terms of the SM horizontal velocity u (k,ω) which, in turn, allows us to express the spectrum of any second-order moment in terms of the SM energy spectrum and of the resolved fields. Such a program, which we discuss in detail in Appendix C, was previously used by the authors to parameterize mesoscales in the ocean interior (CD5, 6) and in the mixed layer  shows that, using the third of Eq. (C6), in the limit (Eq. 1c), u R can be represented as follows: where Ro is the Rossby number and p is the SM pressure field. Thus, when treating mesoscales that are characterized by small Ro, the first relation in Eq. (3d) reduces to the geostrophic relation u R →u g =−ikρ −1 p . On other hand, since SM are characterized by Ro∼1, we must use the complete form of Eq. (C12). This is the reason why we have not called u R the geostrophic component and the divergent component u D (Eq. C11) the a-geostrophic one.

Sub-mesoscale vertical tracer flux
Following the program we have outlined at the end of the previous section, in Appendix C we show how to express w (k, ω) and τ (k, ω) in terms of the SM horizontal velocity u (k, ω) and of the resolved fields; that, in turn, allows us to express the spectrum of F V ≡ w τ in terms of the SM energy spectrum. Finally, integrating the spectra over all wavenumbers, we obtain the following SM vertical tracer flux where: e z is the unit vertical vector and Ro is the Rossby number defined in Eq. (3d). It is worth stressing that in the second relation in Eq. (4a), the second term in the square bracket is a vector: in fact, although e z × u is a pseudo-vector (cross product of the vectors e z and u), f is a pseudo-scalar which is the scalar product of the vector e z and the pseudo-vector 2 and thus the product is a vector. The variable u may be interpreted as the ML baroclinic part of the mean velocity. The parameterization (Eq. 4a,b) is obtained under condition (Eq. 1c) and can be obtained from Eq. (7a,b) of CD9 in the limit (Eq. 1c). In Eq. (4a) the velocity u + S may be interpreted as the submesoscale induced velocity which is a counterpart of the mesoscale induced velocity. As noted earlier, to make the analogy with the mesoscale induced velocity more complete and since in the mixed layer τ z is small due to the strong mixing, one may add to Eq. (4a) the term w + S τ z , where w + S is found from the continuity condition: The only variable in Eq. (4b) that is not yet parameterized is the SM kinetic energy K SM which we study in the next section. Before doing so and for future reference, we next derive the explicit form of the vertical flux itself. Integrating Eq. (4a) over z with the boundary condition F V (0)=0, we account for only the z-dependency of u within the mixed layer. Thus, we obtain: where the submesoscale diffusivity is given by: From Eq. (6b), one observes that at the bottom of the ML, z=−h, we have that: which is a good approximation since SM eddies hardly penetrate the bottom of the mixed layer (Boccaletti et al., 2007). We also have the additional relations:

MS kinetic energy in terms of resolved fields
Assuming that the production of K SM occurs at scales ∼r s and since the eddy kinetic energy equation shows that the vertical buoyancy flux F b V acts as the source of K SM , we employ the following relations: In the case of 3-D turbulence where kinetic energy cascades from large to small scales and a Kolmogorov spectrum sets in, the first relation in Eq. (7a) is simply the statement that production=dissipation with the former is defined in the second of Eq. (7a) while dissipation is represented by a Kolmogorov form. In such a case, C=(3/2)Ko(1+A * ) where A * >0 accounts for the contribution to K SM of the energy spectrum at k<k 0 and Ko is the Kolmogorov constant whose value ranges between 1.4-2.2. With A * =1, Ko=2, we have C=6. However, in the 2-D case of interest here, cascade of K SM to smaller scales does not take place. Instead, K SM transforms into SM potential energy which we denote by W SM . Then, in the quasi-stationary case, the production of K SM given in the second of Eq. (7a), approximately equals the production of W SM . Since the latter cascades to smaller scales where is ultimately dissipated, we have: where C depends on the spectrum of W SM . If we further denote by SM the ratio of the SM kinetic and potential energies, from Eq. (7b) we conclude that the first relation of Eq. (7a) is satisfied with: It is worth noticing that relations analogous to Eq. (7a,b) hold true for mesoscales as well with the proviso that the constants C , and therefore C, are different for SM and M due to their different dynamics since, as already mentioned, SM are trapped in the ML while M form coherent structures throughout the entire water column and thus entail the dynamics of the deep ocean. Using a mesoscale resolving simulation, Eq. (7a) was validated in more than 70 different mesoscale resolving simulations . In the next sections, on the basis of Capet et al. (2008), we estimate that C≈6. Even though at present we have determined C on the basis of only one simulation by C8, we shall show below that the variable of interest to OGCMs, the tracer flux, is only weakly dependent on C. Substituting Eq. (6a,b) and r S =N h/π |f | into Eq. (7d), we obtain the following algebraic equation for K SM : which is however not convenient for the computation of K SM since the latter enters into the right hand side of Eq. (7d) through the variables γ , as one observes from Eq. (4b). From Eqs. (7d) and (4b), one derives the following equation for γ : where: the vector s being the slope of the isopycnals. Equation (7e) is valid under the condition which is the same as condition (Eq. 1c) expressed in terms of the resolved fields. We assess this condition in detail in both strong and weak wind cases.
In summary, for OGCMs that resolve M but not SM, the parameterization of the z-derivative of the vertical SM tracer flux is given by Eqs. (4a, b) and (7e-f).
To illustrate the solutions we have just derived, in the next sections we consider three important cases: (1) strong wind driven flows, (2) wind and buoyancy driven flows, (3) buoyancy-only driven flows.

Wind driven flows
In this section we study flows driven by strong winds when the Ekman velocity exceeds the geostrophic mean velocity and compare with results obtained in the submesoscale resolving simulations of Capet et al. (2008). To obtain results in an analytical form, we further assume that the ML turbulent viscosity ν∼10 −2 m 2 s −1 is z-independent. Under these conditions, the mean velocity field can be decomposed into geostrophic u g and Ekman u E components; with the x axis along the wind direction, we have the relations: where ρu 2 * is the surface stress and δ E is the Ekman layer's depth. Below we analyze flows driven by winds of different directions with respect to the geostrophic component of the mean velocity.

Down-front winds
Down-front winds (blowing in the direction of the surface geostrophic velocity) drive dense water over buoyant one and provide favorable conditions for the generation of submesoscales (Thomas, 2005;Thomas and Lee, 2005). We have: which corresponds to an horizontal buoyancy gradient given by: It is worth noticing that S g is a scalar since e y =e z ×e x is the cross product of the two vectors e z (the unit vector along the Earth radius) and e x (the unit vector along the wind direction) which is a pseudo-vector. To obtain the submesoscale flux and its z-derivative, we need to compute u and u defined in Eqs. (4b), (6b), where u=u E +u g . Assuming that the mixed layer thickness h is considerably larger than the depth of the Ekman layer so that the Ekman number E=δ 2 E h −2 1 and e z 0 /δ E 1, from Eq. (8a,b) we derive that: A strong wind (larger than the geostrophic wind) is represented by the relation: Let us now study the z-derivative of the submesoscale buoyancy flux which we obtain substituting Eq. (8c,e) into Eq. (4a) with τ =b. The result is: As one can check, the z-derivative of this expression is negative which implies that the SM vertical buoyancy flux restratifies the ML. Let us compare the latter effect with that of a down-front mean wind that de-stratifies the ML since "Ekman flow advects dense water over light" (Thomas et al., 2008). In the approximation of strong mixing, in the z-derivative of the mean advection term u·∇ H b, the largest contribution comes from the baroclinic term u·∇ H b (in fact, when we differentiate u·∇ H b w.r.t. z, the mean velocity u=<u>+ u may be substituted by the baroclinic component u since the z-derivative of ∇ H b may be neglected). From Eq. (8c,e) we have: which we compare with Eq. (9a) in the case of a very strong wind: Ro 1, i.e. γ 1 (9c) Under this condition, we may approximate Eq. (9a) as follows: The implication of this relation is that the re-stratification by SM largely compensates the de-stratifying effect of the mean flow, a conclusion in agreement with Mahadevan et al. (2010). To express the eddy kinetic energy in terms of the resolved fields, we first find V from Eq. (7d): Substituting Eq. (10a), together with Eq. (8c), into Eq. (7d), we obtain: We recall that in real flows the mixed layer depth exceeds the Ekman one and thus in Eqs. (10a) we have E<1. Inspection of relations Eqs. (10a,b) shows that both V x,y contribute with the same sign to the SM kinetic energy. In addition, under condition (Eq. 8b) that S g >0, these contributions are positive. Thus, we conclude that the down-front wind generates the most vigorous SM eddies, in agreement with the results of Thomas (2005), Thomas and Lee (2005), Thomas et al. (2008) and Mahadevan et al. (2010) that down-fronts winds provide the most favorable condition for SM generation.
As discussed in the previous section, the variable γ is solution of Eq. (7e) in which A 3,4 are given by (making use of Eqs. 7f and 10a): In the case of a strong down-front wind satisfying condition (Eq. 9c) and with E 1, we have |A 3 |>|A 4 | and the first term in Eq. (7c) may be neglected. Then, we obtain: Next, we discuss condition (Eq. 7g). Using Eq. (8e) for the case of a strong down-front wind, we obtain: Substituting this result into Eq. (7g) together with Eq. (10d), we obtain: which is amply satisfied in the simulations of .

Up-front wind
Up-front wind (blowing in the direction opposite to that of the surface geostrophic velocity). In this case, the contribution of the baroclinic component of the mean flow to Eq. (1a) is given by Eq. (9b) and has the sign opposite to that in the previous case since now S g <0, i.e., in this case u leads to a re-stratification of the ML. It is worth treating strong winds first. In this case, the first term of V x in Eq. (10a) dominates, we have V x >0, V y <0 and thus the expression in the square bracket of Eq. (10b) is negative. Using the relation (1+γ 2 ) −1 = K SM /(K SM + r 2 S f 2 ), we conclude that the only solution of Eq. (10b) is K SM =0 which implies that submesoscales are not generated.
On the other hand, in the case of weak winds when submesoscales are generated, they re-stratify the ML. Thus, the re-stratification effect of the mean flow is strengthened by that due to SM, a conclusion in agreement with the results of

Wind perpendicular to the geostrophic velocity
In such flows we have: v g = v 0 + S g z, u g = 0 (11a) and thus: In contrast with Eqs. (8d-e) and (10a), in this case, the geostrophic component contributes to the y-projections of <u>, u, and V. In particular, we have Substituting this relation into Eq. (7d), we obtain: Inspection of relations Eq. (11c,d) shows that the components V x,y contribute to K SM with the opposite sign. Still, in the case of a strong wind corresponding to a small γ , the first term in Eq. (11d) dominates which yields a positive K SM if S g >0. In this case, the direction of the wind coincides with that of the horizontal buoyancy gradient (up-gradient wind). Under the same condition S g >0, the baroclinic component of u de-stratifies the ML since in this case instead of Eq. (9b), we have: which has a positive z-derivative and thus de-stratifies the mixed layer, while SM tend to re-stratify the mixed layer. As one can see from Eqs. (9b) and (11e), in both cases (Sect. 5.1) and (Sect. 5.3) "Ekman flow advects dense water over light" (Thomas et al., 2008). On the contrary, if ∇ H b has a direction opposite to that of the wind (down-gradient wind), so that S g <0, strong winds tend to re-stratify the mixed layer and do not generate submesoscale eddies.
In conclusion, our analytical results for flows driven by wind and baroclinic instability with different directions of the wind and the buoyancy gradient, are in agreement with the results of eddy resolving simulations by Mahadevan et al. (2010).

Testing the SM parameterization against simulations with wind and buoyancy
In this section, we compare our model results with Capet et al. (2008, C8). We begin with Fig. 12 of C8 which gives ∂ z F T V (z) which we compare with our model (Eq. 4a). To do so, we assume that the direction of ∇ H T coincides with that of ∇ H b and that the wind blows in the down-front direction. Then, in analogy with Eq. (9a), from Eqs. (4a) and (8a) we obtain: e ζ β(ζ )+E 1/2 −γ e ζ α(ζ )+A −1 S g (z+ we have 0.8·10 −5 0 Cm −1 <|∇ H T |<2.7·10 −5 0 Cm −1 , we use |∇ H T |=1.5·10 −5 0 Cm −1 . As for the buoyancy frequency N, we use the characteristic value N=10 −3 s −1 . Next, we need the values of A, S g , h, δ E and K SM defined in the previous section. From the C8 data presented in Fig. 10, we derive the following results: δ E ≈ 10 m, h∼40 m, K SM ≈ 1.8 × 10 −3 m 2 s −2 While for the determination of A and S g we need the profiles of the projections of the mean velocity u(z), in their Fig. 10, C8 present only the absolute value |u(z)|. However, assuming that the wind blows in the down-front direction, one can extract A and S g from their Fig. 10. Finally, from Eqs. (12b), (4b) with r S =(N/π|f |)h and Eq. (10e), we derive that: Using this value of γ in Eq. (7e) together with Eq. (10c), we derive that: Regrettably, the C8 data are the only ones available that allowed us to estimate the parameter C and other data would be welcome to confirm Eq. (12d). Fortunately, as relation Eq. (10d) shows, the dependence of γ on C is rather weak. In addition, under condition (Eq. 9c), which is satisfied by the first result (Eq. 12c), in Eq. (12a,e) the terms linear in γ are small compared to the main terms which do not contain γ . Thus, we conclude that in the case of a strong down-front wind, the effect of C on the SM tracer flux is relatively weak. Next, substituting Eq. (12b,c) into Eq. (1c), one can see that the condition is amply satisfied. Finally, in Fig. 1 we compare the z-profile of −∂ z F T V (z) from Eq. (12a) with that of Fig. 12 of C8 (dashed line). In Fig. 2, we compare the profiles of the fluxes F T V (z) from the present model: with that of C8 which we compute using C8 data for the zderivative of the flux shown in Fig. 12. As one can observe, the profiles of the fluxes are quite close throughout the mixed layer depth. As for the profile of −∂ z F T V (z), they are quite  close in the upper half of the mixed layer but differ somewhat in the lower half. We think this is due to the similarity of the mean velocity profile (Eq. 8a) with that in the C8 at small depths and by an unavoidable difference in the lower part of the mixed layer due to the different profiles of the vertical viscosity used here and in the C8 simulations. While in our analysis we adopted an Ekman profile corresponding to a ν(z) = const., C8 adopted a more realistic model for ν(z). In spite of that difference, the profiles |u(z)| compare well, as seen in Fig. 3.

No-wind case -Buoyancy only
In addition to the realistic case studied by C8, there are data from simulations corresponding to the less realistic case of  no wind (Fox-Kemper et al., 2008, FFH) which can also be used to test of our model predictions. Following FFH, we assume that the mean velocity is in thermal wind balance with the mean buoyancy field and that the mean buoyancy gradient does not vary inside the mixed layer, i.e., u z =f −1 e z ×∇ H b is z independent. Irrespectively of the surface value of u, from Eq. (4b) and the second of Eq. (6b), we derive that: Taking τ =b in Eq. (6a,b), and using Eq. (13a), the buoyancy flux is given by: To complete this parameterization of F V , γ must be obtained from solving Eq. (7e). Substituting Eq. (13a) into the second of Eq. (7d), we obtain: Substituting this result into Eq. (7f), we obtain: The solution of Eq. (7f) is then as follows: γ 2 = C * Ri + Ri 2 + 2Ri/C * , C * = 6 π 2 (2C) 3/2 (14c) If we employ the value C=12 that we have determined from the data of Capet et al. (2008), we can compute the function (Eq. 13b): which we exhibit in Fig. 4. We further notice that in spite of the fact noted before that we have only one set of data to determine C, Fig. 5 shows that different C's have only an overall marginal effect on the buoyancy flux in the interval 1<Ri<100. Next, we find the limits of applicability of the parameterizations (Eqs. 13b and 14c,d). To this end, we find the baroclinic mean kinetic energy averaged over the ML depth K using the definition (Eqs. 7g, 4b and 13a): This result shows that at the surface the mean baroclinic kinetic energy K is twelve times smaller than the mean kinetic energy K M = 1 2 u 2 . To compute the SM kinetic energy, we use the first definition (Eqs. 4b and14c) and derive that: Thus, condition (Eq. 1c) becomes: Ri > π 2 36 3C 3/2 −1 ∼ 2 × 10 −3 (14g) Next, we compare Eq. (13b) with the FFH data and recall that in their Fig. 14e the authors plot the ratio: where: is the parameterization suggested by FFH by fitting their simulation data. To compare the model results with the simulation data, we construct the ratio: Thus, (FFH)=1, by definition (blue line in Fig. 4). To compute Eq. (15d) in our model, we notice that the profile µ(z) in Eq. (15c) has the additional factor (1+5ξ 2 /21) in comparison with the profile (Eq. 13b). The difference does not exceed 25% and we neglect it when substituting Eqs. (13b) and (15b,c) into Eq. (15d). As a result, we obtain: which we show in Fig. 4 (red line). The present parameterization yields a good representation of the FFH simulation data especially the curvature exhibited by the data. The agreement is somewhat worse at Ri 10 3 which is due to the use of Killworth's (2005) approximation to neglect τ z ; in fact, such an approximation becomes questionable when Ri is large, that is, when b z =N 2 is large. Finally, we discuss whether the FFH flux (Eq. 15b,c) without wind can represent the case with wind. To that end, we take the value of S g in the first of Eq. (12b) as determined from C8 simulations and substitute it in Eq. (11b) with the result: Substituting this result in Eq. (15b,c), and using the same mixed layer depth h=40 m as in C8, we obtain: If one compares this value with the value from C8 (Fig. 2): one concludes that the FFH flux formula with no wind underestimates the flux by about an order of magnitude, a conclusion in agreement with Mahadevan and Tandon (2006) who stressed that "winds play a crucial role in inducing submesoscale structure" not to mention the multifaceted and important implications on the mixed layer stratification caused by the down-wind, up-wind vs. buoyancy topology described in Sect. 6.

Conclusions
Recently, there has been a considerable interest in submesoscales which are oceanic structures of O (1 km) size and life times of the order of days. High resolution numerical simulations have been the best source of information to assess the parameterizations of SM fluxes to be used in OGCMs that do not resolve such features. If one considers that the highest resolution of about 1/10 0 in stand-alone OGCMs can represent structures of about 10km size which is 10 times larger than the SM sizes and that OGCMs employed in thousand years runs for climate studies have resolution of about 1 0 (corresponding to structures 100 times larger than SM), it seems clear that a good deal of important physical processes have thus far gone unrepresented in those OGCMs. The present work presents an analytical parameterization of the SM vertical flux of an arbitrary tracer. The main features can be described as follows: (a) no other SM parameterization exists (to the best of our knowledge) that provides an analytical expression for the vertical flux of an arbitrary tracer under arbitrary buoyancy and wind (strength and direction), (b) the SM parameterization by FFH does not include winds and it is limited to the buoyancy field which cannot be used to describe tracers such as the ones needed in the C-cycle models in OGCMs for climate studies, (c) the results of the realistic simulations by C8 and MTF with baroclinic instabilities and winds, are well reproduced by our model, (d) our parameterization given by Eqs. (4a-c), (7e-f), (12d), is to be used in OGCMs that resolve mesoscales but not sub-mesoscales, (e) in a different study , we have derived the parameterization for the tracer vertical flux for OGCMs that do not resolve either submesoscales or mesoscales.
which changes Eq. (C1) to the form: Relation (Eq. C2) implies that the dependence of the submesoscale fields on ω is of the form: and therefore in the (t,k)-space the fields A depend on time as follows: Due to relations (Eqs. C7,C8), after substituting Eq. (3a) in Eq. (2e), the latter may be solved in both (ω, k) and (t,k) representations.

C2 Sub-mesoscale velocity field w
It is convenient to begin by splitting the mesoscale velocity field u into a rotational (divergence free, solenoidal) and a divergent (curl free, potential) components: u D( k) = nu D , n = k/k and thus the third equation in Eq. (2e) becomes: To determine u R,D , we substitute the second relation (Eq. 3a) in the second equation in Eq. (2e) and derive the following expressions: These relations, as well as Eq. (C6), are valid in both (ω, k) and (t,k) representations. Below, we will use them in the (t,k) representation together with Eq. (C8). Due to the third relation of Eq. (3a), the further use of Eqs. (C11) and (C12) is considerably simplified under the condition which allows us to neglect the second terms in the brackets in Eqs. (C11) and (C12). Condition (Eq. C13) coincides with Eq. (1c).

C3 Spectrum of the vertical SM flux
In the dynamical Eq. (1a) one needs ∂ z F V which we shall write as follows: where in the last expression we have neglected the term w τ z in accordance with the adopted approximation τ and because it is of a higher order in z. Substituting Eq. (C11) into Eq. (C10), we obtain the expression for ∂ z w which allows us to compute ∂ z F V using Eq. (C14). The strategy of computing submesoscale fluxes, which are bilinear correlation functions, consists in computing these functions in (t,k)-space which, in the approximation of homogeneous and stationary mean flow, have the form: which, because of relation (Eq. C8), does not depend on t.
The function ReA B * (k) is usually referred to as the density of A B in k-space. The spectrum of the correlation function A B is: i.e., the spectrum of A B is obtained by averaging ReA B * (k) over the directions of k and multiplying the result by π k. Finally, the correlation function A B in physical space is obtained by integrating over the spectrum. Following this procedure, from the first of Eq. (C6) and using Eqs. (C9,C10), we derive the relation: Re w z τ * (k) = (C17) − I m kχ 2 (χ + ik· u) u D u * R (k)n×e z +|u D | 2 (k)n · ∇ H τ Next, from Eqs. (C11,C12) we derive the following relations: Re u D u * R (k) = χf −1 |u R | 2 (k), Substituting Eq. (C18) into Eq. (C17) and averaging over the directions of k, we obtain the spectrum of w z τ (k). Under conditions (Eq. C13), we get: w z τ (k) (C20) = − k 2 0 χ −2 χf −1 |u R | 2 (k) u ×e + |u D | 2 (k) u · ∇ H τ 692 V. M. Canuto and M. S. Dubovikov: Mixed layer sub-mesoscale parameterization -Part 1 where: πk|u R | 2 (k) = 2 1 + χ 2 f −2 −1 E(k), πk|u D | 2 (k) (C21) = 2 1 + χ −2 f 2 −1 E(k), E(k) = 1 2 πk|u | 2 (k) where E(k) is the spectrum of the total (rotational + divergent) eddy kinetic energy. Due to relation (Eq. C14), the left hand side of Eq. (C20) multiplied by πk is the spectrum of the z-derivative of the vertical flux ∂ z F V (k). Thus, multiplying Eq. (C20) by πk, using Eq. (C21), we obtain the following expression for the spectrum of ∂ z F V (k) near the maximum of the energy spectrum: where: Assuming that the spectra F V (k) and E(k) have similar shapes, integration of Eqs. (C22,C23) over k reduces to the substitution of E(k) and F V (k) with the eddy kinetic energy K SM and F V in physical space. Thus, we get Eq. (4a,b).