Derivation and assessment of a mixed layer sub-mesoscale model

Derivation and assessment of a mixed layer sub-mesoscale model V. M. Canuto and M. S. Dubovikov NASA, Goddard Institute for Space Studies, New York, NY, 10025, USA Dept. Appl. Math. and Phys. Columbia University, New York, NY, 10027, USA Center for Clim. Systems Res., Columbia University, New York, NY, 10025, USA Received: 7 August 2009 – Accepted: 5 September 2009 – Published: 17 September 2009 Correspondence to: V. M. Canuto (vcanuto@giss.nasa.gov) Published by Copernicus Publications on behalf of the European Geosciences Union.

ocean cannot be extrapolated to describe sub-mesoscales and a new one must be devised.Most of our knowledge about mixed layer sub-mesoscales comes from high resolution studies carried out by several groups (Levy et al., 2001(Levy et al., , 2009;;Thomas and Lee, 2005;Mahadevan, 2006;Mahadevan and Tandon, 2006;Klein et al., 2008;Thomas et al., 2008;Capet et al., 2008).These studies have revealed many interesting features of sub-mesoscales especially their contribution to the vertical mixing of mass, buoyancy and tracers in the upper ocean.Among the most salient effects of sub-mesoscales on global ocean properties is the well documented tendency to re-stratify the mixed layer (Spall, 1995;Nurser and Zhang, 2000).The effect of sub-mesoscales 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 sub-mesoscale and small scale fluxes leading to the mixed layer re-stratification (e.g., Capet et al., 2008, Fig. 12;Klein et al., 2008, Sect. 4); specifically, Hosegood et al. (2008) estimated that sub-mesoscales 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 sub-mesoscales concerns the location of the WBC that is shifted south by 4 • 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), which found a significant increase (from ∼30 Sv to ∼70 Sv) in the barotropic transport in the Gulf Stream when moving from 1 • to 1/6 • resolution, was recently confirmed and the inclusion of submesoscales further increased the transport by ∼50 Sv.Finally, the structure of the MOC (meridional overturning circulation) was also significantly affected by the presence of sub-mesoscales not so much in its intensity as to its location (Levy et al., 2009, Fig. 12).Introduction

Conclusions References
Tables Figures

