Modelling deep-water formation in the North-West Mediterranean Sea with a new air-sea coupled model: sensitivity to turbulent flux parameterizations

In the north western Mediterranean, the strong, dry, cold winds, the Tramontane and Mistral, produce intense heat and moisture exchange at the interface between the ocean and the atmosphere leading to the formation of deep dense waters, a process that occurs only in certain regions of the world. The purpose of this study is to demonstrate the ability of a new coupled ocean-atmosphere modelling system based on MESONH-SURFEX-SYMPHONIE to simulate a deep-water formation 5 event in real conditions. The study focuses on summer 2012 to spring 2013, a favourable period that is well documented by previous studies and for which many observations are available. Model results are assessed through detailed comparisons with different observation data sets, including measurements from buoys, moorings and floats. The good overall agreement between observations and model results shows that the new coupled system satisfactorily simulates the formation of deep dense water and can be used with confidence to study ocean-atmosphere coupling in the north-western Mediterranean. In addition, to 10 evaluate the uncertainty associated with the representation of turbulent fluxes in strong wind conditions, several simulations were carried out based on different parameterizations of the flux bulk formulas. The results point out that the choice of turbulent flux parameterization strongly influences the simulation of the deep water convection and can modify the calculated volume of the deep water formed by up to one order of magnitude ::::: newly :::::: formed ::::: deep ::::: water :: by : a :::::: factor ::: two.

(1) P3,l16: extra parenthesis Authors' Answer (2): Done (1) p3,l19: an atmospheric model (2): Done (1) p4,l19-20: "In our case, the heat and water fluxes are linearly distributed over the whole mixed layer, the depth of which 10 is given by the depth at which the vertical density gradient becomes negative." meaning not clear to much repetition (depth, which). Try to make shorter sentences.
(2) : Rephrased sentences "In our case, the heat and water fluxes are linearly distributed over the whole mixed layer. The mixed layer is defined by the depth at which the vertical density gradient becomes negative." (1) P5, l26: thee 15 (2): Done (1) p15,l14 "to limit the discrepancies induced differences between by the different parameterizations." I dont understand the meaning of this sentence. Can you reformulate?
(2) : Rephrased sentences "It is likely that these feedback loop effects tend to limit the discrepancies between the different parameterizations." 20

Summary of all changes made in the manuscript
The technical corrections of the referee were taken into account. In the previous version we estimated a factor 11.3 between our different simulations because our reference level was in summer 2012. But looking at the Fig. 12 we think that was not a fair evaluation. Thus, to clarify the Fig. 12 description and help the comparison with other studies the volume of newly formed deep water has been calculated only during increasing phases of this volume. These new numbers have been added in part 5.3.3 25 and in the abstract and conclusion. We apologize for this last minute changes but we think that it is really more representative of the deep water evolution.

