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

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 (Beuvier 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 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, τ:

$$\begin{array}{}\text{(1)}& {\displaystyle }{Q}_E={\mathit{\rho }}_{\mathrm{0}}\mathrm{\Lambda }{C}_E(q_z-q_{\mathrm{0}})\left|\mathrm{\Delta }\mathit{u}\right|,\text{(2)}& {\displaystyle }{Q}_H={\mathit{\rho }}_{\mathrm{0}}{c}_p{C}_H({\mathit{\theta }}_z-\mathrm{SST})\left|\mathrm{\Delta }\mathit{u}\right|,\text{(3)}& {\displaystyle }{Q}_{\mathrm{LW}}=Q_{\mathrm{LW},\mathrm{a}}-\mathit{ϵ}\mathit{\sigma }{\mathrm{SST}}_K^{\mathrm{4}},\text{(4)}& {\displaystyle }\mathit{\tau }={\mathit{\rho }}_{\mathrm{0}}{C}_{\mathrm{D}}\mathrm{\Delta }\mathit{u}\left|\mathrm{\Delta }\mathit{u}\right|.\end{array}$$

where z is the height above the sea surface the atmospheric variables are provided at, with the zero values (subscript “0”) representing the values at the sea surface. u is the horizontal wind vector, with $\mathrm{\Delta }\mathit{u}={\mathit{u}}_z-{\mathit{u}}_{\mathrm{0}}$ as the difference between the wind velocity and sea surface current (assuming a no-slip condition at the ocean surface). q and θ are the specific humidity and potential temperature of air, respectively. Λ and c_p are the latent heat of evaporation and the specific heat of water, respectively. ρ₀ is the air density at the sea surface. SST is the sea surface temperature, with SST_K as the sea surface absolute temperature. ϵ is the sea surface emissivity, σ 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):

$$\begin{array}{}\text{(5)}& Q_{\mathrm{net}}=Q_{\mathrm{SW}}+Q_{\mathrm{LW}}+Q_H+Q_E,\end{array}$$

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 $\mathrm{1}/\mathrm{12}$^∘ (roughly 7 km) and 75 vertical levels (with a variable vertical resolution from 1 m at the surface to 135 m at the bottom). The 3-D temperature and salinity fields are restored towards the ORAS4 global ocean reanalysis (Balmaseda et al., 2013) in the buffer zone. The conservation of volume in the buffer zone is achieved through strong damping of the sea surface height (SSH) towards the ORAS4 reanalysis. The Black Sea, runoff of 33 major rivers, and coastal runoff are represented by climatological data input from Ludwig et al. (2009). A deeper explanation of the configuration and boundary conditions is given by in the following works: Waldman et al. (2018), Hamon et al. (2016), Beuvier et al. (2012), Lebeaupin-Brossier et al. (2011), and Arsouze et al. (2012). The initial conditions were provided by an ocean objective analysis by Estournel et al. (2016b).

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

2.2 Atmospheric forcing

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

For the “control” simulation, the forcings was used as is. For the “seasonal” simulation, the u and v wind components, specific humidity, and potential temperature were filtered (see Fig. 2) over the entire domain shown in Fig. 1a. These variables were chosen as they are the primary variables that affect the surface flux calculations in the bulk formulae (Eqs. 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:

$$\begin{array}{}\text{(6)}& {\mathit{\chi }}_i={\displaystyle \frac{\mathrm{1}}{i+N+\mathrm{1}}}\sum _{j=\mathrm{0}}^{i+N}{x}_j,\end{array}$$

where χ_i is the averaged (filtered) value at index i of a time series of variable x with length n, where $i=\mathrm{0}\to 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 “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.

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.

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

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

The average values for mistral events from 30 August 2012 to 16 February 2013 are $\stackrel{\mathrm{‾}}{\mathrm{\Delta }t}=\mathrm{5.69}$ d and $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathit{\tau }}=\mathrm{10.88}$ d. The standard deviations for the same time frame are σ_Δt=4.22 d and σ_Δτ=4.59 d.

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):

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial z}{\partial t}}={\displaystyle \frac{B\left(t\right)}{N^{\mathrm{2}}\left(z\right)z}},\end{array}$$