Back Close
Full This brief summary of some of the results of very high resolution (1/54 • , ∼2 km) regional studies is sufficient to highlight the importance of sub-mesoscales and thus the question arises as to how much of sub-mesoscale physics is actually accounted for by present day global OGCMs (ocean global circulation models).Even today's highest resolution global ocean models ∼1/10 • (Maltrud and McClean, 2005;Sasaki et al., 2008) are not able to capture the sub-mesoscale field and much less are in the position to do so the OGCMs coupled to an atmospheric model used for climate studies where the resolution is at best 1 • but generally lower.The significant global processes revealed in going from 1 • resolution (∼100 km) and 1/54 • resolution (∼2 km) are presently absent in such global models especially in climate studies.It therefore seems that a reliable parameterization of submesoscales in terms of the resolved fields has become necessary to ensure the physical completeness of mixed layer mixing processes.
When solving the dynamic equations of a coarse resolution OGCM, two key equations are those for T and S (temperature and salinity) but since climate models must also include a carbon cycle, one must consider a general tracer field whose dynamic equation is given by1 : where for completeness we have included the contribution of mesoscales which are however not treated in this work which is restricted to sub-mesoscales.The submesoscales horizontal and vertical fluxes are defined as follows: where a prime denotes the submesoscale fields and an overbar stands for an ensemble average over resolved scales.However, since the horizontal sub-mesoscales flux is smaller than the corresponding one due to mesoscales, we shall concentrate on the parameterization of the vertical sub-mesoscale flux, F V (τ).The methodology we employ to carry out the parameterization contains three steps: first, one solves in Fourier space the sub-mesoscale dynamic equations describing the sub-mesoscale fields w , τ ; second, one forms the averages of the product of the two fields so as to obtain the second-order correlation w τ and third, one integrates over all wavenumbers to obtain the final expression for F v in coordinate space and expressed in terms of the resolved fields, to be used in Eq. ( 1a).The procedure was first worked out for the linear case by Eady and later by Killworth (1997) and for the non-linear case by Canuto andDubovikov (2005, 2006, CD5,6).Though the dynamic equations describing the velocity and temperature fields are formally the same as those describing mesoscales that were discussed in CD5-6, in the present case they must be solved in the regime appropriate to sub-mesoscales represented by Ro=ζ /f =O(1) rather than Ro 1 as in the case of mesoscales which leads to the appearance of terms that were not present in the Ro 1 regime.
The key difficulty in solving the sub-mesoscale dynamic equations is represented by the 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 rather involved, in Appendices A, 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 its use in the present context.Such an assessment was carried out 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 (e.g., Canuto et al., 1996Canuto et al., -1999)).Even so, we consider such assessment nec- essary but not sufficient for the credibility of the parameterization of sub-mesoscale fluxes we derive in Sects.5-6.The additional requirement consists in assessing the model predictions against results from sub-mesoscale resolving simulations of the type cited in the first part of this discussion.The first simulation corresponds to a system forced by baroclinic instabilities and no wind (Fox-Kemper et al., 2008, FFH, 2 km resolution) and the second one consists of a flow under realistic wind and buoyancy forcing (Capet et al., 2008, 0.75 km resolution).We shall present a detailed comparison of our parameterization with these simulation data.
To make the sub-mesoscale model results usable in OGCMs, we looked for analytical solutions of the sub-mesoscale dynamic equations and to achieve that goal, we introduced one approximation that consists in assuming 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 and for an arbitrary Ri, was worth exploring.Following Killworth (2005) suggestion, we adopt the approximation that due the mixed layer strong mixing, one can neglect τ z .Anticipating our main result, the vertical gradient of the vertical flux that enters in the original 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. (1c) the term w + S τ z , where w + S is found from the continuity condition ∂ z w (Killworth, 2005).From Eq. (1c) we shall derive the vertical flux which we write as: where u g,a are the geostrophic and a-geostrophic velocity components, the latter being a key component in the case of strong winds.The structure of the paper is as follows.In Sect. 2 we discuss the dynamic equations for the sub-mesoscale fields in the vicinity of the ocean's surface; in Sect. the turbulence closure model to the non-linear term in the sub-mesoscale tracer equation in Fourier space in the vicinity of the maximum of the energy spectrum and find the solutions for the tracer field τ ; in Sect. 4 we do the same with the sub-mesoscale for the momentum equation and derive an expression for the submesoscale field u .Using these results, in Sect. 5 we compute the spectrum of ∂ z F V in the vicinity of its maximum and then compute ∂ z F V and F v in physical space in terms of resolved fields as well as the sub-mesoscale eddy kinetic energy K E .In Sect.6 we express K E in terms of the resolved fields that, together with the results of the previous section, complete the problem of expressing ∂ z F V in the presence of both frontogenesis and Ekman pumping.In Sect.7 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 wind or of ∇ H b, it tends to de-stratify the mixed layer but at the same time it generates sub-mesoscales that have the tendency to re-stratify the mixed layer.
On the other hand, when the wind blows in opposite directions to the previous ones, it tends to re-stratify the mixed layer and does not generates submesoscale eddies.In Sect.8 we compare the model results with the data from the sub-mesoscale resolving simulations of Capet et al. (2008).In Sect.9 we compare the model results for the no-wind case with the FFH simulation data.In Sect. 10, we present some conclusions.

