Comparison of in situ microstructure measurements to different turbulence closure schemes in a 3-D numerical ocean circulation model

In situ measurements of kinetic energy dissipation rate ε and estimates of eddy viscosity K Z from the Gulf of Lion (NW Mediterranean Sea) are used to assess the ability of − k ɛ and − k ℓ closure schemes to predict microscale turbulence in a 3-D numerical ocean circulation model. Two di ﬀ erent surface boundary conditions are considered in order to investigate their in ﬂ uence on each closure schemes ’ performance. The e ﬀ ect of two types of stability functions and optical schemes on the − k ɛ scheme is also explored. Overall, the 3-D model predictions are much closer to the in situ data in the surface mixed layer as opposed to below it. Above the mixed layer depth, we identify one model ’ s con ﬁ guration that outperforms all the other ones. Such a con ﬁ guration employs a − k ɛ scheme with Canuto A stability functions, surface boundary conditions parameterizing wave breaking and an appropriate photosynthetically available radiation attenuation length. Below the mixed layer depth, reliability is limited by the model ’ s resolution and the speci ﬁ cation of a hard threshold on the minimum turbulent kinetic energy.


Introduction
Turbulence is an essential mechanism for the transport of energy, salinity, and suspended and dissolved matter.Turbulent fluxes of such quantities are the result of correlated, small-scale fluctuations of the velocity field and of the transported quantity itself.The prevalent turbulence production mechanisms in coastal ocean are: mean shear, unstable stratification, Langmuir circulation (Farmer and Li, 1995) and breaking surface waves (Agrawal et al., 1992).For coastal ocean, mean shear is mainly generated by the action of winds and tides, but also by surface waves and baroclinic flows (e.g., Thorpe, 2005), including nonlinear internal waves (Toole and Schmitt, 1987).Unstable stratification results from surface processes such as surface cooling, evaporation or differential advection (e.g., Kantha and Clayson, 2000).Destruction of turbulence occurs by transformation into potential energy during stable stratification or viscous dissipation into heat (e.g., Kantha and Clayson, 2000).The complexity of these processes by themselves and of their interactions requires numerical models to cover a wide range of spatio-temporal scales and Reynolds number (e.g., Burchard et al., 2008).This is especially true in the upper ocean where all the above phenomena concur together to generate turbulence.
Upper ocean connects -through various turbulent mechanisms-the surface forcing from the atmosphere with the quiescent deeper ocean where heat and fresh water are sequestrated and released on longer time and global scales (Ferrari and Wunsh, 2009).Also, upper ocean turbulence plays an important role in biological phenomena by, for example, determining phytoplankton growth rate (Thomas and Gibson, 1990), influencing primary production (Flierl and Davis, 1993) and the onset of blooms (Taylor and Ferrari, 2011).The complexity of modelling such mechanism within ocean circulation numerical models gave rise to several approaches.In particular, many turbulence closure schemes have been proposed.The ones most frequently found in the ocean modelling community's literature are the − kk ℓ by Mellor and Yamada (1982); the − k ɛ by Rodi (1987); the − kk ω by Wilcox (1988); the − k ℓ by Gaspar et al. (1990) and the KPP by Large et al. (1994).Following recent numerical modelling literature (Ilicak et al., 2008;Reffray et al., 2015), in the present study, we consider the − k ɛ and − k ℓ second moments closure (SMC) schemes.Note that other kinds of closure schemes such as the KPP (Large et al., 1994) are not considered here being not as well suited as the other two schemes for a comparison with in situ data of kinetic energy dissipation rate ε.
Additional complexity is added to the modelling by the interplay of the SMC and the choice of boundary conditions.The choice of surface and bottom boundary conditions can also profit from a vast literature (e.g., Craig and Banner, 1994;Stacey and Pond, 1997;Estournel et al., 2001;Warner et al., 2005), aiming at modelling different forcing mechanisms.Furthermore, different stability functions can be chosen in order to include the effect of the parameterized non-local moments and pressure strain correlations in the dynamical equations (e.g., Galperin et al., 1988;Kantha and Clayson, 1994;Canuto et al., 2001).The choice of the optical scheme is particularly important considering the high number of studies coupling Symphonie to biochemical models as it can influences turbulent fluxes and nutrient availability.
Thus, the in situ validation of the closure schemes, boundary conditions, stability functions, optical scheme and their interplay is fundamental for assessing the reliability of numerical models (Warner et al., 2005;Peters and Baumert, 2007;Arneborg et al., 2007;Ilicak et al., 2008).
The current study presents the comparison of kinetic energy dissipation rate ε measurements and vertical eddy viscosity K Z estimates issued from a Self Contained Microstructure Profiler (SCAMP) with the predictions of a 3-D numerical ocean circulation model (Symphonie; Marsaleix et al., 2008) obtained with different model's setup.The aim is to gain some insights on which scheme and/or boundary conditions permit to have the representation of turbulence activity closer to the observations.
The GoL is located in the northwestern Mediterranean Sea and is characterized by a large continental margin (Fig. 1) and complex hydrodynamics (Millot, 1990).Its circulation is strongly influenced by the southwestward along-slope Northern Current.This density current flows in a cyclonic way and constitutes a barrier between the coastal waters of the continental shelf from the open northwestern Mediterranean Sea (Alberola and Millot, 1995;Sammari et al., 1995;Petrenko, 2003).Cross-shore exchanges between the GoL and offshore waters are regulated by wind induced dynamics (Estournel et al., 2003;Hauser et al., 2003;Petrenko et al., 2017) and by processes associated with the Northern Current, such as intrusions into the continental shelf and barotropic and baroclinic instabilities (Conan and Millot, 1992;Flexas et al., 1997;Petrenko et al., 2005;Barrier et al., 2016).The Gulf of Lion is a suitable case study because of the high number of physical (Qiu et al., 2010;Hu et al., 2011), sediment dispersion (Bourrin et al., 2011) and biochemical (Pinazo et al., 2001;Herrmann et al., 2014) numerical studies carried out there.
Symphonie has already been validated on a variety of different aspects like current modelling and eddy generation (Rubio et al., 2009;Hu et al., 2011;Kersalé et al., 2013), river plume dynamics (Reffray et al., 2004;Gatti et al., 2006) and dense water formation (Dufau-Julliand et al., 2004;Estournel et al., 2016).But a study of the different SMC that the user can implement in the Symphonie code has not yet been done.In particular, the modeling of the near-surface physical and biogeochemical processes is sensitive to the choice of SMC and the computed K Z values (Fraysse et al., 2014).
In general, we can regard all modeled large-scale circulation features in an integrated fashion as they result from successive calculation steps and approximations.Hence, a major difficulty in validating numerical models -beside the high number of variables at play-is the possible compensation of different errors between each other.This fact makes difficult to attribute a specific amount of the total error on a certain quantity to a specific step in its calculation, in the present case the turbulence scheme.Here, our goal is to assess the model predictions focusing on turbulence modelling in the most realistic configuration we can achieve: 3-D dynamics with realistic forcing.Indeed, mixing has a primary role in influencing the large-scale circulation motion (Rhines, 1988;Ferrari, 2014).Nevertheless, the 3-D nature of the model brings in play an augmented number of numerical issues, among which spurious numerical diffusion (Marsaleix et al., 2008;Marchesiello et al., 2009;Hu et al., 2009) that are usually neglected in similar, but 1-D, studies (Gaspar et al., 1990;Burchard et al., 2002;Reffray et al., 2015).
The manuscript is organized as follows.In Section 2 we describe the properties of the numerical model, the microstructure in situ data and how we carry out the comparison between them and the numerical data.In Section 3 we report the results of our analysis and we discuss them in Section 4.I nSection 5 we summarize the conclusions of our study.

