Articles | Volume 18, issue 2
Research article
08 Apr 2022
Research article |  | 08 Apr 2022

Untangling the mistral and seasonal atmospheric forcing driving deep convection in the Gulf of Lion: 2012–2013

Douglas Keller Jr., Yonatan Givon, Romain Pennel, Shira Raveh-Rubin, and Philippe Drobinski

Deep convection in the Gulf of Lion is believed to be primarily driven by the mistral winds. However, our findings show that the seasonal atmospheric change provides roughly two-thirds of the buoyancy loss required for deep convection to occur for the year 2012 to 2013, with the mistral supplying the final third. Two NEMOMED12 ocean simulations of the Mediterranean Sea were run from 1 August 2012 to 31 July 2013, forced with two sets of atmospheric-forcing data from a RegIPSL coupled run within the Med-CORDEX framework. One set of atmospheric-forcing data was left unmodified, while the other was filtered to remove the signal of the mistral. The control simulation featured deep convection, while the seasonal simulation did not. A simple model was derived by relating the anomaly scale forcing (the difference between the control and seasonal runs) and the seasonal scale forcing to the ocean response through the stratification index. This simple model revealed that the mistral's effect on buoyancy loss depends more on its strength rather than its frequency or duration. The simple model also revealed that the seasonal cycle of the stratification index is equal to the net surface heat flux over the course of the year, with the stratification maximum and minimum occurring roughly at the fall and spring equinoxes.

1 Introduction

Deep convection, also known as open-ocean convection, is an important ocean circulation process that typically occurs in high-latitude regions (Marshall and Schott1999). Localized events are triggered by the reduction of the stable density gradient through sea surface layer buoyancy loss. One such area of deep convection is the Gulf of Lion (GOL) in the Mediterranean Sea. The deep-convection events that occur in this region aid the general thermohaline circulation of the Mediterranean Sea by forming the Western Mediterranean Dense Water (WDMW) (Robinson et al.2001). After its formation, this dense water spreads out along the northwestern basin among the deeper layers of the Mediterranean Sea (MEDOC1970), with some transported along the northern boundary current towards the Balearic Sea (Send and Testor2017) and some transported to the south within eddies (Beuvier et al.2012; Testor and Gascard2003) into the southern Algerian Basin and towards the Strait of Gibraltar (Béranger et al.2009), completing the cyclonic circulation pattern of the sea. The water column mixing that occurs during a deep-convection event also brings oxygenated water down from the oxygen-rich sea surface layer and injects sea-bottom nutrients upwards towards the surface (Coppola et al.2017; Severin et al.2017), resulting in increased phytoplankton blooms in the following season (Severin et al.2017).