Sub-mesoscales dynamic equations near the surface
Consider an arbitrary tracer field τ.Separating it into a mean and a fluctuating part, τ=τ+τ , the dynamical equations for the sub-mesoscale tracer field τ are 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 and no physical assumptions, we cite only the final result (for the notation, see the footnote): where the function Q's represent the non-linear terms; ∇ H is the horizontal gradient operator, an overbar stands for an ensemble average over resolved scales and the vertical diffusivity K v represents small scale mixing processes.As expected, the average of Eq. ( 2a) yields identically zero.It must be noted that in Eq. ( 2a) no closure has been used for the sub-mesoscale fields.Equation (2a) formally coincides with those describing mixed layer mesoscales tracer fields studied by Killworth (2005).The difference in representing mesoscales and sub-mesoscales lies in the scales over which the averages (represented by an overbar in Eq. 2a), is taken: in the mesoscale description, averages are over scales exceeding mesoscales while in the description 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 submesoscales Ro, Ri∼O(1).Following Killworth (2005), we neglect the terms containing τ z and τ z .Then, the first of Eq. ( 2a) simplifies to: Without the non-linear term, this equation is equivalent to Eq. (2) of Killworth (2005) for the mesoscale buoyancy field in the mixed layer.Within the same approximation, the equation for the horizontal eddy velocity is given by: where e z is the unit vector along z axis.Next, we Fourier transform Eq. (2b,c) in horizontal planes and time.Following Killworth (1997Killworth ( , 2005)) transformed.Thus, we obtain: where we have added the continuity equation that provides the z-derivative of w'.We recall that τ , u and the non-linear 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.
The solution of Eq. ( 2e) provides the necessary ingredients to construct the vertical flux (Eq. 1b).On the other hand, since in the dynamical equation (Eq.1a) we only need ∂ z F V , we shall write it as follows: where in the last expression we have neglected the term w τ z in accordance with the adopted approximation and since it is of a higher order in z.

The sub-mesoscale tracer field τ
As discussed in Appendices A, B, in the vicinity of k=k 0 corresponding to the maximum of the eddy energy spectrum E (k), the non-linear terms Q H have the following forms: where u is here understood to be in physical space and =k −1 0 is the characteristic submesoscale length scale.As it is stressed in the literature on this subject (e.g., review by Thomas et al., 2008;Fox-Kemper and Ferrari, 2008;Boccaletti et al., 2007), is closely related to the deformation radius in the mixed layer, which implies that:

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version Interactive Discussion
where h is the depth of the mixed layer and N is the buoyancy frequency.Equation (3b) can also be derived by solving the eigenvalue problem to which the submesoscale equations reduce and which is analogous to the eigenvalue problem in the case of mesoscale (see CD5).Substituting Eq. (3a) into the tracer equation Eq. ( 2e), we obtain the expression for the submesoscale tracer field: where |k| = k 0 .As in the case of mesoscale discussed in CD5, the frequency ω, obtained from solving the eigenvalue problem mentioned above, yields the following dispersion relation: This relation can be interpreted as the Doppler transformation for the frequency provided that in the system of coordinates moving with the velocity u d , the submesoscale flow is stationary in which case ω=0.Stated in different words, relation Eq. (3d) implies that u d is the eddy drift velocity whose expression in terms of mean fields is analogous to that for the case of mesoscales given in Eq. (4f) of CD6: where L=−∇ H b/N 2 is the slope of isopycnal surfaces and β=∇f .The bracket averaging is defined as follows: Due to the smallness of the scale characterizing submesoscales, the second term in Eq. ( 3e) is negligible and thus:

Conclusions References
Tables Figures

Back Close
Full Relation (Eq.3d) implies that the dependence of the submesoscale fields on ω is of the form: and therefore in (t, k)-space the fields A depend on time as follows: Due to relations Eq. (3i,j), after substituting Eq. (3a) in Eq. ( 2e), the latter may be solved in both (ω, k) and (t,k) representations.