Numerical modelling
The numerical model Symphonie is a 3-D primitive equations, free surface, sigma coordinate ocean model, based on Boussinesq and hydrostatic approximations (Marsaleix et al., 2008;2009;2012).Components of current, temperature and salinity are computed on a C-grid (Arakawa and Lamb, 1977) using a classic finite difference method detailed in Marsaleix et al. (2006);2008).This model has been extensively used in studies of the Mediterranean Sea, mostly at the scale of the continental shelves (Ulses, 2005;Estournel et al., 2003;2005), generally comparing satisfactorily with available in-situ observations of classical hydrological quantities.Symphonie has also been coupled to biochemical models for studies that demonstrated the impact of the turbulence level on determining the vertical flux of nutritive salts, the nutricline depth and -as a consequence-the results given by the biochemical models (Ulses et al., 2016;Herrmann et al., 2013).However, a study of the consequences of choosing a certain model's implementation of the Symphonie code has not yet been done.To fill this gap, we compare the model predictions of ε and K Z with the values measured with the SCAMP profiler.
The model domain we use, shown in Fig. 1, is that of Estournel et al. (2016).Note that all the measurements sites are far from the open boundaries.The horizontal resolution of the model grid is 1/110°(about 1 km).All the numerical experiments we perform cover the whole period in which in situ data are available: from 1 July 2010 to 13 March 2014, plus ten weeks of spin up.
In the vertical the model exploits a generalized sigma coordinate with 50 levels.Surface fluxes are computed using the bulk formulae by (Large and Yaeger, 2004) and the 3-h ECMWF by Estournel et al. (2016).The boundary condition for ε is deduced with a length scale reasoning from the value of the gradient Richardson number (i.e., the ratio between buoyant production of turbulence and shear production of turbulence); see Estournel and Guedalia (1987), Michaud et al. (2012) and Appendix A for more details.In order to simulate the limiting effect of stable stratification, following Galperin et al. (1988), the minimum of ε is linked to the minimum turbulent kinetic energy value k min through:  (Gaspar et al., 1990;Burchard et al., 2002), that is based on the estimate of the internal wave activity.
Similar low-frequency buoyancy conditions are maintained for all numerical experiments using a nudging procedure on temperature and salinity toward the corresponding MERCATOR fields (product PSY2V4R4).The nudging time scale is 30 days, enabling the free development of higher frequencies (including those of the turbulence closure scheme), and, at the same time, ensuring that the different turbulence schemes are tested in similar general conditions of stratification.
We choose the − k ɛ (Burchard and Bolding, 2001) and − k ℓ (Gaspar et al., 1990) closure schemes because, other than being the more exploited by Symphonie's users, they are also widely used in the wider scientific community.Reffray et al. (2015) showed that in a 1-D case the − k ɛ scheme gives more reliable mixing estimate with respect to other schemes widely used in the literature: − kk l , − k ω and − k ℓ.The − k ℓ scheme in Reffray et al. (2015) is based on Gaspar et al. (1990) but simplified for 1-D applications.We want to test if the original scheme by Gaspar et al. (1990) performs better than − k ɛ in a 3-D case.
Moreover, our questioning the numerical results' sensitivity on the value of k min follows from the study by Gaspar et al. (1990).Therein the authors encouraged (but not implemented) the use of a parameterization k min rather than fixing a hard value.Herein we test if a good result can be achieved in a simpler way by specifying a different value of k min .
A study by Burchard and Bolding (2001) showed that, in a 1-D study of temperature and mixed layer depth data of the well-known dataset OWS Papa (northern Pacific), the − k ɛ closure scheme performs better when employing the stability functions proposed by Canuto et al. (2001) -commonly called Canuto A-rather than the ones by Kantha and Clayson (1994), Rodi (1980) and Hossain (1980).On the other hand, Ilicak et al. (2008) showed that -in a 3-D study of the Red Sea outflow-the stability functions of both Canuto et al. (2001) and Kantha and Clayson (1994) perform similarly.We want to further investigate the different performance of these two stability functions in a 3-D case.Note also that this study -as opposed to Burchard and Bolding (2001) and Ilicak et al. (2008)-has a strong focus on microstructure measurements and not only on more standard quantities like mixed layer depth, temperature and salinity.
In the literature there are different formulations for K Z .In particular, it can include the molecular diffusivity D T (so that we always have K Z > D T ; e.g., Burchard and Bolding, 2001) or not (e.g., Han, 2014).Here we want to clarify the differences (if any) between the two approaches.
Not having the necessary computational power to explore all the possible combinations of these factors in a 3-D model (as done for example in a similar study in the 1-D case by Reffray et al. (2015), and in a 3-D case by Ilicak et al. (2008) but on a shorter time span), we restrict the study to a subset of combinations.
In particular, nine different numerical experiments employing different combinations of turbulent closure schemes, boundary conditions, stability functions, values of the minimum of turbulent kinetic energy and optical schemes are analyzed here (see Table 1 for a concise  Gaspar et al. (1990) -hereafter marked by the prefix KL.Details of these nine numerical experiments can be found in Appendix B. We test the effect of two possible surface boundary conditions.The first one (marked by a suffix set) supposes equilibrium between the production and dissipation terms in the dynamic equation for k (see Eq. (B.1)).The second one (marked by a suffix flu), takes into account the effect of breaking waves of all scales, as suggested by Craig and Banner (1994); see Eq. (A.1) and Appendix A for details.The numerical experiments exploiting the KE scheme and using the stability by Canuto et al. (2001) -instead of the ones by Kantha and Clayson (1994)are marked by the suffix CAN.Simulations with a higher minimum TKE -− 10k g m / s 72 instead of − 10k g m / s 8 2 -are marked with a suffix MINk.
All the numerical experiments with a higher threshold on k are such that K Z > D T .Therefore, in the other numerical experiments the total diffusivity could be smaller than the molecular one.
One numerical experiment (marked by the suffix Opt) investigates the effect of the attenuation length of the penetrative solar radiation Q sr .Q sr is parameterized as a two-band exponential scheme (Maraldi et al., 2013): ( 1 ) where the first right-hand term parameterizes the attenuation of red and near-infrared radiation (whose attenuation length is = l 0.35 cm); and the second right-hand term is the one of the visible and ultra-violet radiation; l PAR is the photosynthetic available radiation diffuse attenuation length.Q sr (0) is the fraction of the available penetrative solar radiation at the surface assuming a constant albedo of 6.6% and R = 0.54 determines the fraction in each band of Q sr (0).Following Maraldi et al. (2013), the l PAR default in Symphonie is set to 11 m.However, this value has to be considered as an annual climatological estimate of the photosynthetic available radiation.With this numerical experiment we test the effect of the seasonality of l PAR by setting its value to 23 m coherently with the fact that most of the in situ measurements were acquired in September when we expect a low biological activity in the surface boundary layer of the GoL.With this set of numerical experiments we can answer three principal questions: 1) which SMC between − k ɛ and − k ℓ performs better with respect to our dataset?2) what is the effect of the boundary conditions on the results of the numerical numerical experiments?and 3) what is the effect of the two stability functions?

