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 2/3 of the buoyancy loss required for deep convection to occur, for the 2012 to 2013 year, with the Mistral supplying the final 1/3. Two NEMOMED12 ocean simulations of the Mediterranean Sea were run for the Aug. 1st, 2012 to July 31st, 2013 year, 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 5 was filtered to remove the signal of the Mistral. The Control simulation featured deep convection, while the Seasonal did not. A simple model was derived, 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 10 stratification maximum and minimum occurring roughly at the fall and spring equinoxes. Copyright statement.

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). It's role was determined by running two NEMO ocean simulations 35 of the Med. Sea from Aug. 1st, 2012 to July 31st, 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 was collected during this year in the framework of the HYdrological cycle in the Mediterranean EXperiment (HyMeX) (Estournel et al., 2016b;Drobinski 40 et al., 2014), which provided a solid base of observations to validate the ocean model results.
In particular, our findings quantify: 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. 45 Two additional years were also studied with the same methodology as the 2012-2013 winter: the 1993-1994 and 2004-2005 winter. The 1993-1994 winter does not have a deep convection event, and allows us to compare a deep convecting year versus 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. 50 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/spreading phase during the proceeding summer (MEDOC, 1970;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 used and the methodology is described in the methods section (Sec. 2). Model results and validation are presented where Q SW is the downward shortwave radiation. Snowfall and precipitation are included in the simulation calculations 85 but excluded here for brevity in the following sections, as the buoyancy loss due to the water flux at the surface is essentially negligible ).
The NEMO model was run in the NEMOMED12 configuration using NEMO v3.6. The domain is shown in Fig. 1 (b); it covers the Mediterranean Sea and a buffer zone representing the exchanges with the Atlantic Ocean. This configuration features a horizontal resolution of 1/12 • (roughly 7km) and 75 vertical levels (with a variable vertical resolution from 1m at the surface 90 to 135m 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 is represented by climatological data input from Ludwig et al. (2009). A deeper explanation of the configuration and boundary conditions is given by the works Waldman et al. (2018); Hamon et al. (2016); Beuvier et al. (2012); Lebeaupin-Brossier et al.

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 hind-cast simulation (ERA interim downscaling), 100 performed at 20 km resolution, 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 forcing were 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. 1 (a). These variables 105 were chosen as they are the primary variables that affect the surface flux calculations in the bulk formulae (Eq. set (1)). 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 under a month timescale), effectively removing the Mistral's influence on the ocean response. This creates two separate forcing datasets, one with the anomaly scale forcing included, one with just the seasonal scale forcing (hence the designation of Control and Seasonal).