The 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: 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.(3h), are valid in both (ω,k) and (t,k) representations.Below we will use them in the (t,k) representation together with Eq. (3j).To illustrate the physical content of the second of Eq. (4c), we notice that it can also be written, using the third of Eq. (3h), as: where Ro is the Rossby number.Thus, in the quasi-geostrophic approximation corresponding to small Ro, the first relation in Eq. ( 4d) reduces to the geostrophic relation u R →u g =−i kρ −1 p .However, since in the submesoscale regime typically Ro∼1, one must consider the complete expressions Eq. (4c).That is the reason why we have not called u R the geostrophic component and u D the a-geostrophic one.
Substituting Eq. (4c) into Eq.( 4b), we obtain the expression for ∂ z w which allows us to compute Eq. (2f).NOTE.Readers not interested in the details of the derivation can move directly to the final result Eq. (7a,b).

Sub-mesoscale vertical tracer flux
The strategy used to derive 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: Printer-friendly Version Interactive Discussion 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. (3h) and using Eq.(4a,b), we derive the relation: where we have introduced the notation: Next, from the first Eq.(4d) we derive the following relations: Substituting Eq. (6a) into Eq.( 5a) and averaging over the directions of k, we obtain the spectrum of w z τ (k).Under the condition we obtain: where: the z-derivative of the vertical flux ∂ z F V (k).Thus, multiplying Eq. (6d) by πk, using Eq.(6e), we obtain the following expression for the spectrum of ∂ z F V in the vicinity of the maximum of the energy spectrum: where: With the assumption that the spectra F V (k) and E (k) have similar shapes, integrating Eq. (6g,h) over k reduces to the substitution of E (k) and F V (k) with the eddy kinetic energy K E and F V in physical space.The result is: where: It is worth noticing that in the second relation in Eq. (7a) the term λe z × u 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 (f is the scalar product of the vector e and the pseudo-vector 2Ω) and thus the product is a vector.In Eq. (7a) the velocity u + S may be interpreted as the sub-mesoscale induced velocity which is a counterpart of the mesoscale induced velocity.As Killworth (2005) noticed, to make the analogy with the mesoscale induced velocity more complete, since in the mixed layer τ z is small due to the strong mixing, one may add to Eq. (7a) the term w + S τ z , where w + S is found from the continuity condition:

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion To simplify the use of Eq. (7a,b), we further approximate relation Eq. (3g) by taking the eddy kinetic energy to be constant in the mixed layer and thus we have: Thus, u and K may be interpreted as the ML baroclinic velocity of the mean flow and the ML baroclinic mean kinetic energy.The only variable in Eq. (7a) that is not yet determined is the eddy kinetic energy K E which we shall compute in the next section.
Before doing so and for future purposes we next derive the explicit form of the vertical flux itself.Integrating Eq. (7a) over z with the boundary condition F V (0)=0, we account for only the z-dependency of u within the Ekman layer.Thus, we obtain: where the submesoscale diffusivity is given by: Recall that results (Eq.7a,f) have been obtained under condition (Eq.6c).From relations (Eq.7f), one observes that at the bottom of the ML, z=−h, it follows that: which is a good approximation since submesoscale eddies almost do not penetrate the mixed layer bottom (Boccaletti et al., 2007).This result supports the approximation adopted in Eq. ( 2f) and in deriving the last of Eq. (7d).We also have that:

Conclusions References
Tables Figures

Back Close
Full Assuming that the production of eddy kinetic energy, denoted by P K , occurs at scales , we suggest the relations: Since P K is a power, upon multiplication by the dynamical time scale τ=2K E /ε one obtains an energy; with ε= −1 K 3/2 E , one then easily derives the first of Eq. (7i).A more basic justification can be found in Lesieur's (1990) book on turbulence.Furthermore, since the kinetic energy equation shows that the vertical buoyancy flux F v acts as the source of K E , we further suggest the second relation where h is the depth of the mixed layer.Using a mesoscale resolving simulation (Canuto et al., 2009) we validated Eq. (7i) and in the next section, see after Eq. (17a), we validate Eq. (7i) for submesoscales using the simulation data of FFH and Capet et al. (2008).The coefficient C is related to the Kolmogorov constant Ko by the relation, C=3Ko/2 with Ko=5/3.Substituting Eq. (7e,f) into Eq.(7i), we obtain the following algebraic equation for K E : where the velocity V is defined in terms of the mean velocity by the following relation: We notice that besides a non-trivial solution K E =0, Eq. ( 7j) has the solution K E =0 since the variable η defined in Eq. ( 7b), is proportional K E .The zero solution is realized when the non-trivial one is not positively defined.It means that in such cases sub-mesoscale Introduction

