Biogeosciences Characterization of turbulence from a fine-scale parameterization and microstructure measurements in the Mediterranean Sea during the BOUM experiment

One of the main purposes of the BOUM experiment was to find evidence of the possible impact of submesoscale dynamics on biogeochemical cycles. To this aim physical as well as biogeochemical data were collected along a zonal transect through the western and eastern basins of the Mediterranean Sea. Along this transect 3-day fixed point stations were performed within anticyclonic eddies during which microstructure measurements of the temperature gradient were collected over the top 100 m of the water column. We focus here on the characterization of turbulent mixing. The analysis of microstructure measurements revealed a high level of turbulent kinetic energy (TKE) dissipation rate in the seasonal pycnocline and a moderate level below with mean values of the order of 10 −6 W kg−1 and 10−8 W kg−1, respectively. The Gregg Henyey (Gregg, 1989) fine-scale parameterization of TKE dissipation rate produced by internal wave breaking, and adapted here following Polzin et al. (1995) to take into account the strain to shear ratio, was first compared to these direct measurements with favorable results. The parameterization was then applied to the whole data set. Within the eddies, a significant increase of dissipation at the top and base of eddies associated with strong near-inertial waves is observed. Vertical turbulent diffusivity is increased both in these regions and in the weakly stratified eddy core. The stations collected along the East–West transect provide an overview of parameterized TKE dissipation rates and vertical turbulent diffusivity over a latitudinal section of the Mediterranean Sea. Strong TKE dissipation rates are found within the first 500 m and up to 1500 m above the bottom. Close to the bottom where the stratification is weak, the inferred vertical turbulent diffusivity can reach Kz ' 10−3 m2 s−1 and may therefore have a strong impact on the upward diffusive transport of deep waters masses.


Introduction
During the last two decades, increasing evidence has shown that vertical transport is a key factor controlling biogeochemical fluxes in the ocean.These fluxes need to be accurately quantified in order to adequately characterize biogeochemical processes and their representation in numerical models (Lewis et al., 1986;Denman and Gargett, 1983;Klein and Lapeyre, 2009).Two main processes account for vertical transport: upwelling resulting from divergent Ekman transport and turbulent mixing.Wind shear and convection are the main sources of turbulent mixing in the mixed layer, whereas breaking internal waves are responsible for most of mixing in the stratified ocean (Munk and Wunsch, 1998;Thorpe, 2004).An adequate representation of these processes is notably needed in oligotrophic regions where vertical transport generates an upward nutrient flux that directly sustains primary production in the depleted euphotic zone.
Among oligotrophic environments, anticyclonic eddies have been the subject of several studies since processes intrinsic to the eddy dynamics locally enhance the vertical transport of nutrients in the euphotic layer (McGillicudy et al., 1999;Ledwell et al., 2008).In these regions upward Published by Copernicus Publications on behalf of the European Geosciences Union.
doming of the seasonal thermocline (and downward doming of the permanent thermocline), known as eddy pumping, generates an uplift of a nutrient-enriched deep layer, a process which can be enhanced from the secondary circulation generated by interactions between wind driven surface currents and the eddy motion (Martin and Richards, 2001).Finally, Ledwell et al. (2008) suggest that increased shear and mixing could result from near-inertial wave trapping.Indeed, the negative vorticity ζ of the anticyclonic eddy can shift the effective inertial frequency to a lower value f eff = f + ζ /2 (Kunze, 1985); therefore, near inertial waves which evolve in the frequency band f > f eff will encounter their turning points when propagating away from anticyclonic eddy centers (Bouruet-Aubertot et al., 2005) and remain trapped in the eddy core.
The Mediterranean Sea is an oligotrophic environment in which mesoscale dynamics are important and anticyclonic eddies are predominant (Moutin et al., 2011).Regarding internal wave energy sources, the Mediterranean Sea is also a unique region since internal wave energy is expected to be mainly due to the atmospheric forcing on account of the weak tidal forcing.A large part of the BOUM experiment was dedicated to the precise characterization of biogeochemical processes within three oligotrophic environments, namely anticyclonic eddies (eddies A, B, C) (Moutin and Prieur, 2012).During the experiment special effort was made to determine the physical forcing, and specifically the vertical mixing in these environments, which impacts most of biogeochemical processes studied and modeled within the BOUM project (Bonnet et al., 2011;Mauriac et al., 2011).Another point of interest is the impact of vertical turbulent mixing on the water masses evolution and water massed circulation in the Mediterranean Sea.A number of studies have focused on the complex problem of water masses circulation and transformation in the Mediterranean Sea; see for instance Lascaratos et al. (1999) for a review as well as Moutin and Prieur and Touratier et al. (2012) for the characterization of water masses during BOUM.However the vertical turbulent mixing, a key process for the general circulation of water masses, has never been quantified experimentally at the basin scale from fine-scale parameterization.
Estimates of vertical mixing from in situ measurements are very scarce in the Mediterranean Sea.In a pioneering study, Woods and Wiley (1972) reported some estimates of mixing based on temperature microstructure measurements in the upper 100 m of Malta coastal waters.Very recently Gregg et al. (2012) performed intensive microstructure measurements on the Cycladic plateau.Mixing processes reported in these studies are quite specific to the continental slope.In this paper we focus on the characterization of turbulent dissipation and mixing rates in deeper areas of the Mediterranean Sea.Turbulent dissipation rates are characterized from microstructure measurements in the upper 100 m and are favorably tested against parametrization of energy dissipation based on fine-scale internal wave shear and strain.This pa-rameterization is then applied in order to characterize vertical mixing within the full depth range of the eddies and the effect of possible internal wave trapping is discussed.Finally, vertical mixing is inferred from deep fine-scale measurements all along the Mediterranean East-West transect, thus providing a first insight into the internal wave mixing rates at the basin scale.From this transect, the location of mixing hot spots and their possible impact on the Mediterranean overturning circulation is also briefly discussed.

Hydrographic and current measurements
Conductivity temperature depth (CTD) measurements were performed for each station using a SeaBird SBE911 instrument.Data were averaged over 1 m bins to filter spurious salinity peaks.Simultaneously currents were measured by a 300 kHz lowered broadband acoustic current profiler (LADCP).LADCP data were processed using the Visbeck inversion method and provided vertical profiles of current velocity at 8 m resolution.A total of 30 CTD/LADCP profiles were collected along the BOUM transect, down to the bottom during short duration stations (SD stations hereafter) at a horizontal spatial resolution of 100 km (see Fig. 1).In addition, intensive measurements (every 3 h over 3 days) were realized down to 500 m depth for each of the 3 long duration (LD stations hereafter) stations within anticyclonic eddies A, B and C (Fig. 1).