SCAMP measurements
An in situ estimate of ε can be derived from high-resolution vertical profiles of temperature T. Batchelor (1959) derived the spectral shape of a conserved scalar field that is passively advected by an incompressible turbulent fluid with a molecular Prandtl number = Prν D / greater than 1 (seawater has = Pr 7 at 20 °C), where ν and D are re- spectively the molecular viscosity and the molecular diffusivity of the scalar.In the present case the scalar is the temperature.Gibson and Schwarz (1963) derived the one-dimensional Batchelor spectrum E(K) of temperature gradient as a function of the rate of dissipation of temperature variance χ T , the kinetic energy dissipation rate ε, the molecular diffusivity of temperature D T and the circular wavenumber K: T TB (3) where K B is the inverse of the Batchelor length scale describing the length scales at which fluctuations in scalar concentration (temperature in this case) can still exist before being evened out by molecular diffusion.Therefore, (see Ruddick et al., 2000;Luketina andImberger, 2000 andSteinbuck et al., 2009 for details) once measured the temperature vertical gradient at the millimeter scale and derived χ T and K B by fitting the Batchelor spectrum, a measure of ε follows from (Batchelor, 1959): The temperature gradient profiles were measured with a SCAMP profiler.This instrument is equipped with a 100 Hz FP07 glass rod microthermistors (sensitivity of 0.001 °C).
Our SCAMP was deployed in an upward configuration.After deployment, it sinks to a predetermined depth following an oblique trajectory and then rises up vertically at an approximately constant velocity = − U 10 m/s 1 .This type of trajectory permits the SCAMP to get away from the ship and be free from the influence of the ship's wake when rising up in an undisturbed water column.This allows to have reliable measurements of ε and K Z near the sea surface (Anis, 2006).
The dataset (spanning the period 2010-2014) for this analysis consisted of 126 profiles of variable vertical extent -between 1 m below the surface and 100 m depth-collected in different sites in the Gulf of Lion (Fig. 1).The profiles were collected during various oceanographic campaigns in the GoL conducted by the Mediterranean Institute of Oceanography (MIO -Marseille, France).Because the SCAMP measurements were opportunistically taken in cruises mainly not dedicated to turbulence measurements, repeated casts were often not possible.Data were collected mostly during summer and with meteorological conditions favorable to operations with a small boat, generally with wave heights less than half meter.Precipitation was always absent or negligible.Surface buoyancy flux was positive for 115 profiles indicating a gain of buoyancy by the ocean surface.
Temperature gradient spectra were computed from 128 points ( ≈ 13 cm) windows without overlap.The choice of this segmentation resulted from a sensitivity analysis using different segmentation methods proposed in the literature (see Appendix C).
The vertical eddy viscosity coefficient K Z can be derived on the basis of the turbulence intensity parameter = Reν N ɛ/ , b as follows.Re b expresses the ratio of the destabilizing effect of turbulence to the stabilizing effect of stratification and viscosity.In different Re b regimes, the mixing efficiency, expressing the portion of the energy produced by shear which is dissipated by viscosity, assumes different values and determines different vertical turbulent diffusivity of density K ρ .Here we use a recent field-validated parameterization of K ρ as function of Re b proposed by Bouffard and Boegman (2013) based on a previous parameterization derived by Shih et al. (2005) ), the turbulent regime is regarded as diffusive and K ρ is set equal to the temperature molecular diffusivity ) turbulent mixing tends to be controlled by buoyancy effects with incomplete mixing favoring up-gradient fluxes reducing the mixing efficiency.Here K ρ can be expressed as the mixing efficiency has the classical form derived by Osborn (1980) as we expect in weakly stratified fluids (Osborn, 1980).As in the case of the Mediterranean Sea (e.g., Cuypers et al., 2012), when the density variations are dominated by those of temperature, the density vertical eddy diffusivity coefficient is assumed to be equal to the temperature vertical eddy diffusivity coefficient (Peters et al., 1988).Assuming then a turbulent Prandtl number == PrK K / 1 tZ ρ (Hogg et al., 2001), it is possible to estimate the momentum eddy viscosity (or eddy viscosity) K Z to be equal to K ρ .

Comparison of numerical and in situ data
While the vertical resolution of the in situ profile of ε and K Z is constant (13 cm), the vertical resolution of the model is variable and generally coarser (≈− 12 m ) than the in situ one.In order to compare the numerical data to the measurements, the SCAMP data are grouped in windows centered on each sigma level.Then, the median of each window is calculated and compared to the numerical data at that sigma level.The choice of the median permits to give less importance to outliers and reduce the error due to the fact that the profiles were mainly single casts (Lozovatsky et al., 2005).
For both the in situ and numerical data, we define the surface mixed layer (MLD) as the depth at which the temperature is smaller than the surface value by 0.5 °C (Anis, 2006;Jurado et al., 2012).Then we separate the data in the MLD from the ones below it.The data in the bottom layer are not considered for this analysis because of the low number of in situ profiles near the sea bottom.
Following Burchard et al. (2002),w ed e fine the standard deviation between the logarithmic numerical and in situ values of ε as: where M is the number of points in a given profile, and ε mod and ε obs are the numerical and observed values respectively.Similar definitions hold for the decadal standard deviation of N (σ N ) and K Z (σ KZ ).To calculate these quantities, no division between data above and below the MLD is made because the number of in situ and numerical values in the MLD can differ for a given profile.Hereinafter, the mean values of σ ε , σ N and σ KZ on all the profiles are denoted σ ε , σ N and σ KZ for ease of reading, instead of < σ ε >, <σ N > and < > σ KZ .We also compare the probability density functions (PDFs) of the in situ and numerical values of kinetic energy dissipation rate ε, Brunt-Väisälä frequency N and eddy viscosity K Z .In particular, the PDFs are calculated with a kernel density estimation (Bowman and Azzalini, 1997) on the base of the frequency distributions of values from all the in situ and numerical profiles.
In order to compare a distribution of numerical data g(x) with the distribution of the in situ data f(x), we compute the squared difference between the empirical cumulative distribution functions (F n and G n )o f the two distributions: We define the shift S between two distributions as the shift of G n (x) that permits to minimize Δ 2 .A subscript will tell to which distribution the values of Δ 2 and S refer to: ∆ ɛ 2 and S ε for the PDF of ε; ∆ N 2 and S N for the PDF of N; ∆ K 2 Z and S KZ for the PDF of K Z .An estimate of the error on Δ 2 values is obtained by re-sampling the in situ empirical distribution function (Nerini and Ghattas, 2007) in order to further account for uncertainty in the measurements.