Conclusions References
Tables Figures

Back Close
Full eddies are not generated.Thus, for simulations that resolve mesoscales but not submesoscales, the z-derivative of the vertical tracer flux is given by Eqs.(7a,b) and (7j,k) for K E .To facilitate the solution of Eq. (7j) that yields K E in terms of resolved fields, in the next sections we consider two limiting cases: 1) strong winds, 2) no-wind case.

Wind driven flows
In this section we show that when the Ekman velocity exceeds that of the geostrophic mean velocity, one can derive analytical results that we shall compare with those from the submesoscale resolving simulations of Capet et al. (2008).To obtain results in analytical form, we further take the turbulent viscosity ν∼10 −2 m 2 s −1 to be 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: where ρu 2 * is the surface stress and δ E is the Ekman layer's depth.We begin by considering the case of along-front winds blowing in the direction of the surface geostrophic velocity when the winds drive dense water over buoyant and provide favorable conditions for the generation of submesoscales (e.g., Thomas, 2005;Thomas and Lee, 2005).Thus we have: which corresponds to an horizontal buoyancy gradient given by: To obtain the submesoscale flux and its z-derivative (Eq.7e, a), we need to compute u , u and u defined in Eq. (7d The condition for a strong wind is Let us now study the submesoscale buoyancy flux and its z-derivative which we obtain substituting Eq. (8c,e,f) into Eq.(7a,e) with τ=b.The results are: Let us compare the contribution of the vertical flux Eq. (9a) to mean tracer equation (Eq.1a) (for τ=b) with that of the baroclinic component u•∇ H τ of the mean advection term.From Eq. (8c,e) we have As one can see, this term de-stratifies the ML since its z-derivative is positive.Also, we notice that in the case of a sufficiently strong wind, the parameter λ in Eq. ( 9a) is considerably smaller than unity and thus the first term in the curly bracket in Eq. ( 9a) dominates.Since in this case η≤1, we conclude that:

Conclusions References
Tables Figures

Back Close
Full which means that sub-mesoscales re-stratify the ML.To compute the vertical submesoscale flux, we integrate Eq. ( 9a) over z to obtain: To express the eddy kinetic energy in terms of the resolved fields, we find first V defined in Eq. ( 7k).Use of Eq. (8f) yields: Substituting Eq. ( 10a), together with Eq. ( 8c), into Eq.( 7j), we obtain: Recall that in real flows the mixed layer thickness exceeds the Ekman one and thus in Eq. (10b) 0>ζ −1 0 >−1.Therefore, from Eq. ( 10a) we conclude that V y <0.In addition, under the third condition (Eq.8b) S g >0, we have V x >0.This means that the contributions of both projections of V to K E are positive.In the opposite case when S g <0, since for the strong wind case the first term of V x in Eq. ( 10a) dominates, we have V x >0, V y <0 as in the previous case.Then since η∼K E , as it follows from Eq. ( 7b), the only realizable solution of Eq. ( 10b) is K E =0 and thus submesoscale eddies cannot be generated.It is worth noticing that the contribution of the baroclinic component of the mean flow to Eq. ( 1a) given by Eq. ( 9a), 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.On the other hand, in the case of a weak wind, when the S g term dominates in Eq. ( 10a) and (10b), expression (10b) is positive and thus submesoscales are generated.
Next, we consider the case of a strong wind blowing perpendicularly to the geostrophic flow so that:

Conclusions References
Tables Figures

Back Close
Full and thus: In contrast with Eqs.(8d-f) and (10a), in this case, the geostrophic component contributes to the y-projections of u , u, u and V .In particular, we have that: Substituting this relation into Eq.( 7j), we obtain: In contrast with Eq. ( 10b), in Eq. ( 11d) the terms from the y component of V contribute to K E with the opposite sign compared to the term from the x component of V .This confirms the conclusion of Thomas (2005) and Thomas and Lee (2005) that along fronts winds which drive dense water over buoyant, provide the most favorable condition for submesoscale eddy generation.Still, in the case of a strong wind, when the first term dominates in Eq. (11d), there is a positive solution for K E if S g >0.Under the same condition, the baroclinic component of u de-stratifies the ML since instead of Eq. (9b), we now have: which has a positive z-derivative.At the same time, the submesoscale eddies tend to re-stratify the mixed layer.As one can see from Eqs. ( 9b) and (11e), in both cases "Ekman flow advects dense water over light" (Thomas et al., 2008).When the wind blows in directions opposite to the previous ones, it tends to re-stratify the mixed layer and does not generates mesoscale eddies.Introduction

Conclusions References
Tables Figures

Back Close
Full The only submesoscale resolving simulation data presently available in the literature that include both baroclinic instabilities and wind are those of Capet et al. (2008, C8).We begin with Fig. 12 of C8 showing the z-derivative ∂ z F T V (z) which we compare with our model Eq.(7a).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, from Eqs. (7a) and (7c,e) we obtain: From Fig. 11 of C8 we have 0.8×10 −5 0 Cm −1 <|∇ H T |<2.7×10 −5 0 Cm −1 , and thus we use |∇ H T |=1.5×10 −5 0 Cm −1 .Next, we need the variables A, S g , η, λ, z 0 , Though in their Fig. 10, C8 present only the absolute value |u(z)| and not the components of the mean velocity u that are needed in Eq. ( 7a), assuming that the wind blows in the down-front direction, one can extract the data needed in Eq. (7a, b).Specifically, we derive the following values: Once (12b) are substituted in Eq. (8e), one computes K defined in (Eq.6d) while the kinetic energy K E is computed from Fig. 10 of C8.With K and K , one then computes λ, η defined in Eq. (7b); as for the buoyancy frequency N, we take the characteristic mixed layer value of N=10 −3 s −1 .In Fig. 1 we compare the profile of −∂ z F T V (z) from Eq. (11a) with that of Fig. 12 of C8 (black 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 z-derivative of the flux shown in 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, c) 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.Though in our analysis we adopted the Ekman profile of the mean velocity which correspond to ν(z)=const.while C8 adopted a more realistic model for ν(z), the profiles |u(z)| compare well, as seen in Fig. 3.We carried out a further test the model: we used Eq.(10b) to compute K E averaged over the mixed layer.The results, compared with the value of C8 in their Fig. 10, are: which are quite close.On the other hand, from Eq. ( 12b) we obtain the Ekman kinetic energy (which yields the main contribution to the ML baroclinic mean energy) K ≈1.25•10 −3 m 2 s −2 .Thus condition Eq. ( 6c) of the applicability of the model is satisfied.