Dissipation measurements and K z estimates
For each LD station, repeated profiles with a temperature gradient microstructure profiler, Self-Contained Autonomous Microstructure Profiler (SCAMP; Precision Measurements Engineering, www.PME.com),enabled estimates of the turbulent kinetic energy dissipation rate ( ) and the vertical diffusivity of temperature (K z ) to be made (e.g.Ruddick et al., 2000).SCAMP is limited to 100 m depth and operates at a slow optimal free fall velocity of 0.1-0.2m s −1 .SCAMP was deployed only under calm weather conditions (low wind and swell) and between the CTD/LADCP profiles, therefore the total number of SCAMP profiles was limited to 21. Estimation of the dissipation rate is based on Batchelor curve fitting of the temperature gradient spectrum (Ruddick et al., 2000).The maximum likelihood method of Ruddick et al. (2000) implemented in SCAMP software is used for the curve fitting.However we customize the algorithm by: -including the improvement on the estimation of χ T ( • C 2 s −1 ) (the rate of destruction of temperature variance by molecular diffusion) proposed by Steinbuck et al. (2009).
-switching to the least square fit method of Luketina and Imberger ( 2001) when χ T is smaller than the noise variance.The Ruddick et al. (2000) method includes a model of the noise as a part of the model spectrum that is fitted to the experimental curve.For very low χ T the temperature gradient spectra is dominated by noise, and any small inaccuracies in the noise model results in large errors in the estimation of χ T .In the Luketina et al. (2000) method, the high wavenumber part of the spectrum dominated by noise is simply discarded; this appeared to improve largely the results when noise variance was larger than χ T .

Biogeosciences
We estimate the rate of cross-isopycnal turbulent mixing or diapycnal diffusivity, K turb , following Shih et al. (2005), applying their new parameterization for flows characterized by turbulent intensities larger than 100: ; (1) with the Osborn relation (Osborn, 1980) for turbulent intensities in the range between 7 and 100, and assuming turbulent diffusion is null for turbulent intensities less than 7.Note that a background molecular diffusion κ ρ is always present; since density is mainly set by temperature and since molecular diffusion of heat κ T is much larger than molecular diffusion of salt, we simply set κ ρ = κ T = 1 × 10 −7 m 2 s −1 .The final expression of diapycnal diffusivity including both turbulent and molecular diffusion therefore reads K z = K turb + κ T .