where N² is the Brunt–Väisälä frequency, z is the vertical coordinate along the water column, $\frac{\partial z}{\partial 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

$$\begin{array}{}\text{(8)}& \mathrm{SI}=\underset{\mathrm{0}}{\overset{D}{\int }}{N}^{\mathrm{2}}z\mathrm{d}z,\end{array}$$

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

$$\begin{array}{}\text{(9)}& \mathrm{SI}={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{N}^{\mathrm{2}}.\end{array}$$

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

3 Winter of 2013

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

3.1 Model validation

To validate the model results, data from the HyMeX (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).

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.

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 profiles and their standard deviations are visualized in Figs. 5 and 6. The locations of the CTD profiles are shown in Fig. 4.

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.

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^−2 ∘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.

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, ${\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{CTD}}}$, is marked in red and is present for all plots but may be difficult to see for 7 March and 9 May.

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

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.

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.

3.2 Stratification index and mixed-layer depth

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

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

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_S+δSI, and the seasonal run, SI_S. 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_A through t_D correspond to the timestamps of the plots in Fig. 7: 30 August 2012, 11 December 2012, 13 February 2013, and 1 June 2013. Two definitions of MLD are plotted in (c): one calculated by a vertical change in density less than 0.01 kg m⁻³, denoted by Δρ, and one calculated by a vertical diffusivity less than $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹, denoted by K_z. The MLD denoted by the vertical diffusivity criteria follows the turbocline depth and is taken to represent the mixed-layer depth more accurately, as it matches the deep-convection timing in the stratification index.

To investigate the time series ocean response in more detail, a spatially averaged time series of the SI for both simulations was analyzed at the grid point nearest to 42^∘ N, 5^∘ E. These coordinates were selected as it is the point with the most destratification in Fig. 7 and is the typical center of deep convection in the GOL (Marshall and 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² s⁻² (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² s⁻², 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² s⁻². Additionally, if only the anomaly timescale atmospheric forcing is considered, and therefore δSI is the only stratification change from the initial 1.57 m² s⁻², the roughly −0.6 m² s⁻² 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² s⁻² 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² s⁻² 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^{\mathrm{2}}=\mathit{\delta }{N}^{\mathrm{2}}+N_{\mathrm{S}}^{\mathrm{2}}$, where δ denotes the anomaly terms and “S” denotes the seasonal terms. To determine the mistral's effect, we derive a simple model to describe the SI of the water column in response to atmospheric forcing (for the full derivation, see Appendix A1). We start with the energy equation for incompressible fluids (White, 2011), then multiply the equation by $-g/T_{\mathrm{0}}$ to express the energy equation in terms of buoyancy, assuming that the ocean's density varies negatively proportionally with temperature, $\mathit{\rho }=-\mathit{\beta }T$. We then perform partial differentiation with respect to z to obtain an equation describing the Brunt–Väisälä frequency, N², in response to the atmospheric forcing, given by a forcing function, F(t). Separating by timescale, we arrive at the following partial differential equations:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{D\mathit{\delta }{N}^{\mathrm{2}}}{Dt}}=-\mathit{\delta }F\left(t\right),\text{(11)}& {\displaystyle \frac{D{N}_{\mathrm{S}}^{\mathrm{2}}}{Dt}}=-F_{\mathrm{S}}\left(t\right).\end{array}$$

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_{\mathrm{S}}^{\mathrm{2}}$ is only a function of time, t:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial N_{\mathrm{S}}^{\mathrm{2}}}{\partial t}}=-F_{\mathrm{S}}\left(t\right).\end{array}$$

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_S, and seasonal atmospheric forcing:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}}=-{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{F}_{\mathrm{S}}\left(t\right).\end{array}$$

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

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \mathit{\delta }{N}^{\mathrm{2}}}{\partial t}}+\mathit{\alpha }\left(\mathit{\delta }{N}^{\mathrm{2}}\right)=-\mathit{\delta }F\left(t\right).\end{array}$$

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_S, to the seasonal atmospheric forcing, F_S(t). We have the following definition of F_S(t) from Appendix A1:

$$\begin{array}{}\text{(15)}& F_{\mathrm{S}}\left(t\right)={\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{{\mathit{q}}_{\mathrm{a},\mathrm{S}}g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}\right)={\displaystyle \frac{g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{a},\mathrm{S}}}{\partial z}},\end{array}$$

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

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}}={\displaystyle \frac{g}{\mathrm{2}\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}{Q}_{\mathrm{net},\mathrm{S}},\end{array}$$

where $\frac{g}{\mathrm{2}\mathit{\rho }{c}_p{T}_{\mathrm{0}}}\approx {\mathrm{10}}^{-\mathrm{9}}$ m⁴ J s⁻², which means the derivative of SI_S w.r.t. time, t, multiplied by 10⁹ is on the same order of magnitude as Q_net (with the subscript S now dropped for convenience, as the rest of the subsection discusses seasonal heat fluxes), which is what we see in Fig. 9a for the 2013 winter, with $\frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}\times {\mathrm{10}}^{\mathrm{9}}$ following the curve of Q_net. This relationship means when Q_net crosses zero with a negative derivative, SI_S experiences a maximum and vice versa for a minimum. Additionally, the longer Q_net remains negative, 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. Consequently, the maximum and minimum values for SI_S occur 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⁻², depending on the time of the year. Q_LW remains roughly constant during the year, decreasing Q_net by roughly −100 W m⁻². These results are corroborated by the results of multiple model reanalysis for the region as well (Song and Yu, 2017).

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_S, and its derivative, $\frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}$, comparing it to the seasonal net heat flux, Q_net (the subscript S is dropped for convenience). Panel (b) shows the net heat flux separated into its components: Q_E, Q_H, Q_SW, and Q_LW for latent heat flux, sensible heat flux, shortwave downward flux, and longwave downward flux, respectively (neglecting contributions from precipitation and snowfall). The different line colors correspond to the similarly colored axes.

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

4.3 Anomaly solution and forcing

To solve for the anomaly timescale, described by Eq. (14), we assume δF(t) can be represented by a pulse function shown in Fig. 10. This pulse function assumes the primary forcing at the anomaly timescale is represented by the mistral events. Each Mistral event, k, has a duration, Δt_k, and a period between the start of the current and following event, Δτ_k. δ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². The last assumption is that the restoring coefficient remains constant for each section of the forcing function:

$$\begin{array}{}\text{(17)}& \mathit{\delta }{\mathrm{SI}}_k\left(t\right)=\left\{\begin{array}{l}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)\right]\\ \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }{t}_k)\\ \left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k)\end{array}\right.,\end{array}$$

where α_d and α_a are the restoring coefficients during ($[t_k,t_k+\mathrm{\Delta }{t}_k)$) and after ($[t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k)$) a mistral event, respectively.

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

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

$$\begin{array}{}\text{(18)}& \mathit{\delta }{\mathrm{SI}}_k\left(t\right)=\left\{\begin{array}{l}\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}[\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left(\frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}-\mathrm{1}\right)\\ +\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)]e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }{t}_k)\\ \frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}[\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left(\frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}-\mathrm{1}\right)\\ +\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)]e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k)\end{array}\right.,\end{array}$$

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

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

Figure 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 α_d and a value of 0.021 d⁻¹ for α_a was found. Panel (c) shows the δSI response using the determined restoration coefficients, given an ideal Mistral event with the average values of 5.69 d for the duration and 10.88 d for the period. The average strength of a mistral, $\mathit{\delta }F=\mathrm{4.01}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻² d⁻¹, was taken from values found in Table A1 from Appendix A.

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

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_k\left(t\right)}{\partial t}}=-{\mathit{\alpha }}_{\mathrm{d}}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\right]e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}.\end{array}$$