Results
To quantify the agreement between the numerical and the in situ data we use two methods.One takes into account the deviation of the numerical data from the in situ profiles over the whole water column and expresses it by the decadal standard deviations (see Section 3.1 and Table 2).The other approach looks at the shape of the probability density functions of the numerical values of ε, N and K Z and expresses it by the squared differences and shift values above and below the MLD (see Section 3.2 and Tables 3 and 4).We then also look at the in situ and numerical data in the (ε, N, K Z ) space (see Section 3.3).

In situ and numerical median profiles
Median profiles of dissipation ε for the numerical outputs and observations are shown in Fig. 2 vs. nondimensional depth z/MLD.Note that, the average value of the in situ MLD value is 27.3 m.As reported in Table 2, all the numerical MLDs are systematically lower by ≈− 61 5 % compared to in situ MLDs.Overall, the KE numerical ex- periments predict a slightly deeper MLD, closer to the observed values.
The thick lines in Fig. 2 represent the median dissipation profiles calculated using all profiles, where the shades represent the 95% bootstrap confidence interval.At the surface, all the numerical experiments agree with the in situ data to within one order of magnitude with the in situ data.At this depth, KEflu is the numerical experiment with values the closest to the observations.Also note that, at the surface, the equilibrium boundary condition generates higher dissipation values than the flux boundary condition.
As we descend the water column, beyond = zM L D /1 , the in situ data tend to become more variable due to the lower number of profiles reaching depths greater than the MLD.On the contrary, the numerical data tend to become much less variable.Moreover, there is an evident difference between the SMCs: when the lower threshold on k is applied, the KE numerical experiments always have lower levels of turbulence with respect to the KL numerical experiments.In this case, KEs' median S ε indicates the shift of the numerical distribution of ε in the surface layer respect to the in situ one.∆ ɛ 2 is the squared difference between the numerical and in situ distribution of ε values in the surface layer.Similar definitions apply for different subscript variables.Errors are calculated with a re-sampling procedure.values of ε are one order of magnitude lower than the KLs' ones.When a higher threshold on the kinetic energy k is applied, the difference is slightly reduced.
The in situ data within the mixed layer are characterized by higher mixing close to the surface and lower mixing near the MLD.This is reflected by a change of slope in the ε profiles at = zM L D /0 .5 .K E flu, KLflu, KLsetMINk and KLfluMINk are the numerical experiments that better reflect this behavior.Note also that the in situ data show a clear S ε indicates the shift of the numerical distribution of ε below the surface layer respect to the in situ one.∆ ɛ 2 is the squared difference between the numerical and in situ distribution of ε values below the surface layer.Similar definitions apply for different subscript variables.Errors are calculated with a re-sampling procedure.increment of ε at = zM L D /1 . 1 that is absent in the numerical data.We do not show values for z/MLD > 2 because of the large bootstrap error bars due to the lower number of profiles in that depth range.
Looking at the column as a whole, the numerical experiments with a higher minimum value of kinetic energy (KEfluCANMINk, KEfluCANMINkOpt KLsetMINk and KLfluMINk) have the lower average error with respect to the measurements (see values of σ ε in Table 2) because of the better agreement with in situ data below the MLD.
Median profiles of stratification N for the numerical outputs and observations are shown in Fig. 3 vs.nondimensional depth z/MLD.Near the surface, the in situ data are considerably more variable than the numerical data.However, the numerical data are of the same order of magnitude as the measurements.As we go below the MLD, the median value of the in situ data increases in the depth range 1 < z/MLD < 1.8.Overall, all the numerical experiments do not show such a marked increase in magnitude and are rather similar one to another (as we expect because all of them are nudged to the MERCATOR fields).However, KEfluCANMINkOpt seems to be the numerical experiment that better reproduce the in situ data's trend right below the MLD.
Looking at the column as a whole, KLsetMINk has the lowest average error with respect to the measurements (see values of σ N in Table 2).
The measurements of ε and N make it possible to calculate the turbulence intensity parameter Re b .The probability density functions of the observed Re b are reported in Fig. S1 in the Extra Materials.In particular, in our dataset, 30% of the in situ Re b values fell in the diffusive regime, 15% in the buoyancy-controlled regime, 29% in the transitional regime and 26% in the energetic regime.If only the data above the MLD are considered these numbers respectively become 17%, 10%, 31% and 42%; or 36%, 16%, 28% and 20%, if only the data below the MLD are considered.
Median profiles of eddy diffusivity K Z for the numerical outputs and observations are shown in Fig. 4 vs. nondimensional depth z/MLD.I n the surface layer the numerical values are often more variable than the in situ ones.The opposite is observed below the MLD.
Among the four basic model's configurations (KEset, KEflu, KLset and KLflu), near the surface, KEflu shows the best agreement with the observations among the simulations with the lower threshold on k.KEflu is not as good as KEset, but it is significantly improved by using the Canuto A stability functions.KLset and KLflu perform similarly but have an error in opposite directions.Near the surface, KLset is higher than KLflu.When a higher threshold on k is applied, KLsetMINk is closer to the observations than KLfluMINk.
Below the MLD, the numerical experiments appear to be grouped on the base of the threshold on k:n od i fferences between different SMC and boundary conditions appear obvious.Noticeably, below the MLD, the numerical experiments with a higher threshold on the kinetic energy k are more compatible with the in situ data than the other numerical experiments.The median value of K Z below the MLD appears to have a seasonal behavior.In particular, it is equal to × − 5.90 10 m /s 62 2 in spring, × − 3.02 10 6 in summer, × − 3.27 10 m /s 62 2 in autumn and × − 5.57 10 m /s 72 2 in winter.Overall, the numerical experiments that have the lower average error with respect to the measurements are the four numerical experiments with a high threshold on k, as the values of σ KZ in Table 2 in- dicate.