A fine-scale parameterization
In the absence of microstructure measurements, energy dissipation is commonly inferred from fine-scale parameterization which relates the characteristics of the internal wave field to turbulent kinetic energy dissipation.These parameterizations assume a steady state spectrum of internal waves where wave-wave interactions transfer energy from large to small scale motion.Internal waves eventually break when they reach a critical wavenumber producing turbulent dissipation of the internal wave energy (e.g.Gregg, 1989;Polzin et al., 1995).Away from the boundaries, in the stratified ocean interior, the internal wave continuum is well represented by the steady state Garret and Munk (GM hereafter) spectrum (Garrett and Munk, 1979).Therefore parameterizations of dissipation rate ( ) were developed assuming a Garrett and Munk spectrum for the internal wave field.Henyey et al. (1986) used a ray tracing approach to simulate interactions of small amplitude test waves in a background internal wave field with a Garrett and Munk spectrum.From these simulations they determined a scaling of the dissipation rate on internal wave energy as , where E GM is the canonical Garret and Munk energy level.This scaling has received strong support from experimental observations (Gregg, 1989;Polzin et al., 1995).Gregg (1989) (G89 hereafter) proposed a popular incarnation of the original Henyey et al. (1986) parameterization in the form: where 0 = 6.7.10 −10 W kg −1 is the canonical GM dissipation rate at a latitude of 30 • for N = N 0 , N 0 = 3 cph is the canonical GM buoyancy frequency and j (f/N ) = f a cos h(N/f )/f 30 a cos h(N/f 30 ) is the latidunal dependance equal to unity at N = N 0 and f = f 30 , where f 30 is the coriolis parameter at 30 • .Note that in the final form used by Gregg (1989), j (f/N ) was approximated to unity, following recent works (Kunze et al., 2006;Gregg et al., 2003); we keep the full expression for j (f/N ).In Eq. (3) the original E 2 GM of Henyey et al. (1986) formulation was substituted for by GM in order to take into account possible deviations from the E GM level and to give an expression which is a function of shear (a quantity more conveniently estimated than energy).S 2 10 represents the observed shear variance with vertical wavelength greater than 10 m (or vertical wavenumber smaller than k crit = 0.6 rad m −1 ).This one is estimated from the 10 m velocity difference, with an additional √ 2.11 factor accounting for the attenuation from the first difference filter.S 2 GM is the corresponding shear variance for the canonical GM spectrum integrated up to the wavenumber k crit = 0.6 rad m −1 , which is given by (Cairns and Williams, 1976): with j = 3, b = 1300 m, and E GM = 6.3 × 10 −5 .The wavenumber k crit is assumed to represent the critical wavenumber for which the waves reach a critical Richardson number Ri crit = 0.5 and break.It is noteworthy to mention that this parameterization is able to reproduce the observed levels of dissipation within a factor of two for conditions close to the GM model (Gregg, 1989).
The main drawback of the G89 parameterization is that it fails to take into account possible variations of the shear to buoyancy scaled strain variance ratio, R ω , which is a measure of the internal wave field's aspect ratio (Polzin et al., 1995;Kunze et al., 2006).From our measurements we estimate R ω as: where η z is the strain computed as η z = (N 2 − N 2 )/ N 2 .For LD stations, .denotes a time average over the largest integer number of inertial periods resolved by the duration of the LD station measurements (3 for each long station); for isolated SD stations, .denotes an average over a 150 m depth bin.We apply the correction factor proposed by Polzin et al. (1995) (and applied by Kunze et al., 2006, andGregg et al., 2003, parameterizations): so that the parameterization we apply is given by: Because in the upper part of the water column z < 25 m (where z is the water depth) the shear is not properly measured by LADCP, following Kunze et al. (2006) we estimate shear from strain as S 2 10 = R ω N 2 η 2 z,10 , where η z,10 is a 10 m running average of the strain (equivalent to a 10 m first order finite difference of the isopycnal displacement) and R ω is the shear to buoyancy scaled strain variance ratio averaged over the lower part of the water column  m.Note that near the surface, for z < 25 m the results from the parameterization should be considered with caution due to the lack of shear measurements.Also, the GM spectrum is expected to be an adequate representation of the internal waves field "far" from boundaries and energy sources.Therefore strong deviations from GM are expected above the pycnocline, which has an average location at 15 m depth during BOUM.
Finally, we consider the effect of noise contamination on the LADCP measurements in the computation of shear.We determined the noise contamination by fitting the observed shear spectra with a composite GM and k 2 z shape.Both the energy level of noise and GM shaped spectrum are determined numerically from a fitting process.This allows us to determine the vertical wavenumber k noise at which the noise spectra intersects the fitted GM shape and strongly influences the measured shear.Figure 2a and b shows examples of mean shear spectra computed in the upper (z < 500 m) and lower (z > 500 m) parts of the water columns of deep isolated LADCP profiles performed during short duration stations.Spectra in the upper 500 m are marginally affected by noise in contrast with spectra at larger depths for which we determined k noise /2π = 2 × 10 −2 cpm due to a lower shear level.
Since the parameterization as formulated in Eq. ( 3) is not based on spectral estimates but rather on a local expression of the shear in physical space, we used a numerical Finite Impulse response Filter (FIR) with a cut off wavenumber k c = 0.5k noise to low pass filter the LADCP profiles and get rid of the noise at high wavenumbers.The filtered signal is then used to compute 10 m shear in Eq. (3).Because the high wavenumber energy above k c is removed from the experimental signal in this filtering process, we consider the GM expression for the shear variance for wavenumber smaller than k c by replacing k crit by k c in Eq. (4).A similar approach is used in Kunze et al. (2006), where both experimental and GM shear variance are computed using truncated spectra at k c .The underlying hypothesis in both cases is that the ratio of experimental to GM shear remains constant if both are low pass filtered at the same cut-off wavenumber.
Figure 2c, d and e shows examples of mean shear spectra for the three LD stations in eddies A, B and C covering the upper 500 m of the water column.The shear spectra are characterized by a non-GM shape with a large broad peak around k z /2π = 0.08 cpm.This peak is associated with the presence of strong near inertial internal waves in the eddies.In this case a fit was impossible, and we assumed that the noise level was the same as for SD stations in the upper 500 m.
As a final remark the value of the critical wavenumber k crit = 0.6 rad m −1 was prescribed for the GM spectrum following observations of Gargett et al. (1981).This value is applied in the G89 parameterization, but the actual value of the critical wavenumber depends on the energy of the observed internal waves field and may be significantly smaller than k crit for high energy levels.A criticism of G89 parameterization made by Gargett (1990) is that this one may underestimate the dissipation for high energy level because of an overestimate of the critical wavenumber.Polzin et al. (1995) suggest that significant deviations can occur if the actual critical wavenumber is smaller than 0.5k crit .However in this study wavenumber shear spectra from the first 500 m of the water column (Fig. 2a, c, d and e) do not show evidence of any roll off below 0.5k crit /2π=0.05cpm, so the (modified) G89 scaling can be applied there.For deep profiles shear spectra below 500 m are above the GM and the critical wavenumber is likely reduced.But for noise issues described above, the parameterization is considered for wavenumber smaller than 0.5k noise/2π = 0.01 cpm and no evidence of roll off is observed below this low wavenumber.
3 Observations: direct estimation of dissipation and validation of a fine scale parameterization