Significant deep-convection events occur every few years in the GOL (Somot et al.2016; Houpert et al.2016; Marshall and Schott1999; Mertens and Schott1998), driven by the mistral and tramontane winds. These sister northerly flows bring cool, continental air through the Rhône Valley (mistral) and the Aude Valley (tramontane), leading to large heat transfer events with the warmer ocean surface (Drobinski et al.2017; Flamant2003). These cooling and evaporation events destabilize the water column in the GOL and are widely accepted to be the primary source of buoyancy loss leading to deep convection (Lebeaupin-Brossier et al.2017; Houpert et al.2016; L'Hévéder et al.2012; Lebeaupin-Brossier et al.2012; Herrmann et al.2010; Lebeaupin-Brossier and Drobinski2009; Noh et al.2003; Marshall and Schott1999; Mertens and Schott1998; Madec et al.1996; Schott et al.1996; Madec et al.1991a, b; Gascard1978).

Here, we investigated the mistral's role in deep convection in the GOL (as the mistral and tramontane winds are sister winds, we will refer to them jointly as “mistral” winds). Its role was determined by running two NEMO ocean simulations of the Mediterranean Sea from 1 August 2012 to 31 July 2013, forming a case study of the encapsulated winter. One simulation was forced by unmodified atmospheric-forcing data, while the other was forced by a filtered atmospheric dataset with the signal of the mistral removed from the forcing. Thus, the ocean response due to the mistral events could be separated and examined, revealing the effects of seasonal atmospheric change alone. A multitude of observational data were collected during this year in the framework of the HYdrological cycle in the Mediterranean EXperiment (HyMeX) (Estournel et al.2016b; Drobinski et al.2014), which provided a solid base of observations to validate the ocean model results.

In particular, our findings quantify the following items:

  • the separated and combined effect of the mistral and seasonal atmospheric cycle on deep convection,

  • the dominant attribute of the mistral causing buoyancy loss,

  • the source of the buoyancy loss due to the seasonal atmospheric cycle.

In addition, another 2 years were also studied with the same methodology as the 2012–2013 winter: the 1993–1994 and 2004–2005 winters. The 1993–1994 winter does not have a deep-convection event and allows us to compare a deep-convecting year vs. a non-deep-convecting year. The 2004–2005 winter is a well-studied deep-convecting winter and offers some additional literature to draw analysis upon, as well as an additional deep-convecting year to compare and contrast with, using the same methodology as the 2012–2013 winter.

There are three distinct sections of the deep-convection cycle: the preconditioning phase in the fall, the main large overturning phase in the winter and early spring (when deep convection occurs), and the restratification and spreading phase during the summer (MEDOC1970; The Lab Sea Group1998). The focus of study is on the preconditioning and overturning phase where the mistral is stronger and more frequent (Givon et al.2021) and therefore plays a larger role in the deep-convection cycle.

The model and methodology used are described in Sect. 2. Model results and validation are presented in Sect. 3 for the 2012–2013 winter. Patterns observed in the model results lead to the development of a simple model that describes the effects of the mistral and seasonal cycles. This simple model is presented in Sect. 4. The 2 additional years, 1993–1994 and 2004–2005, are presented in Sect. 5, following the two previous sections that focus on the 2012–2013 winter. Our concluding remarks are presented in Sect. 6.

2 Methodology

In our study we used the NEMO ocean model to run two ocean simulations forced by unmodified and modified atmospheric-forcing data from a coupled WRF/ORCHIDEE simulation. Information on the mistral events, used later when developing the simple model in Sect. 4, was extracted from the unmodified atmospheric-forcing data and from ERA-Interim Reanalysis data. The main metric used in this article to examine the model results and relate them to deep convection is the stratification index (SI). Each of these components are described below in their own subsection.

2.1 NEMO

The Nucleus for European Modeling of the Ocean (NEMO) ocean model (, last access: 24 March 2022) was used in bulk formula configuration to simulate the GOL region with two distinct simulations, both performed from 1 August 2012 to 31 July 2013. In bulk formula configuration, sea surface fluxes are computed from parameterized formulas using atmospheric and oceanic measurable variables as inputs, such as temperature and wind velocity. The following parameterized formulas are used to calculate the latent heat flux, QE, the sensible heat flux, QH, the longwave radiation heat flux, QLW, and the surface shear stress, τ:


where z is the height above the sea surface the atmospheric variables are provided at, with the zero values (subscript “0”) representing the values at the sea surface. u is the horizontal wind vector, with Δu=uz-u0 as the difference between the wind velocity and sea surface current (assuming a no-slip condition at the ocean surface). q and θ are the specific humidity and potential temperature of air, respectively. Λ and cp are the latent heat of evaporation and the specific heat of water, respectively. ρ0 is the air density at the sea surface. SST is the sea surface temperature, with SSTK as the sea surface absolute temperature. ϵ is the sea surface emissivity, σ is the Stefan–Boltzmann constant, and QLW,a is the atmospheric longwave radiation. CE, CH, and CD are the coefficients of latent heat, sensible heat, and drag, respectively, and are defined in Large and Yeager (2004, 2008).

The net downward heat flux, Qnet, is described by the summation of the terms in the following equation (Large and Yeager2004; Estournel et al.2016b):

(5) Q net = Q SW + Q LW + Q H + Q E ,

where QSW is the downward shortwave radiation. Snowfall and precipitation are included in the simulation calculations but excluded here for brevity in the following sections, as the buoyancy loss due to the water flux at the surface is essentially negligible (Somot et al.2016).

The NEMO model was run in the NEMOMED12 configuration using NEMO v3.6. The domain is shown in Fig. 1b; it covers the Mediterranean Sea and a buffer zone representing the exchanges with the Atlantic Ocean. This configuration features a horizontal resolution of 1/12 (roughly 7 km) and 75 vertical levels (with a variable vertical resolution from 1 m at the surface to 135 m at the bottom). The 3-D temperature and salinity fields are restored towards the ORAS4 global ocean reanalysis (Balmaseda et al.2013) in the buffer zone. The conservation of volume in the buffer zone is achieved through strong damping of the sea surface height (SSH) towards the ORAS4 reanalysis. The Black Sea, runoff of 33 major rivers, and coastal runoff are represented by climatological data input from Ludwig et al. (2009). A deeper explanation of the configuration and boundary conditions is given by in the following works: Waldman et al. (2018), Hamon et al. (2016), Beuvier et al. (2012), Lebeaupin-Brossier et al. (2011), and Arsouze et al. (2012). The initial conditions were provided by an ocean objective analysis by Estournel et al. (2016b).

Figure 1The domains of both the WRF domain from the RegIPSL coupled WRF/ORCHIDEE simulation within the Med-CORDEX framework (a) and the NEMOMED12 configuration domain (b). The region of interest, the NW Med., is outlined by the box. This region is later used in Fig. 7. The location used to study the temporal development of deep convection in the GOL is at 42 N, 5 E, and the other location, used in conjunction with the aforementioned point to determine mistral events, is Montélimar, France, at 44.56 N, 4.75 E.

2.2 Atmospheric forcing

The atmospheric forcing used in the simulations were the outputs of RegIPSL, the regional climate model of IPSL (Guion et al.2021), which used the coupling of the Weather Research and Forecasting Model (WRF) (Skamarock et al.2008) and the ORCHIDEE Land Surface Model (Krinner et al.2005). The run was a hindcast simulation (ERA-Interim downscaling) performed at 20 km resolution and spanning the period of 1979–2016 within the Med-CORDEX framework (Ruti et al.2016). The u and v components, specific humidity, potential temperature, shortwave and longwave downward radiation, precipitation, and snowfall were all used to force the ocean simulations.

For the “control” simulation, the forcings was used as is. For the “seasonal” simulation, the u and v wind components, specific humidity, and potential temperature were filtered (see Fig. 2) over the entire domain shown in Fig. 1a. These variables were chosen as they are the primary variables that affect the surface flux calculations in the bulk formulae (Eqs. 14). The variables relating to radiation and precipitation fluxes were left unchanged. The filtering removes the short-term anomaly-scale forcing from the forcing dataset (the phenomena with timescales of under a month), effectively removing the mistral's influence on the ocean response. This creates two separate forcing datasets, one with the anomaly-scale forcing included and one with just the seasonal-scale forcing (hence the designation of control and seasonal).

The filtering process was performed by a moving window average:

(6) χ i = 1 i + N + 1 j = 0 i + N x j ,

where χi is the averaged (filtered) value at index i of a time series of variable x with length n, where i=0n. The window size is equal to 2N+1, which in this case is equal to 31 d. The ends have a reduced window size for averaging and thus show edge effects. The edge effects did not affect the forcing used for the NEMO simulations, as they were before and after the ocean simulation dates, as two full-year atmospheric-forcing datasets were used for the simulations.

The moving window average was applied to each time point per day over a 31 d window (i.e. for 3-hourly data, the time series is split into eight separate series, one for each timestamp per day – 00:00, 03:00, 06:00, etc. – then averaged with a moving window before being recombined). This was done to retain the average intra-day variability yet smooth the intra-monthly patterns, as the diurnal cycle has been shown to retard destratification by temporarily reforming a stratified layer at the sea surface during slight daytime warming. This diurnal restratification has to be overcome first before additional destratification of the water column can continue during the next day (Lebeaupin-Brossier et al.2012, 2011) and is shorter than the typical mistral event length of about 5.69 d (Table 1).

An important note must be made about the filtering process. The mistral primarily acts in the higher-frequency range but at a lower frequency than the diurnal cycle, as mentioned above. However, it also features signal strength in the lower frequencies on the seasonal scale. This is due to the fact that the mistral becomes stronger, longer, and more frequent during the preconditioning phase than during the rest of the annual cycle (Givon et al.2021). The moving window averaging we have applied to filter the mistral out of the atmospheric forcing primarily removes the higher-frequency portion of the mistral's presence. However, it also removes part of the lower-frequency portion as well that other filters, such as the Butterworth filter, struggle with, without removing more of the seasonal signal that is not influenced by the mistral than intended. This reveals a very interesting point about the structure of the mistral and the use of “seasonal” and “anomaly” timescales. Since the Mistral primarily acts in the higher frequencies, the anomaly timescale will refer to both the higher- and lower-frequency portions of the mistral that were filtered out but will be treated mainly as referring to the higher-frequency portion in the discussion. The “seasonal” timescale will then refer to the remaining lower-frequency signal and average diurnal cycle post-filtering.

The result of the filtering is shown in Fig. 2. Temperature and specific humidity were filtered as is, while the wind speed was component (u and v) filtered, preserving the general wind direction (Fig. 2 wind direction polar plot). Due to the slow movement of intermediate and dense water, which is on the order of about a year for newly formed WDMW to move into the southern Algerian Basin (Beuvier et al.2012) and on the order of decades for total circulation (Millot and Taupier-Letage2005), we assume the processes outside the NW Mediterranean subdomain in Fig. 1b that are affected by the filtering have a negligible impact on the GOL processes on the preconditioning phase timescale.

Figure 2An illustration of the filtering (averaging) process described by Eq. (6). Here the variables q, T, u, and v are shown for both the unfiltered (control, black) and filtered (seasonal, blue) datasets at the nearest grid point to 42 N, 5 E. Note how the peaks of the time series are removed and the general wind direction is conserved.


2.3 Mistral events

Mistral events will be used for developing the simple model in Sect. 4 for their role in driving buoyancy loss at the ocean surface. Events were determined from the WRF/ORCHIDEE dataset in combination with the ERA-Interim Reanalysis dataset (Dee et al.2011). Two main criteria were used to define a Mistral event:

  1. northerly flow with a streamwise flow direction ±45 about the south cardinal direction above 2 m s−1 at two locations simultaneously, i.e., at Montélimar, France (45.5569 N, 4.7495 E), and in the GOL (42.6662 N, 4.4372 E),

  2. the presence of a Genoa Low, defined as a closed sea level pressure contour around a minimum in the field, using 0.5 hPa intervals, anywhere in the box defined by the latitudes 38 and 44 N and longitudes 4 and 14 E (a slightly different domain than that of Givon et al.2021).

The events during the preconditioning period, 30 August 2012 to 21 February 2013, were then manually checked and edited to remove single-day gaps to better represent the data according to a visual inspection of the atmospheric-forcing data. For k mistral events, each event's duration, Δtk, and period from the beginning of the event to the next event, Δτk, was determined and is provided in Table 1 for the entire ocean simulation period (for further analysis regarding the selection of these criteria; see Givon et al.2021).

Table 1The start date, duration (Δtk), and period between each event (Δτk) for each Mistral event (k) for the entire NEMO simulation period of 1 August 2012 to 31 July 2013. Superscripts “d” and “a” denote events used as ideal cases for calculating αd and αa, respectively, in Sect. 4.3 and Appendix A2.1 and A2.2.

The average values for mistral events from 30 August 2012 to 16 February 2013 are Δt=5.69 d and Δτ=10.88 d. The standard deviations for the same time frame are σΔt=4.22 d and σΔτ=4.59 d.

Download Print Version | Download XLSX

2.4 Stratification index

A useful metric to quantify the vertical stratification of a column of water is the stratification index (SI; Léger et al.2016; Somot et al.2016; Somot2005; sometimes called the “convection resistance”). It is derived from the non-penetrative growth of the mixed-layer depth (MLD; i.e., without entrainment; Turner1973), which has been shown to be an accurate approximation for open-ocean convection (Marshall and Schott1999):

(7) z t = B ( t ) N 2 ( z ) z ,

where N2 is the Brunt–Väisälä frequency, z is the vertical coordinate along the water column, zt is the growth of the mixed-layer depth, and B is the potential buoyancy loss the water column can endure before removing stratification (in m2 s−3). Separating by variable and integrating results in the equation for SI gives

(8) SI = 0 D N 2 z d z ,

where D is the depth of water column. If N2 is assumed to be constant throughout the water column, the integral simplifies to

(9) SI = D 2 2 N 2 .

SI provides a zero-dimensional index to track stratification and can be easily related to the buoyancy loss experienced by the water column due to the atmosphere. Because of this, in this article SI will be used as the diagnostic to track the atmosphere's impact on the stratification of the GOL waters.

3 Winter of 2013

The results are presented in two parts for the 2012–2013 winter (additionally referred as the 2013 winter): model validation against observational data and the model results of the deep-convection cycle presented from the center of convection, roughly at 42 N, 5 E. This section and Sect. 4 primarily focus on the winter of 2013.

3.1 Model validation

To validate the model results, data from the HyMeX (, last access: 24 March 2022) database were compared to the NEMO control simulation. Sea surface temperature (SST) data from Météo-France's Azur and Lion buoy were compared with the control simulation SST of the nearest grid point in NEMOMED12. Figure 3 shows the comparison. The Azur buoy data were missing SST measurements from 19 January to 10 July 2013, but where the data are available, NEMO corresponds well with the observations. The same is true for the Lion buoy data, which had measurements for the entire time covered by the simulations. This comes as no surprise, as the NEMOMED12 simulations' SST is restored to the observational dataset of Estournel et al. (2016b). However, this also means that the calculated surface sensible heat fluxes should be fairly accurate, as both the sensible heat flux and latent heat flux calculations depend on the SST (Eqs. 14).

Figure 3SST comparison between the NEMO control run and the Azur (a) and Lion (b) buoy SST datasets. Where the data are available, the model results match the buoy data fairly well.


Additionally, the control simulation density and potential temperature profiles were compared to conductivity–temperature–depth (CTD) measurements also procured from the HyMeX database. The CTD measurements were collected during the HyMeX Special Observation Period 2 (Taupier-Letage2013; Estournel et al.2016a; Drobinski et al.2014) mission. The CTD profiles collected at approximately the same time and location were averaged together to adjust for small variances and gaps in the data. The averaged profiles and their standard deviations are visualized in Figs. 5 and 6. The locations of the CTD profiles are shown in Fig. 4.

Figure 4Locations of the CTD and Argo profiles. The red circles represent the CTD locations, and the blue triangles represent the Argo float profile locations. The deep-convection area is marked by the box with a dashed perimeter and 42 N, 5 E is marked by a cross.

Like with the SST comparisons, the profiles from the nearest grid point in the control simulation domain were used for the CTD comparisons. The root-mean-square error (RMSE) and bias (calculated as the difference between the model values and the observation values) for each of the averaged CTD profiles and corresponding control simulation profiles was calculated and is presented in Table 2. Overall, the control simulation and CTD profiles are decently well correlated but not perfect, with low RMSE and bias for both density and potential temperature. The density profiles have an average RMSE that is lower than the average RMSE for the potential temperature profiles: 0.025 kg m−3 and 0.094 C, respectively.

Argo float profiles from the HyMeX database were also compared to the control simulation, again with profiles from the nearest grid point being used. A total of 3118 potential temperature profiles within the box, bounded by the 40 to 44 N latitudes and the 2 to 8 longitudes, representing the GOL area were considered (see Fig. 4). The average RMSE between the Argo profiles and control simulation profiles was 0.43 C, with an average bias of 0.23 C. These values are larger than the values of the comparison with the CTD profiles. However, considering the sheer volume of profiles and the fact that during stratified conditions the temperature can range over a few degrees from the surface to the lower layers, these results are not unexpected.

Temperature differences on the order of 10−2C are potentially all that is required to sustain an ocean convective cycle (Marshall and Schott1999), and density differences for the same order of magnitude (10−2 kg m−3) are used to separate newly formed dense water during deep convection (Houpert et al.2016; Somot et al.2016; Beuvier et al.2012). This means our model results should be studied with a critical eye, as they may not be fully representative of the true ocean response, given the bias and RMSE values from comparing the simulation to CTD and Argo profiles. Additionally, meanders around 40 km in wavelength form due to baroclinic instability along the edge of the convection path (Gascard1978). This could mean the deviations from observations are due to out-of-phase meanders around the convective patch region in the model relative to actuality. Regardless, we believe the simulations are accurate enough to provide interesting results for the transient and regional-scale response of the GOL, which covers the main interest of our study.

Figure 5Comparison of CTD and NEMO control simulation density profiles. The CTD profiles were averaged by combining multiple vertical profiles collected at the date and location into one profile. The standard deviation of this averaging, σρCTD, is marked in red and is present for all plots but may be difficult to see for 7 March and 9 May.


Figure 6The same as Fig. 5 but for potential temperature.


Table 2RMSE and bias between the averaged observed CTD density and potential temperature profiles and the nearest NEMO control grid point profiles for the respective variables.

The average RMSE and bias for the density profiles was 0.025 and −0.009 kg m−3, respectively. The average RMSE and bias for potential temperature was 0.094 and 0.032 C, respectively.

Download Print Version | Download XLSX

3.2 Stratification index and mixed-layer depth

Figure 7 shows the SI calculated over the GOL for both simulations: the top row is for the control simulation, and the bottom row is for the seasonal simulation. An important distinction between the two results is that deep convection is present in the control simulation but not in the seasonal simulation. This is more clearly seen in Fig. 8c (the closest NEMOMED12 grid point to 42 N, 5 E), as the control simulation MLD reaches the sea floor on 13 February 2013, while the seasonal MLD remains close to the sea surface. This confirms that atmospheric forcing with timescales less than a month, e.g., the mistral winds, provide a significant amount of buoyancy loss, as without them deep convection fails to occur. There is, however, still significant loss of stratification at the location of the GOL gyre in the seasonal simulation, which is visible in the bottom row of Fig. 7 for 13 February 2013. This spot of destratification is also present (but less so) in the preceding and proceeding plots of the same row.

Figure 7The stratification index across the GOL (the area marked as NW Med. in Fig. 1b) at different timestamps. The top row displays the values of SI for the control simulation, and the bottom row displays the values of SI for the seasonal simulation. The box denoted by DC indicates the area of deep convection in the GOL that was not seen in the seasonal simulation.


Figure 8The stratification index of the nearest NEMO grid point to 42 N, 5 E and MLD over the year of both simulations. Panel (a) shows the stratification index for the control run, SIS+δSI, and the seasonal run, SIS. Panel (b) shows the difference between the control and seasonal stratification index, δSI. Panel (c) shows the MLD for both simulations. Mistral events are shown in all three panels and are colored green for events during the preconditioning and deep-convection phase and red for events outside of the preconditioning phase. Mistral events with dotted hatching (the blue-colored intervening time between events) are used as ideal destratification (restratification) events to compute the simple model restoration coefficients. The specific timestamps tA through tD correspond to the timestamps of the plots in Fig. 7: 30 August 2012, 11 December 2012, 13 February 2013, and 1 June 2013. Two definitions of MLD are plotted in (c): one calculated by a vertical change in density less than 0.01 kg m−3, denoted by Δρ, and one calculated by a vertical diffusivity less than 5×10-4 m2 s−1, denoted by Kz. The MLD denoted by the vertical diffusivity criteria follows the turbocline depth and is taken to represent the mixed-layer depth more accurately, as it matches the deep-convection timing in the stratification index.


To investigate the time series ocean response in more detail, a spatially averaged time series of the SI for both simulations was analyzed at the grid point nearest to 42 N, 5 E. These coordinates were selected as it is the point with the most destratification in Fig. 7 and is the typical center of deep convection in the GOL (Marshall and Schott1999; MEDOC1970). The spatial averaging involved horizontally averaging the immediately adjacent grid points such that nine grid points in total were averaged, centered around 42 N, 5 E. The stratification index from the control simulation is given as the sum of δSI+SIS, while the stratification index of the seasonal simulation is given as SIS. The difference between the two, δSI, should contain the change in stratification due to shorter timescale atmospheric events, such as the mistral, because of the filtering performed in Sect. 2.2. δSI+SIS, SIS, and δSI are all shown in Fig. 8.

Both the control and seasonal runs start off with an SI value of 1.57 m2 s−2 (beginning of Fig. 8a) and then diverge at the first major mistral event starting on 30 August 2012. After diverging, the two runs remain diverged until the end of the simulation run time, ending with a difference of about −0.22 m2 s−2, which is seen in δSI (shown in Fig. 8b). As noted earlier, the most striking difference between the control and seasonal runs is the occurrence of deep convection in the control run, occurring when δSI+SIS is equal to 0 (also signifying when the MLD reaches the sea floor), and the lack of deep convection in the seasonal run, as SIS only reaches a minimum of 0.43 m2 s−2. Additionally, if only the anomaly timescale atmospheric forcing is considered, and therefore δSI is the only stratification change from the initial 1.57 m2 s−2, the roughly −0.6 m2 s−2 of maximum destratification that the anomaly timescale provides is not enough to overcome the initial stratification. This means that both the intra-monthly and the inter-monthly variability of the buoyancy loss, reflected in δSI+SIS, are required for deep convection to occur.

Another significant result is the timing of the deep convection. Deep convection initially occurs on 13 February 2013, which is before SIS reaches its minimum on 21 February 2013 but after δSI reaches its minimum on 11 December 2012. After δSI reaches its minimum, it stays around −0.43 m2 s−2 until May 2013, where it starts to increase. This means that while the induced destratification from the anomaly-scale forcing would have been able to overcome ∼0.6 m2 s−2 of stratification to form deep convection in December, the seasonal stratification was only low enough in February for both δSI and SIS to have a combined destratification strong enough for the water column to mix. In other words, the seasonal atmospheric-forcing destratified the already preconditioned water column into deep convection along with a simultaneous mistral event. This means buoyancy loss due to the anomaly forcing may not necessarily be the only trigger for deep convection, at least for this year. This can be seen more clearly in the MLD, as the MLD grows over two mistral events preceding it reaching the seafloor.

4 Process analysis

To pick apart how the atmospheric forcing influences the stratification in the Gulf of Lion, a simple model was developed to separate out the individual components of interest for both the seasonal and anomaly timescales.

4.1 Simple model derivation

To connect the mistral to the ocean's response, we make the assumption that the response is a superposition of the seasonal response and the anomaly response. This means the effects of the mistral can be categorized as anomalies affecting the short-term anomaly timescale and studied separately from the seasonal response. In terms of the Brunt–Väisälä frequency, this linear combination is represented by N2=δN2+NS2, where δ denotes the anomaly terms and “S” denotes the seasonal terms. To determine the mistral's effect, we derive a simple model to describe the SI of the water column in response to atmospheric forcing (for the full derivation, see Appendix A1). We start with the energy equation for incompressible fluids (White2011), then multiply the equation by -g/T0 to express the energy equation in terms of buoyancy, assuming that the ocean's density varies negatively proportionally with temperature, ρ=-βT. We then perform partial differentiation with respect to z to obtain an equation describing the Brunt–Väisälä frequency, N2, in response to the atmospheric forcing, given by a forcing function, F(t). Separating by timescale, we arrive at the following partial differential equations:


F(t) is preceded by a minus sign for ease of derivation, as positive quantities of F(t) mean heat (and thus buoyancy) is removed from the water column.

To simplify the seasonal timescale, we assume NS2 is only a function of time, t:

(12) N S 2 t = - F S ( t ) .

Assuming a homogeneous seasonal Brunt–Väisälä frequency over the depth of the water column gives us the relation for the seasonal stratification index, SIS, and seasonal atmospheric forcing:

(13) SI S t = - D 2 2 F S ( t ) .

To simplify the analytical solution for the anomaly timescale, we describe the advection term, V⋅∇(δN2), as a restoring term, R=α(δN2), which relates the overall Brunt–Väisälä frequency to its seasonal component, NS2, using the linear assumption made before. This means the restoring coefficient, α, represents the advective operation. This results in the following differential equation:

(14) δ N 2 t + α ( δ N 2 ) = - δ F ( t ) .

4.2 Seasonal solution and forcing

The solution for the seasonal timescale is relatively straightforward. As already shown, Eq. (13) relates the seasonal stratification, SIS, to the seasonal atmospheric forcing, FS(t). We have the following definition of FS(t) from Appendix A1:

(15) F S ( t ) = z q a , S g ρ c p T 0 = g ρ c p T 0 q a , S z ,

where cp is the specific heat capacity of water, taken as 4184 J kg−1 K−1; g is gravity; ρ is the density of water, taken as 1000 kg m−3; and T0 is the reference temperature, taken as the average seasonal sea surface temperature of 292.4 K. This means SIS can be related to the seasonal volumetric atmospheric heat transfer, qa,S. Setting qa,S=-Qnet,S/D, where Qnet,S is the seasonal net downward heat flux at the ocean surface from Eq. (5), we can calculate SISt from Qnet,S. If we integrate both sides of Eq. (13) using z, after plugging in Eq. (15) and the relationship for Qnet,S, as SIS is constant with respect to (w.r.t.) z, Eq. (13) becomes

(16) SI S t = g 2 ρ c p T 0 Q net , S ,

where g2ρcpT010-9 m4 J s−2, which means the derivative of SIS w.r.t. time, t, multiplied by 109 is on the same order of magnitude as Qnet (with the subscript S now dropped for convenience, as the rest of the subsection discusses seasonal heat fluxes), which is what we see in Fig. 9a for the 2013 winter, with SISt×109 following the curve of Qnet. This relationship means when Qnet crosses zero with a negative derivative, SIS experiences a maximum and vice versa for a minimum. Additionally, the longer Qnet remains negative, the more seasonal destratification is incurred by the ocean. The seasonal variation of Qnet is primarily driven by the solar radiation, QSW, which is evident in Fig. 9b. Consequently, the maximum and minimum values for SIS occur around 21 September and 21 March, the fall and spring equinoxes. The asymmetry in Qnet is mostly caused by the slightly seasonally varying latent heat flux, QE, followed by the sensible heat flux, QH, both of which also decrease the net heat flux by roughly 100 to 200 W m−2, depending on the time of the year. QLW remains roughly constant during the year, decreasing Qnet by roughly −100 W m−2. These results are corroborated by the results of multiple model reanalysis for the region as well (Song and Yu2017).

Figure 9The smoothed (with Eq. 6) seasonal surface heat fluxes over the point 42 N, 5 E for the seasonal simulation. Pane (a) contains the seasonal stratification index, SIS, and its derivative, SISt, comparing it to the seasonal net heat flux, Qnet (the subscript S is dropped for convenience). Panel (b) shows the net heat flux separated into its components: QE, QH, QSW, and QLW for latent heat flux, sensible heat flux, shortwave downward flux, and longwave downward flux, respectively (neglecting contributions from precipitation and snowfall). The different line colors correspond to the similarly colored axes.


Equation (16) and Fig. 9 convey that the seasonal stratification is primarily driven by shortwave downward radiation. The other terms, the longwave, latent heat, and sensible heat fluxes shift the net heat flux negative enough for the ocean to have a destratification–restratification cycle. If the net heat flux was always positive, stratification would continue until the limit of the simple model applicability. This is an important finding because if future years feature less latent and sensible heat exchange due to warming or more humid winters, there will be less seasonal destratification, requiring more destratification from the anomaly timescale to cause deep convection. Consecutive years of decreasing latent and sensible heat fluxes could form a water column that is too stratified to allow deep convection to occur.

4.3 Anomaly solution and forcing

To solve for the anomaly timescale, described by Eq. (14), we assume δF(t) can be represented by a pulse function shown in Fig. 10. This pulse function assumes the primary forcing at the anomaly timescale is represented by the mistral events. Each Mistral event, k, has a duration, Δtk, and a period between the start of the current and following event, Δτk. δFk is the strength of the forcing for each event. Inserting this into Eq. (14) allows us to solve it in a piecewise manner. As with the seasonal timescale, we assume the water column has a homogeneous Brunt–Väisälä frequency, allowing us to make use of Eq. (9). The restoring coefficient then only represents the horizontal advection, as the vertical component becomes zero with our assumption of a homogeneous N2. The last assumption is that the restoring coefficient remains constant for each section of the forcing function:

(17) δ SI k ( t ) = δ SI k - 1 ( t k ) + δ F k α d 1 - e α d ( t - t k ) e - α d ( t - t k ) [ t k , t k + Δ t k ) δ SI k - 1 ( t k ) + δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) [ t k + Δ t k , t k + Δ τ k ) ,