No-wind case
In addition to the case studied by C8, there are also data from simulations corresponding to the less realistic case of no wind (Fox-Kemper et al., 2008, hereafter FFH) which we consider since they serve the purpose of an additional test of our model predictions.
Following FFH, we assume that the mean velocity is in a 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 the second of Eq. ( 7f) and the first of Eq. (7d), we derive that: Taking τ=b in Eq. (7e, f), and using Eq. ( 13a), the buoyancy flux is given by: with: The baroclinic mean kinetic energy K near the surface is obtained from its definition in Eqs. ( 6f) and (13a): As for the dimensionless variable y defined in Eq. ( 7b), it is easy to express it in terms of the Richardson number Ri corresponding to a geostrophic mean velocity ( f π=hN): To solve equation Eq. (7j) for K E , we find V from Eq. (7k) with the result: Substituting this relation into Eq.(7j), using Eq.(7b), we obtain the following relation for x and its solution: where C=2.5.Recall that the model prediction Eq. (7e,f) and, therefore, Eq. ( 13b) are applicable under condition Eq. (6c), i.e.
x ≥ 1 Ri ≥ 1.5 (14d) From the second of Eq. (14c), we further obtain that the limit Ri→∞ corresponds to x max ≈4.Substituting the first of Eq. ( 14c) into the last of Eq. ( 13b) together with the first of Eq. (13a), we obtain: The first of Eq. ( 13b) then becomes: To compare Eq. ( 14f) with the FFH data, we recall that in their Fig.14e the authors plotted the ratio: where: is the parameterization suggested by FFH.Even though the simulation data exhibit a scatter by more than in order of magnitude, FFH interpret the line Λ = 1 representing their model, to be in agreement with the data.To compare our model results with the same simulation data, we construct the ratio: To compute Eq. ( 15d), we notice that the profile µ(z) in Eq. ( 15c) has the additional factor (1+5ξ 2 /21) in comparison with the profile Eq. ( 14f).The difference does not exceed 25% and it is due to the coarse approximation we used in the integration of Eq. (7a) over z as we discussed below Eq. (7d).Neglecting the additional factor, we substitute Eqs. ( 14f) and (15b,c) into Eq.( 15d) and obtain: where Φ(Ri) is given in Eq. (14f, e) which we recall is valid for Ri>1.5.Inspection of Fig. 4 shows that in the Ri interval where the simulation data are available, the FFH and present model are consistent.Finally, we discuss whether the FFH flux formula Eq. (15a,b) 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).The result is: 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 realistic value from C8 (Fig. 2): one concludes that the FFH flux formula underestimates the true flux by about an order of magnitude.Our results confirm the conclusion of Mahadevan and Tandon (2006) that "winds plays a crucial role in inducing submesoscale structure".