The quantity $\frac{\partial \mathit{\delta }{\mathrm{SI}}_k\left(t\right)}{\partial t}$ must be less than zero for destratification to occur, which means if α_d is a positive quantity (refer to Appendix A2 and Fig. 11), $\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}$ must be a positive quantity. If some destratification has already occurred relative to the seasonal stratification, such that $\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)<\mathrm{0}$, then $\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}$ must be larger than $-\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)$ for destratification to occur. Recalling that δ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² s⁻² and ends at −0.41 m² s⁻², a net restratification of 0.07 m² s⁻². This is despite the fact that this event has a positive δF_k value of $\mathrm{3.80}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻² d⁻¹ (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 -\mathrm{0.6}$ m² s⁻² for δSI on 11 December 2012. Proceeding events after this minimum fail to continue to destratify the water column, and restratification instead occurs on the anomaly timescale, even before deep convection occurs. The seasonal stratification, SI_S, and anomaly destratification, δSI, bring the total SI to zero on 13 February 2013, resulting in deep convection.

4.3.2 Dominating mistral attribute

A pertinent question to ask is which attribute of the mistral, the frequency, strength, or duration, is the most important when it drives destratification. Figures 12 and 13 show the results of varying δF, Δt, and Δτ individually (in panels a, b, and c, respectively) in Eq. (18). The other variables are kept at the mean value when not varied. The dashed lines in both figures show the limit of potential destratification per case. What we can see is that stronger mistral events, with an increased value for δF, result in more destratification, with the reverse happening with decreased values. Decreasing the event duration, Δt, results in less destratification; however, increasing event duration causes more destratification up to the limit where the individual events converge into one single long event and the destratification converges to the dashed line limit. After this, there is no additional destratification. Increasing or decreasing the frequency of events (decreasing or increasing the period, Δτ) only minimally changes the accrued destratification, due to the fact that the magnitude of $\frac{\partial \mathit{\delta }\mathrm{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.

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, δF, (b) varies the duration, Δt, and (c) varies the period between events, Δτ. Panel (d) varies the restoration coefficient during the destratification phase, α_d, and (e) varies the restoration coefficient for the restratification phase. Panel (f) varies the number of events.

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

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

$$\begin{array}{}\text{(20)}& \mathrm{d}\mathit{\delta }\mathrm{SI}=\underset{\mathrm{Strength}}{\underbrace{{\displaystyle \frac{\partial \mathit{\delta }\mathrm{SI}}{\partial \mathit{\delta }F}}\mathrm{d}\mathit{\delta }F}}+\underset{\mathrm{Duration}}{\underbrace{{\displaystyle \frac{\partial \mathit{\delta }\mathrm{SI}}{\partial \mathrm{\Delta }t}}\mathrm{d}\mathrm{\Delta }t}}+\underset{\mathrm{Period}}{\underbrace{{\displaystyle \frac{\partial \mathit{\delta }\mathrm{SI}}{\partial \mathrm{\Delta }\mathit{\tau }}}\mathrm{d}\mathrm{\Delta }\mathit{\tau }}},\end{array}$$

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

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_k(t=t_k+\mathrm{\Delta }t)}{\partial k}}={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}{\displaystyle \frac{\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)}{(\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }})}}\\ & \cdot \left(-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}\right)\left(\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau })\end{array},\end{array}$$

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 $-\mathrm{5.93}\times {\mathrm{10}}^{-\mathrm{11}}$ m² s⁻² per event. This confirms the small change in the potential destratification with increasing events. Taking k to infinity and noting that α_d>α_a results in the following equation:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}& =\mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}(t=t_k+\mathrm{\Delta }t)\\ & ={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}}}\right)\end{array},\end{array}$$

We have an equation that describes the potential destratification, δSI_∞, in terms of the mistral attributes, independent of the number of events. Differentiating by the different attributes (see Appendix A5 for the resulting analytical derivations) and plugging in the mean values where appropriate, we arrive at the resulting values: the derivative w.r.t. the strength of the mistrals, $\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}/\partial \mathit{\delta }F$, is equal to a value of $-\mathrm{1.07}\times {\mathrm{10}}^{\mathrm{7}}$ m² d; the derivative w.r.t. the duration, $\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}/\partial \mathrm{\Delta }t$, is equal to $-\mathrm{7.60}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻² d⁻¹; and the derivative w.r.t. the period, $\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}/\partial \mathrm{\Delta }\mathit{\tau }$, is equal to $\mathrm{2.77}\times {\mathrm{10}}^{-\mathrm{3}}$ 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 $-\mathrm{1.28}\times {\mathrm{10}}^{-\mathrm{1}}$ m² s⁻², the duration term has a value of $-\mathrm{3.21}\times {\mathrm{10}}^{-\mathrm{2}}$ m² s⁻², and the period term has a value of $\mathrm{1.27}\times {\mathrm{10}}^{-\mathrm{2}}$ 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 14 The combined effect of Eq. (17) for multiple mistrals with the mistral data from Tables 1 and A1. Panel (a) shows the calculated δSI+SI_S response, while (b) is the calculated simple model δSI vs. the NEMO δSI simulation results. Effects from mistrals after deep convection are included with the dashed blue line and show that mistrals after deep convection can retard the proceeding restratification during the restratification phase.

5 Comparison with additional years

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

5.1 Stratification index

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

Similar to the 2013 winter simulation set, the seasonal run for the 2005 winter does not experience deep convection, again demonstrating the necessity to include forcing on both timescales for deep convection to occur. The 1994 winter reveals something just as interesting: the δSI is more negative for the 1994 winter, at −0.55 m² s⁻², 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 m² s⁻². 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 m² s⁻² for the 2013 winter and 1.73 and 1.79 m² s⁻² for the 1994 and 2005 winters, respectively. This emphasizes the importance of the destratification caused by the seasonal forcing.

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

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

Figure 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 16 The same as Fig. 15 but for the 2005 winter.

5.2 Seasonal forcing

The seasonal sea surface fluxes for both years resemble the fluxes of the 2013 winter (see Fig. 17). The solar radiation component drives the main shape of the SI time series, with the major component contributing to the asymmetry being the latent heat flux, followed by the sensible heat flux. However, in the 2005 winter, the latent heat flux has a larger heat loss value than the 1994 winter, reaching over 300 W m⁻² 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.

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⁹ m⁴ J s⁻². A negative value means heat is leaving the ocean.

5.3 Anomaly forcing

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