Introduction
The North-West Mediterranean Sea (NWMS) is one of the few regions in the world where the deep open ocean convection process is regularly observed (Marshall and Schott, 1999). The strong, dry, cold winds, the Tramontane (north-westerly) and the Mistral (northerly), play a major role in this process. These winds induce intense exchanges between the atmosphere and the sea (Flamant, 2003;Hauser et al., 2003), with a marked loss of surface buoyancy (Schott and Leaman, 1991).The oceanic 5 deep convection can be separated into three phases. In autumn, a cyclonic gyre, bounded to the north by the Northern Current (NC) (Millot, 1999) and to the south by the North Balearic Front (NBF) (Millot and Taupier-Letage, 2005), isolates a weakly stratified water mass whose stratification is progressively eroded by northerly winds (preconditioning phase). In winter, in some years, vertical mixing induced by strong winds leads to the formation of a vertically homogeneous water mass (convective phase) identifiable by its temperature and salinity properties, and generally referred to as new Western Mediter- 10 ranean Deep Water (nWMDW). After the convective phase, the mixed zone undergoes a restratification while the nWMDW is spread throughout the basin by the general circulation (Schott et al., 1996) and submesoscale eddies (Testor and Gascard, 2006) (restratification/spreading phase).
Two components of the MISTRALS programme (http://www.mistrals-home.org/) focused on the study of deep convection in the NWMS and made a major effort in collecting observations during 2012 and 2013. The first component, HyMeX (Drobinski 15 et al., 2013), studied the atmosphere-land-ocean coupled system. In this context, two Special Observation Periods (SOPs) were organized; SOP1 in autumn 2012, during the preconditioning phase (Ducrocq et al., 2013), and SOP2 in winter 2013, during the convection phase . The second component, MERMEX (Sempéré et al., 2010)  Furthermore, the volume of dense water formed by deep convection could be evaluated thanks to the optimal interpolation of the many Conductivity-Temperature-Density (CTD) profiles available (Waldman et al., 2016b). This series of campaigns provided a unique opportunity to test the ability of models to simulate the different phases of the dense water formation process in the western Mediterranean. 25 Several modelling studies of the formation of deep water in the NWMS, have been carried out over different periods and with different models (e.g., Herrmann et al. (2008); Herrmann and Somot (2008); Herrmann et al. (2010); Léger et al. (2016); Estournel et al. (2016a); Waldman et al. (2016a)). Their results show high sensitivity to the initial stratification of the ocean model and to the accuracy of the atmospheric forcing. An attempt to progress on these two issues was proposed by the HYMEX and MER-MEX groups. Concerning the first point, since 2010, oceanographic cruises have been organized in the NWMS each summer by 30 the Mediterranean Ocean Observing System for the Environment (MOOSE; www.moose-network.fr/). These cruises provide a sample of the different water masses of the NWMS based on about 70 CTD profiles. In particular, the observations collected in summer 2012 have been assimilated to provide a more realistic initial state for ocean models Léger et al., 2016) which has been shown to be crucial for the simulation of the winter convection 6 months later. Regarding the atmospheric forcing, the benefit of using a fully coupled system to study air-sea interactions in numerical weather prediction models was already illustrated in previous studies based upon different air-sea coupled systems (e.g., Lebeaupin Brossier and Drobinski (2009);Small et al. (2012); Renault et al. (2012)). These studies have shown that coupled simulations provide a better representation of atmospheric and oceanic surface parameters compared to uncoupled simulations. In particular, during strong wind events coupled simulations capture the rapid SST cooling more accurately, which makes the atmospheric boundary 5 layer more stable and reduces the heat and moistures exchanges. It is likely that this improved representation of the atmospheric forcing could also lead to an improved representation of the deep water formation.
Besides the question related to coupling, there is still significant uncertainty as to the choice of a relevant parameterization to compute the turbulent fluxes for strong wind conditions such as the Mistral and Tramontane. Current parameterizations have been carefully assessed and validated against large data sets. However, due to the limited number of available observations 10 in strong wind conditions, they are known to be inaccurate for wind speeds exceeding 20 ms −1 (e.g., Hauser et al. (2003)).
The sensitivity tests performed by Estournel et al. (2016a) suggest that the uncertainty associated with the turbulent flux computations could have a strong impact on the deep water formation process in the NWMS.
These issues, among others, have motivated the recent development of a new coupling platform (SURFEX OASIS3-MCT) providing better numerical tools to address the scientific and technical questions related to ocean-wave-atmosphere coupling 15 (Voldoire et al., 2017). This coupling platform is based on a multi-surface model SURFEX ( (Masson et al., 2012) and on the OASIS3-MCT (Valcke et al., 2015) code coupler. SURFEX computes the surface-atmosphere fluxes over four surface types (land, town, ocean and inland waters) and can be used in a stand-alone version with prescribed atmospheric forcing or embedded in an atmospheric models ::::: model. The use of OASIS3-MCT allows SURFEX to be linked to various other models including ocean, land, atmosphere, hydrology, waves and sea-ice models. This generic coupling strategy based upon an externalized sur- 20 face model ensures that the surface flux computations are done in a consistent way, independently of the models to be coupled.
As illustrated in Voldoire et al. (2017), this strategy has greatly facilitated the coupling of the different models developed in the French community, including the coupling of the MESONH atmospheric model (Lafore et al., 1997) and the SYMPHONIE ocean model (Marsaleix et al., , 2009.
A first objective of the present study is to show the capacity of the new coupled regional ocean-atmosphere system MESONH-

25
SURFEX-SYMPHONIE to reproduce the formation of deep dense waters during the winter of 2013, in the NWMS. A second objective is to study the sensitivity of the simulations to the parameterization of the turbulent surface fluxes by testing three different parameterizations (Fairall et al., 2003;Andreas et al., 2015;Moon et al., 2007), all based on bulk formulas. This paper is organized as follows. Section 2 describes the coupled modelling system, and the three parameterizations mentioned above that were used for the computation of turbulent fluxes. Section 3 presents the different observation data sets 30 used to evaluate the model results while section 4 details the setup of the numerical experiments. Results are analysed and discussed in section 5. Some conclusions and perspectives are presented in section 6.
The oceanic model uses a horizontal resolution of 1 km. Given the value of the Rossby radius (5-10 km in the NWMS), 1-km grid spacing appears to constitute a reasonable compromise between the computing cost and the necessary resolution. 10 In the vertical, 40 generalized sigma levels are used; 10 of them in the first hundred meters (above the abyssal plain). The resolution just below the sea surface is 1.5 m. The atmospheric model is run with a 10-km horizontal grid spacing and 52 terrain-following vertical levels ranging from 15 m to 15000 m, with 16 of them in the first km. With such resolutions, both atmospheric and oceanic convection must be parameterized. In the case of the ocean, the vertical diffusion is parameterized following Gaspar et al. (1990) with a prognostic equation for the turbulent kinetic energy and a diagnostic relation for the 15 mixing and dissipation lengths. A 1-km resolution is still too coarse to explicitly resolve convective plumes, which thus need to be parameterized. Different parameterizations have been proposed (e.g., Marsland et al. (2003)). The most common and basic one consists in artificially increasing the vertical diffusion coefficient in statically unstable layers (e.g., Waldman et al. (2016a)). In our case, the heat and water fluxes are linearly distributed over the whole mixed layer, the depth of which . :::: The ::::: mixed ::::: layer is given :::::: defined by the depth at which the vertical density gradient becomes negative. By doing so the first level 20 under the surface does not support the entire amount of heat loss by itself, which prevents the development of static instabilities at the surface. Furthermore, this parameterization is consistent with the nearly linear vertical variation of the buoyancy flux in the convective layer (Deardorff et al., 1969). Regarding the atmosphere, shallow and deep convection are parameterized with mass-flux schemes according to Bechtold et al. (2001) and Pergaud et al. (2009), respectively.
In the coupled system, the surface fluxes are computed by SURFEX on the atmospheric model grid. They are sent to the 25 ocean model by the OASIS3-MCT coupler, which also performs their interpolation on the ocean model grid. Conversely, the ocean model computes the sea surface temperature and sends it to SURFEX using the OASIS3-MCT coupler, which takes care of its interpolation on the atmospheric model grid. The coupling frequency is set to 10 minutes and the interpolation between the two model grids is bilinear.
The computational domains used for this study are presented in Fig. 1. The atmospheric grid covers the whole western 30 Mediterranean basin while only part of it is covered by the ocean grid (blue area in Fig. 1) since the Alboran Sea and part of the Tyrrhenian Sea were excluded to avoid straits issues. Outside the ocean grid, the air-sea fluxes are computed using the sea surface temperature provided by the OSTIA data base (Donlon et al., 2012), the horizontal resolution of which is about 6 km.

Turbulent flux parameterizations
Turbulent air fluxes at the air/sea interface are computed from bulk type parameterizations based on the Monin-Obukhov similarity theory (Foken, 2006). These parameterizations compute the turbulent fluxes as |τ | = ρ a u * 2 (1) where τ is the momentum flux, H the sensible heat flux, LE the latent heat flux; ρ a the air surface density and where u * , θ * , q * are scaling parameters for momentum, potential temperature, and humidity, respectively. The momentum scale, u * , is referred to as friction velocity.
Classically, the scale parameters are expressed as a function of the vertical gradients of the mean fields at the air-sea interface, 10 the surface roughness and the atmospheric stability. Although based upon the same formalism, the turbulent flux parameterizations differ in the way they specify the different roughness lengths and the so-called stability functions. In particular, the validity of the Charnok's formulation (Charnock, 1955) which is generally used to relate u * to the dynamic roughness length has often been questioned for strong wind conditions.
In this study, three well-established parameterizations have been used:

15
-The COARE3.0 parameterization (Fairall et al., 2003) is one of the most widely used in the modelling community. This parameterization derives from the COARE 2.6 algorithm (Fairall et al., 1996) originally developed from the observations performed during the TOGA-COARE experiment (Webster and Lukas, 1992) in the North Pacific. An important upgrade in COARE3.0 is a new formulation of the surface (dynamic and scalar) roughness lengths which slightly increases the fluxes for wind speeds exceeding 10 m s −1 . Although COARE3.0 has been validated against a much larger data set 20 (∼ 7000 observations) than the one used for COARE2.6 , COARE3.0 remains mostly reliable for wind speeds below 20 m s −1 due to the limited number of observations available in strong wind conditions. It is worth noting that the influence of waves (available as two possible options in COARE3.0 but not extensively validated) was not activated in our study.
-The ANDREAS parameterization (Andreas et al., 2015) is a novel and more physically-based approach which distinguishes two different contributions to the turbulent heat fluxes: the standard air-sea interfacial fluxes controlled by 25 molecular processes right at thee :: the : air-sea interface on the one hand, and the sea spray fluxes controlled by microphysical processes around sea spray droplets on the other hand. As opposed to the COARE-type algorithm, the friction velocity used to compute the interfacial fluxes is parameterised as a function of the 10 m wind speed at neutral stability, eliminating thereby the uncertainty associated with the definition of the dynamic roughness length and the use of the Charnock (1955) expression. The sea-spray contribution becomes notable only for wind speeds exceeding 13 30 m s −1 . Small droplets are then ejected by surface waves into the atmospheric surface layer. They cool, evaporate and can significantly contribute to the air-sea exchanges of heat and water. The sea spray fluxes are computed using the fast microphysical algorithm described in Andreas (2005). ANDREAS parameterization has been established with a data set of ∼ 4000 observations with wind speeds up to almost 25 m s −1 .
-As opposed to COARE and ANDREAS, the MOON parameterization (Moon et al., 2007) mainly relies upon model results. It has been developed based upon the results of a coupled wave-wind model. The simulations of 10 idealized tropical cyclones have been used to derive a new expression of the dynamic roughness length, which limits the increase 5 of friction velocity with wind for wind speeds exceeding 12.5 m s −1 . This new formulation was indirectly validated using the Geophysical Fluid Dynamics Laboratory coupled hurricane-ocean prediction model (Kurihara et al., 1998).
For 5 hurricanes observed in the Atlantic Ocean, the new formulation lead to better results than the former one (based upon the Charnock's formulation) with a clear improvement of the cyclone intensity and no degradation of its track and central pressure.  Among the three parameterizations, MOON was the one that produced the strongest sensible heat flux over a wide range of wind speeds (from 6 m s −1 to 23 m s −1 ). The latent heat fluxes presented the same hierarchy in intensity as the sensible heat fluxes. However, the difference between COARE3.0 and ANDREAS only occurred from wind speeds greater than 16 m s −1 (as compared to 8 m s −1 for the sensible heat flux). To summarize, the three parameterizations gave significantly different results, particularly MOON, which produced the largest turbulent heat fluxes in a wide range of wind speeds. It should also be 25 noted that the differences between the parameterizations may become noticeable from wind speeds as low as 8 m s −1 and can dramatically increase with the wind speed, especially for the heat fluxes.

Available Observations
Several in situ datasets, collected by the MOOSE, HyMeX and MerMeX programmes, are used in the present study and are briefly described below. More details are available in Estournel et al. (2016b). First, the Lion meteorological buoy provides 30 hourly measurements of the air temperature and humidity at 2 m above the sea surface, the wind speed at 10 m, and the sea surface temperature. A mooring named LION, part of the MOOSE network, and located only 5 km away from the Lion buoy is also available. It provides vertical profiles of temperature and salinity, sampled over 21 levels for temperature and 15 for salinity, between 150 and 2300 m deep (Houpert et al., 2016).
In addition, several oceanographic cruises were carried out in the NWMS during 2012-2013 (Waldman et al., 2016b). In August 2012, the annual cruise organized by MOOSE for the monitoring of the NWMS, provided 69 vertical surface-tobottom profiles of temperature and salinity. The same network of CTD stations was deployed in February 2013 and April 2013, 5 during the DEWEX leg 1 and leg 2 cruises. A smaller network was implemented in autumn 2012 during the DOWEX cruise, with only 41 profiles. Furthermore, 399 temperature and salinity profiles measured by 14 Argo floats, 4 of them specifically deployed by HyMeX with a daily cycle, are also taken into consideration in this study. The above-mentioned cruises allowed the spatial distribution and temporal evolution of ocean stratification to be well described. Ocean stratification is commonly assessed using the Stratification Index (SI Eq. 4). 10 where ρ is the potential density and Z the reference level. SI is expressed in kg m −2 .
Note : that SI(Z) = 0, means that the water column is mixed at least to depth Z or, in other words, that the mixed layer depth is greater than or equal to Z.  -Brossier et al., 2014;Léger et al., 2016;Rainaud et al., 2016). Figure 4 shows the temporal evolution of surface atmospheric and oceanic parameters (10 m wind speed, 2 m temperature and humidity and sea surface temperature) computed at the buoy, together with the corresponding observations. The same surface parameters are also presented in Fig. 5 in the form of scatter plots (simulations versus observations) while the associated statistics (bias, root mean square error and correlation coefficient) are given in Table 2. in autumn leads to a progressive decrease of the sea surface temperature (Fig. 4 d), which reaches its minimum value of 12.9 • C (i.e. the temperature of the deep water) in early winter and then remains nearly constant during the convective period. 30 All three simulations accurately reproduce the time evolution of the wind speed at the buoy throughout the 8-month period, with a correlation of 0.9 and a bias lower than 0.2 m s −1 . In particular, the timing of the strong wind events is well captured.
Moreover, the wind speed maxima are well represented (Fig. 5 a), which is essential to correctly reproduce the intense air-sea exchanges associated with convection (Herrmann and Somot, 2008). Finally, there is no significant difference in wind speed between the different simulations (Figs 4 a, 5 a and Table 2).

5
For the other surface atmospheric parameters (2 m air temperature and relative humidity), slightly larger discrepancies are found from one simulation to another. Air temperature and humidity remain relatively close to observations in terms of correlation (respectively 0.98 and 0.85, Table 2). Bias and root mean square error exhibit larger but still weak differences between simulations. The largest difference is found for humidity. In particular, it is clear from Fig. 4 c that the moisture drops associated with the strong wind episodes are more pronounced in COARE and ANDREAS than in MOON. 10 The calculated sea surface temperature is remarkably well correlated with the observations for all simulations (> 0.98, Table   2). This correlation is mainly due to the representation of the seasonal cycle and to the weak variability of the SST during the winter period when the SST ceases to evolve. The drops of SST associated with the events of Tramontane and Mistral in autumn are well captured by the three simulations. However, there are significant differences between simulations in autumn.
In particular, the COARE SSTs appear overestimated ( Furthermore, during the autumn, as shown in Fig. 4 d, the sea surface temperature evolves differently in each simulation. This is particularly the case for the COARE simulation, which presents significantly warmer SSTs than the other two simulations and than observations. This results in larger turbulent heat fluxes during low wind periods in autumn for the COARE simulation than for ANDREAS and MOON (Fig. 6).
To complement the analysis, the radiative fluxes and precipitation, which also contribute to the air/sea exchanges, have been analysed. The results are in good agreement with the observations and do not reveal any significant difference between the simulations, they are not presented here.

5
In summary, all the simulations generally are in good agreement with the surface parameters observed at the Lion buoy throughout the eight months considered, although significant punctual differences may appear between the different simulations, especially during the Tramontane and Mistral events. These differences mainly affect the heat and moisture exchanges whereas the momentum exchanges are very weakly impacted. Although the differences remain fairly weak and as reflected by the statistical analysis, in our coupled system, the MOON parameterization gives the best agreement with the available 10 observations. However, considering the impossibility of directly validating the air-sea fluxes and also the multiple sources of uncertainty in such a complex modelling system (and their possible compensations) it is not clear whether the MOON flux parameterization is better per se or whether it is simply is the most suitable parameterization for our modelling system.

Buoyancy Mass fluxes and oceanic stratification 15
The air-sea exchanges are now assessed through the Buoyancy Mass Flux (BMF). This flux, directly linked to turbulent fluxes but also to radiative fluxes and precipitation, is formulated as follows: where α is the thermal expansion coefficient in K −1 , Q net is the net heat flux (sum of the net radiative flux, and turbulent heat fluxes) in W m −2 , C p is the specific heat capacity in J kg −1 K −1 , is the saline contraction coefficient in psu −1 , SSS is 20 the surface salinity in psu, ρ 0 is a reference density in kg m −3 , E is the evaporation and P the precipitations in m s −1 . BMF is expressed in kg m −2 s −1 .
The time evolution of the buoyancy mass flux computed at the Lion buoy and for each simulation is shown in Fig. 7  March) with a decrease in the cumulated buoyancy mass flux without significant decrease in SST (since the water column is nearly mixed) and, finally, the start of the restratification period (21 March to 30 April) when the buoyancy mass flux starts to 30 increase. As can be seen in Fig. 7 b, the buoyancy loss during the preconditioning period is nearly equivalent to the buoyancy loss during the convective period. The largest cumulative buoyancy loss at the end of the simulations is obtained with MOON (195 kg m −2 ) followed by COARE (170 kg m −2 ) and ANDREAS (165 kg m −2 ). By the end of the simulations, COARE has produced a slightly larger loss of buoyancy than ANDREAS while turbulent fluxes for the Tramontane and Mistral events are greater in ANDREAS. This is due to the larger buoyancy mass flux during the calm wind periods of the preconditioning phase 5 in COARE.  observed during the autumnal Tramontane and Mistral events at the buoy location (Fig. 7 a), which is situated in the vicinity of the frontal zone. This suggests that, during the preconditioning period, the dynamics of the NBF plays a major role in the loss of surface buoyancy in the deep water zone. Furthermore, the front displacement is modulated by the wind intensity. The increase (decrease) of stratification during the period of weak (strong) winds (Fig. 7 c) is due to the lateral advection of light (heavy) water by the northward (southward) displacement of NBF . Because of these horizontal processes and 25 their feedback on the buoyancy mass flux, the three simulations experience different time evolutions. The large discrepancy seen between the COARE and MOON buoyancy mass fluxes and stratification indices may also result from differences in the NBF progression. In winter, during the convective phase, only the two maxima of the Gulf of Lion remain. In the deep water zone, as the SST has reached a nearly-constant value, there is no significant effect of the SST structures on the buoyancy mass flux. At the Lion buoy, the buoyancy mass flux is not affected by the SST and the local increase in stratification during the 30 period of weak winds (Fig. 7 c) is principally due to the advection of light water into the mixing zone by baroclinic instability (Marshall and Schott, 1999).
The comparison of the three simulations highlights the impact of the surface flux parameterization. In MOON, the buoyancy mass flux losses are stronger than in COARE and affect a much wider area.

Validation of oceanic stratification with oceanographic cruises
To assess the evolution of the oceanic stratification, Stratification Indices collocated in space and time with all CTD and Argo profiles were calculated for each simulation and compared with the corresponding values deduced from the observations.
Results were analysed in terms of bias (observations minus simulations). Figure 9 shows the spatial distribution of the SI (1000 m) bias obtained for each simulation and for the two DEWEX oceanographic cruises (leg 1 in winter 2013 and leg 2 in spring 5 2013) while Table 3 gives the values of the averaged SI for observations and the SI simulations bias and root mean square error However, it is noteworthy that the SI (1000 m) bias is significantly reduced in the MOON simulation (Fig. 9 c).
The good performance of MOON with respect to the stratification is further confirmed by to the south-east and the water column is too mixed while, in the convection area, the bias is fairly low, especially for the MOON simulation (Fig. 9 f). On average (Table 3)  In summary, the evolution of stratification is not simply related to the buoyancy mass flux but results from the complex interaction between buoyancy mass fluxes and advective processes. In this study, the MOON simulation significantly reduces the negative bias of the stratification, in the convection area, during the convective and restratification periods, again indirectly suggesting that the air-sea fluxes are most realistic in the MOON simulation.

Impact on the air-sea exchanges on the deep water formation
In this section observations provided by the Lion buoy and the DEWEX cruises are compared with the results of the different simulations in terms of mixed patch, timing of convective process and volume of deep dense water formed.

Mixed patch
The spatial distribution of the convection is first examined. Figure 10 shows the position of the stratified and mixed CTD and March, destroys this restratification and again leads to the full mixing of the water column and to an increase of the seawater density up to 29.12 kg m −3 . This mixing event is not well represented by COARE and ANDREAS. Overall, despite the tooearly full mixing, the MOON simulation gives the best representation of the deep convection at the LION mooring and correctly captures its three densification events.

Conclusions
This study focused on assessing the ability of a regional ocean atmosphere coupled system, based on the SYMPHONIE, 10 SURFEX and MESONH models, to correctly represent ocean convection and deep water formation in the NWMS. Several realistic simulations, were carried out over a period of 8 months, from summer 2012 to spring 2013, and were used to investigate the sensitivity of the system to the parameterization of turbulent fluxes.
First, this study shows the ability of the air-sea coupled system to reproduce the evolution of the ocean and the atmosphere for several months by relying only on realistic initial and boundary conditions and without resorting to data assimilation or nudging.

15
For all simulations, a good correlation is obtained between the observed and computed surface parameters at the Lion buoy.
During the Tramontane and Mistral events, the turbulent heat fluxes differ significantly from one simulation to another, directly impacting the atmospheric and oceanic surface parameters. In a previous study devoted to the same case study, Rainaud et al. (2016) underlined the difficulty of reproducing air surface temperature and moisture during the Mistral and Tramontane events and advocated the use of an air-sea coupled model and a purposely adjusted turbulent flux parameterization. In addition to 20 air surface temperature and moisture, sea surface temperature is also strongly sensitive to the turbulent flux parameterizations.
Our results also suggest that coupling plays a key role during the autumn storms when the rapid drops of SST reduce the turbulent heat fluxes. In winter, the impact of the coupling is likely to be weaker since the SST does not vary much anymore.
Among the three parameterizations, we found that MOON (i.e. the parameterization yielding the strongest heat turbulent fluxes) significantly reduced the bias between the observed and computed surface parameters. Unfortunately, due to the lack of flux 25 measurements at the buoy, it was not possible to validate the computed turbulent surface fluxes directly.
Then, the buoyancy mass flux was calculated and compared to the evolution of the stratification for each simulation. The stratification evolution is directly impacted by the buoyancy mass flux loss but there is no strict correlation between stratification and buoyancy mass flux. As already shown in Estournel et al. (2016a), this confirms the importance of the advective processes on the evolution of the stratification in the deep water area. Moreover, these advective processes also directly impact the surface 30 buoyancy mass fluxes, particularly during preconditioning period, when the position and dynamics of the North Balearic Front clearly affect these fluxes. This interaction between the buoyancy mass fluxes and the advective processes is clearly an air-sea coupled process, which deserves to be analysed in greater depth. In terms of stratification, the use of MOON also led to a general reduction of the bias between observed and computed parameters.
However, this conclusion regarding the good performance of the MOON flux parameterization needs to be further consoli-10 dated.
First the present results were obtained with a coupled system. They could probably be different with uncoupled simulations.
In air-sea coupled simulations, the interactive evolution of ocean and atmosphere influences the turbulent heat fluxes, which themselves modify the atmospheric and oceanic surface fields involved in the flux calculation. In statically unstable Mistral and Tramontane conditions, if the sensible (or latent) heat flux increases, the vertical temperature (or humidity) gradient is reduced, 15 which in turn limits the increase in the sensible (latent) heat flux. It is likely that these feedback loop effects tend to limit the discrepancies induced differences between by :::::: between : the different parameterizations. The results of partial and preliminary uncoupled simulations (not shown) suggest that these discrepancies could be larger than in the coupled simulations. It would be therefore of great interest to disentangle the effect of the flux formulation from the effect of the air-sea coupling and to check whether the MOON parameterization still improves the results in uncoupled conditions. However, it is not straightforward 20 to isolate the coupling effect in a clean and rigorous way. This requires a series of carefully-designed experiments in which the current coupled system is step-by-step downgraded into an uncoupled system till it exactly mimics the behavior of the atmospheric and ocean models in their stand-alone configuration. In our current system, this type of study is hampered by the fact that the surface fluxes are computed on the atmospheric grid, ie at a coarser resolution that the one used by the ocean model.

25
The differences in resolution between the atmospheric and ocean models (10 and 1 km, respectively), though partly justified by scale considerations, is also a debatable question. A further development will thus investigate the sensitivity to the resolution of the atmospheric model. In the present configuration, the atmospheric model does not have the possibility of representing scales fully adjusted to those of the oceanic model. In particular, with a 10 km resolution, the local maxima and horizontal gradients of the surface parameters are probably too smooth, which may affect the air-sea interactions especially in 30 the vicinity of the oceanic front (Small et al., 2008) and could also modify the response of the coupled system to the different parameterizations.
In addition, the role of the waves necessitates further investigation. In our study, the waves are not considered in COARE and MOON and only indirectly accounted for in ANDREAS. In ANDREAS, the depth of the spray layer is computed as a function of the significant wave height (Andreas et al., 1995). The latter is rather roughly estimated from a simplified parameterization 35 based on wind speed (Andreas and Wang, 2007). Similar crude relationships are used in COARE3.0 for the wave height and wave period when the wave option is turned on. Another envisaged development will couple the current system with a wave model (Michaud et al. (2012) and the references cited therein) and revisit the results obtained with the ANDREAS and COARE3.0 parameterizations.
Acknowledgements. This work is a contribution to the MISTRALS/HyMeX programme through the ASICS-MED (ANR-12-BS06-0003) project funded by the French National Agency for Research (ANR). Data were obtained from the HyMeX programme, sponsored by grants from MISTRALS/HyMeX and Météo-France. The authors acknowledge the international ARGO programme, the LEFE/GMMC programme and the French NAOS project for supporting the deployment of profilers. Argo and CTD data were collected and made freely available by the CORIOLIS project (http://www.coriolis.eu.org) and programmes that contribute to it. We acknowledge the crews of R/V Suroit and Tethys II and the scientists involved in the different cruises mentioned in this paper. Numerical simulations were performed using HPC resources from Bechtold, P., Bazile, E., Guichard, F., Mascart, P., and Richard, E.: A mass-flux convection scheme for regional and global models, Quarterly Ducrocq, V., Braud, I., Davolio, S., Ferretti, R., Flamant, C., Jansa, A., Kalthoff, N., Richard, E., Taupier (Mlawer et al., 1997), (Fouquart and Bonnel, 1980) Microphysics ICE3 (Caniaux et al., 1994), (Pinty and Jabouille, 1998) (Estournel et al., 2009) Mixing Eddy Kinetic Energy (Gaspar et al., 1990) Convection penetrative convection (Deardorff et al., 1969),  River input Lateral condition