Rapid Restratification Processes Control Mixed Layer Turbulence and Phytoplankton Growth in a Deep Convection Region

The Gulf of Lion, Northwestern Mediterranean Sea, is one of few oceanic regions where deep convection occurs. We investigate the restratification following a convection event using measurements from an ocean glider equipped with turbulence microstructure sensors. This unique combination of instruments provides a high‐resolution description of the mixed layer with regard to turbulence, stratification and chlorophyll. We observe a rapid restratification process that proceeds over a timescale of days to one week. We find that restratification exerts a leading order control on surface mixed layer turbulence variability, as abrupt changes in turbulence dissipation rates are associated with the formation of near‐surface stratification. The near‐surface formation of stratification occurs through both the diurnal variability in surface buoyancy fluxes and through lateral advective processes. We conclude that daily near‐surface processes that influence stratification control mixed layer turbulence levels, and thus the phytoplankton response in the critical transition period to spring bloom.


Introduction
Deep convection (DC) is a water mass transformation process that causes ocean surface waters to sink to great depths [O(10 3 ) m] (e.g., Clarke & Gascard, 1983;MEDOC-Group, 1970), and this is mainly achieved by the interplay of the ocean circulation and atmospheric forcing, working together to weaken stratification in winter.DC primarily results in Deep Water formation, which plays a crucial role in regulating global climate by contributing to the thermohaline circulation through three-dimensional heat and salt transports (Marshall & Schott, 1999), and exchanging oxygen and carbon with the deeper ocean (Fröb et al., 2016;Ulses et al., 2021;Wolf et al., 2018).Few oceanic regions are prone to this process; we focus here on the Gulf of Lion, Northwestern Mediterranean Sea.After convection ceases, turbulent convective mixing gives way to a restratification process and, on a timescale of days to months, stratification is restored by surface buoyancy gain and lateral exchange between waters from the DC area and its surroundings (Marshall & Schott, 1999).Simultaneously, this newlyformed, near-surface stratification suppresses turbulence and vertical mixing, and together with a higher input of nutrients in the euphotic layer (d 'Ortenzio et al., 2014;Kessouri et al., 2017) brought by convection, creates conditions conducive for phytoplankton growth the following spring (Frajka-Williams & Rhines, 2010; Lavigne et al., 2013).However, stratification in the Mediterranean Sea is expected to increase due to increased heat and salt contents of its intermediate waters (Margirier et al., 2020) as a response to climate change (Adloff et al., 2015;Coma et al., 2009;Parras-Berrocal et al., 2022).In this context, and due to the significance of DC for global ocean circulation, it is essential to understand restratification and its interaction with ocean turbulence in DC regions.
Despite its relevance to upper ocean dynamics, vertically-resolved turbulence measurements in active DC regions have not, to our knowledge, been performed.This is likely due to the harsh conditions that make ship-based turbulent microstructure profiling difficult to impossible.For this reason many studies rely on models using turbulence parameterizations (André et al., 2005;Giordani et al., 2020;Pal & Chalamalla, 2020) to gain insight into the ocean's turbulent field.Here we overcome these limitations by performing turbulent microstructure measurements from an autonomous underwater glider (Fer et al., 2014;Palmer et al., 2015;Schultze et al., 2017) in the Gulf of Lion DC region.We show that strong variability is present in upper ocean turbulence that can be linked to the variable formation of diurnal stratification, and that relatively rapid time scales are found for the onset of this restratification in the near-surface.Finally, the shut down of turbulence by near-surface stratification provides ideal conditions for a nutrient-rich, newly-formed, shallow mixed layer to favor phytoplankton growth.

Study Area and Methods
Our measurements took place in the Gulf of Lion from 8 to 28 February 2021.They were performed with a Teledyne Webb Research Slocum Electric G2 autonomous underwater glider Comet, equipped with custom Sea-Bird Electronics conductivity, temperature, and depth sensors (Seabird SBE41 CTD) measuring at 1 Hz, as well as various navigation and communications sensors.The glider sampled in a sawtooth pattern at an angle of approximately 26°to the horizontal.It collected a total of 1,452 profiles between ≈2 and 200 m, with an average profile-to-profile horizontal distance of approximately 350 m.The glider surfaced roughly every 3 hr, for about 15 min, obtaining GPS positions and transferring decimated data sets to shore.Figure 1 shows the study area and the glider sampling locations.After the initial transit ( , Figure 1) to the DC region at around 42°N, 4.5°E, Comet was programmed to perform repeated transects around this location every 12 hr ( , DC area, Figure 1).
Each transect has approximately 30 profiles and a horizontal distance of 10 km.We assessed mainly temporal and vertical variations within the study area.A second glider (Dipsy) sampling to a depth of 1,000 m was also deployed and followed a similar path to Comet.Some of these data are used, and further details are found in Text S2 of Supporting Information S1.
In addition to the glider-based sensors above, a Rockland Scientific MicroRider-1,000 for the measurement of turbulent microstructure was mounted on the glider Comet.The turbulence data reported here were obtained from two air-foil shear probes (sampling at 512 Hz), a pressure transducer, electromagnetic current meter, two accelerometers, and an inclinometer measuring at 64 Hz.A time series of turbulent kinetic energy dissipation rate (ϵ, [W kg 1 ]) is obtained from the microstructure shear measurements by accounting for the speed of flow past the sensors, due to the sensitivity of ϵ on this parameter (Merckelbach & Carpenter, 2021;Merckelbach et al., 2019).The procedure used to arrive at ϵ generally follows that of Schultze et al. (2017), and is described in Text S1 of Supporting Information S1.The MicroRider started sampling on 9 February, was turned off in an effort to conserve energy between 18 and 21 February, recorded only downcasts between 21 and 24 February, and malfunctioned toward the end of the experiment, resulting in some data gaps.
To evaluate the surface fluxes of momentum (τ, [N m 2 ]) and buoyancy (B, [m 2 s 3 ]) we used the ERA5 reanalysis provided by the European Centre for Medium-Range Weather Forecasts (ECMWF, Hersbach et al., 2020).The surface heat flux (Q, [W m 2 ]) is calculated as the sum of the sensible and latent heat fluxes and solar radiation.The buoyancy flux is then found through the bulk formula in Grignon et al. (2010).As an additional tool to understand the role of surface fluxes in turbulence variability and stratification, we use the General Ocean Turbulence Model (GOTM, Umlauf & Burchard, 2005) to reproduce our observations in a onedimensional framework (see Text S3 in Supporting Information S1 for a detailed description and results).

The Deep Convection Area
The spatial extent of the Northwestern Mediterranean DC region is often quantified using satellite-derived surface chlorophyll-a (chl-a) levels (Herrmann et al., 2017;Houpert et al., 2016).The intense vertical mixing associated with convection results in reduced growth rates and a dilution of phytoplankton in the mixed layer, hence low surface chl-a visible from satellite imagery (the "blue hole" of Testor et al., 2018).We chose a surface chl-a threshold of 0.25 mg m 3 to define the horizontal extent of the DC area based on previously reported values during periods of large buoyancy losses (Herrmann et al., 2017;Testor et al., 2018).In Figures 1a and 1b chl-a values of 0.25 mg m 3 and smaller, contained by the 0.25 mg m 3 -contour, are found to lie in the center of the DC region (∼42°N, 4.5°E, Houpert et al., 2016), where the glider sampled.
Along the glider transit to the DC region ( , Figure 1), a gyre-scale structure can be seen: a steep, rising mixed layer base (Figures 2c and 2f) along the edge of the cyclonic circulation that corresponds to the strong Northern Current carrying surface fresh Atlantic Water (arrows in Figure 1).In this study, the mixed layer base is defined as the depth at which potential density changes 0.015 kg m 3 relative to 10-m reference depth, the latter chosen such that any contribution from solar heating is avoided.The active mixing layer is the portion of the water column where turbulence is being sustained by surface fluxes.It is found to follow closely the depth at which potential density changes 0.0025 kg m 3 relative to the surface (Figures 2c and 2d), as in Brainerd and Gregg (1995), which we use as a proxy.The observed gyre-scale circulation may result in upwards doming isopycnals that reduce stratification in the DC region, and allow for deep-reaching convection when strong winds blow and trigger large buoyancy losses at the surface (Marshall & Schott, 1999).This can partly be seen in the subsequent deepening of the mixed layer base, and a reduction in the measured stratification strength, as the glider entered the DC region around February 13 (see Figure 2c, and discussion below).We quantify the stratification strength by N 2 ≡ (g/ ρ 0 )∂ρ/∂z, where ρ denotes the potential density with ρ 0 a reference density, z the vertical coordinate (positive upwards), and g the gravitational acceleration.
The time the glider entered the DC region coincided with the period of strongest negative surface buoyancy fluxes during the measurement period (Figure 2b, with a cooling peaking at Q = 460 W m 2 ), resulting in convective turbulence that reached deeper than the 200 m diving depth of the glider on February 14 and 15 (Figure 2d).Nearby observations from the second glider Dipsy show a mixed layer as deep as 300 m (see Text S2 in Supporting Information S1).A significant fraction of the deepening of the mixed layer base can be accounted for by the loss in buoyancy during the convection event, with the remaining presumed to be attributed to existing lateral gradients (GOTM simulation, Text S3 in Supporting Information S1).After the event, the daily-averaged surface buoyancy flux switches to predominantly positive (heating) conditions, but with considerable diurnal variability in the non-averaged time series of surface buoyancy flux, B (gray line, Figure 2b).In addition, the wind stress magnitude in Figure 2a shows that, with the exception of two events on February 22 and 27, winds remained calm following the change in sign of the averaged buoyancy flux, thus allowing for the restratification phase of the DC region (Houpert et al., 2016).

Stratification and Turbulence Variability
The turbulence data show a strong variability that, due to the non-lagrangian nature of the glider measurements, corresponds to a mix of temporal and spatial components (Figure 2f).However, much of this variability can be understood through the interaction between the surface wind and buoyancy forcing, and the formation of diurnal warm layers (DWLs, Price et al., 1986).DWLs can be observed during periods of sufficiently weak winds and strong heating and are associated with the formation of a thin [O(1) m] layer of stratification (triangles, Figure 2c) that is dominated by temperature (Figure 3d).The presence of this stratification shows a rapid shutdown of the turbulence in the mixed layer below (Figures 2d and 3a) with high turbulence confined to the first few meters of the water column when the buoyancy flux is positive (Figures 3a and 3c).This shutdown implies that active mixing happens only in a very thin layer, although the mixed layer is as deep as 100 m (Figures 2c and 2d).
During the cooling phase the diurnal stratification can be mixed away, associated with a rapid increase in ϵ to depths greater than 100 m (e.g., February 10, 14).The 1-D turbulence analysis shows that the variability of turbulence in the mixed layer is predominantly responding to changes in surface forcing (Text S3 in Supporting Information S1).To help understand the rapid increases in the mixed layer turbulence at depth, we look at a characteristic time scale, T, for the penetration of convective turbulence throughout the depth of the mixed layer.This time scale takes the form of T ∼ H/w, where H is the mixed layer depth, and w a characteristic vertical velocity that is dependent on which terms represent the dominant balance in the convection.If buoyancy forcing is balanced by inertia we have w ∼ (BH) 1/3 , and T = O(1) hr, with B = 1 × 10 7 m 2 s 3 and H = 100 m (Marshall & Schott, 1999).On the other hand, if buoyancy forcing is balanced by Coriolis forces we find w ∼ (B/f ) 1/2 , which leads to a similar result of T = O(1) hr (Fernando et al., 1991;Marshall & Schott, 1999).We therefore expect that the mixed layer turbulence responds rapidly over a time scale of ∼1 hr to the mixing caused by convective turbulence following the loss of diurnal stratification, in agreement with the rapid increases in ϵ in Figure 2d.
Figures 3c and 3d shows the time evolution of averaged dissipation and stratification in cases when a clear DWL could be identified (February 9,10,11,13,17,and 23) within two different depth ranges: between 2 and 10 m and between 20 and 40 m, the latter representing the "remnant" layer immediately below the bottom boundary of the DWL (Brainerd & Gregg, 1995).We excluded the remaining DWLs (February 19, 20, 24-26, 28) due to the scarcity (or absence) of dissipation rates.With the onset of the diurnal stratification, the DWL bottom boundary isolates the upper ocean from the remaining mixed layer, creating two distinct turbulence regimes (Brainerd & Gregg, 1993a).In the upper 10 m, mean ϵ levels remain higher than 10 7 W kg 1 , slightly decreasing during heating conditions (red, Figure 3c), when stable stratification immediately forms (red, Figure 3d).In the deeper remnant layer, dissipation rates decrease one order of magnitude from the onset of the diurnal stratification until an hour after the peak of the heat flux (blue, Figure 3c).In other words, once isolated from the surface forcing by stratification, the averaged dissipation in the remnant layer shows a rapid decay with e-folding time scale of 140 min (R 2 = 0.92), and may be associated with a weak restratification occurring within the remnant layer, where increases are seen from near the beginning of the heating period (Figures 3b and 3d).In stratified waters, turbulence is expected to exponentially decay once the background stratification is constant (Brainerd & Gregg, 1993b;Callaghan et al., 2014;Smyth et al., 1997).Our O(1) hr decay shows a rapid turbulence response to the near-surface stratification formation.

Restratification: The Contribution of the Near-Surface
We have seen that upper ocean turbulence is highly variable in this region, and responds rapidly to varying nearsurface stratification over a time scale of roughly O(1) hr.The origins of this observed stratification are now explored and quantified using a buoyancy budget analysis.Recall that, once in the DC area, Comet performed repeated ∼10-km transects ( , Figure 1).Since these sections are much smaller compared to the DC area [O (100) km] and in its evolution [O(1) month], buoyancy variability should be dominated by temporal changes over spatial ones.The Rossby radius being around 10 km in the region, it is however not excluded that mesoscale processes contribute to this variability at smaller temporal scales.
Following the convection event of February 14, a rapid restratification takes place as the daily-mean heat gain becomes positive (Figure 2b).This restratification begins soon after February 16, when the buoyancy profile shows a completely mixed water column in the upper 200 m (①, Figure 4a; see also Figure S2 in Supporting Information S1).We use the buoyancy level of this fully mixed profile as a convenient reference (b mix ) for the mixed state from which buoyancy anomalies b′(z, t) ≡ b(z, t) b mix are calculated.Figures 4a and 4b shows the evolution of the restratification phase through the buoyancy anomaly b′(z, t) as a function of depth and time.
Using b mix , we partition the restratification into positive and negative buoyancy anomalies that we associate with deep and shallow restratification components, respectively.Figures 4a and 4b shows that the buoyancy restratification proceeds over a time scale of T = O(1) week or less, and it is composed of both deep and shallow anomalies.We define the change in buoyancy content of a portion of the water column ( where z b and z t are variable depth limits that yield the partial contribution of different layers to the total restratification.We evaluate the absolute buoyancy anomalies so that both deep and shallow restratification would have a positive contribution to the total. It is helpful to look at various values for the top and bottom integration depths, z t and z b , respectively, and these are plotted in Figure 4c.Namely, the total buoyancy content anomaly ( btot ) results from integration over the whole measured water column, the deep ( bdeep ) and shallow ( bshal ) components come from integrating from the bottom and top to the level z mix where the mixed buoyancy is found, that is, b′(z mix , t) = b mix (red line, Figure 4b).The near-surface contribution we define as bns , the buoyancy anomaly content in the upper 20 m.This partitioning of the buoyancy content results in the sum of the deep and shallow contributions giving the total, that is, btot = bdeep + bshal .Additionally, we can include the contribution of the surface fluxes through where t 0 represents the beginning of the restratification period on 16 February (Figure 4c). Figure 4c shows that from February 16 until the formation of near-surface restratification on February 19, the integrated surface buoyancy flux is negative ( B, Figure 4c) and restratification is mainly through lateral advection, with inputs of buoyancy over the 0-200 m layer (Figures 4b and 4c).By the end of the restratification period shown, the nearsurface (upper 20 m) component accounts for around 19% of the total stratification, with 59% of that value due to the surface buoyancy flux ( B).Interestingly, more than half (60%) of the total buoyancy anomaly content of a 200 m-deep water column is found roughly in the first 100 m, showing the significance of shallow water processes and their interplay with atmospheric forcing.
A remarkable feature of the evolution of the buoyancy content in Figure 4 are the relatively rapid time scales that are present in restratifying the water column after the February 14-15 convection event, with both deep and shallow increases in buoyancy occurring over days to O(1) week.A time scale for the baroclinic mesoscale eddydriven restratification of convection regions has been proposed by Jones and Marshall (1997) as T restrat ≈ 56 r/ (Δbh) 1/2 , with r a radial length scale for a patch of dense fluid with buoyancy difference Δb from the ambient fluid, and a depth h.Taking h = 300 m, Δb = 3 ⋅ 10 4 m s 2 representative of the deep buoyancy anomaly, and r ≈ 50 km (see Figure 1), this results in T restrat = O(100) days.Although this time scale is representative of the full restratification of the region through mesoscale eddies, there are also found to be faster restratification mechanisms operating (e.g., Clemént et al., 2023;Houpert et al., 2016;Testor & Gascard, 2006;Testor et al., 2018), such as the advection of stratified waters by wind-driven Ekman (Bosse et al., 2021) and submesoscale (Boccaletti et al., 2006;Clemént et al., 2023;Testor & Gascard, 2006) processes, and through lateral changes in temperature and salinity being converted to vertical changes (Brainerd & Gregg, 1993b;Caldwell et al., 1997).Using a baroclinic instability model, Testor and Gascard (2006) found a timescale of ∼6 days for the restratification.Couvelard et al. (2015) showed large differences in restratification time scales depending on a model's spatial resolution, hinting at the importance of smaller scale processes; with a 2 km resolution, restratification could occur in less than 10 days.In the buoyancy analysis of Figure 4 we point to the importance of lateral advective processes in both the shallow and near-surface restratification observed.Only roughly 60% of the near-surface buoyancy gain ( bns ) can be accounted for by surface fluxes, with the remaining 40% accounted for presumably through lateral advective processes (see additionally Figure S3 in Supporting Information S1).This fraction becomes much larger (≈80%) for the shallow buoyancy gain, bshal , and has important links to the growth of phytoplankton from restratification processes.

Implications for Phytoplankton Growth
In a classical framework, phytoplankton is expected to bloom when the mixed layer depth shoals beyond a critical value, due to changes in the surface heat flux in spring from cooling to heating (Franks, 2015;Sverdrup, 1953).However, high phytoplankton concentrations may appear before spring, related to a higher input of nutrients from convection events (Figure 1, Kessouri et al., 2018;Severin et al., 2014), the decrease in turbulence levels in the mixed layer (Franks, 2015;Taylor & Ferrari, 2011) and also the feedbacks between phytoplankton and its predators (Behrenfeld & Boss, 2014).The rapid shoaling of the mixed layer and the restratification of the nearsurface that we observe following the convection event (Figures 2 and 4) traps the turbulence in a very thin (<10 m) layer, preventing mixing and allowing phytoplankton to grow (e.g., Kessouri et al., 2018;Ulses et al., 2016).Figure 2e shows the variations in chl-a concentration throughout our time series.Since the sensor was not calibrated prior to deployment, the chl-a estimates should be interpreted with caution, and primarily in qualitative terms.Quenching effects were corrected following Thomalla et al. (2018).One can readily observe higher values of chl-a associated with stratified conditions and lower chl-a levels associated with a mixed state.In particular, it can be seen that increases in chl-a occur when near-surface stratified layers have formed during restratification (e.g., February 19, 25, 26 in Figures 2c and 2e).These layers have turbulence that is also confined only to the near-surface layer, and are unable to mix plankton throughout the generally deeper mixed layer.This situation may be contrasted to the two wind events on February 22 and 27 where chl-a concentrations are found to be homogeneous throughout the mixed layer.
The formation of near-surface stratification largely shuts down turbulence throughout the water column, leading to the growth of phytoplankton within the euphotic zone (i.e., not confined solely to the stratified region near the surface) over a time scale of days.This is a similar mechanism to that suggested by Taylor and Ferrari (2011), who connected the convective forcing associated with wintertime cooling to turbulent mixing, and hypothesized that phytoplankton blooms can be triggered when the net heat flux switches from cooling to heating.We find this switch to occur through the rapid formation of near-surface stratification and its control on upper ocean turbulence.Through such a mechanism, phytoplankton can be trapped in the deeper stratified ocean, as is seen in Figure 2e below the mixed layer base.This could represent a carbon export mechanism associated with restratification (e.g., Lacour et al., 2019).

Summary and Conclusions
Turbulence measurements were collected in the DC region of the Gulf of Lion using an autonomous ocean glider equipped with a turbulent microstructure package.These unique observations provide a detailed characterization of the mixed layer during a 300-m deep convective event, followed by the subsequent restratification of the DC region.The variability in upper ocean turbulence over timescales of hours to days is strongly linked to both changes in surface forcing, and near-surface restratification processes.The near-surface stratification confines the turbulence in the upper meters of the water column, and isolates deeper layers from surface-driven turbulence.
Together with the relatively rapid growth and decay timescales of the turbulence [O(1) hr], these processes are capable of producing the strong variability in upper ocean turbulence in the convection region.Restratification of the upper 200 m of the water column occurs over a timescale of O(1) week through both the diurnal variability in surface buoyancy fluxes and through lateral advective processes.The formation of near-surface stratification and its effects on the turbulence distribution are also found to allow for an effective growth of phytoplankton.Thus we find that complex, and rapidly varying, near-surface processes are crucial to understand restratification and its biophysical feedbacks in the DC region.

Figure 1 .
Figure 1.Chlorophyll-a (chl-a) maps from Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua L3, 4 km, daily products in the Gulf of Lion averaged from (a) 1 to 31 January, (b) 14 and 17 February and (c) 20 and 27 February and the trajectory of the glider showing the transit (yellow) and the deep convection center (green) areas.While at , the glider transected until 28 February.Geostrophic currents (vectors) from Copernicus Marine Environment Monitoring Service (CMEMS, E.U.Copernicus Marine Service Information) are averaged between 8 and 28 February 2021.White contours represent chl-a = 0.25 mg m 3 .Gaps in chl-a were linearly interpolated.

Figure 2 .
Figure 2. (a) Wind stress (τ, black line) and wind magnitude at 10 m above surface (gray vectors), (b) hourly (gray) and daily averaged (blue, red lines) surface buoyancy fluxes (B) from ERA5, (c) squared buoyancy frequency (N 2 ), (d) turbulent kinetic energy dissipation rates (ϵ), and (e) chlorophyll-a (chl-a, uncalibrated) from observations during transit and in the Deep convection (DC) area .Please note the non-linearity of the vertical axis in panels (c) to (f).Orange vertical line marks the transition between transit and the DC area sampling.Black and gold lines in (c) and (d) shows the mixed layer and active mixing layer depths, respectively.(f) Schematic of the main features of interest.Triangles represent occurrences of diurnal warm layers.

Figure 3 .
Figure 3. (a) Turbulence and (b) stratification of the diurnal warm layer on February 17 and the active mixing layer depth (gold) are shown as an example.Time evolution of (c) dissipation rates and (d) stratification, averaged between 2 and 10 m (red), and between 20 and 40 m (blue) from February 9, 10, 11, 13, 17, and 23.In (d), stratification is partitioned into its thermal (dot-dashed) and haline (dotted) contributions, highlighting the dominance of temperature in stratification in the upper 10 m.Vertical lines mark the change in sign of the surface heat flux, reaching its maximum at noon.Hatched areas correspond to the 95% confidence interval around the mean values.Black thick line in (c) is the linear regression fit estimate of the exponential decay of ϵ, that is, ϵ ∼ e ( t/γ) , within the 20 40 m layer.

Figure 4 .
Figure 4. (a) Profiles of buoyancy anomalies at ① midnight February 16, ② midnight February 19, and ③ midday February 27, (b) buoyancy anomaly in the upper 200 m and (c) the time series of the integrated buoyancy in different layers.Red curve in (b) represents the depth (z mix ) where the buoyancy equals b mix .The definitions of bns , bshal , bdeep , and btot are drawn in (a) for density profile-③.B is defined in Equation 2.
Funding for this work has been provided through the Helmholtz Association PoF-IV program, as well as being a contribution to project T2 (Ocean Surface Layer Energetics) of the Collaborative Research Center TRR 181 "Energy transfers in Atmosphere and Ocean" funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Projektnummer 274762653.The glider data analysis received support from GROOM-RI funded by the European Union's Horizon 2020 research and innovation program under Grant agreement no.951842.Open Access funding enabled and organized by Projekt DEAL.