The different components of the mistral (strength, δF, duration, Δt, and period, Δτ) are separated in the same manner as for the 2013 winter to determine the main factor of the mistral leading to destratification according to these additional years. For the 1994 winter, the contributions due to strength, duration, and period are equal to $-\mathrm{8.71}\times {\mathrm{10}}^{-\mathrm{2}}$, $-\mathrm{6.39}\times {\mathrm{10}}^{-\mathrm{2}}$, and $\mathrm{5.83}\times {\mathrm{10}}^{-\mathrm{2}}$ m² s⁻² (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 $-\mathrm{1.51}\times {\mathrm{10}}^{-\mathrm{1}}$, $-\mathrm{5.23}\times {\mathrm{10}}^{-\mathrm{2}}$, and $\mathrm{4.89}\times {\mathrm{10}}^{-\mathrm{2}}$ m² s⁻², respectively. The 2005 winter results have the same order of magnitude as the 2013 winter results, with all three winters having the same order of importance for the mistral attributes: first strength, then duration, and finally the length of the period. Only for the 1994 winter was the strength term found to not be as dominant as in the other years, with the term having the same order of magnitude as the other terms. However, as the order of importance was still the same for all 3 years, this aids the conclusion that the strength term of the mistral is generally its most important factor driving destratification.

6 Conclusions

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

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

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

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

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

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

Another question brought about by the simulation results is what effect the maximal stratification at the beginning of the preconditioning phase has on the ocean's ability to experience deep convection for a given year. The maximal stratification must be overcome by the mistral and seasonal forcing for deep convection to occur, with the ability of the forcing to do so varying per year. For our results, the 2013 winter had a larger maximum control SI than the other two winters investigated and is still deeply convected. Similarly, the 2005 winter had a larger maximum stratification than the winter of 1994: 1.79 vs. 1.73 m² s⁻². However, all three winters had a maximal stratification within 0.1 m² s⁻² 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² s⁻² occurs for the 2012–2013 winter. By December, the SI has already dropped to ∼1 m² s⁻², which means almost half of the destratification has already occurred, with their calculation missing about half of the preconditioning phase. The same is true when looking at the 1994 and 2005 winters. A maximum SI calculated near the beginning of September may be more representative of the water column's ability to deeply convect and needs to be investigated. A related question is how the accumulation (or reduction) of stratification being transferred to the proceeding year affects deep convection in following years. For the winter of 2005, where very few of the preceding 15 years deeply convected due to milder winters, the result was warmer and saltier WDMW production (Herrmann et al., 2010).

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

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

Appendix A: Equations

A1 Simple stratification index model

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

$$\begin{array}{}\text{(A1)}& \mathit{\rho }{c}_p{\displaystyle \frac{DT}{Dt}}={\displaystyle \frac{D\mathit{q}}{Dt}},\end{array}$$

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. $\frac{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:

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{D\mathit{q}}{Dt}}=-{\mathit{q}}_{\mathrm{a}},\end{array}$$

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

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{DT}{Dt}}=-{\displaystyle \frac{{\mathit{q}}_{\mathrm{a}}}{\mathit{\rho }{c}_p}}.\end{array}$$

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

$$\begin{array}{}\text{(A4)}& N^{\mathrm{2}}={\displaystyle \frac{\partial b}{\partial z}}=-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial z}}.\end{array}$$

Assuming a fluid whose density varies negatively proportionally to the temperature, $\mathit{\rho }=-\mathit{\beta }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:

$$\begin{array}{}\text{(A5)}& N^{\mathrm{2}}={\displaystyle \frac{\partial b}{\partial z}}=-{\displaystyle \frac{g}{T_{\mathrm{0}}}}{\displaystyle \frac{\partial T}{\partial z}}.\end{array}$$

Introducing buoyancy as $b=-\frac{g}{T_{\mathrm{0}}}T$, we can rearrange the energy equation into terms of buoyancy:

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{g}{T_{\mathrm{0}}}}\left[{\displaystyle \frac{DT}{Dt}}=-{\displaystyle \frac{{\mathit{q}}_{\mathrm{a}}}{\mathit{\rho }{c}_p}}\right]\Rightarrow {\displaystyle \frac{Db}{Dt}}=-{\displaystyle \frac{{\mathit{q}}_{\mathrm{a}}g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}.\end{array}$$

If we then differentiate by $\frac{\partial }{\partial z}$, we can reorganize the equation in terms of the Brunt–Väisälä frequency, N²:

$$\begin{array}{}\text{(A7)}& {\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{Db}{Dt}}=-{\displaystyle \frac{{\mathit{q}}_{\mathrm{a}}g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}\right)\Rightarrow {\displaystyle \frac{D{N}^{\mathrm{2}}}{Dt}}=-{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{{\mathit{q}}_{\mathrm{a}}g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}\right),\end{array}$$

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:

$$\begin{array}{}\text{(A8)}& F\left(t\right)={\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{{\mathit{q}}_{\mathrm{a}}g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}\right).\end{array}$$

This brings us to

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{D{N}^{\mathrm{2}}}{Dt}}=-F\left(t\right).\end{array}$$

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:

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}& N^{\mathrm{2}}=N_{\mathrm{S}}^{\mathrm{2}}+\mathit{\delta }{N}^{\mathrm{2}},\\ & F\left(t\right)=\mathit{\delta }F\left(t\right)+F_{\mathrm{S}}\left(t\right),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A11)}& {\displaystyle \frac{D\mathit{\delta }{N}^{\mathrm{2}}}{Dt}}=-\mathit{\delta }F\left(t\right),\end{array}$$

and the response of N² on the seasonal timescale, $N_{\mathrm{S}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{D{N}_{\mathrm{S}}^{\mathrm{2}}}{Dt}}=-F_{\mathrm{S}}\left(t\right).\end{array}$$

For the seasonal response, we further make the assumption that $N_{\mathrm{S}}^{\mathrm{2}}$ negligibly depends on the x, y, and z coordinate directions:

$$\begin{array}{}\text{(A13)}& {\displaystyle \frac{\partial N_{\mathrm{S}}^{\mathrm{2}}}{\partial t}}=-F_{\mathrm{S}}\left(t\right).\end{array}$$

