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

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

► To cite this version:

Douglas Keller, Yonatan Givon, Romain Pennel, Shira Raveh-Rubin, Philippe Drobinski. Untangling the mistral and seasonal atmospheric forcing driving deep convection in the Gulf of Lion: 2012-2013. Ocean Science, 2022, 18, pp.483-510. 10.5194/os-18-483-2022. insu-03726902

HAL Id: insu-03726902

https://insu.hal.science/insu-03726902

Submitted on 22 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

Douglas Keller Jr.<fnref n="1" />, Yonatan Givon<fnref n="2" />, Romain Pennel<fnref n="1" />, Shira Raveh-Rubin<fnref n="2" />, and Philippe Drobinski<fnref n="1" />

**Correspondence:** Douglas Keller Jr. (dg.kllr.jr@gmail.com)

Received: 22 July 2021 – Discussion started: 9 August 2021

Revised: 14 January 2022 – Accepted: 2 February 2022 – Published: 8 April 2022

**Abstract.** 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 Schott, 1999). 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 (MEDOC, 1970), with some transported along the northern boundary current towards the Balearic Sea (Send and Testor, 2017) and some transported to the south within eddies (Beuviér et al., 2012; Testor and Gascard, 2003) 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 Schott, 1999; Mertens and Schott, 1998), 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; Flamant, 2003). These cooling and evaporation events destabilize the water column in the GOL and are widely accepted to be the primary source of buoyancy loss

<footnote n="1">¹LMD/IPSL, École Polytechnique, Institut Polytechnique de Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Palaiseau, France</footnote>

<footnote n="2">²Department of Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot, Israel</footnote>
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 Drobinski, 2009; Noh et al., 2003; Marshall and Schott, 1999; Mertens and Schott, 1998; Madec et al., 1996; Schott et al., 1996; Madec et al., 1991a, b; Gascard, 1978).

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 (MEDOC, 1970; The Lab Sea Group, 1998). 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 (https://www.nemo-ocean.eu/, 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, $Q_E$, the sensible heat flux, $Q_H$, the longwave radiation heat flux, $Q_{LW}$, and the surface shear stress, $\tau$:

$$Q_E = \rho_0 \Delta C_E (q_z - q_0) |\Delta \mathbf{u}|, \quad (1)$$

$$Q_H = \rho_0 c_p C_H (\theta_z - \text{SST}) |\Delta \mathbf{u}|, \quad (2)$$

$$Q_{LW} = Q_{LW,a} - \epsilon \sigma \text{SST}_K^4, \quad (3)$$

$$\tau = \rho_0 C_D \Delta \mathbf{u} |\Delta \mathbf{u}|. \quad (4)$$

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. $\mathbf{u}$ is the horizontal wind vector, with $\Delta \mathbf{u} = \mathbf{u}_z - \mathbf{u}_0$ as the difference between the wind velocity and sea surface current (assuming a no-slip condition at the ocean surface). $q$ and $\theta$ are the specific humidity and potential temperature of air, respectively. $\Lambda$ and $c_p$ are the latent heat of evaporation and the specific heat of water, respectively. $\rho_0$ is the air density at the sea surface. SST is the sea surface temperature, with SST$_K$ as the sea surface absolute temperature. $\epsilon$ is the sea surface emissivity, $\sigma$ is the Stefan–Boltzmann constant, and $Q_{LW,a}$ is the atmospheric longwave radiation. $C_E$, $C_H$, and $C_D$ 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, $Q_{net}$, is described by the summation of the terms in the following equation (Large and Yeager, 2004; Estournel et al., 2016b):

$$Q_{net} = Q_{SW} + Q_{LW} + Q_{H} + Q_{E}, \qquad (5)$$

where $Q_{SW}$ 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).

## 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 long-wave 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. 1–4). 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:

$$\chi_i = \frac{1}{i + N + 1} \sum_{j=0}^{i + N} x_j, \qquad (6)$$

where $\chi_i$ is the averaged (filtered) value at index $i$ of a time series of variable $x$ with length $n$, where $i = 0 \rightarrow n$. 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
<figure><img src="image_4.png" /><figcaption>Figure 1. The 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.</figcaption></figure>

“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-Letage, 2005), 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.

## 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⁻¹ 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, Δt_k, 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).

## 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; Somot, 2005; sometimes called the “convection resistance”). It is derived from the non-penetrative growth of the mixed-layer depth (MLD; i.e., without entrainment; Turner, 1973), which has been shown to be an accurate approximation for open-ocean convection (Marshall and Schott, 1999):

$$ \frac{\partial z}{\partial t} = \frac{B(t)}{N^2(z)z}, \qquad (7) $$

where N² is the Brunt-Väisälä frequency, z is the vertical coordinate along the water column, ∂z/∂t is the growth of the mixed-layer depth, and B is the potential buoyancy loss the water column can endure before removing stratification (in m² s⁻³). Separating by variable and integrating results in the equation for SI gives

$$ SI = \int_{0}^{D} N^{2} zdz, \qquad (8) $$

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

$$ SI = \frac{D^2}{2} N^2. \qquad (9) $$

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.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** An 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.</figcaption></figure>

<table><caption>**Table 1.** The start date, duration (Δ*t*<sub>k</sub>), and period between each event (Δτ<sub>k</sub>) 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 α<sub>d</sub> and α<sub>a</sub>, respectively, in Sect. 4.3 and Appendix A2.1 and A2.2.</caption><thead><tr><th>Start date<br>(yyyy-mm-dd)</th><th>Δt<sub>k</sub><br>(d)</th><th>Δτ<sub>k</sub><br>(d)</th><th>Start date<br>(yyyy-mm-dd)</th><th>Δt<sub>k</sub><br>(d)</th><th>Δτ<sub>k</sub><br>(d)</th></tr></thead><tbody><tr><td>2012-08-03</td><td>1</td><td>3</td><td>2012-12-26</td><td>5</td><td>7</td></tr><tr><td>2012-08-06</td><td>1</td><td>2</td><td>2013-01-02</td><td>17</td><td>21</td></tr><tr><td>2012-08-08</td><td>1</td><td>5</td><td>2013-01-23</td><td>6</td><td>10</td></tr><tr><td>2012-08-13</td><td>1</td><td>12</td><td>2013-02-02<sup>d</sup></td><td>15</td><td>18</td></tr><tr><td>2012-08-25</td><td>2</td><td>5</td><td>2013-02-20</td><td>7</td><td>10</td></tr><tr><td>2012-08-30<sup>d,a</sup></td><td>8</td><td>13</td><td>2013-03-02</td><td>1</td><td>11</td></tr><tr><td>2012-09-12<sup>d,a</sup></td><td>4</td><td>7</td><td>2013-03-13</td><td>3</td><td>7</td></tr><tr><td>2012-09-19<sup>d,a</sup></td><td>2</td><td>9</td><td>2013-03-20</td><td>1</td><td>6</td></tr><tr><td>2012-09-28<sup>d,a</sup></td><td>5</td><td>14</td><td>2013-03-26</td><td>1</td><td>5</td></tr><tr><td>2012-10-12<sup>d,a</sup></td><td>4</td><td>15</td><td>2013-03-31</td><td>1</td><td>6</td></tr><tr><td>2012-10-27<sup>d</sup></td><td>5</td><td>15</td><td>2013-04-06</td><td>2</td><td>13</td></tr><tr><td>2012-11-11<sup>d,a</sup></td><td>3</td><td>8</td><td>2013-04-19</td><td>4</td><td>8</td></tr><tr><td>2012-11-19</td><td>2</td><td>8</td><td>2013-04-27</td><td>1</td><td>25</td></tr><tr><td>2012-11-27<sup>d</sup></td><td>6</td><td>11</td><td>2013-05-22</td><td>2</td><td>10</td></tr><tr><td>2012-12-08<sup>d</sup></td><td>4</td><td>9</td><td>2013-06-01</td><td>2</td><td>23</td></tr><tr><td>2012-12-17<sup>d</sup></td><td>3</td><td>4</td><td>2013-06-24</td><td>1</td><td>4</td></tr><tr><td>2012-12-21</td><td>2</td><td>5</td><td>2013-06-28</td><td>1</td><td>41</td></tr></tbody></table>

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 σ<sub>Δt</sub> = 4.22 d and σ<sub>Δτ</sub> = 4.59 d.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** SST 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.</figcaption></figure>

# 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 (https://mistrals.sedoo.fr/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. 1–4).

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-Letage, 2013; 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

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Locations 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.</figcaption></figure>