110
The filtering process was performed by a moving window average: where χ i is the averaged (filtered) value at index i of a time series of variable x with length n, where i = 0 → n. The window size is equal to 2N + 1, which, in this case, is equal to 31 days. 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 115 the ocean simulation dates, as two, full year atmospheric forcing data were used for the simulations.
The moving window average was applied to each time point per day over a 31 day window (i.e. for 3 hourly data, the time series is split into 8 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 120 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 , 2011 and is shorter than the typical Mistral event length of about 5.69 days (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 125 on the seasonal scale. This is due to the fact 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. But 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 than intended that isn't influenced by the Mistral. This reveals a very interesting point 130 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 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   and on the order of decades for total circulation (Millot and Taupier-Letage, 2005), we assume the processes outside the NW Med. subdomain in Fig. 1 (b), that are affected by the filtering have a negligible impact on the GOL processes on the preconditioning phase timescale.

Mistral events
Mistral events will be used for developing the simple model in the process analysis section (Sec. 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:  Givon et al. (2021)).
The events during the preconditioning period, Aug. 30th, 2012 to Feb. 21st, 2013, were then manually checked and edited 150 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  Turner (1973)), which has been shown to be an accurate approximation for open ocean convection (Marshall and Schott, 1999): where N 2 is the Brunt-Väisälä frequency, z is the vertical coordinate along the water column, ∂z ∂t is the growth of the mixed layer depth, and B is the potential buoyancy loss the water column can endure before removing stratification (in units of m 2 /s 3 ). Separating by variables and integrating results in the equation for SI: where D is the depth of water column. If N 2 is assumed to be constant throughout the water column, the integral simplifies 165 to: SI provides a 0 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.

175
To validate the model results, data from the HyMeX (https://www.hymex.org/) database was 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 was missing SST measurements from Jan. 19th, 2013 to July 10th, 2013, but where the data is available, NEMO corresponds well to the observations. The same is true for the Lion buoy data, which had measurements for the entire time covered by the simulations.

180
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 (Eq. (1)).
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 The CTD profiles collected at approximately the same time and location were averaged together to adjust for small variances  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 less than the average RMSE for the potential temperature profiles: 0.025 kg/m 3 and 0.094 • C, respectively.

195
Argo float profiles from the HyMeX database were also compared to the Control simulation, again with profiles from the nearest grid point being used. 3118 potential temperature profiles within the box bounded by the 40 to 44 • N latitudes and the 2 to 8 • longitudes, to represent 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 shear volume of profiles and, during stratified conditions the 200 temperature can range a few degrees from the surface to the lower layers, these results aren't 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 −2 kg/m 3 , are used to separate newly formed dense water during deep convection 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, 205 given the bias and RMSE values from comparing the simulation to CTD and Argo profiles. Additionally, meanders around 40km in wavelength form due to baroclinic instability along edge of the convection patch (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.    To investigate the time series ocean response in more detail, a spatially averaged time series of the SI for both simulations 220 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 9 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 , The average RMSE and bias for the density profiles was 0.025 kg/m 3 and -0.009 kg/m 3 , respectively. The average RMSE and bias for potential temperature was 0.094 • C and 0.032 • C, respectively. while the Stratification Index of the Seasonal simulation is given as SI S . The difference between the two, δSI, should contain 225 the change in stratification due to shorter timescale atmospheric events, such as the Mistral, because of the filtering performed in Sec. 2.2. δSI + SI S , SI S , and δSI are all shown in Fig. 8.
Both the Control and Seasonal runs start off with an SI value of 1.57 m 2 /s 2 (beginning of Fig. 8 (a)), then diverge at the first major Mistral event starting on August 30th, 2012. After diverging, the two runs remain diverged until the end of the simulation run time, ending with a difference of about -0.22 m 2 /s 2 , which is seen in δSI (shown in Fig. 8 (b)). As commented earlier, 230 the most striking difference between the Control and Seasonal run is the occurrence of deep convection in the Control run, occurring when δSI + SI S is equal to 0 (signified also when the MLD reaches the sea floor), and the lack of deep convection in the Seasonal run, as SI S only reaches a minimum of 0.43 m 2 /s 2 . Additionally, if only the anomaly timescale atmospheric Kz. The MLD denoted by the vertical diffusivity criteria follows the turbocline depth and is taken to represent the mixed layer depth more accurately, as it matches the deep convection timing in the Stratification Index.
forcing is considered, hence δSI is the only stratification change from the initial 1.57 m 2 /s 2 , the roughly -0.6 m 2 /s 2 of maximum destratification that the anomaly timescale provides is not enough to overcome the initial stratification. This means 235 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 Feb. 13th, 2013, which is before SI S reaches its minimum on Feb. 21th, 2013, but after δSI reaches its minimum on Dec. 11th, 2012. After δSI reaches its minimum, it stays around -0.43 m 2 /s 2 until May 2013, where it starts to increase. This means that while the 240 induced destratification from the anomaly scale forcing would have been able to overcome 0.6 m 2 /s 2 of stratification to form deep convection in Dec., the seasonal stratification was only low enough in Feb. 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 245 more clearly in the MLD, as the MLD grows over two Mistral events preceding it reaching the seafloor.

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 time scales.

Simple model derivation 250
To connect the Mistral to the ocean's response, we make the assumption that the response is a superposition of the seasonal response and the anomaly response. This means the effects of the Mistral can be categorized as anomalies affecting the short term anomaly timescale and studied separately from the seasonal response. In terms of the Brunt-Väisälä frequency, this linear combination is represented by N 2 = δN 2 + N 2 S , 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 255 forcing (note: the full derivation is found in App. A1). We start with the energy equation for incompressible fluids (White, 2011), then multiply the equation by −g/T 0 to express the energy equation in terms of buoyancy, assuming that the ocean's density varies negatively proportionally with temperature, ρ = −βT . We then perform partial differentiation with respect to z to obtain an equation describing the Brunt-Väisälä frequency, N 2 , in response to the atmospheric forcing, given by a forcing function, F (t). Separating by timescale, we arrive at the following partial differential equations: F (t) is preceded by a minus sign for ease in derivation, as positive quantities of F (t) mean heat, hence buoyancy, is removed from the water column.
To simplify the seasonal time scale, we assume N 2 S is only a function of time, t: Assuming a homogeneous seasonal Brunt-Väisälä frequency over the depth of the water column gives us the relation for the seasonal Stratification Index, SI S , and seasonal atmospheric forcing: To simplify the analytical solution for the anomaly timescale, we describe the advection term, V · ∇(δN 2 ), as a restoring term, R = α(δN 2 ), which relates the overall Brunt-Väisälä frequency to its seasonal component, N 2 S , using the linear assump-270 tion made before. This means the restoring coefficient, α, represents the advective operation. This results in the following differential equation:

Seasonal solution and forcing
The solution for the seasonal timescale is relatively straight forward. As shown before, Eq. (9) relates the seasonal stratification,

275
SI S , to the seasonal atmospheric forcing, F S (t). We have the following definition of F S (t) from App. A1: where c p is the specific heat capacity of water, taken as 4184 Jkg −1 K −1 , g is gravity, ρ is the density of water, taken as 1000 kgm −3 , and T 0 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 q a,S = −Q net,S /D, where Q net,S is the 280 seasonal net downward heat flux at the ocean surface from Eq.
(2), we can calculate dSI S dt from Q net,S . If we integrate both sides of Eq. (9) by z, after plugging in Eq. (11) and the relationship for Q net,S , as SI S is constant with respect to (w.r.t.) z, Eq.
(9) becomes: g 2ρcpT0 ≈ 10 −9 m 4 /Js 2 , which means the derivative of SI S w.r.t. time, t, multiplied by 10 9 is on the same order of magnitude 285 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. 9 (a) for the 2013 winter, with dSI S dt × 10 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. 9 (b). Consequently, the maximum and minimum values 290 for SI S occur around Sept. 21st and March 21st, 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 W/m 2 to 200 W/m 2 , depending on the time of the year. Q LW remains roughly constant during the year, decreasing Q net by roughly -100 W/m 2 . These results are corroborated by the results of multiple model reanalysis for the region as well (Song and Yu, 2017).

295
Equation (12) 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, 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, as, 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 300 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.

Anomaly solution and forcing
To solve for the anomaly timescale, described by Eq. (10), we assume δF (t) can be represented by a pulse function shown in  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. (10) allows us to solve it in a piecewise manner. Like what we did for the seasonal time scale, we assume the water column has a homogeneous Brunt-Väisälä frequency, allowing us to make use of Eq.
(6). The restoring coefficient then only represents the horizontal advection, as the vertical component becomes zero with our where α d and α a are the restoring coefficients during ([t k , t k + ∆t k )) and after ([t k + ∆t k , t k + ∆τ k )) a Mistral event, respectively.
Further assuming δF k = δF , ∆t k = ∆t, and ∆τ k = ∆τ for all k, which results in a periodic pulse function with constant 315 amplitude and period, we can simplify Eq. (13) using the sum of a finite geometric series. At the beginning of the preconditioning period, destratification hasn't yet begun, therefore the initial δSI is zero, resulting in the following equation set: 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.

320
To determine the value of the restoring coefficients, a normalized function was derived for each section of a Mistral event (derivation shown in App. A2.1 for during an event and App. 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 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  Table 1 and δSI values from the NEMO results in Fig. 8. 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.

330
The strength of each Mistral event, δF k , was found in a similar way by solving for δF k after noting the initial value of Then 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 the appendix for the resulting values).

Mistral strength and destratification
Mistral events do not always lead to destratification. Some events in Fig. 8 The quantity ∂δSI k (t) ∂t must be less than zero for destratification to occur, which means if α d is a positive quantity (refer to 340 App. A2 or Fig. 11), δSI k−1 (t k ) + D 2 2 δF k α d must be a positive quantity. If some destratification has already occurred relative to the seasonal stratification, such that δSI k−1 (t k ) < 0, then D 2 2 δF k α d must be larger than −δSI k−1 (t k ) 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 even restratification occurs. An example of this can be seen with the Mistral event starting on Jan. 2nd, 2013, that 345 lasts for 17 days in Fig. 8 (b). The event starts off with an initial destratification of -0.48 m 2 /s 2 and ends at -0.41 m 2 /s 2 , a net restratification of 0.07 m 2 /s 2 . This is despite the fact this event has a positive δF k value of 3.80 × 10 −8 s −2 day −1 (from Table A1).
The combined overall effect of this result can be seen in Fig. 8 (b), as the consecutive Mistral events during the preconditioning phase cause destratification to a minimum of -0.6 m 2 /s 2 for δSI on Dec. 11th, 2012. Proceeding events after this 350 minimum fail to continue to destratify the water column and, instead, restratification occurs on the anomaly time scale, even before deep convection occurs. The seasonal stratification, SI S , and along with the anomaly destratification, δSI, brings the total SI to zero on Feb. 13, 2013, resulting in deep convection.

Dominating Mistral attribute
A pertinent question to ask is which attribute of the Mistral, the frequency, strength, or duration, is the most important when 355 it drives destratification. Figure 13 and 12 show the results of varying δF , ∆t, and ∆τ individually (in subplots (a), (b), and (c)), respectively) in Eq. (14). 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 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 360 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 ∂δSI ∂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 365 and the event strength, leading to destratification that almost reaches the same level as the case with more frequent events.
Increasing the frequency has a similar effect to increasing the duration; when the period is zero, the forcing becomes one large event, converging the resulting destratification to the dashed line.
To more accurately quantify the effect of each attribute, we separate δSI into its total derivative in terms of the Mistral attributes: 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 Fig. 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 notation to make is the character of the potential destratifi-375 cation limit: it approaches some asymptotic value as k approaches infinity. We can take advantage of this by differentiating the destratification phase of Eq. set (14) with respect to k, taking t = ∆τ , at the end of the phase, where the destratification equals the potential destratification: Figure 13. Same as Fig. 12, however, SIS is added to the results from Eq. (14).
Plugging in the mean values of ∆t, ∆τ , and ∆F , and taking k = 16, for the 16 events that occurred during the preconditioning phase, the above derivative equates a very small value of −5.93 × 10 −11 m 2 /s 2 per event. This confirms the small change in the potential destratification with increasing events. Taking k to infinity and noting that α d > α a results in the following: term is equal to −1.28 × 10 −1 m 2 /s 2 , the duration term has a value of −3.21 × 10 −2 m 2 /s 2 , and the period term has a value of 1.27 × 10 −2 m 2 /s 2 . With the strength term an order of magnitude larger than the other two terms, according to this simple model, the strength of the Mistral event is the most sensitive attribute when it comes to the effect of the Mistral on destratification, followed by its duration.

395
A complete and average Mistral destratification and restratification event according to Eq. (14) is given in Fig. 11 c, which took the average Mistral values from Table 1 and A1, and the restoring coefficients from App. 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. (13) with the duration and period information from    (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 in June rather than July was an arbitrary decision and is not believed to significantly affect the results or comparisons.

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. Fig. 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's too stratified to deep convect (0.22 m 2 /s 2 ), even though it's more destratified than the 1994 winter (0.36 m 2 /s 2 ). The MLD, however, clarifies that the 2005 winter experiences deep convection in the simulation results and the 1994 winter does not (not shown), which is consistent with other findings Herrmann et al., 2010;Somot et al., 2016).  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.

Anomaly forcing
The simple model for the anomaly scale was calculated for both the 1994 and 2005 winters, following the same steps in Sec.
4.3. The value of the restoration coefficients, α d and α a , were carried over from the 2013 winter analysis, with the Mistral 465 dates determined through the same process outlined in Sec. 2.3 and were manually adjusted to fit visual data (again with the first strength, then duration, followed by the length of the period. Only for the 1994 winter was the strength term found not to 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 three years, this aids the conclusion that, in general, the strength term of the

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

485
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 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 500 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 on a order of magnitude smaller than changing the Mistral strength.
Two additional years were also studied with the same method of running a Control simulation and a filtered atmospheric forcing Seasonal simulation: the year of 1993-1994 (1994 winter)  milder winters, the result was warmer and saltier WDMW production (Herrmann et al., 2010).

545
These questions are outside the scope of the current article, as they rely on investigating a large number of multiple years to evaluate the inter-annual variability of the atmospheric forcing. The 2013 winter featured 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; 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 of 36%, 61%, and 54% of the preconditioning days 550 having a Mistral event for the winters of 1994, 2005, and 2013, respectively. But 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 isn't clearly discernible with just 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 larger set of multiple years and we are preparing additional works to address them. We hope 555 these works 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): where ρ is density, c p is the specific heat of the fluid with constant pressure, T is temperature, t is time, and q is energy per volume from heat. D Dt is the material derivative. In this model we're assuming the heat transfer term is equal to the heat removed by the atmosphere: where q a is the volumetric heat forcing from the atmosphere. This leaves us with the following equation: The Brunt-Väisälä frequency is defined as: Assuming a fluid whose density varies negatively proportionally to the temperature, ρ = −βT , which is an acceptable 570 approximation as the density only varies only a few tenths of a kg/m 3 and temperature only varies about 10 degrees Celsius, we can describe N 2 in terms of temperature: Introducing buoyancy as b = − g T0 T , we can rearrange the energy equation into terms of buoyancy: If we then differentiate by ∂ ∂z , we can reorganize the equation in terms of the Brunt-Väisälä frequency, N 2 : 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: This brings us to: The main assumption we make is that the ocean column is a linear system and responds to the large and short/anomaly time scale atmospheric forcing independently: which describes the response of N 2 on the anomaly timescale, δN 2 : and the response of N 2 on the seasonal timescale, N 2 S : For the seasonal response, we further make the assumption that N 2 S negligibly depends on the x, y, and z coordinate direc-590 tions: If we want to connect the overall Brunt-Väisälä frequency, N 2 , to the seasonal one, N 2 S , we can formulate a restoring term, R, in terms of T , or in terms of N 2 following the steps mentioned above: Or, with δN 2 = N 2 − N 2 S : Where α is the restoring term coefficient. Separating the material derivative into its time and advective components for Eq. (A11): we will replace the advective component, V · ∇(δN 2 ), with R, which essentially swallows the advective operation into the restoring coefficient, α. This results in the partial differential equation that we will study further:

A1.1 Solution for seasonal SI
To solve the response of N 2 for the seasonal timescale, given by Eq. (A13), we will assume N 2 S is vertically homogeneous,

605
giving us the stratification index response, or Eq. (9), through the use of Eq. (6). We can then separate back out F S (t) into its components: Dividing q a,S by D gives us the atmospheric cooling in terms of a surface flux, −Q net,S . If we plug this relationship back into Eq. (9), we get: Integrating this equation by z gives us: And, therefore, we have SI S expressed in terms of Q net,S .

A1.2 Solution with Mistral forcing function 615
Focusing on just the anomaly time scale, 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 , 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, during and after, the advective components, hence α, 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 + ∆t k ), we get: With the result: After the event, [t k + ∆t k , t k+1 ): With the result: Or, to have the results more succinctly put: is a recursive initial condition, as its initial condition is the event before it, and so on: Therefore, δN 2 k−1 (t k ) can be simplified in expression by combining the initial conditions: If m = k and δN 2 0 = 0: Assuming δF k = δF , ∆t k = ∆t, and ∆τ k = ∆τ for all k, or a periodic pulse function, then δN 2 k−1 can be expressed as: Taking the sum of a finite geometric series: we then we get: where (α a − α d )∆t − α a ∆τ ̸ = 0. Plugging Eq. (A32) into Eq. (A22) and (A24) results in the equation for the response of the Brunt-Väisälä frequency forced by a periodic pulse function: For the anomaly response, Eq. (A25) and (A33), assuming a vertically homogeneous δN 2 leads to the stratification index through Eq. (6), leads us to δSI being expressed as: And: 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 + ∆t k ), with δSI given as: We first reference the time, t, to the starting time of event k as t ′ = t − t k , giving us: Next, we normalize δSI k to the value of zero at t ′ = 0, resulting in δSI k,N I : Then the height or magnitude of destratification for each event is normalized to 1, resulting in δSI k,N H : The extremum value for δSI k,N I (t ′ ) is when t ′ = ∆t k , or at the end of the Mistral event. This simplifies δSI k,N H to: Then, to normalize the length of the event duration from 0 to 1, we divide t ′ by the event length, ∆t k , which results in t ′′ : Plugging t ′′ into δSI k,N H (t ′ ) returns δSI k,N T : This final equation, δSI k,N T , 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 + ∆t k , t k + ∆τ k ), with δSI given as: Referencing the time, t, to the end of the event, t ′′′ = t − (t k + ∆t k ) and plugging t ′′′ into Eq. (A43), we get: Normalizing the vertical intercept of δSI k (t ′′′ ) results in δSI k,N I : δSI k,N I = δSI k (t ′′′ ) − δSI k (t ′′′ = 0) = δSI 2 k−1 (t k ) + 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 ′′′′ : which leads to δSI k,N T : This leaves us with an equation of α a , which can be fitted against NEMO model data, if ∆t k and ∆τ k are provided.
The average duration and period, ∆t and ∆τ , of events during the preconditioning period in Table 1 are used for ∆t and ∆τ in these normalized equations, Eq. (A42) and (A48). The result of the fitting is shown in Fig. 11, with α d having a fitted valued 695 of 0.235 day −1 and α a having a fitted value of 0.021 day −1 .