In situ and numerical probability density functions
We now focus on the comparison of the in situ and numerical PDFs of ε, N and K Z values.
Above the MLD, the in situ dissipation rate data (in black in Fig. 5ab), show a clear bimodal behavior with one peak at ≈× − ɛ 51 0m / s 92 3 and one at ≈× − ɛ 51 0m / s 62 3 .The numerical distributions show this feature in a less marked way and are generally biased low with respect to the in situ data's PDF.As values of S ε and ∆ ɛ 2 in Table 3 show, above the MLD, the distribution more aligned to the in situ data one are KE-fluMINk and KEfluMINkOpt, while the distribution with the shape more resembling the in situ data one is KEfluMINkOpt.
Below the MLD (Fig. 5c-d), all the numerical PDFs have mean values similar to the ones of the corresponding PDF in the mixed layer, but they all are much more narrow than the in situ data.KL numerical experiments have higher mean values than KE ones and the numerical experiments with a threshold on k have higher mean values than the corresponding numerical experiments without threshold.As we see in Table 4, below the MLD, KLsetMINk and KLfluMINk have the same shift   with respect to the in situ distribution, while KEfluCANMINkOpt and KLfluCANMINkOpt have the lower value of ∆ ɛ 2 .But we note that the differences between the values of ∆ ɛ 2 below the MLD are due to small differences in the low right-end tails that are unlikely to be significant.
The probability density functions of the in situ stratification data (in black) and of the numerical data are shown in Fig. 6.Both above and below the MLD, there is no drastic difference between the numerical distributions in both shape and peak.The similar behavior of all the numerical experiments can be due to the low frequency nudging to the MERCATOR fields.Overall, above the MLD (Fig. 6a-b), the numerical experiments overestimate the stratification.Note that this is not evident by just looking at the median profiles in Fig. 3.As values of ∆ N 2 in Table 3 indicate, the distribution more resembling the in situ one is KEfluCANMINkOpt.
Below the MLD (Fig. 6c-d) the peaks of the numerical PDFs are less aligned than above the MLD.In this layer, the numerical distributions show a peak at ≈× −− N 51 0 s 3.3 1 that is not present in the measurements.As values of ∆ N 2 in Table 4 show, KEfluCANMINk is the distribution with the shape more resembling the in situ one above the MLD, while KEfluCAN is the more resembling below the MLD.However, we note that, when the secondary peak is removed, the numerical experiments with the better PDF of N are KEfluCANMINk and KEfluC-ANMINkOpt.
The probability density functions of the in situ eddy viscosity data (in black) and of the numerical data are shown in Fig. 7. Above the MLD (Fig. 7a-b), the in situ values have three peaks with a higher peak at ≈× − K 51 0m / s  3 indicate that KEfluCANMINkOpt is the distribution with the shape most similar to the observations above the MLD.The effect of raising k min is similar on the two SMC families: it shifts the distribution towards higher values and modifies the peaks' height.However, KLsetMINk and KLfluMINk show a prominent first peak and thus perform less well than KEfluCANMINk and KEfluCANMINkOpt.
Below the MLD (Fig. 7c-d 4 indicate that KEfluCAN, KEfluCANMINk, KEfluC-ANMINkOpt and KLfluMINk are the distributions with the shape more resembling the in situ one.But, as we already noted for the distributions of ε below the MLD (Fig. 5), in this layer there is not a real difference between the numerical experiments.

(ε,N ,K Z ) space
All the measurements of ε, N and K Z and the corresponding numerical data are shown in Fig. 8.The point clouds reveal that the in situ values are much less disperse than the model ones (not evident from the picture; see Supplementary Data).In Fig. 8a-b we can see six main planes corresponding to: KEfluCAN, KEset and KEflu, KLset, KLflu and some in situ data at low K Z .Indeed, these last points lie on a plane with constant = kk min as prescribed by the parameterization we used (see Methods and Discussion).Data belonging to the numerical experiments with a higher threshold on k (KEfluCANMINk, KEfluCANMINkOpt, KLsetMINk and KLfluMINk) are first found at  differentiation between above and below the MLD is not shown in Fig. 8).