Stratification and dynamics within the three anticyclonic eddies
The three-day duration of the LD station sampling corresponds to 3.75 inertial period at station A (38 • N), and 3.37 inertial periods at stations B and C (34 • N).This sampling of the eddies allows us to characterize the background subinertial state as well as a large part of the internal wave band We define here the background state as the subinertial currents and the time-mean stratification averaged over 3 inertial periods.
The mean surface stratification presents a sharp pycnocline at ∼ 15 m depth (Fig. 3g-i Time-depth sections also illustrate temporal variability in the internal waves band and highlight a dominance of variability at the inertial frequency (e.g.Fig. 3).Oscillations with sloping iso-phases indicate baroclinic waves (as opposed to the barotropic signal which is constant with depth).Interestingly these waves are localized at the top and base of the eddy (e.g.Fig. 3a and d), leading to significant shear in these well stratified regions.
In order to characterize the internal wave spectrum, we have computed frequency shear spectra (Fig. 4).Note that only part of the internal wave range is resolved by our 3 h sampling profiles since the maximum frequency of these waves, the buoyancy frequency N, reaches values up to 0.05 rad s −1 (0.04 h period, also the spectral resolution is limited by the duration of the stations, it is of ±0.15 f for eddies B and C and of ±0.13 f for eddy A).A main peak around the inertial frequency is observed at stations A and C, which is consistent with the time depth plots of the currents described in Fig. 3.The peak is shifted to 0.8 f for eddy A in the first 100 m, which strongly suggests a modification of the effective inertial frequency by the negative eddy vorticity at this location, as will be discussed in Sect. 5.
A small peak around the semi-diurnal tidal frequency is also observed below 100 m depth for these stations.The shape of the spectra (Fig. 4) corresponds fairly well to that predicted by the GM model, the reference internal wave spectrum, although spectra at station B show a flatter slope.The spectra level is below the GM level in the upper layer for the three stations (red curves in Fig. 4) and slightly above the GM level at the base of eddies C and A (black curves in Fig. 4) where a strong near inertial signal is observed.Overall the experimental spectra are comparable to GM, and we shall test in the following section a fine-scale parameterization of energy dissipation that have been developed in this context of weakly nonlinear internal waves.

Dissipation measurements in the upper oceanic layer and comparison with a fine-scale parameterization
Figure 5 shows the first 50 m of individual profiles of dissipation recorded by the SCAMP during each LD station A, B, C and the background strain η z (see Sect. 2.3 for definition), as obtained from fine-scale measurements.
The strongly intermittent nature of is clearly apparent on these profiles, with values spanning several orders of magnitude [10 −11 , 5 × 10 −6 ] W kg −1 .A strong increase in dissipation rate [10 −8 , 5 × 10 −6 ] W kg −1 is observed between 10 and 20 m depth.This depth range typically corresponds to the variation of the pycnocline location due to internal waves heaving for the three stations (Fig. 5), and it will be referred to hereafter as the pycnocline region.The few values recorded above 10 m in the mixed layer were comparable to pycnocline values but were not considered further in the analysis because of the specific physics of the mixed layer (out of the scope of this paper) and possible ship contamination.Below 20 m depth, low values of dissipation (< 10 −9 W kg −1 ) are recorded with some sporadic events of high dissipations reaching 10 −6 W kg −1 .In Fig. 5, the strain appears clearly related to internal wave induced isopycnal displacement; this is most obvious for station C where the isopycnal displacement shows a dominant period at ∼0.85 day (considering 2 oscillations between day 178.7 and day 181.4) in close agreement with the inertial period of 0.89 day (Fig. 4).As observed for the dissipation rate, the strain values are generally maximum in the pycnocline region, which suggests internal wave strain importance in breaking processes, as already noted by several authors (Alford and Pinkel, 2000;Alford, 2010).Still, there is no systematic maximum of the dissipation rate associated with maximum strain.A similar situation was observed by Alford (2010) for tidal and near inertial internal waves in the Mendocino escapement, which suggests that dissipation results in this case from a cascading process as assumed by fine-scale parameterization of dissipation (Sect.2.3) rather than through direct breaking events of the dominant internal waves.
The question of the appropriate way to average such an intermittent variable as or K z using experimentally limited number of samples has long been debated (Baker and Gibson, 1987;Gregg et al., 1993;Davis, 1996;Gargett, 1999).We use here three estimates of the mean: (1) a simple arithmetic mean as suggested by Gregg et al. (1993), Davis (1996); (2) a geometric mean which is sometimes used in order to reduce the dispersion of dissipation rate data (Gargett, 1999;Smyth et al., 1997)  estimate (MLE, Priestley, 1981) of the mean dissipation following Baker and Gibson (1987).In this last case a lognormal distribution for is assumed and the mean dissipation rate is given by: where µ = log( ) and σ = std(log( )).
The statistical distribution of dissipation for all data and each LD station separately is represented as a probability density function (PDF) of in Fig. 6.The region within [10,20] m (the pycnocline region) and below 20 m depth were considered separately.The PDF were truncated below 10 −11 W kg −1 and above 10 −5 W kg −1 which represent the upper and lower bound of the SCAMP resolution.A log-normal distribution truncated at these resolution bounds was fitted to each PDF by maximum likelihood.An estimate of the mean was then obtained from the fit by Eq. ( 8).This MLE of the mean as well as the arithmetic and geometric mean and their confidence intervals are given in Table 1.
For all stations, the PDFs show two dynamical regions.For data in the pycnocline region, the most probable value (mode value) is ∼ 2×10 −7 W kg −1 , which characterizes highly dissipative processes, whereas it is close to a fairly low dissipation rate value of 10 −10 W kg −1 below.Below the pycnocline, the PDF is rather close to a log-normal distribution (a probability distribution of a variable whose logarithm is normally distributed) for station C, but is more skewed and spiky for stations B and A, which may result partially from a lack of convergence of observed PDFs.The lack of statistics in the   A possible log-normal MLE fit and consistent estimate of the mean was only obtained in this region when the whole data set was considered (Fig. 6a).
Overall the arithmetic mean and MLE fit mean are quite close below the pycnocline for all stations with values ∼ 10 −8 W kg −1 (Table 1).These mean values are almost two orders of magnitude larger than the most probable value, which illustrates the large intermittency of the data.The geometric mean largely underestimates the dissipation rate with a value ∼ 10 −10 W kg −1 , which is closer to the mode value.
In the pycnocline region, the mean dissipation rate is almost two orders of magnitude higher.The arithmetic mean of dissipation rate is lower at station C, where it is of the order of ∼ 10 −7 W kg −1 , than at stations B and A , where it is of the order of ∼ 10 −6 W kg −1 .
Figure 7 shows PDFs of K z truncated in the range [10 −7 , 10 −3 ] m 2 s −1 .The PDFs also show two dynamical regions: (1) a flat distribution with a mode value of K z ∼ 3 × 10 −7 m 2 s −1 is observed below 20 m; and (2) a highly spiked distribution with mode value of ∼ 5 × 10 −5 m 2 s −1 is found in the pycnocline region.Below the pycnocline some reasonable agreement is found between the MLE log-normal fit and www.biogeosciences.net/9/3131/2012/Biogeosciences, 9, 3131-3149, 2012  K z distribution for station C; elsewhere no log-normal behavior could be observed.The arithmetic mean of K z below the pycnocline is ∼ 1 × 10 −5 m 2 s −1 for all stations, whereas the MLE estimate at station C reaches 1.9 × 10 −5 m 2 s −1 .In the pycnocline region depth, arithmetic mean K z values are higher by a factor ∼5. Geometric means largely underestimate K z by nearly 2 order of magnitude below the pycnocline, but agree within a factor of 3 above 20 m.The increase of K z in the pycnocline region observed here is unusual because the larger stratification generally prevents the increase of mixing.It results here from a high mean dissipation rate reaching ∼ 6.6 × 10 −6 W kg −1 in the pycnocline region and its dramatic decrease below the seasonal pycnocline.
From this analysis we also find that MLE of from a log-normal distribution when applicable and arithmetic mean gave similar results; therefore, following the advice of Davis (1996), we will simply consider arithmetic mean in the following.
We next look at the ensemble averaged vertical profiles of and compare them with the parameterization proposed in Sect.2.3.In order to reduce dispersion of and K z values and to allow better comparison with parameterization based on 10 m scale shear, and K z profiles were first smoothed using a 10 m running average.Depth average profiles of and K z were then computed (Figs. 8 and 9).A total 95 % confidence intervals were computed from bootstrap percentiles (Efron and Tibshirani, 1994), except for station A and B where the number of samples was not sufficient.When averaged over all profiles, the dissipation rate decreases from high values of 10 −6 W kg −1 to moderate values of 10 −8 W kg −1 between 10 m and 40 m depth.Below dissipation rates are approximately constant with a value of 10 −8 W kg −1 .This dissipation rate level is comparable with the GM reference level GM that mostly falls within the 95 % confidence interval below 40 m, but overcomes GM by more than one order of magnitude in the pycnocline [10 m, 20 m].
Station average profiles show larger variability at depth that likely results from the lack of statistics, as illustrated by the large confidence intervals.Station B, however, shows a clear decrease of the dissipation rate 10 −10 W kg −1 below GM between 50 and 70 m depth, which is one order of magnitude smaller than GM level.
The parameterized dissipation rate param shows a good agreement with SCAMP measurements.When the average of the whole set of profiles is considered, param falls within the 95 % interval of the SCAMP measurements over 85 % of the profile depth range.The agreement is also good when the average is computed independently for each station; the overall shape of the SCAMP average profile is well reproduced by the parameterization, notably the decrease of dissipation rate in the first thirty meters and the lower dissipation rate at station B around 55 m depth.Large discrepancies exceeding one order of magnitude are observed, but those occur mostly for stations B and A where a very small number of profiles is available.This good agreement suggests that the dissipation rate observed is mostly due to breaking internal waves, and that the wave-field is in equilibrium such that the rate of downscale energy cascade at larger scales is proportional to the rate of dissipation.K z depth averaged profiles are shown in Fig. 9: the average of all LD station profiles shows decreasing values from 5 × 10 −5 m 2 s −1 to 10 −6 m 2 s −1 between 10 m and 40 m depth and then slowly increasing values up to to 5 × 10 −5 m 2 s −1 at 95 m depth.Top and bottom values are significantly higher than the nearly constant GM value of 5 × 10 −6 m 2 s −1 , whereas the local minimum at 40 m depth is lower.Individual station averages evolve in the same range with a noticeable minimum of 10 −6 m 2 s −1 for station B between 45 and 70 m depth.For station A, K z remains highly variable in the range 10 −6 m 2 s −1 to 5 × 10 −5 m 2 s −1 , likely due to the lack of observations.
The overall average values of parameterized K z,param are close to experimental values and fall within the experimental confidence interval over 90 % of the profile depth range.However the local minima around 40 m depth is not reproduced.
The proportion of the different diffusion regimes found according to the Shih et al. ( 2005) classification (Sect.2.2) is also shown in Figs. 8 and 9. Intermediate and strong turbulence regimes dominate in the pycnocline, whereas the     molecular diffusion dominates below 25 m depth; the minimum of turbulent diffusion around 40 m depth corresponds to a region where molecular diffusion regime is found for 80 % of the samples.Finally, we compared the dependance of and param with the fine-scale variability as represented by the 10m squared shear S 2 10 and the strain scaled by buoyancy N 2 η 2 z .To this end we have averaged observed and param into logarithmically spaced bins of S 2 10 (below z = 25 m where shear measurements are available) or N 2 η 2 z .The 95 % confidence interval were computed from the bootstrap method when at least 10 values were averaged in a bin.Both param and show a clear increase trend with increasing values of the buoyancy scaled strain (Fig. 10b).For high values of the buoyancy scaled strain (N 2 η 2 z > 10 −5 s −2 ), a good agreement is found between param and .The values of param mostly fall within the 95 % confidence interval.This suggests that the dependance of the dissipation rate on the buoyancy scaled strain is well reproduced by the parameterization for the higher range of N 2 η 2 z values.The parametrization seems to slightly overestimate the dissipation rate for lower values of the buoyancy scaled strain (N 2 η 2 z < 10 −5 s −2 ).However fewer data are available in this region (less than 10 values for each bin average).Figure 10a shows the dependence of and param with S 2 10 .param and show a comparable increase with S 2 10 values, but the measured dissipation are largely dispersed around the param (S 2 10 ) curve.

Dissipation rate and turbulent mixing inferred from fine-scale parameterization
The modified G89 parameterization, which favorably compared with SCAMP measurements, was applied to the full data set.Because abyssal waters are weakly stratified, N values can be noisy at depth; N was therefore smoothed with a 10 m running average in the computation of param and K z,param .Time depth plots of param and K z,param are displayed for the three LD stations in Fig. 11, whereas the stations averaged profiles are shown in Fig. 12.The range of variation is large: within [10 −12 , 10 −5 ] W kg −1 for param and within [10 −7 , 10 −3 ] m 2 s −1 for K z,param .The eddy cores are characterized by weak values of the dissipation rate (10 −10 W kg −1 ); as shown in Fig. 11a to c, this is most apparent for eddy C, which has the largest vertical extension.In contrast, the highest values of the dissipation rates are observed at the base (eddy C and eddy A) and at the top of the eddies (C, B and A).This increase in param is related to the high shear values at the boundaries of the eddies.Most of the shear at the base of eddies A and C results from the strong near-inertial internal waves that generate velocity oscillations with a vertical wavelength of the order of ∼100-150 m (Fig. 3).The shear induced by these near inertial internal waves has a clear signature on the shear wavenumber spectra of Fig. 2c and e, which show a strong peak at 8 × 10 −3 cpm (125 m wavelength).This increase in shear results in an increase of param by a factor five at the base of eddy C [400-500] m compared to the low param values in We also present the shear to buoyancy scaled strain variance ratio since it is used here to correct the G89 parameterization, following Polzin et al. (1995), and it gives information on the characteristics of the wave field (see Sect. 2.3).For eddy A, R ω remains close to the GM value R ω,GM = 3.A strong increase of R ω is observed for eddy C below 100 m with an average value of 15 above the GM value.R ω also shows high values for eddy B below 300 m depth with an average value of 12.5.In these regions where R ω largely overcomes the GM value, the original formulation of G89 could have overestimated the dissipation by a factor 3 and the shear to strain ratio correction is crucial.
The spatial distribution of vertical diffusivities K z,param (Fig. 11) differs from that of param due to the impact of stratification: regions of the weakest stratification, typically the eddy cores (see eddies C and B), are characterized by relatively large values of K z,param , 10 −4 m 2 s −1 , which are of the same order as those at the base of the eddy where strong near inertial internal waves are observed.
The parameterization was also applied to the deep short duration (SD) stations performed all along the BOUM transect, thus providing a snapshot of parameterized dissipation and mixing rates (Fig. 14).Values of the dissipation rate were, however, not computed in the upper 25 m where direct measures of shear was not available.Figure 13 also shows the evolution of shear (based on 10 m difference of the velocity) and stratification along the transect.In addition, the main water masses of the Mediterranean Sea were represented in order to discuss the possible impact of turbulent mixing with respect to the Mediterranean thermohaline circulation.Indeed, vertical turbulent mixing is a key process for the global thermohaline circulation, which has never been quantified experimentally at the basin scale in the Mediterranean Sea.We have identified three water masses: (1) the Levantine intermediate waters (LIW), (2) the Western Deep Mediterranean Water (WDMW), and (3) the Eastern Deep Mediterranean Water (EDMW).The Levantine intermediate waters (LIW) is an intermediate water mass (typically lying within 200-500 m) formed in the permanent Rhodes cyclonic gyre in the northwestern part of the Levantine Basin (Lascaratos, 1993).We display here the depth of its core, which can be tracked as the local maximum of salinity.Deep waters (WDMW and EDMW) are formed during winter as a result of cooling and evaporation.The main site of formation is the Gulf of Lion for the WDMW and the Adriatic shelf and the Agean Sea for the EDMW (Lascaratos et al., 1999).Following Touratier et al. (2012), the EDMW and WDMW can be tracked as water with potential density anomaly larger than 29.18 kg m −3 and 29.106 kg m −3 , respectively.
The signature of eddies A, B and C is also seen in Fig. 13 as the depression of upper isopycnals as well as local region of minimum stratification.The same features are observed for the Ierapetra anticyclonic eddy in the south of Crete that was also sampled during a SD station.The highest shear and dissipation rates were found in the upper 500 m and up to 1500 m above the bottom, and were generally at a minimum in a region between 500 m and 1500 m depth.As for the LD stations, high shear and dissipation rates are found at the base of eddies A and C (at x = 433, and 3130 km), but some enhancement is also seen around 600 m at the base of eddy B (x = 1810 km).Similarly high shear and dissipation rates are found at the base of Ierapetra eddy (x = 2478 km, z 170 m).Another noticeable feature is the large enhancement of shear and dissipation rate at x = 1478 km in the strait of Sicily, which is possibly associated with topographic induced effects such as internal tides (Garrett and Kunze, 2007) or increased bottom shear associated with the Sicily strait overflow (Beranger et al., 2004;Stansfield et al., 2001).
K z,param show maximum values in the first 1500 m above the bottom where enhancement of the dissipation rate and very weak stratification combine to give vertical turbulent diffusivity reaching locally 10 −3 m 2 s −1 .It is interesting to notice that this region of strong mixing corresponds to the limit of the EDMW, suggesting that deep mixing may be an important factor controlling the diffusive upward transport of this water mass.High values are also found in a region above 500 m where dissipation rate is relatively high and below the 28.5 kg m −3 isopycnal where stratification is already much weaker than in the upper sea (Fig. 13).This region of relatively strong mixing corresponds to the core of the LIW, which suggests a specific impact of turbulent mixing on this water mass all along its path from the Levantine Basin to the strait of Gibraltar.
Despite a relatively strong shear and dissipation rate, a region of low K z,param is found in the upper well stratified sea above the 28.5 kg m −3 isopycnal as a result of a strong stratification.Note that the parameterization was not applied above 25 m depth for SD stations.A second region of weak K z,param is found between ∼500 m and ∼1000 m and is associated with the corresponding region of minimum dissipation rate.(red), GM (dashed magenta) and R ω (black) for the three LD stations (A, B, C), the red dotted line indicates the R ω,GM = 3 value.Second row: arithmetic mean of K z,param (blue), K z GM profiles (dashed magenta) and R ω (black) for the three LD stations (A, B, C), the red dotted line indicates the R ω,GM = 3 value.

Discussion
Lacking dedicated physical measurements of dissipation rate, previous biogeochemistry oriented studies have considered rough estimates of dissipation rate as a constant value in the computation of vertical turbulent diffusion, ignoring large variations resulting from the fine-scale internal wave field (Moutin and Raimbault, 2002;Copin-Montegut, 2000).For instance, Moutin and Raimbault (2002) considered a constant value of = 7 × 10 −10 W kg −1 to estimate vertical diffusion and nutrient fluxes at the nitracline (located below the seasonal pycnocline and generally above 100 m depth) during the MINOS cruise along the Mediterranean Sea.This value corresponds to a GM level dissipation rate at a reference stratification N 0 = 3 cph (a typical value at the base of the pycnocline) in the G89 parameterization.As noted by Moutin and Raimbault (2002), this value is nearly two orders of magnitude smaller than the constant value of 5 × 10 −8 W kg −1 (derived form Denman and Gargett, 1983) considered by Copin-Montegut (2000) in the northwestern Mediterranean Sea.Clearly fluxes estimations could drastically change, depending on the chosen value for .The adapted G89 parametrization used here will significantly improve estimation of mixing compared to these previous rough estimates.Parameterized param estimated here from SD LADCP/CTD profiles along the Mediterranean transect show a mean value 1.5 × 10 −9 W kg −1 between the base of the seasonal pycnocline and 100 m depth, which is slightly higher than Moutin and Raimbault (2002) value.Higher dissipation rates were, however, found within eddies where the average dissipation rate directly observed by SCAMP measurements and estimated from the fine scale parameterization reaches 8.5 × 10 −9 W kg −1 .These values are in between the Copin-Montegut (2000) and Moutin and Raimbault (2002) values.At greater depth the parameterization also shows that there is a strong increase of the dissipation rate in region of high near inertial shear at the base of eddies A and C.
These results suggest a strong influence of anticyclonic eddies on near inertial waves dynamics and mixing.Indeed, anticyclonic eddies induce a negative background vorticity which influences inertial waves propagation.Moutin and Prieur (2012) estimated the eddies vorticity from an analysis of drifting mooring trajectories deployed during BOUM.They find the strongest negative vorticity for eddy A, which reaches ζ = −0.397f; a slightly weaker vorticity for eddy C, reaching ζ = −0.32f; and a weaker vorticity of −0.297f for eddy B. Negative background vorticity can result in a trapping of near inertial waves and explain enhanced near inertial shear.Indeed, as shown by the theoretical work of Kunze (1985), anticyclonic mesoscale vorticity ζ induces a slight decrease of the effective inertial frequency f eff = f + ζ /2 locally; therefore, near inertial waves which evolve in the frequency band f > f eff will encounter their turning points when propagating away from anticyclonic eddy centers (Bouruet-Aubertot et al., 2005) and remained trapped in the eddy core.Numerical studies by Lee and Niiler (1998) have also shown some increase of near inertial shear resulting from the interaction of near inertial waves with frontal structures or eddies, validating partially the mechanism proposed by Kunze (1985).Evidence of a subinertial peak is indeed found at f eff = 0.8f at station A in the first 100 m, which is in agreement with a shift of 0.5ζ .Although spectral analysis (Sect.3.1) did not reveal a subinertial peak at stations B and C, slightly subinertial waves cannot be ruled out for these stations because of the coarse frequency resolution of ± 0.15f and the weaker vorticity of these eddies.The large increase of near inertial shear at the base of eddies A and C may result more specifically from the vertical trapping of near inertial waves specific to baroclinic anticyclonic structures where vorticity and thus f eff increases with depth (Kunze, 1985).This mechanism was observed in a warm core ring (anticyclonic) of the Gulf Stream (Kunze, 1995) together with a (10-100) increase of dissipation rate.Another possible mechanism is the radiation of near inertial waves from the baroclinic adjustment of the eddy.Further work is needed to explore the possibility of this near inertial waves generation at the eddy base.
Estimates of vertical mixing below the pycnocline and above the first 100 m of the eddies A, B, C can be compared with previous studies of mixing within eddies based on tracer release experiments.Such tracer experiments integrate vertical transport processes over a large spatial scale (typically a region below the seasonal pycnocline and above 100 m depth) and temporal scale (a few days).This is an advantage for estimating turbulent mixing because it avoids problems related to the under sampling of highly intermittent turbulent mixing processes, but it also precludes a dynamical characterization of intermittent internal wave breaking as done here.Kim et al. (2005) and Ledwell et al. (2008) found K z 3 × 10 −5 m 2 s −1 between the base of a shallow seasonal mixed layer and 100 m depth in North Atlantic anticyclonic eddies (46 • N) and in the Sargasso Sea (31 • N), respectively.One order of magnitude higher values were found by Law et al. (2001) also in a North Atlantic anticyclonic eddy (59 • N).In BOUM experiment the overall averaged K z found within eddies A, B and C pycnocline base and 100 m depth is 10 −5 m 2 s −1 , which is two times the GM value but is still three times smaller than Kim et al. (2005) and Ledwell et al. (2008) estimates.However wind forcing was relatively weak during BOUM, whereas all the experiments cited above were affected by the passage of storms (Law et al., 2001;Ledwell et al., 2008) or strong wind gusts (Kim et al., 2005) that likely increased internal waves energy and induced dissipation rate.It should also be noted that Greenan (2008) provides smaller K z estimates than Ledwell et al. (2008) from the G89 parameterization of dissipation rate for the same experiment (EDDIES).

Concluding remarks
The microstructure estimates of the dissipation rate showed good agreement with a slightly adapted G89 parameterization (following Polzin et al., 1995) over the first 100 m depth of three anticyclonic eddies.This parameterization was used to infer the dissipation rate param and the turbulent vertical diffusivity K z,param from fine scale CTD/LADCP measurements.The parameterized dissipation rate found from deep SD stations along the transect and outside of eddies is relatively weak in a region below 500 m depth and away from the bottom, whereas it increases toward the bottom and in the upper 500 m.K z,param show maximum values in the range [10 −5 , 10 −3 ] m 2 s −1 in the first 1500 m above the bottom where the stratification is very weak and param is strong.This region corresponds to the location of the EDMW, which suggests that vertical mixing is an important process for the vertical diffusive transport of this deep and dense water mass.These high vertical turbulent diffusivities above the bottom may result from interaction between abyssal flows and the bottom topography.Such a process can generate internal lee waves and can result in enhanced dissipation and mixing, a mechanism proposed by Nikurashin and Ferrari (2010a, b) for the Southern Ocean.A relatively strong bottom circulation was observed in numerical simulations of the Mediterranean Sea circulation by Zavatarelli and Mellor (1994).As well Bouche et al. (2009) show the presence of important bottom mesoscale currents in the Ionian Sea (36 • 19 N 16 • 05 E, 3050 m).Our deep LADCP measurements also reveal strong velocities (up to 0.5 m s −1 not shown) and shear in the first 1500 m above the bottom, particulary in the Ionian Sea at the east of the strait of Sicily.However deep microstructure profiles would be needed to confirm the parameterized dissipation rates since it has been shown that fine scale parameterization can overestimate the dissipation rate near the bottom (Waterman et al., 2012).A second region of enhanced vertical eddy diffusivity is observed above 500 m and below the 28.5 kg m −3 isopycnal where relatively weak stratification N 2 < 5 × 10 −5 rad s −1 and relatively high param coexist.There the vertical mixing could affect the LIW properties as it flows from the Levantine Basin into the Atlantic.Outside of these two regions, K z,param is generally comparable or smaller than the GM canonical value K z,GM = 5 × 10 −6 m 2 s −1 .
This picture is modified in eddies where large near inertial shear at the base of the eddies is associated with dissipation rates exceeding the canonical GM level.Turbulent vertical diffusion increases in these regions of high shear and dissipation, and also within the eddy core because of a much weaker stratification there.The spectacular increase of near inertial shear found in eddies may result from trapping or channeling of near inertial energy input at the surface, a mechanism highlighted in several numerical and experimental studies (Kunze, 1985(Kunze, , 1995;;Lee and Niiler, 1998;Bouruet-Aubertot et al., 2005).Both microstructure estimates of and parameterized param also show strong dissipation rates ( 10 −6 W kg −1 ) in the pycnocline of the three anticyclonic eddies despite a relatively weak wind forcing during these three long stations (<10 m s −1 , data not shown).These strong dissipation rates may result from the proximity of the surface forcing since the average pycnocline location was only 15 m depth.Near inertial energy may also be enhanced within the eddies because of the possible near inertial waves trapping.These strong dissipations rates result in a large vertical turbulent diffusivity in the pycnocline reaching 10 −4 m 2 s −1 .Such vertical turbulent diffusivity associated with strong gradients of temperature and salinity will likely influence the mixed layer heat and salinity budget.Future Y. Cuypers et al.: Turbulence in the Mediterranean sea work including estimations of advective lateral transport will allow us to establish the mixed layer heat and salinity budget and assess the importance of the vertical turbulent fluxes.
Further studies using BOUM observations and numerical models will also allow us a thorough characterization of the impact of mesoscale eddies on biogeochemical processes.The statistical distribution of vertical diffusion may notably be used to reproduce the impact of the strong intermittency of turbulence in one-dimensional biogeochemical models available already (Mauriac et al., 2011).
Regarding the impact on nutrients fuxes, however, Bonnet et al. (2011) have shown that the vertical nitrogen turbulent fuxes determined from K z values obtained in eddies still balance only a small fraction of the nitrogen fluxes resulting from primary production, suggesting a main contribution of regenerated production.Vertical advection was not considered in this study and may also provide significant vertical transport of nutrients, as suggested in previous studies (Ledwell, 2008).

Figure 1 :Fig. 1 .
Figure 1: Bathymetric map of the Mediterranean sea, short duration stations (SD) proles are marked by black dots, long duration LD stations (A,B,C) at the eddy centers are marked by red dots Fig. 1.Bathymetric map of the Mediterranean Sea; short duration stations (SD) profiles are marked by black dots, long duration LD stations (A, B, C) at the eddy centers are marked by red dots.

Figure 2 :Fig. 2 .
Figure2: Upper row: wavenumber vertical shear spectrum for deep SD stations LADCP proles (a) 0 < z < 500 m (b) z > 500m and lower row ensemble averaged wavenumber vertical shear spectra for LD stations LADCP proles at (c) Eddy C (d) Eddy B and (e) Eddy A. In blue, raw data vertical shear spectrum, in magenta dashed line GM level shear spectrum, in black dashed line, the tted GM shear spectrum, in cyan dashed line, the noise tted spectrum, and in red the composite of noise and GM tted spectra.The vertical bars indicate the 95% condence intervals.For spectra in eddies (c, d, e) no t was performed and the noise spectrum is assumed to be identical to the one in (a) (see text)

Figure 5 :Fig. 5 .
Figure5: Strain in gray scale (black high strain, white weak strain) , 0.03kg/m 3 iso-density contours in the mixed layer in red lines, 0.2kg/m 3 iso-density contours below the mixed layer in black and TKE dissipation rates proles in Log10(W.kg−1 ) in colored square marks, station C upper panel, station B middle panel, station A lower panel

Figure 6 :Fig. 6 .
Figure 6: Experimental PDF of the dissipation rate Log 10 (ϵ(W.kg.−1 )) in green for 10m<z<20m in blue for 20m<z<100m, in red MLE t of a log-normal PDF, (a) for the complete LD station data set (stations A,B,C), (b) for station C, (c) for station B and (d) for station A.
does not allow us to state clearly whether distributions of dissipation rate are log-normal there.

Figure 7 :
Figure 7: Experimental PDF of Log 10 (K z (m 2 .s. −1 )) in blue for 10m<z<20m in blue for 20m<z<100m in red MLE t of a log-normal PDF, (a) for the complete LD station data set (stations A,B,C), (b) for station C, (c) for station B and (d) for station A. 44

Fig. 7 .
Fig. 7. Experimental PDF of Log 10 (K z (m 2 s. −1 )) in blue for 10 m < z < 20 m in blue for 20 m < z < 100 m in red MLE fit of a log-normal PDF, (a) for the complete LD station data set (stations A, B, C), (b) for station C, (c) for station B and (d) for station A.

Figure 8 :
Figure 8: Panel 1 proportion of the various diusion regimes (strong, intermediate, molecular) found from SCAMP dissipation rate measurements.Panel 2,3,4,5 Overall and station averaged proles of dissipation rate from SCAMP in blue plain line with 95% condence intervals in gray shading, parameterization in black, and reference GM level in magenta dashed lines

Fig. 8 .Figure 9 :Fig. 9 .
Fig. 8. Panel 1: proportion of the various diffusion regimes (strong, intermediate, molecular) found from SCAMP dissipation rate measurements.Panel 2, 3, 4, 5: overall and station averaged profiles of dissipation rate from SCAMP in blue plain line with 95 % confidence intervals in gray shading, parameterization in black, and reference GM level in magenta dashed lines.

Figure 10 :Fig. 10 .
Figure10: Averaged dissipation rate ϵ (blue circles) and parameterized dissipation rate ϵ param (red line) calculated in bins of shear S 2 10 and buoyancy scaled strain N 2 ζ 2 , the gray shading represents the 95% condence interval calculated when at least 10 data point were averaged in a bin.

Figure 11 :Fig. 11 .
Figure 11: Parameterized dissipation rate ϵ param (a) and diapycnal diusion coecient K z,param (b) as a function of time and depth for the three stations.

Figure 12 :Fig. 12 .
Figure12: First row: arithmetic mean of ϵ param (red), ϵ GM (dashed magenta) and R ω (black) for the three LD stations (A, B, C), the red dotted line indicates the R ω,GM = 3 value.Second row: arithmetic mean of K z,param (blue), K z GM proles (dashed magenta) and R ω (black) for the three LD stations (A, B, C), the red dotted line indicates the R ω,GM = 3 value

Figure 13 :Fig. 13 .
Figure 13: Buoyancy frequency squared (N 2 )(upper panel), shear squared computed from 10m dierence (lower panel) along the BOUM transect.The abscissae x is expressed as the distance along the transect from the st station in the Rhone river mouth (not represented).Black cross indicate approximate positions of the base of eddy A, B, C and Ierapetra eddy.Black lines represent isopycnal 28.5 and 29.1 kg/m −3 , black dashed lines represent the core of the LIW, and the limits of EDMW and WDMW

Figure 14 :Fig. 14 .
Figure 14: parameterized dissipation rate (upper panel) and parameterized diapycnal diusivity (lower panel) along the BOUM transect.The abscissae x is expressed as the distance along the transect from the st station in the Rhone river mouth (not represented).Black cross indicate approximate positions of the base of eddy A, B, C and Ierapetra eddy.Black lines represent isopycnal 28.5 and 29.1 kg/m −3 , black dashed lines represent the core of the LIW, and the limits of EDMW and WDMW

Table 1 .
Mean SCAMP estimations (W kg −1 ) from various methods: arithmetic mean, MLE of the log-normal distribution mean, geometric mean.

Table 2 .
Mean K z estimations (m 2 s −1 ) from various methods, arithmetic mean, MLE of the log-normal distribution mean, geometric mean.