profiles and their standard deviations are visualized in Figs. 5 and 6. The locations of the CTD profiles are shown in Fig. 4.

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⁻³ 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⁻² °C are potentially all that is required to sustain an ocean convective cycle (Marshall and Schott, 1999), and density differences for the same order of magnitude (10⁻² kg m⁻³) 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 (Gascard, 1978). 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.

## 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.

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 Schott, 1999; MEDOC, 1970). 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 + SI_S$, while the stratification index of the seasonal simulation is given as $SI_S$. 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+SI_S$, $SI_S$, and $δSI$ are all shown in Fig. 8.

Both the control and seasonal runs start off with an SI value of $1.57 m^2 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 $m^2 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 + SI_S$ is equal to 0 (also signifying when the MLD reaches the sea floor), and the lack of deep convection in the seasonal run, as $SI_S$ only reaches a minimum of $0.43 m^2 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 m^2 s^{-2}$, the roughly -0.6 $m^2 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 + SI_S$, 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 $SI_S$ 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 $m^2 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 $m^2 s^{-2}$ of stratification to form deep convection in December, the seasonal stratification was only low enough in February for both $δSI$ and $SI_S$ 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 $N^2 = δN^2 + N_S^2$, where $δ$ denotes the anomaly terms and “S” denotes the seasonal terms. To determine the mistral's effect,
<figure><img src="image_8.png" /><figcaption>Figure 5. Comparison 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, $\sigma_{\rho_{CTD}}$, is marked in red and is present for all plots but may be difficult to see for 7 March and 9 May.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6. The same as Fig. 5 but for potential temperature.</figcaption></figure>

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 (White, 2011), then multiply the equation by $-g/T_0$ to express the energy equation in terms of buoyancy, assuming that the ocean's density varies negatively proportionally with temperature, $\rho = -\beta T$. We then perform partial differentiation with respect to $z$ to obtain an equation describing the Brunt-Väisälä frequency, $N^2$, in response to the atmospheric forcing, given by a forcing function, $F(t)$. Separating by timescale, we arrive at the following partial differential equations:

$$ \frac{D\delta N^2}{Dt} = -\delta F(t), \qquad (10) $$

$$ \frac{DN_S^2}{Dt} = -F_S(t). \qquad (11) $$

$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 $N_S^2$ is only a function of time, $t$:

$$ \frac{\partial N_S^2}{\partial t} = -F_S(t). \qquad (12) $$
<table><caption>**Table 2.** RMSE and bias between the averaged observed CTD density and potential temperature profiles and the nearest NEMO control grid point profiles for the respective variables.</caption><thead><tr><th>Date</th><th>Lat<br>(degree)</th><th>Long<br>(degree)</th><th>RMSE<sub>ρ</sub><br>(kg m<sup>-3</sup>)</th><th>RMSE<sub>θ</sub><br>(°C)</th><th>Bias<sub>ρ</sub><br>(kg m<sup>-3</sup>)</th><th>Bias<sub>θ</sub><br>(°C)</th></tr></thead><tbody><tr><td>29 January 2013, 18:49</td><td>42.126</td><td>5.061</td><td>0.004</td><td>0.041</td><td>0.0032</td><td>-0.0265</td></tr><tr><td>30 January 2013, 11:49</td><td>42.558</td><td>5.277</td><td>0.030</td><td>0.055</td><td>-0.0093</td><td>0.0348</td></tr><tr><td>20 February 2013, 14:29</td><td>42.167</td><td>6.161</td><td>0.004</td><td>0.050</td><td>-0.0026</td><td>-0.0059</td></tr><tr><td>21 February 2013, 17:11</td><td>41.376</td><td>5.001</td><td>0.003</td><td>0.033</td><td>-0.0019</td><td>-0.0120</td></tr><tr><td>7 March 2013, 16:14</td><td>42.588</td><td>5.636</td><td>0.077</td><td>0.233</td><td>-0.0377</td><td>0.1751</td></tr><tr><td>9 March 2013, 14:29</td><td>41.645</td><td>4.227</td><td>0.005</td><td>0.043</td><td>0.0036</td><td>-0.0414</td></tr><tr><td>10 March 2013, 11:52</td><td>42.506</td><td>4.990</td><td>0.058</td><td>0.224</td><td>-0.0320</td><td>0.1719</td></tr><tr><td>9 May 2013, 14:54</td><td>42.017</td><td>4.727</td><td>0.018</td><td>0.077</td><td>0.0046</td><td>-0.0367</td></tr></tbody></table>

The average RMSE and bias for the density profiles was 0.025 and -0.009 kg m⁻³, respectively. The average RMSE and bias for potential temperature was 0.094 and 0.032 °C, respectively.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** The 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.</figcaption></figure>

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, SI<sub>S</sub>, and seasonal atmospheric forcing:

$$ \frac{\partial SI_S}{\partial t} = -\frac{D^2}{2} F_S(t). \qquad (13) $$

To simplify the analytical solution for the anomaly timescale, we describe the advection term, $V \cdot \nabla(\delta N^2)$, as a restoring term, $R = \alpha(\delta N^2)$, which relates the overall Brunt–Väisälä frequency to its seasonal component, $N_S^2$, using the linear assumption made before. This means the restoring coefficient, $\alpha$, represents the advective operation. This results in the following differential equation:

$$ \frac{\partial \delta N^2}{\partial t} + \alpha(\delta N^2) = -\delta F(t). \qquad (14) $$

## 4.2 Seasonal solution and forcing

The solution for the seasonal timescale is relatively straightforward. As already shown, Eq. (13) relates the seasonal stratification, SI<sub>S</sub>, to the seasonal atmospheric forcing,

$F_S(t)$. We have the following definition of $F_S(t)$ from Appendix A1:

$$ F_S(t) = \frac{\partial}{\partial z} \left( \frac{\mathbf{q}_{a,S} g}{\rho c_p T_0} \right) = \frac{g}{\rho c_p T_0} \frac{\partial \mathbf{q}_{a,S}}{\partial z}, \qquad (15) $$

where $c_p$ is the specific heat capacity of water, taken as 4184 J kg⁻¹ K⁻¹; $g$ is gravity; $\rho$ is the density of water, taken as 1000 kg m⁻³; and $T_0$ is the reference temperature, taken as the average seasonal sea surface temperature of 292.4 K. This means SI<sub>S</sub> can be related to the seasonal volumetric atmospheric heat transfer, $\mathbf{q}_{a,S}$. Setting $\mathbf{q}_{a,S} = -Q_{\text{net,S}}/D$, where $Q_{\text{net,S}}$ is the seasonal net downward heat flux at the ocean surface from Eq. (5), we can calculate $\frac{\partial SI_S}{\partial t}$ from $Q_{\text{net,S}}$. If we integrate both sides of Eq. (13) using $z$, after plugging in Eq. (15) and the relationship for $Q_{\text{net,S}}$, as SI<sub>S</sub> is constant with respect to (w.r.t.) $z$, Eq. (13) becomes

$$ \frac{\partial SI_S}{\partial t} = \frac{g}{2\rho c_p T_0} Q_{\text{net,S}}, \qquad (16) $$

where $\frac{g}{2\rho c_p T_0} \approx 10^{-9} \text{ m}^4 \text{ J s}^{-2}$, which means the derivative of SI<sub>S</sub> w.r.t. time, $t$, multiplied by $10^9$ is on the same order of magnitude as $Q_{\text{net}}$ (with the subscript S now dropped
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The 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, SI<sub>S</sub> + δSI, and the seasonal run, SI<sub>S</sub>. 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 t<sub>A</sub> through t<sub>D</sub> 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 (e): one calculated by a vertical change in density less than 0.01 kg m<sup>-3</sup>, denoted by Δρ, and one calculated by a vertical diffusivity less than 5 × 10<sup>-4</sup> m<sup>2</sup> s<sup>-1</sup>, denoted by K<sub>Z</sub>. 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.</figcaption></figure>

for convenience, as the rest of the subsection discusses sea-
sonal heat fluxes), which is what we see in Fig. 9a for the
2013 winter, with $\frac{\partial SI_s}{\partial t} \times 10^9$ following the curve of $Q_{net}$.
This relationship means when $Q_{net}$ crosses zero with a neg-
ative derivative, SI_s experiences a maximum and vice versa
for a minimum. Additionally, the longer $Q_{net}$ remains neg-
ative, the more seasonal destratification is incurred by the
ocean. The seasonal variation of $Q_{net}$ is primarily driven by
the solar radiation, $Q_{SW}$, which is evident in Fig. 9b. Con-
sequently, the maximum and minimum values for SI_s oc-
cur around 21 September and 21 March, the fall and spring
equinoxes. The asymmetry in $Q_{net}$ is mostly caused by the
slightly seasonally varying latent heat flux, $Q_E$, followed by
the sensible heat flux, $Q_H$, both of which also decrease the
net heat flux by roughly 100 to 200 W m$^{-2}$, depending on the
time of the year. $Q_{LW}$ remains roughly constant during the
year, decreasing $Q_{net}$ by roughly -100 W m$^{-2}$. These results