Discussion
By considering the median profiles, the values of logarithmic standard deviation and the shape of the probability density functions, we can progress towards the identification of an optimal Symphonie's configuration with respect to our dataset.
From visual inspection of the median profiles it is apparent that the default threshold on k (used for KEset, KEflu, KLset, KLflu) underestimates the turbulent activity (Fig. 5a-b).By using Eq. ( 1) and a value of ≈ − ɛ 10 m /s 42 3 that we can infer from the in situ data below the MLD (Fig. 2), we estimated the value of a more appropriate threshold ( ≈ − k 10 m /s 72 2 ) that was used for the numerical experiments with the higher threshold (KEfluMINk, KEfluMINkOpt, KLsetMINk and KLflu-MINk).As Figs. 2 and 6 and the values of σ ε and σ KZ in Table 2 indicate, this new threshold permitted to get better estimates of both ε and K Z .
The decadal standard deviation values are suited for whole-column comparisons between numerical experiments and in situ data, when it is interesting to have the general behavior of the numerical experiment (e.g., Burchard et al., 2002).However, we are also interested to details of the model's configuration that likely have a major effect only in the MLD and cannot be appreciated by this metric.
The comparison of the shape of the numerical and in situ probability density functions is more appropriate to examine the model's behavior above the MLD.In particular, it permits to investigate its response to the different configuration choices on the full range of values of ε, N and K Z .For example, Fig. 5a-b shows that the effect of raising the threshold value of k is not just simply a shift, as is suggested by the σ ε values and the median profiles (Fig. 2).In fact, the numerical experiments with a higher threshold (KEfluCANMINk, KEfluCANMINkOpt, KLsetMINk, KLfluMINk) are not merely shifted towards higher values but also show significant differences in the peaks' heights.This is due to two facts: (i) with higher values of kinetic energy, it is generally more probable to observe higher values of dissipation, and (ii) the highest values in each of the ε PDF' peaks cannot be arbitrarily large because the left-hand peak mainly contains values below the MLD and viceversa for the right-hand peak (data not shown).
Regarding the average value of the MLD estimated by the different numerical experiments, we note that it is systematically lower than the observed one.Generally, the numerical MLD coincides with the depth at which the kinetic energy reaches its minimum value (Fig. S2 in the Extra Materials).Possible reasons for the mismatch could be some missing or not well parameterized surface turbulent phenomena (e.g., Langmuir turbulence), as well as the uncertainty on the drag coefficients.
The analysis of ε data above the MLD allows examination of the effect of boundary conditions, stability functions and optical scheme.The effect of the surface boundary conditions is readily seen by visually inspecting the ε median profiles (Fig. 2).The KEset, KEflu, KLset and KLflu profiles highlight that, at the surface, the equilibrium boundary condition determines ε values lower than both the observations and the values of the corresponding numerical experiments employing a flux boundary condition.Note that the values of ε close to the surface that one obtains using a given boundary condition also depend on the employed SMC.In fact, the − k ɛ and the − k ℓ schemes have different values of the constants in the surface boundary condition for the kinetic energy (see Appendix A).
As a consequence, an influence of the boundary conditions is also found in the PDFs of ε (Fig. 5).In fact, the value of ε in the upper tail of the PDF is different for the schemes employing the equilibrium boundary condition and the ones employing the flux boundary condition as one would expect as the higher values near the surface are predominantly affected by the surface boundary conditions rather than by the threshold on k.Indeed, the right peak at higher ε values is mainly due to values from the depth range 0 ≤ z/MLD ≤ 0.5.
Regarding the average values of the MLD, we find that all the numerical experiments predict a shallower mixed layer depth.A reason for this behavior could be the uncertainty in the drag coefficient and/or to the presence of some not well-parameterized mixing phenomenon.
Looking at the profiles of turbulent production terms P and B in Eqs.(B.1) and (B.2) (Figs.S3 and S4 in the Extra Materials), we find that turbulence is mainly due to shear.In particular, we find that shear is enhanced at the mixed layer depth.This indicates that the increased levels of turbulence right below the MLD that we observe in situ (Fig. 2) are due to internal waves s or intermittent Kelvin-Helmholtz instability within the stratified layer (Woods, 1968;Grant and Belcher, 2011).The fact that the model does not reproduce this enhanced turbulence is due to the too coarse resolution that does not permit resolving such processes.
Secondly, the effect of the stability functions is made clear by both the median profiles and the PDFs of ε.I nFig. 2 it appears that the Canuto A stability functions improve the agreement with respect to the observations (cf.KEflu vs. KEfluCAN), especially near the surface.The same information is given by the PDFs (Fig. 5; Table 3): KEfluCAN's PDF has just two peaks as the in situ PDF and, thus, it is in better agreement with the observations.This result did not change when we raised the threshold on the kinetic energy.
Thirdly, the effect of the optical scheme is best seen by looking at the ε PDF (Fig. 5).Indeed, we see that the modified optical scheme enables KEfluCANMINkOpt to reproduce remarkably well the in situ distribution.We attribute this significant improvement to the better estimate of the temperature and salinity profiles (see Figs. S5 and S6 in the Extra Materials) that, in turn, improves the stratification.Indeed, the profile of N (Fig. 6)ofKEfluCANMINkOpt is slightly better than the one of KEfluCANMINk.Even if this is not captured by the values of S N (Table 3), we note that, right below the MLD, KEfluCANMINkOpt better reproduces the peak that we observe in the in situ stratification.However, it is not straightforward to say how this can influence the dynamics above the MLD.(indeed the KEfluCANMINkOpt's N PDFs are totally comparable with the ones of KEfluCANMINk).The only obvious difference that we observe above the MLD between the KEfluC-ANMINk's and KEfluCANMINkOpt's stratifications is that the median profile of the buoyancy turbulent production term B of KEfluC-ANMINkOpt is significantly higher than the one of KEfluCANMINk (Fig. S4 in the Extra Materials).
The analysis of the PDFs of N indicates that all the numerical experiments overestimate the stratification above in the surface layer (Fig. 6a-b).On the other hand the stratification is well reproduced by all the numerical experiments below the MLD.The only absent feature in the numerical data of N is a bump of high values right below the MLD where the entrainment takes place and where we also observe higher levels of turbulence (see Fig. 2).Overall, below the MLD the numerical PDFs of ε and K Z markedly differ from the in situ ones: they are more peaked than the in situ PDFs.A possible cause for this behavior could be the too coarse resolution of the model that does not permit to adequately resolve turbulent processes below the MLD or to the different natures of mixing mechanisms.
Also the analysis of the K Z probability density functions leads to the conclusion that KEfluCANMINkOpt is the better configuration.Even though the values of ∆ K 2 Z in Table 3 indicate that there is not a clear difference between KEfluCANMINk and KEfluCANMINkOpt, we note that KEfluCANMINkOpt K Z PDF has three peaks as the in situ PDF whereas KEfluCANMINk's PDF has only two.In general, the − k ɛ scheme performs significantly better in reproducing the in situ PDF of K Z .This is likely due to the advantage of − k ɛ in estimating ε with respect to − k ℓ.The inclusion of the molecular diffusivity D T in the formulation of K Z did not cause significant differences, probably because the turbulence level is too high to allow to appreciate the difference.Here we did not show a numerical experiment in which we change only this aspect of the model's implementation.However, one example is given in the Extra Materials (Fig. S7).
The analysis of the data below the MLD shows a sharp disagreement between the numerical and the in situ values of ε and K Z .However, we can observe three main features: i) the typical value of numerical ε depends on the SMC that one employs; ii) the model values of ε and K Z are much less variable than the measurements; iii) setting a hard threshold on the kinetic energy limits the reliability of the modeling in this layer.
The first point is made clear by Fig. 2. In fact, we can see that when the same threshold on k is applied, the ε values obtained with the − k ɛ scheme are one order of magnitude smaller than the ones obtained with the − k ℓ scheme.This is partially explained by the fact that the two SMCs use different values for the constants in Eqs. .However, this fact alone cannot explain the one-order-of-magnitude difference that we observe below the MLD.
The second point is illustrated by Figs.5c-d and 6 c-d.Therein, the numerical PDFs in the lower layer are much less spread out than the in situ ones.A possible cause for this could be the resolution of the model that does not permit to adequately resolve the turbulent processes in this layer.
Also, a possible cause for the sharp discrepancy between numerical and in situ PDFs below the MLD could be the use of a hard threshold for k.In fact, Fig. 8 suggests that this choice leads to a non-physical result, with many low-energy points on a single plane in the space (K Z ,N , ε).This behavior results directly from Eqs. for the KE scheme (where c k ,c ε and c 0 are constants).This is exactly the behavior we observe in Fig. 8b, with the planes following the line . Thus, the points on the planes in Fig. 8 are values with a minimal kinetic energy.In situ data do not show this behavior, suggesting that specifying a hard value for k min is likely not the best choice.A more suitable approach would probably be a parameterization of k min on the base of the different turbulent processes at play.For example, Gaspar et al. (1990) hypothesized the use of a parameterization on the internal wave field activity when dealing with data in a depth range similar to the one in our analysis.Note here that such a solution is likely to depend on the amount of numerical mixing and to the processes that are resolved by each model.