where αd and αa are the restoring coefficients during ([tk,tk+Δtk)) and after ([tk+Δtk,tk+Δτk)) a mistral event, respectively.

Figure 10The mistral forcing as a pulse function used to solve Eq. (14). k corresponds to the event, and δFk corresponds to the forcing strength of the mistral event. Δtk corresponds to the duration of the of the Mistral event, and Δτk corresponds to the period between events, with tk denoting the start of event k.


Further assuming δFk=δF, Δtkt, and Δτkτ for all k, which results in a periodic pulse function with constant amplitude and period, we can simplify Eq. (17) using the sum of a finite geometric series. At the beginning of the preconditioning period, destratification has not yet begun, and therefore the initial δSI is zero, resulting in the following equation set:

(18) δ SI k ( t ) = D 2 2 δ F α d [ 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 + 1 - e α d ( t - t k ) ] e - α d ( t - t k ) [ t k , t k + Δ t k ) D 2 2 δ F α d [ 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 + 1 - e α d Δ t k ] e ( α a - α d ) Δ t k - α a ( t - t k ) [ t k + Δ t k , t k + Δ τ k ) ,

This final equation set allows us to describe the integrated effect of consecutive mistrals and to easily pick apart the effects of the mistral's different attributes, including the frequency of the events.

To determine the value of the restoring coefficients, a normalized function was derived for each section of a mistral event (derivation shown in Appendix A2.1 for during an event and Appendix A2.2 for after an event). The resulting normalized functions were fitted against the NEMO δSI results in Fig. 8 for the denoted ideal events in Table 1 (denoted as “d” for the dates with ideal destratification taking place during the event and “a” for the dates with ideal restratification taking place after the event) and given the average event values of Δt=5.69 d and Δτ=10.88 d. The result of the fitting is shown in Fig. 11, with αd having a fitted value of 0.235 d−1 and αa having a fitted value of 0.021 d−1. If we recall the meaning of αd and αa from the derivation of the simple model in Sect. 4.1, this means the advective term in Eq. (14) has a larger role in the destratification phase of the mistral event than in the restratification phase, as it is an order of magnitude larger. This result suggests horizontal mixing occurs between events, as a smaller value for the restoration coefficient during the restratification phase means the existence of weaker horizontal gradients than during the preceding destratification phase.

Figure 11The normalized theoretical solutions (Eqs. A42 and A48) for during (a) and after (b) a destratification event fitted to the ideal mistral events from Table 1 and δSI values from the NEMO results in Fig. 8. A value of 0.235 d−1 for αd and a value of 0.021 d−1 for αa was found. Panel (c) shows the δSI response using the determined restoration coefficients, given an ideal Mistral event with the average values of 5.69 d for the duration and 10.88 d for the period. The average strength of a mistral, δF=4.01×10-8 s−2 d−1, was taken from values found in Table A1 from Appendix A.