If we want to connect the overall Brunt–Väisälä frequency, N², to the seasonal one, $N_{\mathrm{S}}^{\mathrm{2}}$, we can formulate a restoring term, R, in terms of T or in terms of N², following the steps mentioned above:

$$\begin{array}{}\text{(A14)}& R={\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{g}{T_{\mathrm{0}}}}\mathit{\alpha }\mathit{\rho }{c}_p(T-T_{\mathrm{S}})\right)\Rightarrow \mathit{\alpha }(N^{\mathrm{2}}-N_{\mathrm{S}}^{\mathrm{2}}),\end{array}$$

or with $\mathit{\delta }{N}^{\mathrm{2}}=N^{\mathrm{2}}-N_{\mathrm{S}}^{\mathrm{2}}$ as follows:

$$\begin{array}{}\text{(A15)}& R=\mathit{\alpha }\mathit{\delta }{N}^{\mathrm{2}},\end{array}$$

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

$$\begin{array}{}\text{(A16)}& {\displaystyle \frac{\partial \mathit{\delta }{N}^{\mathrm{2}}}{\partial t}}+\mathit{V}\cdot \mathrm{\nabla }\left(\mathit{\delta }{N}^{\mathrm{2}}\right)=-\mathit{\delta }F\left(t\right).\end{array}$$

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

$$\begin{array}{}\text{(A17)}& {\displaystyle \frac{\partial \mathit{\delta }{N}^{\mathrm{2}}}{\partial t}}+\mathit{\alpha }\mathit{\delta }{N}^{\mathrm{2}}=-\mathit{\delta }F\left(t\right).\end{array}$$

A1.1 Solution for seasonal SI

To solve the response of N² for the seasonal timescale, given by Eq. (A13), we will assume $N_{\mathrm{S}}^{\mathrm{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:

$$\begin{array}{}\text{(A18)}& F_{\mathrm{S}}\left(t\right)={\displaystyle \frac{g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{a},\mathrm{S}}}{\partial z}}.\end{array}$$

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

$$\begin{array}{}\text{(A19)}& {\displaystyle \frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}}={\displaystyle \frac{D}{\mathrm{2}}}{\displaystyle \frac{g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}{\displaystyle \frac{\partial Q_{\mathrm{net},\mathrm{S}}}{\partial z}},\end{array}$$

Integrating this equation by z gives us

$$\begin{array}{}\text{(A20)}& \begin{array}{rl}& \underset{\mathrm{0}}{\overset{D}{\int }}{\displaystyle \frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}}\partial z={\displaystyle \frac{D}{\mathrm{2}}}{\displaystyle \frac{g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}\underset{\mathrm{0}}{\overset{D}{\int }}{\displaystyle \frac{\partial Q_{\mathrm{net},\mathrm{S}}}{\partial z}}\partial z\\ & {\displaystyle \frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}}D={\displaystyle \frac{D}{\mathrm{2}}}{\displaystyle \frac{g}{\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}{Q}_{\mathrm{net},\mathrm{S}}\\ & {\displaystyle \frac{\partial {\mathrm{SI}}_{\mathrm{S}}}{\partial t}}={\displaystyle \frac{g}{\mathrm{2}\mathit{\rho }{c}_p{T}_{\mathrm{0}}}}{Q}_{\mathrm{net},\mathrm{S}},\end{array}\end{array}$$

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, δF(t), as a series of k pulse functions of magnitude δF_k over a duration of Δt_k and with a period of Δτ_k, as visualized in Fig. 10.

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

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\delta }{N}^{\mathrm{2}}}{\partial t}}+{\mathit{\alpha }}_{\mathrm{d}}\left(\mathit{\delta }{N}^{\mathrm{2}}\right)=-\mathit{\delta }F\left(t\right)\\ & \mathit{\delta }{N}_k^{\mathrm{2}}\left(t\right)=-{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}+c_{\mathrm{0}}{e}^{-{\mathit{\alpha }}_{\mathrm{d}}t}\\ & \mathit{\delta }{N}_k^{\mathrm{2}}\left(t_k\right)=-{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}+c_{\mathrm{0}}{e}^{-{\mathit{\alpha }}_{\mathrm{d}}{t}_k}=\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)\\ & c_{\mathrm{0}}=\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\right]e^{{\mathit{\alpha }}_{\mathrm{d}}{t}_k},\end{array}\end{array}$$

with the following result:

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}\mathit{\delta }{N}_k^{\mathrm{2}}\left(t\right)& =\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)\right]\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\delta }{N}^{\mathrm{2}}}{\partial t}}+{\mathit{\alpha }}_{\mathrm{a}}\left(\mathit{\delta }{N}^{\mathrm{2}}\right)=\mathrm{0}\\ & \mathit{\delta }{N}_k^{\mathrm{2}}\left(t\right)=c_{\mathrm{1}}{e}^{-{\mathit{\alpha }}_{\mathrm{a}}t}\\ & \mathit{\delta }{N}_k^{\mathrm{2}}(t_k+\mathrm{\Delta }{t}_k)=c_{\mathrm{1}}{e}^{-{\mathit{\alpha }}_{\mathrm{a}}(t_k+\mathrm{\Delta }{t}_k)}\\ & =\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\\ & c_{\mathrm{1}}=\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k+{\mathit{\alpha }}_{\mathrm{a}}{t}_k},\end{array}\end{array}$$

With the following result:

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}\mathit{\delta }{N}_k^{\mathrm{2}}\left(t\right)& =\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)},\end{array}\end{array}$$

or to have the results put more succinctly:

$$\begin{array}{}\text{(A25)}& \mathit{\delta }{N}_k^{\mathrm{2}}\left(t\right)=\left\{\begin{array}{l}\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)\right]\\ \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }{t}_k)\\ \left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k).\end{array}\right.\end{array}$$

A1.3 $\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}$ initial condition

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