Conclusions
Vertical mixing in the surface layer of the ocean plays an essential role in both physical and biochemical phenomena.Therefore, it must be correctly estimated by numerical models.Different turbulence closure schemes have different performances in predicting mixing.Moreover, the problem is made more complicated by the interplay of the SMC with other aspects of ocean dynamical numerical modelling like boundary conditions, stability functions and optical scheme.The influence of the minimum value of kinetic energy that is allowed has also been investigated.
We performed nine numerical experiments and compared their estimates of the turbulent quantities ε, N and K Z to in situ microstructure measurements in the Gulf of Lion.In particular, two SMCs were considered: a − k ɛ scheme proposed by Burchard and Bolding (2001)  (herein KE); and a − k ℓ based on Gaspar et al. (1990) (herein KL).We considered two surface boundary conditions: one supposing equilibrium between the production and dissipation terms in the dynamic equation for k (herein set); and one taking into account the effect of breaking waves of all scales based on Craig and Banner (1994) (herein flu).We considered the stability functions proposed by Kantha and Clayson (1994) and Canuto et al. (2001).In addition, we also considered two different attenuation lengths for the photosynthetically available radiation in the model's optical scheme.The combinations of these factors that we explored are resumed in Table 1.
A recent study by Reffray et al. (2015) showed that, in a 1-D situation, the KE closure scheme gives better mixing estimates than other widely used SMCs.However, that study did not compare KE to a mixing length scheme based on Gaspar et al. (1990).Our study fills this gap.Moreover, our study exploits a 3-D model that, by definition, considers more terms in the dynamical equations than what 1-D models do.Additionally, it compares the modeled turbulence activity directly to microstructure turbulence measurements rather than with derived quantities like the mixed layer depth.
The two SMCs do not show relevant differences when estimating K Z both above and below the MLD.The fact that the two SMC perform similarly in estimating K Z when KL has an advantage in estimating ε,is not explained by differences in estimating N and need further investigation in the future.
Previous studies, in 1-D and 3-D numerical simulations, gave conflicting indications on the effect of the stability functions proposed by Kantha and Clayson (1994) and Canuto et al. (2001) on the KE scheme.We find that the effect of the Canuto A stability functions on the KE is to improve the performance of the closure scheme when the PDFs of ε and K Z are considered.The value of the kinetic energy threshold plays a pivotal role in approaching the observations below the MLD (see Figs. 2 and 4).Nonetheless, we found a non physical behavior of the numerical experiments for low kinetic energy levels.This supports the idea that the minimum of kinetic energy should rather be parameterized as function of different turbulent mechanisms rather than being a hard threshold (as already hypothesized by Gaspar et al., 1990).Such a parameterization should also depend on the SMC that is employed.Indeed, we found that the two SMC predict different turbulence levels below the MLD even when they employ the same threshold on the turbulent kinetic energy.
We found that the attenuation length of the photosynthetically available radiation plays an important role in determining the stratification and, as a consequence, the performance of the model in predicting ε and K Z .This result highlights the importance of biological activity in influencing physical processes.
Our study shows that the comparison with in situ microstructure data can effectively help in setting up the implementation of a SMC for an ocean numerical model.In the future, we should expect more utilization of microstructure data as new automated instrumentation becomes available.In particular, microstructure and wave height probes mounted on autonomous platforms, such as drifting profilers or gliders, will permit new and bigger datasets, especially during adverse meteorological situations.This kind of data will provide essential data to developers of future ocean-wave coupled models, as these are expected to significantly advance the ocean modelling state of the art by reducing the need of general surface boundary conditions.Herein, S M is obtained by supposing a null stratification (i.e., = G 0 h in Eq. ( 24) in Galperin et al., 1988).Under this hypothesis, S M becomes constant.In the general case, the expression for S M are given by Kantha and Clayson (1994) and Canuto et al. (2001).The value of the stability coefficient c 0 for the − k ɛ scheme is issued from Warner et al. (2005).The value of c 0 for the − k ℓ scheme is issued from Gaspar et al. (1990).

A1.2. Flux
Alternatively, the boundary conditions can be specified as surface flux conditions, namely: where the surface flux can be computed as = Fτ ρ 100( / ) 0 3/2 (Craig and Banner, 1994) or directly prescribed from the 'wave to ocean' turbulence flux computed by a wave model when available.

A2. Bottom and surface boundary conditions for ε
The ε surface and bottom conditions are computed on the first level under the surface and above.Let z 1 denotes the distance between this level and the considered boundary.Boundary conditions for ε are obtained from k and Eq.(B.2), using the latter with some appropriate hypothesis for l B a boundary length scale value.A simple formulation (Warner et al., 2005) is eventually given by =+ lz z 0.4( ) B 10 , where z 0 is a length scale representing the roughness of the bottom boundaries.Unfortunately, this formulation potentially leads to unrealistic high values when the underlying hypothesis of neutral stratification is no longer valid (a situation that is more likely to occur in deep zone where the vertical grid resolution near the bottom is generally coarse).One way to solve this problem is to introduce a dependency on the gradient Richardson number (Estournel and Guedalia, 1987;Michaud et al., 2012) Following Burchard and Bolding (2001), the equations describing the dynamics for k and ε are: are the production terms due to shear and buoyancy respectively.Parameters values issued from Warner et al. (2005).
The eddy viscosity K Z and the temperature eddy diffusivity K T used in Eqs.(B.3) and (B.4) are given by: The turbulent length l is related to k, TKE and ε according to S M and S H , the quasi-equilibrium stability functions defined by Kantha and Clayson (1994), depend on the Richardson flux number.