The strength of each Mistral event, δFk, was found in a similar way by solving for δFk after noting that the initial value of δSIk(tk) is equal to δSIk−1(tk) (derivation found in Appendix A3). Following this, the values of δSI from the NEMO results in Fig. 8 were plugged in to determine the values of δFk (see Table A1 in Appendix A for the resulting values).

4.3.1 Mistral strength and destratification

Mistral events do not always lead to destratification. Some events in Fig. 8 fail to create further destratification and actually continue to restratify the water column. The simple model can describe this phenomena. To determine which events lead to destratification vs. those that do not, we take the derivative with respect to time of Eq. (17) for during an event. This results in the following equation:

(19) δ SI k ( t ) t = - α d δ SI k - 1 ( t k ) + D 2 2 δ F k α d e - α d ( t - t k ) .

The quantity δSIk(t)t must be less than zero for destratification to occur, which means if αd is a positive quantity (refer to Appendix A2 and Fig. 11), δSIk-1(tk)+D22δFkαd must be a positive quantity. If some destratification has already occurred relative to the seasonal stratification, such that δSIk-1(tk)<0, then D22δFkαd must be larger than -δSIk-1(tk) for destratification to occur. Recalling that δFk is positive when heat is removed from the water column, this means that additional mistral events must overcome the current amount of destratification to further destratify the water column. Otherwise, no destratification occurs or restratification may even occur. An example of this can be seen with the mistral event starting on 2 January 2013, which lasts for 17 d in Fig. 8b. The event starts off with an initial destratification of −0.48 m2 s−2 and ends at −0.41 m2 s−2, a net restratification of 0.07 m2 s−2. This is despite the fact that this event has a positive δFk value of 3.80×10-8 s−2 d−1 (from Table A1).

The combined overall effect of this result can be seen in Fig. 8b, as the consecutive mistral events during the preconditioning phase cause destratification to a minimum of -0.6 m2 s−2 for δSI on 11 December 2012. Proceeding events after this minimum fail to continue to destratify the water column, and restratification instead occurs on the anomaly timescale, even before deep convection occurs. The seasonal stratification, SIS, and anomaly destratification, δSI, bring the total SI to zero on 13 February 2013, resulting in deep convection.

4.3.2 Dominating mistral attribute

A pertinent question to ask is which attribute of the mistral, the frequency, strength, or duration, is the most important when it drives destratification. Figures 12 and 13 show the results of varying δF, Δt, and Δτ individually (in panels a, b, and c, respectively) in Eq. (18). The other variables are kept at the mean value when not varied. The dashed lines in both figures show the limit of potential destratification per case. What we can see is that stronger mistral events, with an increased value for δF, result in more destratification, with the reverse happening with decreased values. Decreasing the event duration, Δt, results in less destratification; however, increasing event duration causes more destratification up to the limit where the individual events converge into one single long event and the destratification converges to the dashed line limit. After this, there is no additional destratification. Increasing or decreasing the frequency of events (decreasing or increasing the period, Δτ) only minimally changes the accrued destratification, due to the fact that the magnitude of δSIt is dependent on the strength of the current mistral event and the already achieved destratification. Decreasing the frequency (increasing the period), allows for more restratification to occur after an event, but the proceeding event has a larger difference between current destratification and the event strength, leading to destratification that almost reaches the same level as the case with more frequent events. Increasing the frequency has a similar effect to increasing the duration; when the period is zero, the forcing becomes one large event, converging the resulting destratification to the dashed line.

Figure 12Equation (18) plotted with one variable varying in each plot with the other variables held constant at the mean value. Panel (a) varies the strength of the mistral, δF, (b) varies the duration, Δt, and (c) varies the period between events, Δτ. Panel (d) varies the restoration coefficient during the destratification phase, αd, and (e) varies the restoration coefficient for the restratification phase. Panel (f) varies the number of events.


Figure 13The same as Fig. 12, but SIS is added to the results from Eq. (18).


To more accurately quantify the effect of each attribute, we separate δSI into its total derivative in terms of the mistral attributes:

(20) d δ SI = δ SI δ F d δ F Strength + δ SI Δ t d Δ t Duration + δ SI Δ τ d Δ τ Period ,

Due to the lack of available total derivatives for δF, Δt, and Δτ, we approximate them with their respective standard deviation: σx≈dx. Before we determine the partial derivatives for each attribute, note that in Figs. 12 and 13 subplot f that changing the number events, k, does not change the potential destratification limit (the dashed line). This means the potential destratification does not change with the number of events. Another important factor is the character of the potential destratification limit: it approaches some asymptotic value as k approaches infinity. We can take advantage of this by differentiating the destratification phase of equation set (18) with respect to k, taking tτ at the end of the phase, where the destratification equals the potential destratification:

(21) δ SI k ( t = t k + Δ t ) k = D 2 2 δ F α d e - α d Δ t - 1 ( 1 - e ( α a - α d ) Δ t - α a Δ τ ) - e [ ( α a - α d ) Δ t - α a Δ τ ] k ( ( α a - α d ) Δ t - α a Δ τ ) ,

Plugging in the mean values of Δt, Δτ, and δF, and taking k=16, for the 16 events that occurred during the preconditioning phase, the above derivative equates a very small value of -5.93×10-11 m2 s−2 per event. This confirms the small change in the potential destratification with increasing events. Taking k to infinity and noting that αd>αa results in the following equation:

(22) δ SI = δ SI ( t = t k + Δ t ) = D 2 2 δ F α d e - α d Δ t - 1 1 1 - e ( α a - α d ) Δ t - α a Δ τ ,

We have an equation that describes the potential destratification, δSI, in terms of the mistral attributes, independent of the number of events. Differentiating by the different attributes (see Appendix A5 for the resulting analytical derivations) and plugging in the mean values where appropriate, we arrive at the resulting values: the derivative w.r.t. the strength of the mistrals, δSI/δF, is equal to a value of -1.07×107 m2 d; the derivative w.r.t. the duration, δSI/Δt, is equal to -7.60×10-3 m2 s−2 d−1; and the derivative w.r.t. the period, δSI/Δτ, is equal to 2.77×10-3 m2 s−2 d−1 (larger periods mean less frequent mistral events and thus less destratification). Replacing δSI with δSI in Eq. (20), we can now multiply the partial derivatives with the standard deviations to determine which attribute leads to the most potential destratification. The strength term is equal to -1.28×10-1 m2 s−2, the duration term has a value of -3.21×10-2 m2 s−2, and the period term has a value of 1.27×10-2 m2 s−2. With the strength term an order of magnitude larger than the other two terms according to this simple model, the strength of the mistral event is the most sensitive attribute when it comes to the effect of the mistral on destratification, followed by its duration.

4.4 Simple model results

A complete and average mistral destratification and restratification event according to Eq. (18) is given in Fig. 11c, which took the average mistral values from Tables 1 and A1 and the restoring coefficients from Appendix A2. During the event, marked in green, the mistral causes destratification. After the event, marked in blue, the ocean column restratifies until another event occurs (denoted by the dashed line). This is the same behavior we see in Fig. 8.

If we put together Eq. (17) with the duration and period information from Table 1 and mistral strength information from Table A1, we can create a time series of δSI to compare the integrated response of the simple model to the NEMO model results. This comparison is presented in Fig. 14. The simple model results resemble the NEMO simulation results quite well, which is expected as the fitted values for the restoring coefficients and the values for the Mistral event strengths are extracted from the NEMO model results. However, this means that a series of variable pulse like mistral events can recreate the patterns that we see in the NEMO results for δSI with decent accuracy. This essentially confirms that the mistral events are the primary driving component of the heat loss at the anomaly timescale that leads to destratification.

Figure 14The combined effect of Eq. (17) for multiple mistrals with the mistral data from Tables 1 and A1. Panel (a) shows the calculated δSI+SIS response, while (b) is the calculated simple model δSI vs. the NEMO δSI simulation results. Effects from mistrals after deep convection are included with the dashed blue line and show that mistrals after deep convection can retard the proceeding restratification during the restratification phase.


5 Comparison with additional years

To understand the results of the 2013 deep-convection year in a more generalized context, 2 additional years were simulated and analyzed in a similar fashion: the winters of 1994 and 2005 (the 1 June 1993 to 31 May 1994 year and the 1 June 2004 to 31 May 2005 year, respectively). The winter of 2005 featured a deep-convection event (Beuvier et al.2012; Herrmann et al.2010), whereas the 1994 winter did not (Somot et al.2016). These years were chosen for the sake of having an additional deep-convection year and a year without deep convection to see if there are any significant differences for non-deep-convection years and other deep-convecting years. Simulations for the additional years were run in the same manner as the 2012 to 2013 year and with the same model and configuration. The seasonal run similarly had its atmospheric forcing filtered with the same method as in Sect. 2.2, with the control run left unmodified. The only differences between these additional simulations and the 2012–2013 simulations are the initial conditions, restoration data, and start dates of the simulations. For the additional years, the NEMO simulations were initialized with and restored to the MEDRYS reanalysis (Hamon et al.2016). This was done as the initial conditions and restoration data for the 2012–2013 simulations were only available for that year. The starting time beginning in June rather than July was an arbitrary decision and is not believed to significantly affect the results or comparisons.

5.1 Stratification index

As we are comparing separate years together, the time series simulation results were spatially averaged over a larger area for the additional years: from 42 to 42.5 N and 4.25 to 5 E. Figures 15 and 16 show the SI time series for the 1994 and 2005 winters, respectively. The increased spatial averaging reduces the extent at which the SI destratifies due to surrounding stratified water being averaged in, which makes the 2005 winter appear as though it is too stratified to deeply convect (0.22 m2 s−2), even though it is more destratified than the 1994 winter (0.36 m2 s−2). The MLD, however, clarifies that the 2005 winter experiences deep convection in the simulation results and the 1994 winter does not (not shown), which is consistent with other findings (Beuvier et al.2012; Herrmann et al.2010; Somot et al.2016).

Similar to the 2013 winter simulation set, the seasonal run for the 2005 winter does not experience deep convection, again demonstrating the necessity to include forcing on both timescales for deep convection to occur. The 1994 winter reveals something just as interesting: the δSI is more negative for the 1994 winter, at −0.55 m2 s−2, than for the 2013 winter at the time of minimum control SI (with deep convection in the latter but not the former), at −0.43 m2 s−2. This means even a larger anomaly-driven destratification is not able to overcome the residual stratification in a non-convecting year, despite the fact that both the 1994 and 2005 winters each featured a lower maximum control SI than the 2013 winter: 1.83 m2 s−2 for the 2013 winter and 1.73 and 1.79 m2 s−2 for the 1994 and 2005 winters, respectively. This emphasizes the importance of the destratification caused by the seasonal forcing.

A note of interest for the 2005 and 2013 winters is the occurrence of the δSI minimum. In the 2013 winter, the minimum occurs significantly before deep convection (in December). For the 2005 winter, the minimum occurs roughly about the time of deep convection (around the beginning of March; seen more clearly in the MLD and not shown here), and also during a mistral event, much like the 2013 winter. However, unlike the 2013 winter, the seasonal destratification is less active at the time of deep convection, whereas the anomaly destratification drops almost −0.62 m2 s−2, most of it occurring during a larger mistral event, to start deep convection. This suggests that the mistral event occurring during this time triggers the deep-convection event. While deep convection does not occur in the 1994 winter, the δSI minimum is also at about the same time as the minimum in the control SI, with a small mistral event occurring at that date and with seasonal destratification remaining roughly constant. This suggests that if this year had further seasonal destratification or less initial destratification, the mistral may have been the main trigger for deep convection, along with the larger mistral event preceding it.

A note of interest for all three winters is the location of the seasonal SI minimum and the time of deep convection (or minimum of control SI for 1994). For the 2013 winter, the minimum is at roughly the time of deep convection but is almost a month later in the 2005 winter. And for the 1994 winter, the minimum occurs before the Control minimum. This begs the question of whether the location of the Seasonal SI minimum relative to the Control SI minimum is important, and, if this is the case, how important it is in terms of deep convection occurring or not occurring.