$$\begin{array}{}\text{(A26)}& \begin{array}{rl}& \mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)=\left[\mathit{\delta }{N}_{k-\mathrm{2}}^{\mathrm{2}}\left(t_{k-\mathrm{1}}\right)+{\displaystyle \frac{\mathit{\delta }{F}_{k-\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{d}}}}(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_{k-\mathrm{1}}})\right]\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_{k-\mathrm{1}}-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }{\mathit{\tau }}_{k-\mathrm{1}}}\\ & \mathit{\delta }{N}_{k-\mathrm{2}}^{\mathrm{2}}\left(t_{k-\mathrm{1}}\right)=\left[\mathit{\delta }{N}_{k-\mathrm{3}}^{\mathrm{2}}\left(t_{k-\mathrm{2}}\right)+{\displaystyle \frac{\mathit{\delta }{F}_{k-\mathrm{2}}}{{\mathit{\alpha }}_{\mathrm{d}}}}(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_{k-\mathrm{2}}})\right]\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_{k-\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }{\mathit{\tau }}_{k-\mathrm{2}}}.\end{array}\end{array}$$

Therefore, $\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)$ can be simplified in expression by combining the initial conditions:

$$\begin{array}{}\text{(A27)}& \begin{array}{rl}\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)& =\mathit{\delta }{N}_{k-m}^{\mathrm{2}}\left(t_{k-(m-\mathrm{1})}\right)e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}}){\sum }_{i=\mathrm{1}}^{m-\mathrm{1}}\mathrm{\Delta }{t}_{k-i}}\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{a}}{\sum }_{i=\mathrm{1}}^{m-\mathrm{1}}\mathrm{\Delta }{\mathit{\tau }}_{k-i}}+\sum _{j=\mathrm{1}}^{m-\mathrm{1}}{\displaystyle \frac{\mathit{\delta }{F}_{k_j}}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_{k-j}}\right)\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}}){\sum }_{i=\mathrm{1}}^j\mathrm{\Delta }{t}_{k-i}}{e}^{-{\mathit{\alpha }}_{\mathrm{a}}{\sum }_{i=\mathrm{1}}^j\mathrm{\Delta }{\mathit{\tau }}_{k-i}}.\end{array}\end{array}$$

If m=k and $\mathit{\delta }{N}_{\mathrm{0}}^{\mathrm{2}}=\mathrm{0}$ then

$$\begin{array}{}\text{(A28)}& \begin{array}{rl}\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)& =\sum _{j=\mathrm{1}}^{k-\mathrm{1}}{\displaystyle \frac{\mathit{\delta }{F}_{k_j}}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_{k-j}}\right)\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}}){\sum }_{i=\mathrm{1}}^j\mathrm{\Delta }{t}_{k-i}}{e}^{-{\mathit{\alpha }}_{\mathrm{a}}{\sum }_{i=\mathrm{1}}^j\mathrm{\Delta }{\mathit{\tau }}_{k-i}}.\end{array}\end{array}$$

Assuming δF_k=δF, Δt_k=Δt, and Δτ_k=Δτ for all k, or a periodic pulse function, then $\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}$ can be expressed as follows:

$$\begin{array}{}\text{(A29)}& \mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)={\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\sum _{j=\mathrm{1}}^{k-\mathrm{1}}{e}^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]j}.\end{array}$$

Taking the sum of a finite geometric series as follows:

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}\sum _{n=\mathrm{0}}^{m-\mathrm{1}}{r}^n& =\left({\displaystyle \frac{\mathrm{1}-r^m}{\mathrm{1}-r}}\right),\\ \sum _{n=\mathrm{0}}^{m-\mathrm{1}}{r}^n& =\sum _{n=\mathrm{1}}^{m-\mathrm{1}}{r}^n+\mathrm{1},\\ \sum _{n=\mathrm{1}}^{m-\mathrm{1}}{r}^n& =\sum _{n=\mathrm{0}}^{m-\mathrm{1}}{r}^n-\mathrm{1},\end{array}\end{array}$$

where r≠1. If $r=e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}$, then

$$\begin{array}{}\text{(A31)}& \sum _{j=\mathrm{1}}^{k-\mathrm{1}}{e}^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]j}=\left({\displaystyle \frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}}-\mathrm{1}\right),\end{array}$$

and we then we get

$$\begin{array}{}\text{(A32)}& \mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)={\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left({\displaystyle \frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}}-\mathrm{1}\right),\end{array}$$

where $({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }\ne \mathrm{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:

$$\begin{array}{}\text{(A33)}& \mathit{\delta }{N}^{\mathrm{2}}\left(t\right)=\left\{\begin{array}{l}\frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}[\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left(\frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}-\mathrm{1}\right)\\ +\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)]e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }t)\\ \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}[\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left(\frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}-\mathrm{1}\right)\\ +\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)]e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }t,t_k+\mathrm{\Delta }\mathit{\tau }).\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(A34)}& \mathit{\delta }{\mathrm{SI}}_k\left(t\right)=\left\{\begin{array}{l}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)\right]\\ \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }{t}_k)\\ \left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k),\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(A35)}& \mathit{\delta }{\mathrm{SI}}_k\left(t\right)=\left\{\begin{array}{l}\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}[\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left(\frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}-\mathrm{1}\right)\\ +\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)]e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }t)\\ \frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}[\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)\left(\frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]}}-\mathrm{1}\right)\\ +\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}\right)]e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }t,t_k+\mathrm{\Delta }\mathit{\tau }),\end{array}\right.\end{array}$$

for the period pulse function case.

A2 Restoring coefficients α_d and α_a

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

A2.1 Restoration coefficient α_d during a mistral

To solve for α_d, we normalize Eq. (A34) for during a mistral event, $[t_k,t_k+\mathrm{\Delta }{t}_k)$, with δSI given as follows:

$$\begin{array}{}\text{(A36)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_k\left(t\right)& =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\right)\right]\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)},\end{array}\end{array}$$

We first reference the time, t, to the starting time of event k as $t^{\prime }=t-t_k$, giving us

$$\begin{array}{}\text{(A37)}& \mathit{\delta }{\mathrm{SI}}_k\left(t^{\prime }\right)=\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{\mathit{\alpha }}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }}\right)\right]e^{-{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }},\end{array}$$