B2. − k ℓ scheme
We used the − k ℓ closure scheme proposed by Gaspar et al. (1990).Therein the authors assume ε to be (following Kolmogorov, 1942): where ρ is the water density.
As highlighted by Gaspar et al. (1990), these length scales have a straightforward physical interpretation: they are the distances traveled upward/ downward by a fluid particle by converting all of its kinetic energy into potential energy.
Hypothesizing the turbulence to be stationary and homogeneous Gaspar et al. (1990) show that, in stably stratified regions, the model parameterizations yields: In order to apply the Batchelor's theory to the temperature gradient one has to divide the data into segments in which turbulence can be considered homogeneous, viz. the turbulent motion can be regarded as a random motion which average properties are independent of position in the fluid.Given the fact that the first aim is to fit the data to a theoretical Batchelor's spectrum, arguably the best segmentation method is the one giving the fits with the best quality.
To determine the fit quality we use a Maximum Likelihood Estimation technique.Following Ruddick et al. (2000), we calculated the goodness of a fit by maximizing the joint probability of the measured spectrum with respect to the theoretical Batchelor spectrum.Hence, the higher the joint probability, the higher the fit goodness.In the following we call this joint probability JP (corresponding to C11 in Ruddick et al., 2000).
To make the fit goodness criteria more rigorous, following Sanchez et al. (2011) we also require that: i) the mean absolute deviation (MAD) of the ratio between the in situ and theoretical spectra within the fitting domain be lower than 1.1; ii) the signal to noise ratio (SNR) be lower than 1.3; iii) the likelihood ratio (LHR) -which quantifies if the measured spectrum fits the Batchelor's spectra or a power-law spectrum better-be lower than 2.
We tested which of the following segmentation methods permitted to have the better fits in our dataset: a segmentation based on an eight order AR model by Imberger and Ivey (1991)  TIn Fig. 9 the histogram of the JP values is plotted with the cumulative distribution of the latter.Note that each histogram is normalized to the total number of segments obtained with the corresponding segmentation method.As it is evident by eye, the M12 is the segmentation method that gives the highest percentage of good profiles, as highlighted by the higher histogram bars on the right end of the distribution with respect to the other cases.Moreover this is the method that has the vertical line lying at the highest value, meaning that a higher percentage of segments have a better fit.Furthermore cases b), e) and d) show that progressively reducing the segmentation window gives better results.Case f) also shows, when compared to c), that adding overlap, though permitting to have a higher resolution, does not augment the percentage of good fits.These evidences suggest to apply a 128 points ( ≈ 12.5 cm) segmentation with no overlap.Anyway, not needing such a high resolution in order to compare the in situ data with the numerical data, we avoided the overlap.The distribution of JP values of this segmentation method did not significantly differ from the one of M12 (data not shown).

Fig. 1 .
Fig. 1.Numerical model domain.The color code represents the water depth.The Gulf of Lion is magnified in the smaller box where the measurements sites are represented by red dots.Note that many profiles were taken at the same location over time.The black lines in the smaller box represent the 0, 50, 100, 500, 1000 and 1500 m isobaths.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) ρb .For high Re b (Re b > 400), turbulence is in an energetic regime, where = K νR e 4 ρb .Therefore, for high Re b the buoyancy flux K ρ N tends to vanish together with stratifi-

Fig. 2 .
Fig. 2. The thick lines represent median values of ε estimated from in situ data (in black) and numerical experiments.The shades indicate 95% bootstrap confidence intervals.Water depth z is non-dimensionalised with respect to the mixed layer depth MLD for each profile.

Fig. 3 .
Fig. 3.The thick lines represent median values of N estimated from in situ data (in black) and numerical experiments.The shades indicate 95% bootstrap confidence intervals.Water depth z is non-dimensionalised with respect to the mixed layer depth MLD for each profile.

Fig. 4 .
Fig. 4. The thick lines represent median values of K Z estimated from in situ data (in black) and numerical experiments.The shades indicate 95% bootstrap confidence intervals.Water depth z is non-dimensionalised with respect to the mixed layer depth MLD for each profile.

Fig. 5 .
Fig. 5. a) and b): ε probability density functions issued from all the in situ (black shade) and numerical data (thick lines) in the surface mixed layer for the KE closure scheme and the KL closure scheme respectively.c) and d): the same but below the surface mixed layer.

Fig. 6 .
Fig. 6. a) and b): N probability density functions issued from all the in situ (black shade) and numerical data (thick lines) in the surface mixed layer for the KE closure scheme and the KL closure scheme respectively.c) and d): the same but below the surface mixed layer.
), similarly to what seen for ε, the numerical distributions are clearly different from the in situ one.∆ K 2 Z values in Table planes with an average value of ε that is lower than the average value of ε of the KL numerical experiments with = − k 10 m /s min 82 2 .These last ones have lower average values of ε than the numerical experiments with = − k 10 m /s min 72 2 .This closely mirrors what we already noted in Fig.5c-d.But it turns out that the points lying on the planes in Fig.8do not belong exclusively to depths greater than the MLD (the

Fig. 7 .
Fig. 7. a) and b): K Z probability density functions issued from all the in situ (black shade) and numerical data (thick lines) in the surface mixed layer for the KE closure scheme and the KL closure scheme respectively.c) and d): the same but below the surface mixed layer.Note the difference in the range of the y axes in a) and b).

Fig. 8 .
Fig. 8. a) All in situ (in black) and numerical data of K Z ,Nand ε. b) All in situ (in black) and numerical data of K Z and ε.The color code is the same as in the other figures.
(B.1)  and (B.2) that relate ε to the kinetic energy.Indeed, − k ɛ uses a lower value ( wind stress and ρ 0 is the reference density of sea water.By inserting this last equation in the expression for the shear production term (Eq.(B.1)) and using the expression for ε proposed byGalperin et al. (1988) (Eq.(B.1)), we obtain the surface boundary condition for k: and Lacarrere, 1989).The eddy momentum diffusivity is related to TKE according to: has to be determined.The dissipation and mixing length scales l k and l ε are the ones proposed byBougeault and Lacarrere (1989): (II91); a constant segmentation of 1024 data points in which the values of temperature gradient variance in 7 sub-segments are in the same order of magnitude proposed by Sanchez et al. (2011) (S11), a constant segmentation of 1024 data points with no overlap suggested by Cuypers et al. (2012) (C12); a constant segmentation of 128 data points with 50% overlap suggested by Moniz et al. (2012) (M12); a constant 512 data points segmentation with no overlap; and a constant 1024 data points segmentation with 30% overlap.In practice, we fitted the segmented data to the Batchelor's spectrum following the fitting procedure of Steinbuck et al. (2009) and we compared the results searching for the distribution of JP values with the greater proportion of high values.The results are depicted in Fig. 9.

Fig. 9 .
Fig. 9. a) JP values distribution obtained with the II91 segmentation method, b) same as a) for S11, c) same as a) for C12, d) same as a) for M12.e) same as a) for a constant 512 points segmentation.f) same as a) for a 1024 points segmentation with 30% overlap.The red curve is the cumulative distribution of JP values for each case and the vertical line marks the point where the cumulative distribution equals 0.6.The database for this analysis comprehended 126 profiles with an average depth of 50 m in various meteorological conditions in the Gulf of Lion.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1
Set of numerical experiments performed for this study.

Table 2
MLD ratio indicates the average ratio between the numerical prediction of the MLD and the in situ value.σ ε indicates the mean value of the decadal standard deviation of the ε profiles.σ N indicates the mean value of the decadal standard deviation of the N profiles.σ K Z indicates the mean value of the decadal standard deviation of the K Z profiles.

Table 4
None of the numerical experiments reproduces this behavior.As values of σ KZ in Table 2 suggest, above the MLD the distribution most aligned to the in situ one is KLsetMINk.Besides, ∆ K Z 42 2 .