Conclusions
Recently, there has been a considerable interest in sub-mesoscales which are oceanic structures of O (1 km) size and a life time of the order of days.Langevin equation in k-space: in which the non-linear (NL) term of NSE is represented by the two terms: the turbulent forcing f t i (k, t) which is due to the infrared (small k) part of the NL interactions and ultraviolet part which is represented by the enhanced k-dependent dynamical viscosity reduces the two NL terms in Eq. ( A1) to the second one only which, in the notation of Eq. (2d,e), implies that: Contrary to the kinematic viscosity ν which does not depend on the size of the eddies, the turbulent viscosity ν t (k) which is due to the NL interactions, depends on the eddy size and its sum to ν is called the dynamical viscosity, ν d (k)=ν + ν t (k).The search for a suitable expression for ν t (k) dates back many decades and the first explicit expression is the heuristic one proposed by Heisenberg as discussed in Batchelor's book (1970, Sect. 6.6, Eq. 6.6.13),Equation (B2) has several interesting features worth discussing.First, it says that an eddy of arbitrary size (∼k −1 ) feels the effects of all the eddies smaller than itself, as the integral begins at k and accounts for all the wavenumbers from k to infinity.Equation (B2) naturally reduces to the kinematic viscosity ν when the size of the eddy is very small and k tends to infinity.Due to the presence of the kinematic viscosity, Eq. ( B2) is valid for arbitrary Reynolds number since it can be rewritten as: If one employs the Kolmogorov spectrum E (p)=Koε 2/3 k −5/3 , one obtains in the large Re regime: Re 1 : ν t (k) = (ν 2 + 3Ko 16 ) 1/2 ≈ ε 1/3 4/3 (B4) which is the well-known Richardson 4/3 law diffusivity ∼ 4/3 .Finally, relation Eq. (B2) shows that there is no such a thing as a unique turbulent viscosity since each eddy feels its own turbulent viscosity.In Eqs.(2e) and (3a) we are interested in the function ν t (k)≈ν d (k) in the vicinity of the maximum of the energy spectrum k=k 0 .Assuming that most of the energy is contained in that region, from Eq. (B2) we get ν d ∼k Thus, from Eq. (A7) it follows which is the closure form in Eq. (3a).The closure for the tracer field is analogous.
characterized by a small Rossby number Ro=ζ /f 1 (where ζ and f are the relative and planetary vorticities respectively), sub-mesoscales are characterized by Ro∼1 and Ri∼Ro −1/2 ∼O(1), where Ri is the Richardson number (Thomas et al., Sect.2).For these reason, mesoscale parameterizations devised for the deep , we keep the same notation u , τ for the submesoscale fields in the k−ω space and assume that the mean fields u and ∇ H τ are constant in time and horizontal coordinates when Eq. (2b,c) are Fourier Screen / EscPrinter-friendly Version Interactive Discussion which changes Eq. (3c) to the form: and, because of relation Eq. (3j), relation Eq. (4e) does not depend on t.The function Re[A B * (k)] is usually referred to as the density of A B in k-space.The spectrum of the correlation function A B is: is the spectrum of the total (rotational+divergent) eddy kinetic energy.Due to relation Eq. (2f), the left hand side of Eq. (6d) multiplied by πk is the spectrum of determine the eddy kinetic energy K E defined in the last relation in Eq. (3a).
Fig. 12.As one can observe, the profiles of the fluxes are quite close throughout 2177 2 E (p) 1/2 d p, γ = O(1) (B1)where E (k) defined in Eq. (A2) is the kinetic energy spectrum whose integral over all wavenumbers yields the eddy kinetic energy K E .As discussed by Batchelor, Eq. (B1) was successfully used to derive the Kolmogorov spectrum.A non heuristic derivation of ν t (k) has however been lacking until recently with the advent of methods to treat the Navier-Stokes equation borrowed to a large extent from quantum field theory.A full presentation was made by the present authors in 1997 with the final result:ν t (k) = (ν 2 + 1
Though the full picture Introduction , numerical simulations of increasingly higher resolution have been a rich source of information that serves as a test bed for parameterizations of the submesoscale fluxes to be then used in low resolution OGCMs to represent structures they cannot resolve.If one considers that the highest resolution of about 1/10 • stand alone OGCMs can resolve structures of about 10 km size which is 10 times larger than submesoscale sizes and that OGCMs employed in thousand years runs for climate studies can hardly afford even a 1• resolution which corresponds to structures 100 times larger than sub-mesoscales, one realizes that a good deal of important physical processes have thus far gone unrepresented in low resolution OGCMs.The present paper presents a model of the key unresolved quantity in the dynamic equations, namely the vertical flux of an arbitrary tracer which we express in terms of the resolved fields.The results can be applied to the equations for the mean T (temperature), S (salinity) and C (passive scalar such as Carbon).The key relations are given by Eq. (7a, b, j, k), is to be used in OGCMs that resolve mesoscales but not sub-mesoscales.The model predictions have been tested against the results of sub-mesoscale resolving simulations.In future work we shall derive expressions for the tracer vertical flux to be used in OGCMs that do not resolve either submesoscales or mesoscales.Introduction where ν is the kinematic viscosity while ν t (k) is a turbulent viscosity discussed in Appendix B. As discussed in the references cited above, the dynamic equation for the energy spectrum E (k) is obtained by multiplying (A1) by u Introduction Thomas, L. N. and Lee, C. M.: Intensification of ocean fronts by down-front winds, J. Phys.Oceanogr., 35, 1086-1102, 2005.Thomas, L. N., Tandon, A., and Mahadevan, A.: Submesoscale processes and dynamics.In: Eddy Resolving Ocean Modeling, edited by: Hecht, M. and Hasumi, H., (AGU Monograph), Introduction