Next, we normalize δSI_k to the value of zero at $t^{\prime }=\mathrm{0}$, resulting in δSI_k,NI:

$$\begin{array}{}\text{(A38)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NI}}\left(t^{\prime }\right)& =\mathit{\delta }{\mathrm{SI}}_k\left(t^{\prime }\right)-\mathit{\delta }{\mathrm{SI}}_k(t^{\prime }=\mathrm{0})\\ & =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }}\right)\right]\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }}-\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)\\ & =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\right]\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }}-\mathrm{1}\right),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A39)}& \mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NH}}\left(t^{\prime }\right)={\displaystyle \frac{\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NI}}\left(t^{\prime }\right)}{\mathrm{extremum}\left(\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NI}}\left(t^{\prime }\right)\right)}}.\end{array}$$

The extremum value for $\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NI}}\left(t^{\prime }\right)$ is when $t^{\prime }=\mathrm{\Delta }{t}_k$ or at the end of the mistral event. This simplifies δSI_k,NH to

$$\begin{array}{}\text{(A40)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NH}}\left(t^{\prime }\right)& ={\displaystyle \frac{\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\right]\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }}-\mathrm{1}\right)}{\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\right]\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}-\mathrm{1}\right)}}\\ & ={\displaystyle \frac{\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime }}-\mathrm{1}\right)}{\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}-\mathrm{1}\right)}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A41)}& t^{\prime \prime }={\displaystyle \frac{t^{\prime }}{\mathrm{\Delta }{t}_k}}\Rightarrow t^{\prime }=t^{\prime \prime }\mathrm{\Delta }{t}_k.\end{array}$$

Plugging t^′′ into $\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NH}}\left(t^{\prime }\right)$ returns δSI_k,NT:

$$\begin{array}{}\text{(A42)}& \mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NT}}\left(t^{\prime \prime }\right)={\displaystyle \frac{e^{-{\mathit{\alpha }}_{\mathrm{d}}{t}^{\prime \prime }\mathrm{\Delta }{t}_k}-\mathrm{1}}{e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}-\mathrm{1}}}.\end{array}$$

This final equation, δSI_k,NT, can be used with a fitting function to solve for α_d if Δt_k is supplied.

A2.2 Restoration coefficient α_a, after a mistral

To solve for α_a, we normalize Eq. (A34) for after a mistral event, $[t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k)$, with δSI given as

$$\begin{array}{}\text{(A43)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_k\left(t\right)& =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ & \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)},\end{array}\end{array}$$

Referencing the time, t, to the end of the event, $t^{\prime \prime \prime }=t-(t_k+\mathrm{\Delta }{t}_k)$ and plugging t^′′′ into Eq. (A43), we get

$$\begin{array}{}\text{(A44)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_k\left(t^{\prime \prime \prime }\right)& =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}{e}^{-{\mathit{\alpha }}_{\mathrm{a}}{t}^{\prime \prime \prime }},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A45)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NI}}& =\mathit{\delta }{\mathrm{SI}}_k\left(t^{\prime \prime \prime }\right)-\mathit{\delta }{\mathrm{SI}}_k(t^{\prime \prime \prime }=\mathrm{0})\\ & =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}{e}^{-{\mathit{\alpha }}_{\mathrm{a}}{t}^{\prime \prime \prime }}\\ & -\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\\ & =\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ & \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\left(e^{-{\mathit{\alpha }}_{\mathrm{a}}{t}^{\prime \prime \prime }}-\mathrm{1}\right),\end{array}\end{array}$$

Each post-event restratification is normalized to the height of 1 by dividing δSI_k,NI by $\left(\mathit{\delta }{\mathrm{SI}}_k\right(t^{\prime \prime \prime }=\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k)-\mathit{\delta }{\mathrm{SI}}_k(t^{\prime \prime \prime }=\mathrm{0}\left)\right)$:

$$\begin{array}{}\text{(A46)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NH}}& ={\displaystyle \frac{\mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NI}}}{\mathit{\delta }{\mathrm{SI}}_k(t^{\prime \prime \prime }=\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k)-\mathit{\delta }{\mathrm{SI}}_k(t^{\prime \prime \prime }=\mathrm{0})}}\\ & \\ & ={\displaystyle \frac{\begin{array}{c}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\left(e^{-{\mathit{\alpha }}_{\mathrm{a}}{t}^{\prime \prime \prime }}-\mathrm{1}\right)\end{array}}{\begin{array}{c}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\left(e^{-{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k)}-\mathrm{1}\right)\end{array}}}\\ & ={\displaystyle \frac{e^{-{\mathit{\alpha }}_{\mathrm{a}}{t}^{\prime \prime \prime }}-\mathrm{1}}{e^{-{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k)}-\mathrm{1}}},\end{array}\end{array}$$

which gives us δSI_k,NH.

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

$$\begin{array}{}\text{(A47)}& t^{\prime \prime \prime \prime }={\displaystyle \frac{t^{\prime \prime \prime }}{\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k}}\Rightarrow t^{\prime \prime \prime }=t^{\prime \prime \prime \prime }(\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k),\end{array}$$

which leads to δSI_k,NT:

$$\begin{array}{}\text{(A48)}& \mathit{\delta }{\mathrm{SI}}_{k,\mathrm{NT}}={\displaystyle \frac{e^{-{\mathit{\alpha }}_{\mathrm{a}}{t}^{\prime \prime \prime \prime }(\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k)}-\mathrm{1}}{e^{-{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{\Delta }{\mathit{\tau }}_k-\mathrm{\Delta }{t}_k)}-\mathrm{1}}}.\end{array}$$

This leaves us with an equation of α_a that can be fitted against NEMO model data if Δt_k and Δτ_k are provided.

The average duration and period of events during the preconditioning period in Table 1, $\stackrel{\mathrm{‾}}{\mathrm{\Delta }t}$ and $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathit{\tau }}$, respectively, are used for Δt and Δτ in these normalized equations, Eqs. (A42) and (A48). The result of the fitting is shown in Fig. 11, with α_d having a fitted valued of 0.235 d⁻¹ and α_a having a fitted value of 0.021 d⁻¹.