are corroborated by the results of multiple model reanalysis
for the region as well (Song and Yu, 2017).

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.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The 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, SI<sub>S</sub>, and its derivative, $\frac{\partial SI_S}{\partial t}$, comparing it to the seasonal net heat flux, Q<sub>net</sub> (the subscript S is dropped for convenience). Panel (b) shows the net heat flux separated into its components: Q<sub>E</sub>, Q<sub>H</sub>, Q<sub>SW</sub>, and Q<sub>LW</sub> 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.</figcaption></figure>

## 4.3 Anomaly solution and forcing

To solve for the anomaly timescale, described by Eq. (14), we assume $\delta 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, $\Delta t_k$, and a period between the start of the current and following event, $\Delta \tau_k$. $\delta F_k$ 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 $N^2$. The last assumption is that the restoring coefficient remains constant for each section of the forcing function:

$$ \delta \mathrm{SI}_k(t) = \begin{cases} \left[ \delta \mathrm{SI}_{k-1}(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d(t-t_k)}) \right] e^{-\alpha_d(t-t_k)} & \\ [t_k, t_k + \Delta t_k) & \\ \left[ \delta \mathrm{SI}_{k-1}(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] e^{(\alpha_a - \alpha_d) \Delta t_k - \alpha_a(t-t_k)} & \\ [t_k + \Delta t_k, t_k + \Delta \tau_k) & \end{cases}, \quad (17) $$

where $\alpha_d$ and $\alpha_a$ are the restoring coefficients during ($[t_k, t_k + \Delta t_k)$) and after ($[t_k + \Delta t_k, t_k + \Delta \tau_k)$) a mistral event, respectively.

Further assuming $\delta F_k = \delta F$, $\Delta t_k = \Delta t$, and $\Delta \tau_k = \Delta \tau$ for all $k$, which results in a periodic pulse function with constant amplitude and period, we can simplify Eq. (17) using the sum

<figure><img src="image_13.png" /><figcaption>**Figure 10.** The mistral forcing as a pulse function used to solve Eq. (14). $k$ corresponds to the event, and $\delta F_k$ corresponds to the forcing strength of the mistral event. $\Delta t_k$ corresponds to the duration of the event of the Mistral event, and $\Delta \tau_k$ corresponds to the period between events, with $t_k$ denoting the start of event $k$.</figcaption></figure>

of a finite geometric series. At the beginning of the preconditioning period, destratification has not yet begun, and therefore the initial $\delta SI$ is zero, resulting in the following equation set:

$$ \delta \mathrm{SI}_k(t) = \begin{cases} \frac{D^2}{2} \frac{\delta F}{\alpha_d} \left[ (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{\left[ (\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau \right] k}}{1 - e^{\left[ (\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau \right]}} - 1 \right) + (1 - e^{\alpha_d (t-t_k)}) \right] e^{-\alpha_d (t-t_k)} & \\ [t_k, t_k + \Delta t_k) & \\ \frac{D^2}{2} \frac{\delta F}{\alpha_d} \left[ (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{\left[ (\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau \right] k}}{1 - e^{\left[ (\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau \right]}} - 1 \right) + (1 - e^{\alpha_d \Delta t_k}) \right] e^{(\alpha_a - \alpha_d) \Delta t_k - \alpha_a (t-t_k)} & \\ [t_k + \Delta t_k, t_k + \Delta \tau_k) & \end{cases}, \quad (18) $$

This final equation set allows us to describe the integrated effect of consecutive mistral events 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 $\Delta t = 5.69$ d and $\Delta \tau = 10.88$ d. The result of the fitting is shown in Fig. 11, with $\alpha_d$ having a fitted value of $0.235$ d$^{-1}$ and $\alpha_a$ having a fitted value of $0.021$ d$^{-1}$. If we recall the meaning of $\alpha_d$ and $\alpha_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.

The strength of each Mistral event, $\delta F_k$, was found in a similar way by solving for $\delta F_k$ after noting that the initial value of $\delta SI_k(t_k)$ is equal to $\delta SI_{k-1}(t_k)$ (derivation found in Appendix A3). Following this, the values of $\delta SI$ from the NEMO results in Fig. 8 were plugged in to determine the values of $\delta F_k$ (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:

$$ \frac{\partial \delta SI_k(t)}{\partial t} = -\alpha_d \left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} \right] e^{-\alpha_d(t-t_k)}. \quad (19) $$

The quantity $\frac{\partial \delta SI_k(t)}{\partial t}$ must be less than zero for destratification to occur, which means if $\alpha_d$ is a positive quantity (refer to Appendix A2 and Fig. 11), $\delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d}$ must be a positive quantity. If some destratification has already occurred relative to the seasonal stratification, such that $\delta SI_{k-1}(t_k) < 0$, then $\frac{D^2}{2} \frac{\delta F_k}{\alpha_d}$ must be larger than $-\delta SI_{k-1}(t_k)$ for destratification to occur. Recalling that $\delta F_k$ 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$ m$^2$ s$^{-2}$ and ends at $-0.41$ m$^2$ s$^{-2}$, a net restratification of $0.07$ m$^2$ s$^{-2}$. This is despite the fact that this event has a positive $\delta F_k$ value of $3.80 \times 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 $\sim -0.6$ m$^2$ s$^{-2}$ for $\delta 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, $SI_S$, and anomaly destratification, $\delta 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 $\delta F$, $\Delta t$, and $\Delta \tau$ 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 $\delta F$, result in more destratification, with the reverse happening with decreased values. Decreasing the event duration, $\Delta 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, $\Delta \tau$) only minimally changes the accrued destratification, due to the fact that the magnitude of $\frac{\partial \delta SI}{\partial t}$ 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.

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

$$ d\delta SI = \underbrace{\frac{\partial \delta SI}{\partial \delta F}}_{\text{Strength}} d\delta F + \underbrace{\frac{\partial \delta SI}{\partial \Delta t}}_{\text{Duration}} d\Delta t + \underbrace{\frac{\partial \delta SI}{\partial \Delta \tau}}_{\text{Period}} d\Delta \tau, \quad (20) $$

Due to the lack of available total derivatives for $\delta F$, $\Delta t$, and $\Delta \tau$, we approximate them with their respective standard de-
<figure><img src="image_14.png" /><figcaption>Figure 11. The 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⁻¹ for α₉ and a value of 0.021 d⁻¹ for αₐ 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⁻⁸ s⁻² d⁻¹, was taken from values found in Table A1 from Appendix A.</figcaption></figure>

viation: σₓ ≈ 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:

$$ \frac{\partial \delta \text{SI}_k(t = t_k + \Delta t)}{\partial k} = \frac{D^2}{2} \frac{\delta F}{\alpha_d} \frac{(e^{-\alpha_d \Delta t} - 1)}{(1 - e^{(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau})} \cdot (-e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]k}) ((\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau) , \quad (21) $$

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⁻¹¹ m² s⁻² per event. This confirms the small change in the potential destratification with increasing events. Taking k to infinity and noting that α₉ > αₐ results in the following equation:

$$ \begin{aligned} \delta \text{SI}_{\infty} &= \delta \text{SI}_{\infty}(t = t_k + \Delta t) \\ &= \frac{D^2}{2} \frac{\delta F}{\alpha_d} (e^{-\alpha_d \Delta t} - 1) \left( \frac{1}{1 - e^{(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau}} \right), \quad (22) \end{aligned} $$

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 mistral, ∂δSI_∞/∂δF, is equal to a value of −1.07 × 10⁷ m² d; the derivative w.r.t. the duration, ∂δSI_∞/∂Δt, is equal to −7.60 × 10⁻³ m² s⁻² d⁻¹; and the derivative w.r.t. the period, ∂δSI_∞/∂Δτ, is equal to 2.77 × 10⁻³ m² s⁻² d⁻¹ (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⁻¹ m² s⁻², the duration term has a value of −3.21 × 10⁻² m² s⁻², and the period term has a value of 1.27 × 10⁻² m² s⁻². 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><img src="image_15.png" /></figure>

**Figure 12.** Equation (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, $\delta F$, (b) varies the duration, $\Delta t$, and (c) varies the period between events, $\Delta\tau$. Panel (d) varies the restoration coefficient during the destratification phase, $\alpha_d$, and (e) varies the restoration coefficient for the restratification phase. Panel (f) varies the number of events.

## 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
<figure><img src="image_16.png" /><figcaption>Figure 13. The same as Fig. 12, but SI<sub>S</sub> is added to the results from Eq. (18).</figcaption></figure>

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, restora-
tion data, and start dates of the simulations. For the additional
years, the NEMO simulations were initialized with and re-
stored 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 arbi-
<figure><img src="image_17.png" /></figure>

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

trary 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 restratifies 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 \text{ m}^2 \text{ s}^{-2}$), even though it is more destratified than the 1994 winter ($0.36 \text{ m}^2 \text{ 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 $\delta$SI is more negative for the 1994 winter, at $-0.55 \text{ m}^2 \text{ 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 \text{ m}^2 \text{ 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 \text{ m}^2 \text{ s}^{-2}$ for the 2013 winter and 1.73 and $1.79 \text{ m}^2 \text{ 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 $\delta$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 \text{ m}^2 \text{ 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 $\delta$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.

## 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⁻² vs. under 200 W m⁻² 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.

## 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, $\alpha_d$ and $\alpha_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 violating 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, $\delta F$, duration, $\Delta t$, and period, $\Delta \tau$) 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 \times 10^{-2}$, $-6.39 \times 10^{-2}$, and $5.83 \times 10^{-2} \text{ m}^2 \text{ 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 \times 10^{-1}$, $-5.23 \times 10^{-2}$, and $4.89 \times 10^{-2} \text{ m}^2 \text{ 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 WR-F/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.
<figure><img src="image_18.png" /></figure>

**Figure 15.** The 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><img src="image_19.png" /></figure>

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

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 con-
<figure><img src="image_20.png" /><figcaption>**Figure 17.** The 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 $10^9 \text{ m}^4 \text{ J s}^{-2}$. A negative value means heat is leaving the ocean.</figcaption></figure>

vection. 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 m$^2$ s$^{-2}$. However, all three winters had a maximal stratification within 0.1 m$^2$ 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 m$^2$ s$^{-2}$ occurs for the 2012–2013 winter. By December, the SI has already dropped to ~1 m$^2$ 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 forc-
ing. 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 (White, 2011):

$$ \rho c_p \frac{DT}{Dt} = \frac{Dq}{Dt}, \qquad (A1) $$

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

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

$$ \frac{Dq}{Dt} = -q_a, \qquad (A2) $$

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

$$ \frac{DT}{Dt} = -\frac{q_a}{\rho c_p}. \qquad (A3) $$

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

$$ N^2 = \frac{\partial b}{\partial z} = -\frac{g}{\rho_0} \frac{\partial \rho}{\partial z}. \qquad (A4) $$

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 N² in terms of temperature as follows:

$$ N^2 = \frac{\partial b}{\partial z} = -\frac{g}{T_0} \frac{\partial T}{\partial z}. \qquad (A5) $$

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

$$ \frac{g}{T_0} \left[ \frac{DT}{Dt} - \frac{q_a}{\rho c_p} \right] \Rightarrow \frac{Db}{Dt} = - \frac{q_a g}{\rho c_p T_0}. \qquad (A6) $$

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

$$ \frac{\partial}{\partial z} \left( \frac{Db}{Dt} - \frac{q_a g}{\rho c_p T_0} \right) \Rightarrow \frac{DN^2}{Dt} = - \frac{\partial}{\partial z} \left( \frac{q_a g}{\rho c_p T_0} \right), \qquad (A7) $$

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:

$$ F(t) = \frac{\partial}{\partial z} \left( \frac{q_a g}{\rho c_p T_0} \right). \qquad (A8) $$

This brings us to

$$ \frac{DN^2}{Dt} = -F(t). \qquad (A9) $$

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:

$$ N^2 = N_S^2 + \delta N^2, \\ F(t) = \delta F(t) + F_S(t), \qquad (A10) $$

which describes the response of N² on the anomaly timescale, δN²:

$$ \frac{D\delta N^2}{Dt} = -\delta F(t), \qquad (A11) $$

and the response of N² on the seasonal timescale, N_S²:

$$ \frac{DN_S^2}{Dt} = -F_S(t). \qquad (A12) $$

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

$$ \frac{\partial N_S^2}{\partial t} = -F_S(t). \qquad (A13) $$

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

$$R = \frac{\partial}{\partial z} \left( \frac{g}{T_0} \alpha \rho c_p (T - T_S) \right) \Rightarrow \alpha(N^2 - N_S^2), \quad (A14)$$

or with $\delta N^2 = N^2 - N_S^2$ as follows:

$$R = \alpha \delta N^2, \quad (A15)$$

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

$$\frac{\partial \delta N^2}{\partial t} + \mathbf{V} \cdot \nabla (\delta N^2) = -\delta F(t). \quad (A16)$$

We will replace the advective component, $\mathbf{V} \cdot \nabla(\delta N^2)$, with R, which essentially swallows the advective operation into the restoring coefficient, $\alpha$. This results in the partial differential equation that we will study further:

$$\frac{\partial \delta N^2}{\partial t} + \alpha \delta N^2 = -\delta F(t). \quad (A17)$$

## A1.1 Solution for seasonal SI

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

$$F_S(t) = \frac{g}{\rho c_p T_0} \frac{\partial \mathbf{q}_{a,S}}{\partial z}. \quad (A18)$$

Dividing $\mathbf{q}_{a,S}$ by D gives us the atmospheric cooling in terms of a surface flux, $-Q_{net,S}$. If we plug this relationship back into Eq. (13), we get

$$\frac{\partial \mathrm{SI}_S}{\partial t} = \frac{D}{2} \frac{g}{\rho c_p T_0} \frac{\partial Q_{\mathrm{net},S}}{\partial z}, \quad (A19)$$

Integrating this equation by z gives us

$$\int_{0}^{D} \frac{\partial \mathrm{SI}_S}{\partial t} \partial z = \frac{D}{2} \frac{g}{\rho c_p T_0} \int_{0}^{D} \frac{\partial Q_{\mathrm{net},S}}{\partial z} \partial z \\ \frac{\partial \mathrm{SI}_S}{\partial t} D = \frac{D}{2} \frac{g}{\rho c_p T_0} Q_{\mathrm{net},S} \\ \frac{\partial \mathrm{SI}_S}{\partial t} = \frac{g}{2 \rho c_p T_0} Q_{\mathrm{net},S}, \quad (A20)$$

Therefore, we have SI$_S$ expressed in terms of $Q_{net,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, $\delta F(t)$, as a series of k pulse functions of magnitude $\delta F_k$ over a duration of $\Delta t_k$ and with a period of $\Delta \tau_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 ($\alpha$) remain constant with respect to time. This leads to $\alpha_d$ and $\alpha_a$ representing the advective components during and after an event, respectively. During a Mistral event, [$t_k, t_k + \Delta t_k$), we get

$$\frac{\partial \delta N^2}{\partial t} + \alpha_d(\delta N^2) = -\delta F(t) \\ \delta N_k^2(t) = -\frac{\delta F_k}{\alpha_d} + c_0 e^{-\alpha_d t} \\ \delta N_k^2(t_k) = -\frac{\delta F_k}{\alpha_d} + c_0 e^{-\alpha_d t_k} = \delta N_{k-1}^2(t_k) \\ c_0 = \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} \right] e^{\alpha_d t_k}, \quad (A21)$$

with the following result:

$$\delta N_k^2(t) = \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{-\alpha_d(t-t_k)}) \right] \cdot e^{-\alpha_d(t-t_k)}, \quad (A22)$$

After the event, [$t_k + \Delta t_k, t_{k+1}$):

$$\frac{\partial \delta N^2}{\partial t} + \alpha_a(\delta N^2) = 0 \\ \delta N_k^2(t) = c_1 e^{-\alpha_a t} \\ \delta N_k^2(t_k + \Delta t_k) = c_1 e^{-\alpha_a(t_k + \Delta t_k)} = \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{-\alpha_d \Delta t_k}) \right] e^{-\alpha_d \Delta t_k} \\ c_1 = \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{-\alpha_d \Delta t_k}) \right] \cdot e^{(\alpha_a - \alpha_d) \Delta t_k + \alpha_a t_k}, \quad (A23)$$

With the following result:

$$\delta N_k^2(t) = \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{-\alpha_d \Delta t_k}) \right] \cdot e^{(\alpha_a - \alpha_d) \Delta t_k - \alpha_a(t - t_k)}, \quad (A24)$$
or to have the results put more succinctly:

$$ \delta N_k^2(t) = \begin{cases} \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d(t-t_k)}) \right] \\ \cdot e^{-\alpha_d(t-t_k)} \\ [t_k, t_k + \Delta t_k] \\ \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \\ \cdot e^{(\alpha_a - \alpha_d) \Delta t_k - \alpha_a(t-t_k)} \\ [t_k + \Delta t_k, t_k + \Delta \tau_k]. \end{cases} \quad (A25) $$

### A1.3 $\delta N_{k-1}^2$ initial condition

$\delta N_{k-1}^2(t_k)$ is a recursive initial condition, as its initial condition is the event before it (and so on):

$$ \begin{aligned} \delta N_{k-1}^2(t_k) &= \left[ \delta N_{k-2}^2(t_{k-1}) + \frac{\delta F_{k-1}}{\alpha_d} (1 - e^{\alpha_d \Delta t_{k-1}}) \right] \\ & \quad \cdot e^{(\alpha_a - \alpha_d) \Delta t_{k-1} - \alpha_a \Delta \tau_{k-1}} \\ \delta N_{k-2}^2(t_{k-1}) &= \left[ \delta N_{k-3}^2(t_{k-2}) + \frac{\delta F_{k-2}}{\alpha_d} (1 - e^{\alpha_d \Delta t_{k-2}}) \right] \\ & \quad \cdot e^{(\alpha_a - \alpha_d) \Delta t_{k-2} - \alpha_a \Delta \tau_{k-2}}. \end{aligned} \quad (A26) $$

Therefore, $\delta N_{k-1}^2(t_k)$ can be simplified in expression by combining the initial conditions:

$$ \begin{aligned} \delta N_{k-1}^2(t_k) &= \delta N_{k-m}^2(t_{k-(m-1)}) e^{(\alpha_a - \alpha_d) \sum_{i=1}^{m-1} \Delta t_{k-i}} \\ & \quad \cdot e^{-\alpha_a \sum_{i=1}^{m-1} \Delta \tau_{k-i}} + \sum_{j=1}^{m-1} \frac{\delta F_{k_j}}{\alpha_d} (1 - e^{\alpha_d \Delta t_{k-j}}) \\ & \quad \cdot e^{(\alpha_a - \alpha_d) \sum_{i=1}^{j} \Delta t_{k-i}} e^{-\alpha_a \sum_{i=1}^{j} \Delta \tau_{k-i}}. \end{aligned} \quad (A27) $$

If $m=k$ and $\delta N_0^2 = 0$ then

$$ \begin{aligned} \delta N_{k-1}^2(t_k) &= \sum_{j=1}^{k-1} \frac{\delta F_{k_j}}{\alpha_d} (1 - e^{\alpha_d \Delta t_{k-j}}) \\ & \quad \cdot e^{(\alpha_a - \alpha_d) \sum_{i=1}^{j} \Delta t_{k-i}} e^{-\alpha_a \sum_{i=1}^{j} \Delta \tau_{k-i}}. \end{aligned} \quad (A28) $$

Assuming $\delta F_k = \delta F$, $\Delta t_k = \Delta t$, and $\Delta \tau_k = \Delta \tau$ for all $k$, or a periodic pulse function, then $\delta N_{k-1}^2$ can be expressed as follows:

$$ \delta N_{k-1}^2(t_k) = \frac{\delta F}{\alpha_d} (1 - e^{\alpha_d \Delta t}) \sum_{j=1}^{k-1} e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] j}. \quad (A29) $$

Taking the sum of a finite geometric series as follows:

$$ \begin{aligned} \sum_{n=0}^{m-1} r^n &= \left( \frac{1-r^m}{1-r} \right), \\ \sum_{n=0}^{m-1} r^n &= \sum_{n=1}^{m-1} r^n + 1, \\ \sum_{n=1}^{m-1} r^n &= \sum_{n=0}^{m-1} r^n - 1, \end{aligned} \quad (A30) $$

where $r \neq 1$. If $r = e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau}$, then

$$ \sum_{j=1}^{k-1} e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] j} = \left( \frac{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] k}}{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]}} - 1 \right), \quad (A31) $$

and we then get

$$ \delta N_{k-1}^2(t_k) = \frac{\delta F}{\alpha_d} (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] k}}{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]}} - 1 \right), \quad (A32) $$

where $(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau \neq 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:

$$ \delta N^2(t) = \begin{cases} \frac{\delta F}{\alpha_d} \left[ (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] k}}{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]}} - 1 \right) + (1 - e^{\alpha_d (t - t_k)}) \right] e^{-\alpha_d (t - t_k)} \\ [t_k, t_k + \Delta t) \\ \frac{\delta F}{\alpha_d} \left[ (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] k}}{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]}} - 1 \right) + (1 - e^{\alpha_d \Delta t}) \right] e^{(\alpha_a - \alpha_d) \Delta t - \alpha_a (t - t_k)} \\ [t_k + \Delta t, t_k + \Delta \tau). \end{cases} \quad (A33) $$

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

$$ \delta SI_k(t) = \begin{cases} \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d(t-t_k)}) \right] \\ \cdot e^{-\alpha_d(t-t_k)} \\ [t_k, t_k + \Delta t_k) \\ \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \\ \cdot e^{(\alpha_a - \alpha_d) \Delta t_k - \alpha_a(t-t_k)} \\ [t_k + \Delta t_k, t_k + \Delta \tau_k), \end{cases} \quad (A34) $$

and

$$ \delta SI_k(t) = \begin{cases} \frac{D^2}{2} \frac{\delta F}{\alpha_d} \left[ (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] k}}{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]}} - 1 \right) + (1 - e^{\alpha_d (t - t_k)}) \right] e^{-\alpha_d (t - t_k)} \\ [t_k, t_k + \Delta t) \\ \frac{D^2}{2} \frac{\delta F}{\alpha_d} \left[ (1 - e^{\alpha_d \Delta t}) \left( \frac{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau] k}}{1 - e^{[(\alpha_a - \alpha_d) \Delta t - \alpha_a \Delta \tau]}} - 1 \right) + (1 - e^{\alpha_d \Delta t}) \right] e^{(\alpha_a - \alpha_d) \Delta t - \alpha_a (t - t_k)} \\ [t_k + \Delta t, t_k + \Delta \tau), \end{cases} \quad (A35) $$

for the period pulse function case.

### A2 Restoring coefficients $\alpha_d$ and $\alpha_a$

The restoration coefficients, $\alpha_d$ and $\alpha_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 stratification and restratification cases that are highlighted in Table 1 and Fig. 8 to retrieve the values of the restoration coefficients.
## A2.1 Restoration coefficient $\alpha_d$ during a mistral

To solve for $\alpha_d$, we normalize Eq. (A34) for during a mistral event, $[t_k, t_k + \Delta t_k)$, with $\delta SI$ given as follows:

$$ \delta SI_k(t) = \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d(t-t_k)}) \right] \cdot e^{-\alpha_d(t-t_k)}, \quad (A36) $$

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

$$ \delta SI_k(t') = \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha} (1 - e^{\alpha_d t'}) \right] e^{-\alpha_d t'}, \quad (A37) $$

Next, we normalize $\delta SI_k$ to the value of zero at $t' = 0$, resulting in $\delta SI_{k,NI}$:

$$ \begin{aligned} \delta SI_{k,NI}(t') &= \delta SI_k(t') - \delta SI_k(t' = 0) \\ &= \left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d t'}) \right] \cdot e^{-\alpha_d t'} - \delta SI_{k-1}(t_k) \\ &= \left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} \right] (e^{-\alpha_d t'} - 1), \quad (A38) \end{aligned} $$

Following this, the height or magnitude of destratification for each event is normalized to 1, resulting in $\delta SI_{k,NH}$:

$$ \delta SI_{k,NH}(t') = \frac{\delta SI_{k,NI}(t')}{\text{extremum}(\delta SI_{k,NI}(t'))}. \quad (A39) $$

The extremum value for $\delta SI_{k,NI}(t')$ is when $t' = \Delta t_k$ or at the end of the mistral event. This simplifies $\delta SI_{k,NH}$ to

$$ \begin{aligned} \delta SI_{k,NH}(t') &= \frac{\left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} \right] (e^{-\alpha_d t'} - 1)}{\left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} \right] (e^{-\alpha_d \Delta t_k} - 1)} \\ &= \frac{(e^{-\alpha_d t'} - 1)}{(e^{-\alpha_d \Delta t_k} - 1)}, \quad (A40) \end{aligned} $$

To normalize the length of the event duration from 0 to 1, we then divide $t'$ by the event length, $\Delta t_k$, which results in $t''$:

$$ t'' = \frac{t'}{\Delta t_k} \Rightarrow t' = t'' \Delta t_k. \quad (A41) $$

Plugging $t''$ into $\delta SI_{k,NH}(t')$ returns $\delta SI_{k,NT}$:

$$ \delta SI_{k,NT}(t'') = \frac{e^{-\alpha_d t'' \Delta t_k} - 1}{e^{-\alpha_d \Delta t_k} - 1}. \quad (A42) $$

This final equation, $\delta SI_{k,NT}$, can be used with a fitting function to solve for $\alpha_d$ if $\Delta t_k$ is supplied.

## A2.2 Restoration coefficient $\alpha_a$, after a mistral

To solve for $\alpha_a$, we normalize Eq. (A34) for after a mistral event, $[t_k + \Delta t_k, t_k + \Delta \tau_k)$, with $\delta SI$ given as

$$ \delta SI_k(t) = \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \cdot e^{(\alpha_a - \alpha_d) \Delta t_k - \alpha_a(t - t_k)}, \quad (A43) $$

Referencing the time, $t$, to the end of the event, $t'' = t - (t_k + \Delta t_k)$ and plugging $t''$ into Eq. (A43), we get

$$ \begin{aligned} \delta SI_k(t'') &= \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \cdot e^{-\alpha_d \Delta t_k} e^{-\alpha_a t''}, \quad (A44) \end{aligned} $$

Normalizing the vertical intercept of $\delta SI_k(t'')$ results in $\delta SI_{k,NI}$:

$$ \begin{aligned} \delta SI_{k,NI} &= \delta SI_k(t'') - \delta SI_k(t'' = 0) \\ &= \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \cdot e^{-\alpha_d \Delta t_k} e^{-\alpha_a t''} \\ &\quad - \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] e^{-\alpha_d \Delta t_k} \\ &= \left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \cdot e^{-\alpha_d \Delta t_k} (e^{-\alpha_a t''} - 1), \quad (A45) \end{aligned} $$

Each post-event restratification is normalized to the height of 1 by dividing $\delta SI_{k,NI}$ by $(\delta SI_k(t'' = \Delta \tau_k - \Delta t_k) - \delta SI_k(t''' = 0))$:

$$ \begin{aligned} \delta SI_{k,NH} &= \frac{\delta SI_{k,NI}}{\delta SI_k(t'''' = \Delta \tau_k - \Delta t_k) - \delta SI_k(t''' = 0)} \\ &= \frac{\left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \cdot e^{-\alpha_d \Delta t_k} (e^{-\alpha_a t''''} - 1)}{\left[ \delta SI_{k-1}^2(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (1 - e^{\alpha_d \Delta t_k}) \right] \cdot e^{-\alpha_d \Delta t_k} (e^{-\alpha_a (\Delta \tau_k - \Delta t_k)} - 1)} \\ &= \frac{e^{-\alpha_a t''''} - 1}{e^{-\alpha_a (\Delta \tau_k - \Delta t_k)} - 1}, \quad (A46) \end{aligned} $$

which gives us $\delta SI_{k,NH}$.

To normalize the length of time of post-event restratification, we will divide $t''''$ by the post-event time length, $\Delta \tau_k - \Delta t_k$, resulting in $t''''$:

$$ t'''' = \frac{t'''}{\Delta \tau_k - \Delta t_k} \Rightarrow t''' = t'''' (\Delta \tau_k - \Delta t_k), \quad (A47) $$
<table><caption>**Table A1.** The mistral strengths, $\delta F_k$, for each of the preconditioning phase events, using $\alpha_d$ and $\alpha_a$ from Sect. A2 and the rest of the preconditioning period mistral characteristics from Table 1 plugged into Eq. (A50).</caption><thead><tr><th>Date<br>(yyyy-mm-dd) &times; 10<sup>-8</sup> s<sup>-2</sup> d<sup>-1</sup>)</th><th>&delta;F<sub>k</sub></th><th>Date<br>(yyyy-mm-dd) &times; 10<sup>-8</sup> s<sup>-2</sup> d<sup>-1</sup>)</th><th>&delta;F<sub>k</sub></th></tr></thead><tbody><tr><td>2012-08-30</td><td>1.81</td><td>2012-11-27</td><td>5.46</td></tr><tr><td>2012-09-12</td><td>2.67</td><td>2012-12-08</td><td>6.37</td></tr><tr><td>2012-09-19</td><td>2.66</td><td>2012-12-17</td><td>5.30</td></tr><tr><td>2012-09-28</td><td>2.73</td><td>2012-12-21</td><td>5.03</td></tr><tr><td>2012-10-12</td><td>3.04</td><td>2012-12-26</td><td>4.39</td></tr><tr><td>2012-10-27</td><td>4.59</td><td>2013-01-02</td><td>3.80</td></tr><tr><td>2012-11-11</td><td>4.84</td><td>2013-01-23</td><td>3.53</td></tr><tr><td>2012-11-19</td><td>3.92</td><td></td><td></td></tr></tbody></table>

$\delta F_k = 4.01 \times 10^{-8} s^{-2} d^{-1}$, and $\sigma \delta F_k = 1.196 \times 10^{-8} s^{-2} d^{-1}$.

which leads to $\delta SI_{k,NT}$:

$$ \delta SI_{k,NT} = \frac{e^{-\alpha_a t''''(\Delta\tau_k - \Delta t_k)} - 1}{e^{-\alpha_a(\Delta\tau_k - \Delta t_k)} - 1}. \quad (A48) $$

This leaves us with an equation of $\alpha_a$ that can be fitted against NEMO model data if $\Delta t_k$ and $\Delta \tau_k$ are provided.

The average duration and period of events during the preconditioning period in Table 1, $\overline{\Delta t}$ and $\overline{\Delta \tau}$, respectively, are used for $\Delta t$ and $\Delta \tau$ in these normalized equations, Eqs. (A42) and (A48). The result of the fitting is shown in Fig. 11, with $\alpha_d$ having a fitted value of 0.235 d$^{-1}$ and $\alpha_a$ having a fitted value of 0.021 d$^{-1}$.

## A3 Determining $\delta F_k$

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, $\delta F_k$, can be determined. If we take Eq. (A36) and note that the value of $\delta SI_{k-1}$ is the same as $\delta SI_k(t_k)$ at the beginning of an event, we can simplify the equation in the following steps:

$$ \begin{aligned} \delta SI_k(t_k) &= \delta SI_{k-1}, \\ \delta SI_k(t_k + \Delta t_k) &= \delta SI_{k-1}e^{-\alpha_d \Delta t_k} + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} (e^{-\alpha_d \Delta t_k} - 1), \end{aligned} \quad (A49) $$

and then solve for $\delta F_k$:

$$ \delta F_k = \frac{2(\delta SI_k(t_k + \Delta t_k) - \delta SI_k(t_k)e^{-\alpha_d \Delta t_k})\alpha_d}{(e^{-\alpha_d \Delta t_k} - 1)D^2}. \quad (A50) $$

The results of $\delta F_k$ for each event in the preconditioning phase are given in Table A1, along with mean and standard deviation.

## A4 Time derivative of $\delta N^2$

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

$$ \frac{\partial \delta N_k^2(t)}{\partial t} = \begin{cases} -\alpha_d \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} \right] e^{-\alpha_d(t-t_k)} \\ \qquad \left[t_k, t_k + \Delta t_k\right) \\ -\alpha_a \left[ \delta N_{k-1}^2(t_k) + \frac{\delta F_k}{\alpha_d} \left(1 - e^{\alpha_d \Delta t_k}\right) \right] \\ \qquad \cdot e^{(\alpha_a - \alpha_d)\Delta t_k - \alpha_a(t-t_k)} \\ \qquad \left[t_k + \Delta t_k, t_k + \Delta \tau_k\right), \end{cases} \quad (A51) $$

and

$$ \frac{\partial \delta SI_k(t)}{\partial t} = \begin{cases} -\alpha_d \left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} \right] e^{-\alpha_d(t-t_k)} \\ \qquad \left[t_k, t_k + \Delta t_k\right) \\ -\alpha_a \left[ \delta SI_{k-1}(t_k) + \frac{D^2}{2} \frac{\delta F_k}{\alpha_d} \left(1 - e^{\alpha_d \Delta t_k}\right) \right] \\ \qquad \cdot e^{(\alpha_a - \alpha_d)\Delta t_k - \alpha_a(t-t_k)} \\ \qquad \left[t_k + \Delta t_k, t_k + \Delta \tau_k\right). \end{cases} \quad (A52) $$

## A5 Asymptotic destratification

The following sections differentiate Eq. (A35) at $t = t_k + \Delta t_k$ or at the end of a mistral event, where the destratification is the largest, using $k$ or using the other components, i.e., $\delta F$, $\Delta t$, and $\Delta \tau$, once $k \to \infty$.

$$ A5.1 \quad \frac{\partial \delta SI_k}{\partial k} $$

Equation (A35), at $t = t_k + \Delta t$ results in

$$ \delta SI_k(t_k + \Delta t) = \frac{D^2}{2} \frac{\delta F}{\alpha_d} (e^{-\alpha_d \Delta t} - 1) \cdot \left( \frac{1 - e^{[(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau]k}}{1 - e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau}} \right), \quad (A53) $$

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

$$ \frac{\partial \delta SI_k}{\partial k} = \frac{D^2}{2} \frac{\delta F}{\alpha_d} \frac{(e^{-\alpha_d \Delta t} - 1)}{(1 - e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau})} \cdot \left( -e^{\left[(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau\right]k} \right) \left[ (\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau \right], \quad (A54) $$

As $k \to \infty$, with $\alpha_d > \alpha_a$ Eq. (A53) goes to

$$ \delta SI_{\infty} = \frac{D^2}{2} \frac{\delta F}{\alpha_d} (e^{-\alpha_d \Delta t} - 1) \left( \frac{1}{1 - e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau}} \right). \quad (A55) $$

$$ A5.2 \quad \frac{\partial \delta SI_{\infty}}{\partial \delta F} $$

The derivative of Eq. (A55) w.r.t. $\delta F$ is

$$ \frac{\partial \delta SI_{\infty}}{\partial \delta F} = \frac{D^2}{2} \frac{1}{\alpha_d} (e^{-\alpha_d \Delta t} - 1) \left( \frac{1}{1 - e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau}} \right). \quad (A56) $$
A5.3 $\frac{\partial \delta SI_{\infty}}{\partial \Delta t}$

The derivative of Eq. (A55) w.r.t. $\Delta t$ is

$$ \frac{\partial \delta SI_{\infty}}{\partial \Delta t} = \frac{D^2}{2} \frac{\delta F}{\alpha_d} \frac{1}{(1 - e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau})} \cdot \left[ -\alpha_d e^{-\alpha_d \Delta t} + (e^{-\alpha_d \Delta t} - 1) \cdot \left( \frac{(\alpha_a - \alpha_d)e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau}}{1 - e^{(\alpha_a - \alpha_d)\Delta t - \alpha_a \Delta \tau}} \right) \right]. \quad (A57) $$

A5.4 $\frac{\partial \delta SI_{\infty}}{\partial \Delta \tau}$

The derivative of Eq. (A55) w.r.t. $\Delta \tau$ is

$$ \frac{\partial \delta SI_{\infty}}{\partial \Delta \tau} = \frac{D^2}{2} \frac{\delta F}{\alpha_d} (e^{-\alpha_d \Delta t} - 1) \frac{(-\alpha_a e^{(\alpha_a - \alpha_d) \Delta t} - \alpha_a \Delta \tau)}{(1 - e^{(\alpha_a - \alpha_d) \Delta t} - \alpha_a \Delta \tau)^2}. \quad (A58) $$

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><caption>**Table B1.** The start date, duration ($\Delta t_k$, in d), and period between each event ($\Delta \tau_k$, in d) for each Mistral event ($k$) and event strength, ($\delta F_k$, in $s^{-2} d^{-1}$) for the preconditioning phase of the NEMO simulation period of 1 June 1993 to 31 May 1994.</caption><thead><tr><th>Date<br>(yyyy-mm-dd)</th><th>Δt<sub>k</sub></th><th>Δτ<sub>k</sub></th><th>δF<sub>k</sub></th><th>Date<br>(yyyy-mm-dd)</th><th>Δt<sub>k</sub></th><th>Δτ<sub>k</sub></th><th>δF<sub>k</sub></th></tr></thead><tbody><tr><td>1993-08-28</td><td>4</td><td>7</td><td>9.23 × 10<sup>-9</sup></td><td>1993-12-05</td><td>1</td><td>21</td><td>2.57 × 10<sup>-8</sup></td></tr><tr><td>1993-09-04</td><td>2</td><td>21</td><td>8.39 × 10<sup>-9</sup></td><td>1993-12-26</td><td>3</td><td>13</td><td>2.52 × 10<sup>-8</sup></td></tr><tr><td>1993-09-25</td><td>4</td><td>25</td><td>1.94 × 10<sup>-8</sup></td><td>1994-01-08</td><td>1</td><td>3</td><td>2.17 × 10<sup>-8</sup></td></tr><tr><td>1993-10-20</td><td>6</td><td>16</td><td>2.92 × 10<sup>-8</sup></td><td>1994-01-11</td><td>2</td><td>6</td><td>2.07 × 10<sup>-8</sup></td></tr><tr><td>1993-11-05</td><td>1</td><td>3</td><td>2.96 × 10<sup>-8</sup></td><td>1994-01-17</td><td>15</td><td>19</td><td>2.91 × 10<sup>-8</sup></td></tr><tr><td>1993-11-08</td><td>2</td><td>4</td><td>2.54 × 10<sup>-8</sup></td><td>1994-02-05</td><td>7</td><td>14</td><td>3.62 × 10<sup>-8</sup></td></tr><tr><td>1993-11-12</td><td>5</td><td>9</td><td>2.42 × 10<sup>-8</sup></td><td>1994-02-19</td><td>1</td><td>15</td><td>4.19 × 10<sup>-8</sup></td></tr><tr><td>1993-11-21</td><td>2</td><td>5</td><td>2.42 × 10<sup>-8</sup></td><td>1994-03-06</td><td>3</td><td>7</td><td>3.88 × 10<sup>-8</sup></td></tr><tr><td>1993-11-26</td><td>7</td><td>9</td><td>2.57 × 10<sup>-8</sup></td><td>1994-03-13</td><td>6</td><td>7</td><td>3.74 × 10<sup>-8</sup></td></tr><tr><td>1993-12-05</td><td>1</td><td>21</td><td>2.57 × 10<sup>-8</sup></td><td>1994-03-20</td><td>3</td><td>6</td><td>3.72 × 10<sup>-8</sup></td></tr></tbody></table>

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

<table><caption>**Table B2.** The same as Table B1 but for the NEMO simulation period from 1 June 2004 to 31 May 2005.</caption><thead><tr><th>Date<br>(yyyy-mm-dd)</th><th>Δt<sub>k</sub></th><th>Δτ<sub>k</sub></th><th>δF<sub>k</sub></th><th>Date<br>(yyyy-mm-dd)</th><th>Δt<sub>k</sub></th><th>Δτ<sub>k</sub></th><th>δF<sub>k</sub></th></tr></thead><tbody><tr><td>2004-09-15</td><td>2</td><td>4</td><td>-2.61 × 10<sup>-10</sup></td><td>2004-11-23</td><td>2</td><td>4</td><td>1.93 × 10<sup>-8</sup></td></tr><tr><td>2004-09-19</td><td>10</td><td>37</td><td>6.93 × 10<sup>-9</sup></td><td>2004-11-27</td><td>1</td><td>7</td><td>1.88 × 10<sup>-8</sup></td></tr><tr><td>2004-10-26</td><td>1</td><td>6</td><td>7.86 × 10<sup>-9</sup></td><td>2004-12-04</td><td>7</td><td>14</td><td>2.30 × 10<sup>-8</sup></td></tr><tr><td>2004-11-01</td><td>2</td><td>4</td><td>6.13 × 10<sup>-9</sup></td><td>2004-12-18</td><td>17</td><td>31</td><td>2.30 × 10<sup>-8</sup></td></tr><tr><td>2004-11-05</td><td>16</td><td>18</td><td>2.06 × 10<sup>-8</sup></td><td>2005-01-18</td><td>17</td><td>24</td><td>4.10 × 10<sup>-8</sup></td></tr><tr><td>2004-11-23</td><td>2</td><td>4</td><td>1.93 × 10<sup>-8</sup></td><td>2005-02-11</td><td>18</td><td>21</td><td>4.43 × 10<sup>-8</sup></td></tr></tbody></table>

The preconditioning phase for this year is considered to range from 15 September 2004 to 8 April 2005 (170 d).
**Code availability.** The RegIPSL model and the NEMO model and its configuration used in this paper can be accessed at the following Git repository: https://gitlab.in2p3.fr/ipsl/lmd/intro/regipsl/regipsl (Drobinski et al., 2012). Code used to generate figures and implement the simple model can be found at the following Git repository: https://gitlab.com/dkllrjr/untangling-deep Conv_2012_2013_code (Keller, 2022).

**Data availability.** The ERA-Interim Reanalysis can be found at https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim (Dee et al., 2011). HyMeX data are available at https://mistrals.sedoo.fr/HyMeX/ with registration (last access: 24 March 2022) (Drobinski et al., 2014; Taupier-Letage, 2013). 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.

**Disclaimer.** 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/G-MD/HESS/NHESS/OS inter-journal SI)". It is not associated with a conference.

**Acknowledgements.** 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.

## References

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, http://www.umr-cnrm.fr/spip.php?article1197&lang=fr (last access: 24 March 2022), 2012.

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, https://doi.org/10.1002/grl.50382, 2013.

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.

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, https://doi.org/10.1029/2011jc007679, 2012.

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-O₂ floats in the north-western Mediterranean Sea in 2013, J. Geophys. Res.-Oceans, 122, 6325–6341, https://doi.org/10.1002/2016jc012594, 2017.

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, https://doi.org/10.1002/qj.828, 2011 (data available at: https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim, last access: 24 March 2022).

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, https://doi.org/10.1175/bamsd-12-00242.1, 2014 (data available at: https://mistrals.sedoo.fr/HyMeX/, last access: 24 March 2022, registration required).

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., Fougols, M.-A., Khvorostyanov, D., Madec, G., Mancip, M., Mason, 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 stud-
ies in vulnerable regions, Environ. Modell. Softw., 35, 1–18,
https://doi.org/10.1016/j.envsoft.2012.01.017, 2012 (code avail-
able at: https://gitlab.in2p3.fr/ipsl/lmd/intro/regipsl/regipsl, last
access: 24 March 2022).

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, https://doi.org/10.1002/2016jd025530, 2017.

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, https://doi.org/10.5670/oceanog.2016.94, 2016a.

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, https://doi.org/10.1002/2016jc011935, 2016b.

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

Gascard, J. C.: Mediterranean deep water formation baroclinic in-
stability and oceanic eddies, Oceanol. Acta, 1, 315–330, 1978.

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, https://doi.org/10.5194/wcd-2-609-2021, 2021.

Grignon, L., Smeed, D. A., Bryden, H. L., and Schroeder, K.: Im-
portance of the variability of hydrographic preconditioning for
deep convection in the Gulf of Lion, NW Mediterranean, Ocean
Sci., 6, 573–586, https://doi.org/10.5194/os-6-573-2010, 2010.

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., https://doi.org/10.1007/s00382-021-05938-y, online first, 2021.

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, https://doi.org/10.5194/os-12-577-2016, 2016.

Herrmann, M., Sevault, F., Beuvier, J., and Somot, S.: What in-
duced the exceptional 2005 convection event in the northwestern
Mediterranean basin? Answers from a modeling study, J. Geo-
phys. Res., 115, C12051, https://doi.org/10.1029/2010jc006162,
2010.

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, https://doi.org/10.1002/2016jc011857, 2016.

Keller Jr., D.: Code to Analyze and Plot Data for Untangling the
mistral and seasonal atmospheric forcing driving deep convec-
tion in the Gulf of Lion: 2012–2013, GitLab [code], https://
gitlab.com/dkllrjr/untangling_deepoj conv_2012_2013_code, last
access: 5 April 2022.

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogee, 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, https://doi.org/10.1029/2003gb002199, 2005.

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,
https://doi.org/10.5065/D6KK98Q6, 2004.

Large, W. G. and Yeager, S. G.: The global climatology of an inter-
annually varying air-sea flux data set, Clim. Dynam., 33, 341–
364, https://doi.org/10.1007/s00382-008-0441-3, 2008.

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, https://doi.org/10.1029/2008jd011601, 2009.

Lebeaupin-Brossier, C., Béranger, K., Deltel, C., and Drobin-
ski, P.: The Mediterranean response to different space-
time resolution atmospheric forcings using perpetual
mode sensitivity simulations, Ocean Model., 36, 1–25,
https://doi.org/10.1016/j.ocemod.2010.10.008, 2011.

Lebeaupin-Brossier, C., Béranger, K., and Drobinski, P.: Sensi-
tivity of the northwestern Mediterranean Sea coastal and ther-
mohaline circulations simulated by the 1/12°-resolution ocean
model NEMO-MED12 to the spatial and temporal resolu-
tion of atmospheric forcing, Ocean Model., 43-44, 94–107,
https://doi.org/10.1016/j.ocemod.2011.12.007, 2012.

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, https://doi.org/10.1002/2016jc012526, 2017.

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, https://doi.org/10.1002/2015jc011542, 2016.

L'Hévéder, B., Li, L., Sevault, F., and Somot, S.: Interannual vari-
ability of deep convection in the Northwestern Mediterranean
simulated with a coupled AORCM, Clim. Dynam., 41, 937–960,
https://doi.org/10.1007/s00382-012-1527-5, 2012.

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 dur-
ing past and future decades?, Prog. Oceanogr., 80, 199-217,
https://doi.org/10.1016/j.pocean.2009.02.001, 2009.

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, https://doi.org/10.1016/0377-0265(91)90024-a, 1991a.

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, https://doi.org/10.1175/1520-0485(1991)021<1349:atdnso>2.0.co;2, 1991b.

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, https://doi.org/10.1175/1520-0485(1996)026<1393:LSPODW>2.0.CO;2, 1996.

Marshall, J. and Schott, F.: Open-ocean convection: Observations, theory, and models, Rev. Geophys., 37, 1–64, https://doi.org/10.1029/98rg02739, 1999.

MEDOC: Observation of Formation of Deep Water in the Mediterranean Sea, 1969, Nature, 227, 1037–1040, https://doi.org/10.1038/2271037a0, 1970.

Mertens, C. and Schott, F.: Interannual Variability of Deep-Water Formation in the Northwestern Mediterranean, J. Phys. Oceanogr., 28, 1410–1424, https://doi.org/10.1175/1520-0485(1998)028<1410:ivodwf>2.0.co;2, 1998.

Millot, C. and Taupier-Letage, I.: Circulation in the Mediterranean Sea, in: The Mediterranean Sea, edited by: Saliot, A., 29–66, Springer Berlin Heidelberg, https://doi.org/10.1007/b107143, 2005.

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, https://doi.org/10.1175/1520-0485(2003)033<1145:tropit>2.0.co;2, 2003.

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, https://doi.org/10.1006/rwos.2001.0376, 2001.

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., Drobin-ski, 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, https://doi.org/10.1175/bams-d-14-00176.1, 2016.

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, https://doi.org/10.1175/1520-0485(1996)026<0505:oodcit>2.0.co;2, 1996.

Send, U. and Testor, P.: Direct Observations Reveal the Deep Circulation of the Western Mediterranean Sea, J. Geophys. Res.-Oceans, 122, 10091–10098, https://doi.org/10.1002/2016jc012679, 2017.

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, https://doi.org/10.1002/2016jc012664, 2017.

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, https://doi.org/10.5065/D68S4MVH, 2008.

Somot, S.: Modélisation Climatique du Bassin Méditerranéen: Variabilité et Scénarios de Changement Climatique, PhD thesis, Université Toulouse III, https://tel.archives-ouvertes.fr/tel-00165252 (last access: 24 March 2022), 2005.

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, https://doi.org/10.1007/s00382-016-3295-0, 2016.

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, https://doi.org/10.1002/2016jc012254, 2017.

Taupier-Letage, I.: CTD SOP2, Provence - Tethys 2, SE-DOO OMP [data set], https://doi.org/10.6096/MISTRALS-HYMEX.950, 2013 (https://mistrals.sedoo.fr/HyMeX/, last access: 24 March 2022, registration required).

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, https://doi.org/10.1175/1520-0485(2003)033<0075:1ssodw>2.0.co;2, 2003.

The Lab Sea Group: The Labrador Sea Deep Convection Experiment, B. Am. Meteorol. Soc., 79, 2033–2058, https://doi.org/10.1175/1520-0477(1998)079<2033:tlsdce>2.0.co;2, 1998.

Turner, J. S.: Buoyancy Effects in Fluids, Cambridge University Press, https://doi.org/10.1017/CBO9780511608827, 1973.

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, https://doi.org/10.1029/2018g1078502, 2018.

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