A3 Determining δF k
With the restoring coefficients determined in the prior section, Sec. 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 the value of δSI k−1 to be the same as δSI k (t k ) at the beginning of an event, we can simplify the equation in the following steps: 700 δSI k (t k ) = δSI k−1 δSI k (t k + ∆t k ) = δSI k−1 e −α d ∆t k + D 2 2 δF k α d e −α d ∆t k − 1 And then solve for δF k : The results of δF k for each event in the preconditioning phase is given in Table A1, along with mean and standard deviation. Taking the derivative with respect to time of Eq. (A25) and (A34) results in: −α a δN 2 k−1 (t k ) + δF k α d 1 − e α d ∆t k e (αa−α d )∆t k −αa(t−t k ) [t k + ∆t k , t k + ∆τ k ) A5 Asymptotic destratification 710 The following sections under Sec. A5 differentiate Eq. (A35) at t = t k + ∆t k , or at the end of a Mistral event, where the destratification is the largest, by k, and by the other components, δF , ∆t, and ∆τ , once k → ∞.

Appendix B: 1994 and 2005
The Mistral event data for the additional two years of 1993-1994 and 2004-2005 are found below in Tab. B1 and B2.
INSU-MISTRALS support and the Med-CORDEX program (COordinated Regional climate Downscaling EXperiment -Mediterranean region). It was also supported by a joint CNRS -Weizmann Institute of Science collaborative project. The authors acknowledge Météo-France for supplying the data and the HyMeX database teams (ESPRI/IPSL and SEDOO/Observatoire Midi-Pyrénées) for their help in accessing the data. Table B1. The start date of, duration of, ∆t k (days), and period between each event, ∆τ k (days), for each Mistral event, k, and event strength, δF k (s −2 days −1 ), for the preconditioning phase of the NEMO simulation period of Jun. 1st, 1993 to May 31st, 1994.