A3 Determining δ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, δF_k, can be determined. If we take Eq. (A36) and note that the value of δSI_k−1 is the same as δSI_k(t_k) at the beginning of an event, we can simplify the equation in the following steps:

$$\begin{array}{}\text{(A49)}& \begin{array}{rl}& \mathit{\delta }{\mathrm{SI}}_k\left(t_k\right)=\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}},\\ & \mathit{\delta }{\mathrm{SI}}_k(t_k+\mathrm{\Delta }{t}_k)=\\ & \mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}{e}^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}+{\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}-\mathrm{1}\right),\end{array}\end{array}$$

and then solve for δF_k:

$$\begin{array}{}\text{(A50)}& \mathit{\delta }{F}_k={\displaystyle \frac{\mathrm{2}\left(\mathit{\delta }{\mathrm{SI}}_k(t_k+\mathrm{\Delta }{t}_k)-\mathit{\delta }{\mathrm{SI}}_k\left(t_k\right)e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right){\mathit{\alpha }}_{\mathrm{d}}}{\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}-\mathrm{1}\right)D^{\mathrm{2}}}}.\end{array}$$

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

Table A1 The mistral strengths, δF_k, for each of the preconditioning phase events, using α_d and α_a from Sect. A2 and the rest of the preconditioning period mistral characteristics from Table 1 plugged into Eq. (A50).

$\stackrel{\mathrm{‾}}{\mathit{\delta }{F}_k}=\mathrm{4.01}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻² d⁻¹, and ${\mathit{\sigma }}_{\mathit{\delta }{F}_k}=\mathrm{1.196}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻² d⁻¹.

A4 Time derivative of δN²

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

$$\begin{array}{}\text{(A51)}& {\displaystyle \frac{\partial \mathit{\delta }{N}_k^{\mathrm{2}}\left(t\right)}{\partial t}}=\left\{\begin{array}{l}-{\mathit{\alpha }}_{\mathrm{d}}\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\right]e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }{t}_k)\\ -{\mathit{\alpha }}_{\mathrm{a}}\left[\mathit{\delta }{N}_{k-\mathrm{1}}^{\mathrm{2}}\left(t_k\right)+\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k),\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(A52)}& {\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_k\left(t\right)}{\partial t}}=\left\{\begin{array}{l}-{\mathit{\alpha }}_{\mathrm{d}}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\right]e^{-{\mathit{\alpha }}_{\mathrm{d}}(t-t_k)}\\ [t_k,t_k+\mathrm{\Delta }{t}_k)\\ -{\mathit{\alpha }}_{\mathrm{a}}\left[\mathit{\delta }{\mathrm{SI}}_{k-\mathrm{1}}\left(t_k\right)+\frac{D^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\delta }{F}_k}{{\mathit{\alpha }}_{\mathrm{d}}}\left(\mathrm{1}-e^{{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }{t}_k}\right)\right]\\ \cdot e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }{t}_k-{\mathit{\alpha }}_{\mathrm{a}}(t-t_k)}\\ [t_k+\mathrm{\Delta }{t}_k,t_k+\mathrm{\Delta }{\mathit{\tau }}_k).\end{array}\right.\end{array}$$

A5 Asymptotic destratification

The following sections differentiate Eq. (A35) at $t=t_k+\mathrm{\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., δF, Δt, and Δτ, once k→∞.

A5.1  $\frac{\partial \mathit{\delta }{\mathrm{SI}}_k}{\partial k}$

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

$$\begin{array}{}\text{(A53)}& \begin{array}{rl}\mathit{\delta }{\mathrm{SI}}_k(t_k+\mathrm{\Delta }t)& ={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)\\ & \cdot \left({\displaystyle \frac{\mathrm{1}-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}}{\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}}}\right),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A54)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_k}{\partial k}}& ={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}{\displaystyle \frac{\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)}{(\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }})}}\\ & \cdot \left(-e^{\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }]k}\right)\left[\right({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }],\end{array}\end{array}$$

As k→∞, with α_d>α_a Eq. (A53) goes to

$$\begin{array}{}\text{(A55)}& \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}}}\right).\end{array}$$

A5.2  $\frac{\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}}{\partial \mathit{\delta }F}$

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

$$\begin{array}{}\text{(A56)}& {\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}}{\partial \mathit{\delta }F}}={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}}}\right).\end{array}$$

A5.3  $\frac{\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}}{\partial \mathrm{\Delta }t}$

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

$$\begin{array}{}\text{(A57)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}}{\partial \mathrm{\Delta }t}}& ={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }})}}\\ & \cdot [-{\mathit{\alpha }}_{\mathrm{d}}{e}^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}+\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right)\\ & \cdot \left({\displaystyle \frac{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}}{\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}}}\right)].\end{array}\end{array}$$

A5.4  $\frac{\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}}{\partial \mathrm{\Delta }\mathit{\tau }}$

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

$$\begin{array}{}\text{(A58)}& {\displaystyle \frac{\partial \mathit{\delta }{\mathrm{SI}}_{\mathrm{\infty }}}{\partial \mathrm{\Delta }\mathit{\tau }}}={\displaystyle \frac{D^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\delta }F}{{\mathit{\alpha }}_{\mathrm{d}}}}\left(e^{-{\mathit{\alpha }}_{\mathrm{d}}\mathrm{\Delta }t}-\mathrm{1}\right){\displaystyle \frac{\left(-{\mathit{\alpha }}_{\mathrm{a}}{e}^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}\right)}{{\left(\mathrm{1}-e^{({\mathit{\alpha }}_{\mathrm{a}}-{\mathit{\alpha }}_{\mathrm{d}})\mathrm{\Delta }t-{\mathit{\alpha }}_{\mathrm{a}}\mathrm{\Delta }\mathit{\tau }}\right)}^{\mathrm{2}}}}.\end{array}$$

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 B1 The start date, duration (Δt_k, in d), and period between each event (Δτ_k, in d) for each Mistral event (k) and event strength, (δF_k, in s⁻² d⁻¹) for the preconditioning phase of the NEMO simulation period of 1 June 1993 to 31 May 1994.

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

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

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