Figure 15The stratification index for the 1994 winter for both the control and seasonal runs are shown in (a) and are spatially averaged over the area of 42 to 42.5 N and 4.25 to 5 E. The simplified model anomaly solution added to the seasonal SI is denoted by the dashed line. Panel (b) shows the NEMO-determined δSI and the δSI calculated from the anomaly solution of the simple model.


Figure 16The same as Fig. 15 but for the 2005 winter.


5.2 Seasonal forcing

The seasonal sea surface fluxes for both years resemble the fluxes of the 2013 winter (see Fig. 17). The solar radiation component drives the main shape of the SI time series, with the major component contributing to the asymmetry being the latent heat flux, followed by the sensible heat flux. However, in the 2005 winter, the latent heat flux has a larger heat loss value than the 1994 winter, reaching over 300 W m−2 vs. under 200 W m−2 in the latter, driving the net heat flux to be more negative and causing more destratification according to the simple model, resulting in deep convection encompassing a few days on either side of the beginning of March.

The simple model for the seasonal SI is fairly accurate for the 1994 winter, similar to the 2013 winter, but not quite as accurate for the 2005 winter (see Fig. 17). For the 2005 winter, the simple model deviates during the deepest part of the winter. This could be due to more advective behavior captured by the destratification with the larger spatial averaging, which is neglected in the seasonal component of the simple model.

Figure 17The smoothed seasonal sea surface fluxes spatially averaged over the 42 to 42.5 N and 4.25 to 5 E area, with the dashed green line denoting the estimated derivative w.r.t. time of the seasonal SI from the NEMO results, multiplied by 109 m4 J s−2. A negative value means heat is leaving the ocean.


5.3 Anomaly forcing

The simple model for the anomaly scale was calculated for both the 1994 and 2005 winters, following the same steps as in Sect. 4.3. The value of the restoration coefficients, αd and αa, were carried over from the 2013 winter analysis, with the mistral dates determined through the same process outlined in Sect. 2.3 and manually adjusted to fit the visual data (again with the same method described in Sect. 2.3). The mistral strengths of the events during the preconditioning phase were determined through the same process as in Sect. 4.3. The mistral dates and strengths are presented in Appendix B. The results are shown in Figs. 15 and 16, and the simple model follows the NEMO simulation results quite closely, only deviating majorly at extreme peaks and troughs, despite utilizing the restoration coefficients from the 2013 winter. This reinforces the importance of the mistral as a dominating factor for destratification on the anomaly timescale.

The different components of the mistral (strength, δF, duration, Δt, and period, Δτ) are separated in the same manner as for the 2013 winter to determine the main factor of the mistral leading to destratification according to these additional years. For the 1994 winter, the contributions due to strength, duration, and period are equal to -8.71×10-2, -6.39×10-2, and 5.83×10-2 m2 s−2 (recall that a larger period means less frequent Mistral events), respectively. For the 2005 winter, the contribution due to strength, duration, and period are equal to -1.51×10-1, -5.23×10-2, and 4.89×10-2 m2 s−2, respectively. The 2005 winter results have the same order of magnitude as the 2013 winter results, with all three winters having the same order of importance for the mistral attributes: first strength, then duration, and finally the length of the period. Only for the 1994 winter was the strength term found to not be as dominant as in the other years, with the term having the same order of magnitude as the other terms. However, as the order of importance was still the same for all 3 years, this aids the conclusion that the strength term of the mistral is generally its most important factor driving destratification.

6 Conclusions

The 2012–2013 deep-convection year (2013 winter) in the Gulf of Lion was investigated to determine the effect the Mistral winds have on deep convection. Two NEMO ocean simulations were run, one forced with unmodified WRF/ORCHIDEE atmospheric forcing (control) and one forced with atmospheric fields filtered to remove the mistral signature (seasonal). Separating the atmospheric forcing into the long-term and anomaly timescales revealed that the mistral winds do not act alone to destabilize the northwestern Mediterranean Sea. Both the seasonal atmospheric change, reflected in the long-term timescales, and the mistral winds, reflected in the anomaly timescales, combine to destabilize and destratify the water columns in the GOL in roughly equal amounts (favoring the seasonal change).

When the NEMO simulation results were probed further by developing a simple model, the simple model conveyed the underlying drivers of the long-term or seasonal timescale. The evolution of the seasonal stratification index is proportional to the net heat flux leaving the ocean. As the net heat flux follows the shape of the incoming solar radiation, the maximum and minimum values for the seasonal stratification index occur around 21 September and 21 March, respectively, i.e., the fall and spring equinoxes. Shifted negative by the latent, sensible, and longwave radiation heat fluxes, the net heat flux allows for a seasonal cycle of destratification during the winter and restratification during the summer. If any of the three negative shifting components are unable to cool the ocean surface enough, deep convection may fail to appear, unless the contribution of the mistral winds is able to compensate.

The simple model results go on to confirm the hypothesis that the mistral acts on the anomaly timescale to destratify the water column and is the primary driver in this timescale. These results further conveyed that additional mistral events need to be stronger in terms of heat transfer than previous events to create further destratification. Otherwise, no destratification, or even restratification, occurs. The simple model then goes on to reveal, after some additional derivation, that the most important part of a Mistral event is its strength, in regards to potential destratification. Changing the duration or frequency has an effect, but this effect is an order of magnitude smaller than changing the mistral strength.

Another 2 years were also studied with the same method of running a control simulation and a filtered atmospheric-forcing seasonal simulation: the year 1993–1994 (1994 winter) and the year 2004–2005 (2005 winter). These years were then also studied with the simple model framework. The 2005 winter featured a deep convection event like the 2013 winter, but the 1994 winter did not, allowing for the comparison between deep-convecting and non-deep-convecting years. The conclusions determined in the 2013 winter are largely supported by the results of the two additional years. The seasonal change in SI accounts for a larger part of the destratification, while the 2005 winter still required destratification from the anomaly scale to deeply convect. The solar radiation component of the seasonal forcing was also found to be the component giving the cyclical structure to the stratification index, with the latent and sensible heat fluxes creating the asymmetry. For the anomaly-scale forcing, the mistral strength was again found to be the dominating component leading to destratification, although its magnitude is slightly less pronounced in the 1994 winter.

However, some of the NEMO simulation results bring further questions. Section 3 noted that the seasonal change in stratification brought the preconditioned water column in the GOL to the point of deep convection simultaneously with a mistral event for the winter of 2012–2013, as the mistral-induced preconditioning had already passed its minimum destratification beforehand, with both the seasonal change and a mistral event acting to destratify at the moment of deep convection. In the 2005 winter, the maximum seasonal destratification had not yet occurred during the time of deep convection, but the mistral-induced destratification brought the stratification down to point of deep convection, triggering it with a mistral event. Despite deep convection not occurring, a similar structure appears in the SI and δSI for the winter of 1994. Additionally, the date of the seasonal SI minimum was different relative to the date of deep convection (control SI minimum for the 1994 winter) for all three winters.

This brings up two questions. The first is whether the mistral truly triggers deep convection for all deep-convection events or the change in the seasonal destratification at the time of deep convection is instead a more prominent factor. The second is what the importance of the location of the seasonal SI minimum is and whether it makes a difference in regard to the possibility of deep convection. Is it possible for the mistral-induced destratification to cause deep convection after the minimum has passed, despite being hindered by the restratifying seasonal SI.

Another question brought about by the simulation results is what effect the maximal stratification at the beginning of the preconditioning phase has on the ocean's ability to experience deep convection for a given year. The maximal stratification must be overcome by the mistral and seasonal forcing for deep convection to occur, with the ability of the forcing to do so varying per year. For our results, the 2013 winter had a larger maximum control SI than the other two winters investigated and is still deeply convected. Similarly, the 2005 winter had a larger maximum stratification than the winter of 1994: 1.79 vs. 1.73 m2 s−2. However, all three winters had a maximal stratification within 0.1 m2 s−2 of each other, which is only about 6 % of the maximum.

For example, for the 2004–2005 winter the atmospheric forcing was more than enough to overcome the initial stratification, and a milder winter would have lead to deep convection as well, according to other works (Grignon et al.2010). Somot et al. (2016) investigated initial stratification and over a longer time period than Grignon et al. (2010) (1995–2005 vs. 1980–2013); however, they calculated the stratification index at the beginning of December rather than at the beginning of September, where the maximum stratification of ∼1.83 m2 s−2 occurs for the 2012–2013 winter. By December, the SI has already dropped to ∼1 m2 s−2, which means almost half of the destratification has already occurred, with their calculation missing about half of the preconditioning phase. The same is true when looking at the 1994 and 2005 winters. A maximum SI calculated near the beginning of September may be more representative of the water column's ability to deeply convect and needs to be investigated. A related question is how the accumulation (or reduction) of stratification being transferred to the proceeding year affects deep convection in following years. For the winter of 2005, where very few of the preceding 15 years deeply convected due to milder winters, the result was warmer and saltier WDMW production (Herrmann et al.2010).

These questions are outside the scope of the current article, as they rely on investigating a large number of years to evaluate the inter-annual variability of the atmospheric forcing. The 2013 winter was an above-average year in terms of destratification, leading to deep convection, while multiple years in the 1990s saw minimal MLD growth (including the 1994 winter; see Somot et al. (2016), and references therein). This may have been due to an above-average number of mistral (stormy) days, as suggested by Somot et al. (2016), which coincides with our results, which show that 36 %, 61 %, and 54 % of the preconditioning days had a mistral event for the winters of 1994, 2005, and 2013, respectively. However, it may have also been due to a larger than average contribution from the seasonal forcing, as the seasonal SI saw more destratification than the anomaly stratification index, δSI, which is not clearly discernible using only three winters.

We believe the approach of separating the atmospheric forcing into the seasonal and anomaly components will reveal more answers to these questions over a larger set of years, and we are preparing additional work to address this. We hope this future work will provide us with more information on how the Gulf of Lion deep-convection system will evolve in the future.

Appendix A: Equations

A1 Simple stratification index model

The purpose of this simple model is to separate the seasonal-scale atmospheric forcing from the anomaly-scale forcing. We start the derivation with the energy equation for incompressible flow (White2011):

(A1) ρ c p D T D t = D q D t ,

where ρ is density, cp is the specific heat of the fluid with constant pressure, T is temperature, t is time, and q is energy per volume from heat. DDt is the material derivative.

In this model we are assuming that the heat transfer term is equal to the heat removed by the atmosphere:

(A2) D q D t = - q a ,

where qa is the volumetric heat forcing from the atmosphere. This leaves us with the following equation:

(A3) D T D t = - q a ρ c p .

The Brunt–Väisälä frequency is defined as follows:

(A4) N 2 = b z = - g ρ 0 ρ z .

Assuming a fluid whose density varies negatively proportionally to the temperature, ρ=-βT, which is an acceptable approximation as the density only varies only a few tenths of a kilogram per cubic meter and temperature only varies by about 10 C, we can describe N2 in terms of temperature as follows:

(A5) N 2 = b z = - g T 0 T z .

Introducing buoyancy as b=-gT0T, we can rearrange the energy equation into terms of buoyancy:

(A6) g T 0 D T D t = - q a ρ c p D b D t = - q a g ρ c p T 0 .

If we then differentiate by z, we can reorganize the equation in terms of the Brunt–Väisälä frequency, N2:

(A7) z D b D t = - q a g ρ c p T 0 D N 2 D t = - z q a g ρ c p T 0 ,

By renaming the atmospheric forcing term on the right-hand side of Eq. (A7) to F(t), we can make this equation easier to follow:

(A8) F ( t ) = z q a g ρ c p T 0 .

This brings us to

(A9) D N 2 D t = - F ( t ) .

The main assumption we make is that the ocean column is a linear system and responds to the long-term and short-term anomaly timescale atmospheric forcing independently:

(A10) N 2 = N S 2 + δ N 2 , F ( t ) = δ F ( t ) + F S ( t ) ,

which describes the response of N2 on the anomaly timescale, δN2:

(A11) D δ N 2 D t = - δ F ( t ) ,

and the response of N2 on the seasonal timescale, NS2:

(A12) D N S 2 D t = - F S ( t ) .

For the seasonal response, we further make the assumption that NS2 negligibly depends on the x, y, and z coordinate directions:

(A13) N S 2 t = - F S ( t ) .

If we want to connect the overall Brunt–Väisälä frequency, N2, to the seasonal one, NS2, we can formulate a restoring term, R, in terms of T or in terms of N2, following the steps mentioned above:

(A14) R = z g T 0 α ρ c p ( T - T S ) α ( N 2 - N S 2 ) ,

or with δN2=N2-NS2 as follows:

(A15) R = α δ N 2 ,

where α is the restoring term coefficient. Separating the material derivative into its time and advective components for Eq. (A11) gives

(A16) δ N 2 t + V ( δ N 2 ) = - δ F ( t ) .

We will replace the advective component, V⋅∇(δN2), with R, which essentially swallows the advective operation into the restoring coefficient, α. This results in the partial differential equation that we will study further:

(A17) δ N 2 t + α δ N 2 = - δ F ( t ) .

A1.1 Solution for seasonal SI

To solve the response of N2 for the seasonal timescale, given by Eq. (A13), we will assume NS2 is vertically homogeneous, giving us the stratification index response, or Eq. (13), through the use of Eq. (9). We can then separate FS(t) back out into its components:

(A18) F S ( t ) = g ρ c p T 0 q a , S z .

Dividing qa,S by D gives us the atmospheric cooling in terms of a surface flux, -Qnet,S. If we plug this relationship back into Eq. (13), we get

(A19) SI S t = D 2 g ρ c p T 0 Q net , S z ,

Integrating this equation by z gives us

(A20) 0 D SI S t z = D 2 g ρ c p T 0 0 D Q net , S z z SI S t D = D 2 g ρ c p T 0 Q net , S SI S t = g 2 ρ c p T 0 Q net , S ,

Therefore, we have SIS expressed in terms of Qnet,S.

A1.2 Solution with mistral forcing function

Focusing on just the anomaly timescale, we will assume the mistral is the primary source of forcing. To model the atmospheric cooling of the mistral, we will model the forcing function, δF(t), as a series of k pulse functions of magnitude δFk over a duration of Δtk and with a period of Δτk, as visualized in Fig. 10.

To solve the Brunt–Väisälä frequency response with the mistral pulse forcing function, we solve Eq. (A17) in a piecewise manner, with a solution for each section of the pulse function. We will also make the assumption that for each portion of the Mistral event, both during and after, the advective components (α) remain constant with respect to time. This leads to αd and αa representing the advective components during and after an event, respectively. During a Mistral event, [tk,tk+Δtk), we get

(A21) δ N 2 t + α d ( δ N 2 ) = - δ F ( t ) δ N k 2 ( t ) = - δ F k α d + c 0 e - α d t δ N k 2 ( t k ) = - δ F k α d + c 0 e - α d t k = δ N k - 1 2 ( t k ) c 0 = δ N k - 1 2 ( t k ) + δ F k α d e α d t k ,

with the following result:

(A22) δ N k 2 ( t ) = δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d ( t - t k ) e - α d ( t - t k ) ,

After the event, [tk+Δtk,tk+1):

(A23) δ N 2 t + α a ( δ N 2 ) = 0 δ N k 2 ( t ) = c 1 e - α a t δ N k 2 ( t k + Δ t k ) = c 1 e - α a ( t k + Δ t k ) = δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d Δ t k e - α d Δ t k c 1 = δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k + α a t k ,

With the following result:

(A24) δ N k 2 ( t ) = δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) ,

or to have the results put more succinctly:

(A25) δ N k 2 ( t ) = δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d ( t - t k ) e - α d ( t - t k ) [ t k , t k + Δ t k ) δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) [ t k + Δ t k , t k + Δ τ k ) .

A1.3δNk-12 initial condition

δNk-12(tk) is a recursive initial condition, as its initial condition is the event before it (and so on):

(A26) δ N k - 1 2 ( t k ) = δ N k - 2 2 ( t k - 1 ) + δ F k - 1 α d ( 1 - e α d Δ t k - 1 ) e ( α a - α d ) Δ t k - 1 - α a Δ τ k - 1 δ N k - 2 2 ( t k - 1 ) = δ N k - 3 2 ( t k - 2 ) + δ F k - 2 α d ( 1 - e α d Δ t k - 2 ) e ( α a - α d ) Δ t k - 2 - α a Δ τ k - 2 .

Therefore, δNk-12(tk) can be simplified in expression by combining the initial conditions:

(A27) δ N k - 1 2 ( t k ) = δ N k - m 2 ( t k - ( m - 1 ) ) e ( α a - α d ) i = 1 m - 1 Δ t k - i e - α a i = 1 m - 1 Δ τ k - i + j = 1 m - 1 δ F k j α d 1 - e α d Δ t k - j e ( α a - α d ) i = 1 j Δ t k - i e - α a i = 1 j Δ τ k - i .

If m=k and δN02=0 then

(A28) δ N k - 1 2 ( t k ) = j = 1 k - 1 δ F k j α d 1 - e α d Δ t k - j e ( α a - α d ) i = 1 j Δ t k - i e - α a i = 1 j Δ τ k - i .

Assuming δFk=δF, Δtkt, and Δτkτ for all k, or a periodic pulse function, then δNk-12 can be expressed as follows:

(A29) δ N k - 1 2 ( t k ) = δ F α d 1 - e α d Δ t j = 1 k - 1 e [ ( α a - α d ) Δ t - α a Δ τ ] j .

Taking the sum of a finite geometric series as follows:

(A30) n = 0 m - 1 r n = 1 - r m 1 - r , n = 0 m - 1 r n = n = 1 m - 1 r n + 1 , n = 1 m - 1 r n = n = 0 m - 1 r n - 1 ,

where r≠1. If r=e(αa-αd)Δt-αaΔτ, then

(A31) j = 1 k - 1 e [ ( α a - α d ) Δ t - α a Δ τ ] j = 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 ,

and we then we get

(A32) δ N k - 1 2 ( t k ) = δ F α d 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 ,

where (αa-αd)Δt-αaΔτ0. Plugging Eq. (A32) into Eqs. (A22) and (A24) results in the equation for the response of the Brunt–Väisälä frequency forced by a periodic pulse function:

(A33) δ N 2 ( t ) = δ F α d [ 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 + 1 - e α d ( t - t k ) ] e - α d ( t - t k ) [ t k , t k + Δ t ) δ F α d [ 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 + 1 - e α d Δ t ] e ( α a - α d ) Δ t - α a ( t - t k ) [ t k + Δ t , t k + Δ τ ) .

For the anomaly response, Eqs. (A25) and (A33), which assume a vertically homogeneous δN2 leads to the stratification index through Eq. (9), lead us to δSI being expressed as follows:

(A34) δ SI k ( t ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d ( t - t k ) e - α d ( t - t k ) [ t k , t k + Δ t k ) δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) [ t k + Δ t k , t k + Δ τ k ) ,


(A35) δ SI k ( t ) = D 2 2 δ F α d [ 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 + 1 - e α d ( t - t k ) ] e - α d ( t - t k ) [ t k , t k + Δ t ) D 2 2 δ F α d [ 1 - e α d Δ t 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] - 1 + 1 - e α d Δ t ] e ( α a - α d ) Δ t - α a ( t - t k ) [ t k + Δ t , t k + Δ τ ) ,

for the period pulse function case.

A2 Restoring coefficients αd and αa

The restoration coefficients, αd and αa, can be solved for the separate phases of a mistral event in Eq. (A34) by normalizing the equations during their respective phases. These normalized equations are then fitted to selected ideal mistral destratification and restratification cases that are highlighted in Table 1 and Fig. 8 to retrieve the values of the restoration coefficients.

A2.1 Restoration coefficient αd during a mistral

To solve for αd, we normalize Eq. (A34) for during a mistral event, [tk,tk+Δtk), with δSI given as follows:

(A36) δ SI k ( t ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d ( t - t k ) e - α d ( t - t k ) ,

We first reference the time, t, to the starting time of event k as t=t-tk, giving us

(A37) δ SI k ( t ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α 1 - e α d t e - α d t ,

Next, we normalize δSIk to the value of zero at t=0, resulting in δSIk,NI:

(A38) δ SI k , NI ( t ) = δ SI k ( t ) - δ SI k ( t = 0 ) = δ SI k - 1 ( t k ) + D 2 2 δ F k α d 1 - e α d t e - α d t - δ SI k - 1 ( t k ) = δ SI k - 1 ( t k ) + D 2 2 δ F k α d e - α d t - 1 ,

Following this, the height or magnitude of destratification for each event is normalized to 1, resulting in δSIk,NH:

(A39) δ SI k , NH ( t ) = δ SI k , NI ( t ) extremum δ SI k , NI ( t ) .

The extremum value for δSIk,NI(t) is when t=Δtk or at the end of the mistral event. This simplifies δSIk,NH to

(A40) δ SI k , NH ( t ) = δ SI k - 1 ( t k ) + D 2 2 δ F k α d e - α d t - 1 δ SI k - 1 ( t k ) + D 2 2 δ F k α d e - α d Δ t k - 1 = e - α d t - 1 e - α d Δ t k - 1 ,

To normalize the length of the event duration from 0 to 1, we then divide t by the event length, Δtk, which results in t′′:

(A41) t ′′ = t Δ t k t = t ′′ Δ t k .

Plugging t′′ into δSIk,NH(t) returns δSIk,NT:

(A42) δ SI k , NT ( t ′′ ) = e - α d t ′′ Δ t k - 1 e - α d Δ t k - 1 .

This final equation, δSIk,NT, can be used with a fitting function to solve for αd if Δtk is supplied.

A2.2 Restoration coefficient αa, after a mistral

To solve for αa, we normalize Eq. (A34) for after a mistral event, [tk+Δtk,tk+Δτk), with δSI given as

(A43) δ SI k ( t ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) ,

Referencing the time, t, to the end of the event, t′′′=t-(tk+Δtk) and plugging t′′′ into Eq. (A43), we get

(A44) δ SI k ( t ′′′ ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e - α d Δ t k e - α a t ′′′ ,

Normalizing the vertical intercept of δSIk(t′′′) results in δSIk,NI:

(A45) δ SI k , NI = δ SI k ( t ′′′ ) - δ SI k ( t ′′′ = 0 ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e - α d Δ t k e - α a t ′′′ - δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e - α d Δ t k = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e - α d Δ t k e - α a t ′′′ - 1 ,

Each post-event restratification is normalized to the height of 1 by dividing δSIk,NI by (δSIk(t′′′=Δτk-Δtk)-δSIk(t′′′=0)):

(A46) δ SI k , NH = δ SI k , NI δ SI k ( t ′′′ = Δ τ k - Δ t k ) - δ SI k ( t ′′′ = 0 ) = δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e - α d Δ t k e - α a t ′′′ - 1 δ SI k - 1 2 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e - α d Δ t k e - α a ( Δ τ k - Δ t k ) - 1 = e - α a t ′′′ - 1 e - α a ( Δ τ k - Δ t k ) - 1 ,

which gives us δSIk,NH.

To normalize the length of time of post-event restratification, we will divide t′′′ by the post-event time length, Δτk−Δtk, resulting in t′′′′:

(A47) t ′′′′ = t ′′′ Δ τ k - Δ t k t ′′′ = t ′′′′ ( Δ τ k - Δ t k ) ,

which leads to δSIk,NT:

(A48) δ SI k , NT = e - α a t ′′′′ ( Δ τ k - Δ t k ) - 1 e - α a ( Δ τ k - Δ t k ) - 1 .

This leaves us with an equation of αa that can be fitted against NEMO model data if Δtk and Δτk are provided.

The average duration and period of events during the preconditioning period in Table 1, Δt and Δτ, respectively, are used for Δt and Δτ in these normalized equations, Eqs. (A42) and (A48). The result of the fitting is shown in Fig. 11, with αd having a fitted valued of 0.235 d−1 and αa having a fitted value of 0.021 d−1.

A3 Determining δFk

With the restoring coefficients determined in Sect. A2 and the duration and period of each event available in Table 1, the strength of each mistral event, δFk, can be determined. If we take Eq. (A36) and note that the value of δSIk−1 is the same as δSIk(tk) at the beginning of an event, we can simplify the equation in the following steps:

(A49) δ SI k ( t k ) = δ SI k - 1 , δ SI k ( t k + Δ t k ) = δ SI k - 1 e - α d Δ t k + D 2 2 δ F k α d e - α d Δ t k - 1 ,

and then solve for δFk:

(A50) δ F k = 2 δ SI k ( t k + Δ t k ) - δ SI k ( t k ) e - α d Δ t k α d e - α d Δ t k - 1 D 2 .

The results of δFk for each event in the preconditioning phase are given in Table A1, along with mean and standard deviation.

Table A1The mistral strengths, δFk, for each of the preconditioning phase events, using αd and αa from Sect. A2 and the rest of the preconditioning period mistral characteristics from Table 1 plugged into Eq. (A50).

δFk=4.01×10-8 s−2 d−1, and σδFk=1.196×10-8 s−2 d−1.

Download Print Version | Download XLSX

A4 Time derivative of δN2

Taking the derivative with respect to time from Eqs. (A25) and (A34) results in

(A51) δ N k 2 ( t ) t = - α d δ N k - 1 2 ( t k ) + δ F k α d e - α d ( t - t k ) [ t k , t k + Δ t k ) - α a δ N k - 1 2 ( t k ) + δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) [ t k + Δ t k , t k + Δ τ k ) ,


(A52) δ SI k ( t ) t = - α d δ SI k - 1 ( t k ) + D 2 2 δ F k α d e - α d ( t - t k ) [ t k , t k + Δ t k ) - α a δ SI k - 1 ( t k ) + D 2 2 δ F k α d 1 - e α d Δ t k e ( α a - α d ) Δ t k - α a ( t - t k ) [ t k + Δ t k , t k + Δ τ k ) .

A5 Asymptotic destratification

The following sections differentiate Eq. (A35) at t=tk+Δtk or at the end of a mistral event, where the destratification is the largest, using k or using the other components, i.e., δF, Δt, and Δτ, once k→∞.


Equation (A35), at t=tk+Δt results in

(A53) δ SI k ( t k + Δ t ) = D 2 2 δ F α d e - α d Δ t - 1 1 - e [ ( α a - α d ) Δ t - α a Δ τ ] k 1 - e ( α a - α d ) Δ t - α a Δ τ ,

The derivative of Eq. (A53) with respect to (w.r.t.) k is

(A54) δ SI k k = D 2 2 δ F α d e - α d Δ t - 1 ( 1 - e ( α a - α d ) Δ t - α a Δ τ ) - e [ ( α a - α d ) Δ t - α a Δ τ ] k [ ( α a - α d ) Δ t - α a Δ τ ] ,

As k→∞, with αd>αa Eq. (A53) goes to

(A55) δ SI = D 2 2 δ F α d e - α d Δ t - 1 1 1 - e ( α a - α d ) Δ t - α a Δ τ .


The derivative of Eq. (A55) w.r.t. δF is

(A56) δ SI δ F = D 2 2 1 α d e - α d Δ t - 1 1 1 - e ( α a - α d ) Δ t - α a Δ τ .


The derivative of Eq. (A55) w.r.t. Δt is

(A57) δ SI Δ t = D 2 2 δ F α d 1 ( 1 - e ( α a - α d ) Δ t - α a Δ τ ) [ - α d e - α d Δ t + e - α d Δ t - 1 ( α a - α d ) e ( α a - α d ) Δ - α a Δ τ 1 - e ( α a - α d ) Δ t - α a Δ τ ] .


The derivative of Eq. (A55) w.r.t. Δτ is

(A58) δ SI Δ τ = D 2 2 δ F α d e - α d Δ t - 1 - α a e ( α a - α d ) Δ t - α a Δ τ 1 - e ( α a - α d ) Δ t - α a Δ τ 2 .
Appendix B: The 1994 and 2005 data

The mistral event data for the additional 2 years of 1993–1994 and 2004–2005 are found below in Tables B1 and B2.

Table B1The start date, duration (Δtk, in d), and period between each event (Δτk, in d) for each Mistral event (k) and event strength, (δFk, in s−2 d−1) for the preconditioning phase of the NEMO simulation period of 1 June 1993 to 31 May 1994.

The preconditioning phase for this year is considered to range from 28 August 1993 to 3 April 1994 (210 d).

Download Print Version | Download XLSX

Table B2The same as Table B1 but for the NEMO simulation period from 1 June 2004 to 31 May 2005.

The preconditioning phase for this year is considered to range from 15 September 2004 to 8 April 2005 (170 d).

Download Print Version | Download XLSX

Code availability

The RegIPSL model and the NEMO model and its configuration used in this paper can be accessed at the following Git repository: (Drobinski et al.2012). Code used to generate figures and implement the simple model can be found at the following Git repository: (Keller2022).

Data availability

The ERA-Interim Reanalysis can be found at (Dee et al.2011). HyMeX data are available at with registration (last access: 24 March 2022) (Drobinski et al.2014; Taupier-Letage2013). The RegIPSL and NEMO dataset can be made available by the corresponding author upon request.

Author contributions

DK Jr. performed the majority of the analysis presented. YG provided analysis with respect to mistral events. RP provided assistance and guidance with running the ocean simulations. SRR and PD provided analytical assistance and guidance.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Hydrological cycle in the Mediterranean (ACP/AMT/GMD/HESS/NHESS/OS inter-journal SI)”. It is not associated with a conference.


This work is a contribution to the HyMeX program (HYdrological cycle in the Mediterranean EXperiment) through INSU-MISTRALS support and the Med-CORDEX program (COordinated Regional climate Downscaling EXperiment – Mediterranean region). It was also supported by a joint CNRS–Weizmann Institute of Science collaborative project. The authors acknowledge Météo-France for supplying the data and the HyMeX database teams (ESPRI/IPSL and SEDOO/Observatoire Midi-Pyrénées) for their help in accessing the data.

Financial support

This research has been supported by the Centre National de la Recherche Scientifique (CNRS) in the framework of a collaborative project shared between the CNRS and the Weizmann Institute of Science.

Review statement

This paper was edited by Markus Meier and reviewed by two anonymous referees.


Arsouze, T., Beuvier, J., Béranger, K., Bourdallé-Badie, R., Deltel, C., Drillet, Y., Drobinski, P., Ferry, N., Lebeaupin-Brossier, C., Lyard, F., Sevault, F., and Somot, S.: Release note of the high-resolution oceanic model in the Mediterranean Sea NEMO-MED12 based on NEMO 3.2 version, Tech. rep., Centre National de Recherches Météorologiques, (last access: 24 March 2022), 2012. a

Balmaseda, M. A., Trenberth, K. E., and Källén, E.: Distinctive climate signals in reanalysis of global ocean heat content, Geophys. Res. Lett., 40, 1754–1759,, 2013. a

Béranger, K., Testor, P., and Crépon, M.: Modelling water mass formation in the Gulf of Lions (Mediterranean Sea), CIESM Workshop Monographs, Workshop, 27–30 May 2009, Malta, 2009. a

Beuvier, J., Béranger, K., Brossier, C. L., Somot, S., Sevault, F., Drillet, Y., Bourdallé-Badie, R., Ferry, N., and Lyard, F.: Spreading of the Western Mediterranean Deep Water after winter 2005: Time scales and deep cyclone transport, J. Geophys. Res.-Oceans, 117, C07022,, 2012. a, b, c, d, e, f

Coppola, L., Prieur, L., Taupier-Letage, I., Estournel, C., Testor, P., Lefevre, D., Belamari, S., LeReste, S., and Taillandier, V.: Observation of oxygen ventilation into deep waters through targeted deployment of multiple Argo-O2 floats in the north-western Mediterranean Sea in 2013, J. Geophys. Res.-Oceans, 122, 6325–6341,, 2017. a

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597,, 2011 (data available at:, last access: 24 March 2022). a, b

Drobinski, P., Ducrocq, V., Alpert, P., Anagnostou, E., Béranger, K., Borga, M., Braud, I., Chanzy, A., Davolio, S., Delrieu, G., Estournel, C., Boubrahmi, N. F., Font, J., Grubišić, V., Gualdi, S., Homar, V., Ivančan-Picek, B., Kottmeier, C., Kotroni, V., Lagouvardos, K., Lionello, P., Llasat, M. C., Ludwig, W., Lutoff, C., Mariotti, A., Richard, E., Romero, R., Rotunno, R., Roussot, O., Ruin, I., Somot, S., Taupier-Letage, I., Tintore, J., Uijlenhoet, R., and Wernli, H.: HyMeX: A 10-Year Multidisciplinary Program on the Mediterranean Water Cycle, B. Am. Meteorol. Soc., 95, 1063–1082,, 2014 (data available at:, last access: 24 March 2022, registration required). a, b, c

Drobinski, P., Anav, A., Lebeaupin Brossier, C., Samson, G., Stéfanon, M., Bastin, S., Baklouti, M., Béranger, K., Beuvier, J., Bourdallé-Badie, R., Coquart, L., D'Andrea, F., de Noblet-Ducoudré, N., Diaz, F., Dutay, J.-C., Ethe, C., Foujols, M.-A., Khvorostyanov, D., Madec, G., Mancip, M., Masson, S., Menut, L., Palmieri, J., Polcher, J., Turquety, S., Valcke, S., and Viovy, N.: Model of the Regional Coupled Earth system (MORCE): Application to process and climate studies in vulnerable regions, Environ. Modell. Softw., 35, 1–18,, 2012 (code available at:, last access: 24 March 2022). a

Drobinski, P., Alonzo, B., Basdevant, C., Cocquerez, P., Doerenbecher, A., Fourrié, N., and Nuret, M.: Lagrangian dynamics of the mistral during the HyMeX SOP2, J. Geophys. Res.-Atmos., 122, 1387–1402,, 2017. a

Estournel, C., , Testor, P., Taupier-Letage, I., Bouin, M.-N., Coppola, L., Durand, P., Conan, P., Bosse, A., Brilouet, P.-E., Beguery, L., Belamari, S., Béranger, K., Beuvier, J., Bourras, D., Canut, G., Doerenbecher, A., de Madron, X. D., D'Ortenzio, F., Drobinski, P., Ducrocq, V., Fourrié, N., Giordani, H., Houpert, L., Labatut, L., Brossier, C. L., Nuret, M., Prieur, L., Roussot, O., Seyfried, L., and Somot, S.: HyMeX-SOP2: The Field Campaign Dedicated to Dense Water Formation in the Northwestern Mediterranean, Oceanography, 29, 196–206,, 2016a. a

Estournel, C., Testor, P., Damien, P., D'Ortenzio, F., Marsaleix, P., Conan, P., Kessouri, F., de Madron, X. D., Coppola, L., Lellouche, J.-M., Belamari, S., Mortier, L., Ulses, C., Bouin, M.-N., and Prieur, L.: High resolution modeling of dense water formation in the north-western Mediterranean during winter 2012–2013: Processes and budget, J. Geophys. Res.-Oceans, 121, 5367–5392,, 2016b. a, b, c, d

Flamant, C.: Alpine lee cyclogenesis influence on air-sea heat exchanges and marine atmospheric boundary layer thermodynamics over the western Mediterranean during a Tramontane/Mistral event, J. Geophys. Res.-Oceans, 108, 8057,, 2003. a

Gascard, J. C.: Mediterranean deep water formation baroclinic instability and oceanic eddies, Oceanol. Acta, 1, 315–330, 1978. a, b

Givon, Y., Keller Jr., D., Silverman, V., Pennel, R., Drobinski, P., and Raveh-Rubin, S.: Large-scale drivers of the mistral wind: link to Rossby wave life cycles and seasonal variability, Weather Clim. Dynam., 2, 609–630,, 2021. a, b, c, d

Grignon, L., Smeed, D. A., Bryden, H. L., and Schroeder, K.: Importance of the variability of hydrographic preconditioning for deep convection in the Gulf of Lion, NW Mediterranean, Ocean Sci., 6, 573–586,, 2010. a, b

Guion, A., Turquety, S., Polcher, J., Pennel, R., Bastin, S., and Arsouze, T.: Droughts and heatwaves in the Western Mediterranean: impact on vegetation and wildfires using the coupled WRF-ORCHIDEE regional model (RegIPSL), Clim. Dynam.,, online first, 2021. a

Hamon, M., Beuvier, J., Somot, S., Lellouche, J.-M., Greiner, E., Jordà, G., Bouin, M.-N., Arsouze, T., Béranger, K., Sevault, F., Dubois, C., Drevillon, M., and Drillet, Y.: Design and validation of MEDRYS, a Mediterranean Sea reanalysis over the period 1992–2013, Ocean Sci., 12, 577–599,, 2016. a, b

Herrmann, M., Sevault, F., Beuvier, J., and Somot, S.: What induced the exceptional 2005 convection event in the northwestern Mediterranean basin? Answers from a modeling study, J. Geophys. Res., 115, C12051,, 2010. a, b, c, d

Houpert, L., de Madron, X. D., Testor, P., Bosse, A., D'Ortenzio, F., Bouin, M. N., Dausse, D., Goff, H. L., Kunesch, S., Labaste, M., Coppola, L., Mortier, L., and Raimbault, P.: Observations of open-ocean deep convection in the northwestern Mediterranean Sea: Seasonal and interannual variability of mixing and deep water masses for the 2007–2013 Period, J. Geophys. Res.-Oceans, 121, 8139–8171,, 2016. a, b, c

Keller Jr., D.: Code to Analyze and Plot Data for Untangling the mistral and seasonal atmospheric forcing driving deep convection in the Gulf of Lion: 2012–2013, GitLab [code],, last access: 5 April 2022. a

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., and Prentice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Global Biogeochem. Cy., 19, GB1015,, 2005. a

Large, W. and Yeager, S.: Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies, Tech. rep., University Corporation for Atmospheric Research,, 2004. a, b

Large, W. G. and Yeager, S. G.: The global climatology of an interannually varying air–sea flux data set, Clim. Dynam., 33, 341–364,, 2008. a

Lebeaupin-Brossier, C. and Drobinski, P.: Numerical high-resolution air-sea coupling over the Gulf of Lions during two tramontane/mistral events, J. Geophys. Res., 114, D10110,, 2009. a

Lebeaupin-Brossier, C., Béranger, K., Deltel, C., and Drobinski, P.: The Mediterranean response to different space–time resolution atmospheric forcings using perpetual mode sensitivity simulations, Ocean Model., 36, 1–25,, 2011. a, b

Lebeaupin-Brossier, C., Béranger, K., and Drobinski, P.: Sensitivity of the northwestern Mediterranean Sea coastal and thermohaline circulations simulated by the 1/12-resolution ocean model NEMO-MED12 to the spatial and temporal resolution of atmospheric forcing, Ocean Model., 43-44, 94–107,, 2012. a, b

Lebeaupin-Brossier, C., Léger, F., Giordani, H., Beuvier, J., Bouin, M.-N., Ducrocq, V., and Fourrié, N.: Dense water formation in the north-western M editerranean area during HyMeX-SOP2 in 1/36 ocean simulations: Ocean-atmosphere coupling impact, J. Geophys. Res.-Oceans, 122, 5749–5773,, 2017. a

Léger, F., Brossier, C. L., Giordani, H., Arsouze, T., Beuvier, J., Bouin, M.-N., Bresson, É., Ducrocq, V., Fourrié, N., and Nuret, M.: Dense water formation in the north-western Mediterranean area during HyMeX-SOP2 in 1/36 ocean simulations: Sensitivity to initial conditions, J. Geophys. Res.-Oceans, 121, 5549–5569,, 2016. a

L'Hévéder, B., Li, L., Sevault, F., and Somot, S.: Interannual variability of deep convection in the Northwestern Mediterranean simulated with a coupled AORCM, Clim. Dynam., 41, 937–960,, 2012. a

Ludwig, W., Dumont, E., Meybeck, M., and Heussner, S.: River discharges of water and nutrients to the Mediterranean and Black Sea: Major drivers for ecosystem changes during past and future decades?, Prog. Oceanogr., 80, 199–217,, 2009. a

Madec, G., Chartier, M., and Crépon, M.: The effect of thermohaline forcing variability on deep water formation in the western Mediterranean Sea: a high-resolution three-dimensional numerical study, Dynam. Atmos. Oceans, 15, 301–332,, 1991a. a

Madec, G., Delecluse, P., Crepon, M., and Chartier, M.: A Three-Dimensional Numerical Study of Deep-Water Formation in the Northwestern Mediterranean Sea, J. Phys. Oceanogr., 21, 1349–1371,<1349:atdnso>;2, 1991b. a

Madec, G., Lott, F., Delecluse, P., and Crépon, M.: Large-Scale Preconditioning of Deep-Water Formation in the Northwestern Mediterranean Sea, J. Phys. Oceanogr., 26, 1393–1408,<1393:LSPODW>2.0.CO;2, 1996. a

Marshall, J. and Schott, F.: Open-ocean convection: Observations, theory, and models, Rev. Geophys., 37, 1–64,, 1999. a, b, c, d, e, f

MEDOC: Observation of Formation of Deep Water in the Mediterranean Sea, 1969, Nature, 227, 1037–1040,, 1970. a, b, c

Mertens, C. and Schott, F.: Interannual Variability of Deep-Water Formation in the Northwestern Mediterranean, J. Phys. Oceanogr., 28, 1410–1424,<1410:ivodwf>;2, 1998. a, b

Millot, C. and Taupier-Letage, I.: Circulation in the Mediterranean Sea, in: The Mediterranean Sea, edited by: Saliot, A., 29–66, Springer Berlin Heidelberg,, 2005. a

Noh, Y., Cheon, W. G., and Raasch, S.: The Role of Preconditioning in the Evolution of Open-Ocean Deep Convection, J. Phys. Oceanogr., 33, 1145–1166,<1145:tropit>;2, 2003. a

Robinson, A., Leslie, W., Theocharis, A., and Lascaratos, A.: Mediterranean Sea Circulation, in: Encyclopedia of Ocean Sciences, edited by: Steele, J. H., 1689–1705, Elsevier,, 2001. a

Ruti, P. M., Somot, S., Giorgi, F., Dubois, C., Flaounas, E., Obermann, A., Dell'Aquila, A., Pisacane, G., Harzallah, A., Lombardi, E., Ahrens, B., Akhtar, N., Alias, A., Arsouze, T., Aznar, R., Bastin, S., Bartholy, J., Béranger, K., Beuvier, J., Bouffies-Cloché, S., Brauch, J., Cabos, W., Calmanti, S., Calvet, J.-C., Carillo, A., Conte, D., Coppola, E., Djurdjevic, V., Drobinski, P., Elizalde-Arellano, A., Gaertner, M., Galàn, P., Gallardo, C., Gualdi, S., Goncalves, M., Jorba, O., Jordà, G., L'Heveder, B., Lebeaupin-Brossier, C., Li, L., Liguori, G., Lionello, P., Maciàs, D., Nabat, P., Önol, B., Raikovic, B., Ramage, K., Sevault, F., Sannino, G., Struglia, M. V., Sanna, A., Torma, C., and Vervatis, V.: Med-CORDEX Initiative for Mediterranean Climate Studies, B. Am. Meteorol. Soc., 97, 1187–1208,, 2016. a

Schott, F., Visbeck, M., Send, U., Fischer, J., Stramma, L., and Desaubies, Y.: Observations of Deep Convection in the Gulf of Lions, Northern Mediterranean, during the Winter of 1991/92, J. Phys. Oceanogr., 26, 505–524,<0505:oodcit>;2, 1996. a

Send, U. and Testor, P.: Direct Observations Reveal the Deep Circulation of the Western Mediterranean Sea, J. Geophys. Res.-Oceans, 122, 10091–10098,, 2017. a

Severin, T., Kessouri, F., Rembauville, M., Sánchez-Pérez, E. D., Oriol, L., Caparros, J., Pujo-Pay, M., Ghiglione, J.-F., D'Ortenzio, F., Taillandier, V., Mayot, N., Madron, X. D. D., Ulses, C., Estournel, C., and Conan, P.: Open-ocean convection process: A driver of the winter nutrient supply and the spring phytoplankton distribution in the Northwestern Mediterranean Sea, J. Geophys. Res.-Oceans, 122, 4587–4601,, 2017. a, b

Skamarock, W., Klemp, J., Dudhia, J., Gill, D., Barker, D., Wang, W., Huang, X.-Y., and Duda, M.: A Description of the Advanced Research WRF Version 3, Tech. rep., University Corporation for Atmospheric Research,, 2008. a

Somot, S.: Modélisation Climatique du Bassin Méditerranéen: Variabilité et Scénarios de Changement Climatique, PhD thesis, Université Toulouse III, (last access: 24 March 2022), 2005. a

Somot, S., Houpert, L., Sevault, F., Testor, P., Bosse, A., Taupier-Letage, I., Bouin, M.-N., Waldman, R., Cassou, C., Sanchez-Gomez, E., de Madron, X. D., Adloff, F., Nabat, P., and Herrmann, M.: Characterizing, modelling and understanding the climate variability of the deep water formation in the North-Western Mediterranean Sea, Clim. Dynam., 51, 1179–1210,, 2016. a, b, c, d, e, f, g, h, i

Song, X. and Yu, L.: Air-sea heat flux climatologies in the Mediterranean Sea: Surface energy balance and its consistency with ocean heat storage, J. Geophys. Res.-Oceans, 122, 4068–4087,, 2017. a

Taupier-Letage, I.: CTD SOP2, Provence – Tethys 2, SEDOO OMP [data set],, 2013 (, last access: 24 March 2022, registration required). a, b

Testor, P. and Gascard, J.-C.: Large-Scale Spreading of Deep Waters in the Western Mediterranean Sea by Submesoscale Coherent Eddies, J. Phys. Oceanogr., 33, 75–87,<0075:lssodw>;2, 2003. a

The Lab Sea Group: The Labrador Sea Deep Convection Experiment, B. Am. Meteorol. Soc., 79, 2033–2058,<2033:tlsdce>;2, 1998. a

Turner, J. S.: Buoyancy Effects in Fluids, Cambridge University Press,, 1973. a

Waldman, R., Brüggemann, N., Bosse, A., Spall, M., Somot, S., and Sevault, F.: Overturning the Mediterranean Thermohaline Circulation, Geophys. Res. Lett., 45, 8407–8415,, 2018. a

White, F.: Viscous fluid flow, McGraw-Hill Higher Education, New York, NY, ISBN-13 978-0072402315, 2011. a, b

Short summary
The mistral winds are believed to be the primary source of cooling of the Gulf of Lion, leading to deep convection in the region, a process that mixes the ocean column from the seafloor to the sea surface. However, we have found that seasonal atmospheric changes also significantly cool the Gulf of Lion waters to cause deep convection, rather than mistral winds being the sole source, contributing roughly two-thirds of the required cooling, with the mistral winds contributing the final third.