Inhibited vertical mixing and seasonal persistence of a thin cyanobacterial layer in a stratified lake

Harmful blooms of the filamentous cyanobacteria Planktothrix rubescens have become common in many lakes as they have recovered from eutrophication over the last decades. These cyanobacteria, capable of regulating their vertical position, often flourish at the thermocline to form a deep chlorophyll maximum. In Lake Zurich (Switzerland), they accumulate during stratified season (May–October) as a persistent metalimnetic thin layer (~2 m wide). This study investigated the role of turbulent mixing in springtime layer formation, its persistence over the summer, and its breakdown in autumn. We characterised seasonal variation of turbulence in Lake Zurich with four surveys conducted in April, July and October of 2018 and September of 2019. Surveys included microstructure profiles and high-resolution mooring measurements. In July and October, the thin layer occurred within a strong thermocline (N≳0.05\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \gtrsim 0.05$$\end{document} s-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}) and withstood significant turbulence, observed as turbulent kinetic energy dissipation rates (ε≈10-8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \approx 10^{-8}$$\end{document} W kg-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}). Vertical turbulent overturns –monitored by the Thorpe scale– went mostly undetected and on average fell below those estimated by the Ozmidov scale (LO≈1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_O \approx 1$$\end{document} cm). Consistently, vertical diffusivity was close to molecular values, indicating negligible turbulent fluxes. This reduced metalimnetic mixing explains the persistence of the thin layer, which disappears with the deepening of the surface mixed layer in autumn. Bi-weekly temperature profiles in 2018 and a nighttime microstructure sampling in September 2019 showed that nighttime convection serves as the main mechanism driving the breakdown of the cyanobacterial layer in autumn. These results highlight the importance of light winds and convective mixing in the seasonal cycling of P. rubescens communities within a strongly stratified medium-sized lake.


Introduction
Temperate lakes are density-stratified throughout most of their annual cycle (Boehrer and Schultze 2008). During the stratified period, the water column typically assumes a vertical density structure consisting of three layers. These include a warmer surface layer (epilimnion) and a colder, weakly stratified deep layer (hypolimnion) separated by a region of continuous density stratification (metalimnion or thermocline). Stable conditions and the balance between nutrient and light availability, which exhibit opposing vertical distributions, make the thermocline an ecological niche for certain planktonic microbial organisms (Abbott et al. 1984;Sharples et al. 2001). These conditions often give rise to a sub-surface chlorophyll maximum of variable thickness. These may span less than a few meters and are termed thin Aquatic Sciences * Bieito Fernández Castro bieito.fernandezcastro@epfl.ch; b.fernandez-castro@soton.ac.uk layers in oceanographic literature (Dekshenieks et al. 2001;Sullivan et al. 2010;Durham and Stocker 2012). Thin layers have received growing attention in the scientific community due to their complex roles as ecological niches (Benoit-Bird et al. 2009), and because they often host harmful microbial communities (Nielsen et al. 1990;Broullón et al. 2020). Efforts to reestablish oligotrophic conditions in some lakes has led to the proliferation of the harmful filamentous cyanobacterium Planktothrix rubescens, an unintended side-effect which poses a risk to human health (Jacquet et al. 2005;Ernst et al. 2009;Marti et al. 2016). In Lake Zurich (Switzerland), P. rubescens has become the dominant phytoplanktonic species over the last 40 years as the lake has recovered from strong eutrophication during the 1970s (Posch et al. 2012). Planktothrix rubescens' adaptation to low light and their capacity to control their vertical position via buoyancy regulation with gas vesicles (Walsby and Schanz 2002;Walsby 2005) allow them to proliferate in a metalimnetic thin layer that persists throughout the entire stratified season (May-October; Posch et al. 2012). This layer differs from other deep chlorophyll maxima in lakes due to the former's reduced thickness ( < 2 m, Leach et al. 2018) and from similar structures in the ocean due to its temporal persistence (~6 months). Marine thin layers typically last for a few hours to a few days (Dekshenieks et al. 2001;Durham and Stocker 2012), by which time advection or turbulent mixing typically disrupt the layers Shroyer et al. 2014).
The few turbulence measurements available for marine thin layers indicate that turbulent mixing plays a critical role in regulating their formation and persistence (McManus et al. 2003;Cheriton et al. 2009;Steinbuck et al. 2009). Due to the scarcity of turbulence measurements in the presence of deep chlorophyll maxima in lakes however, the behaviour and scales of lacustrine thin layers is often inferred only from stratification (Leach et al. 2018). In Lake Zurich, for example, long-term persistence of P. rubescens is attributed to the strengthening and deepening of the metalimnion and to the increased duration of seasonal stratification (Yankova et al. 2016). Both of these sets of factors arise from warming experienced over the last decades (Livingstone 2003;Schmid and Köster 2016). However, different stratification regimes alone cannot explain the high abundance of P. rubescens in certain lakes (Jacquet et al. 2005). An adequate characterisation of turbulent mixing and its influence on P. rubescens communities can further constrain physical factors influencing their residence in Lake Zurich and similar systems.
Mixing in the metalimnion of lakes depends on the dissipation of internal waves (Lorke et al. 2005;Preusse et al. 2010) and turbulence originating at the surface, either by wind-driven shear (Imberger 1985;Pernica et al. 2014) or cooling-induced convection (Serra et al. 2007;Tedford et al. 2014). Previous turbulence measurements in the interior of lakes during the stratified season revealed weak levels of internal mixing away from lake boundaries (Saggio and Imberger 2001;Etemad-Shahidi and Imberger 2001). This factor may explain the long-lasting layer persistence but turbulence is highly episodic and modulated by highly variable wind forcing (Saggio and Imberger 2001). Seasonal and higher-frequency variability in turbulent mixing may regulate formation of the P. rubescens thin layer in spring, its persistence over summer and breakdown in autumn, however, evidence for these mechanisms remains lacking.
To fill this gap in understanding, we conducted four field surveys in Lake Zurich in April, July and October of 2018, and in September of 2019. Surveys included microstructure profiles and high-resolution measurements of velocity and temperature fluctuations to characterise the seasonal and short-term variability of turbulent mixing in the metalimnion and its influence on the P. rubescens thin layer. Biophysical data collected during bi-weekly monitoring was used to interpret the high-resolution turbulence measurements. Comparison of Lake Zurich results with those from a similar lake (Lake Geneva) identified key physical factors that contribute to the proliferation of P. rubescens in some lakes but not in others.

Study site and sampling overview
Lower Lake Zurich is a long (28 km), narrow (3.8 km), deep (max. depth 136 m), and medium-sized (65 km 2 ) perialpine lake. The lake is monomictic meaning that it experiences one seasonal deep convective mixing occasionally reaching the lake bottom in February-March, and subsequent thermal stratification during summer. Lake Zurich is exposed to light winds mainly from the northwest (mean ± SD of 2.6 ± 1.8 m s −1 for 2018), flowing along the main lake axis (Baracchini et al. 2019;Fig. S1).
Three surveys were carried out in order to investigate the seasonal variations in mixing conditions. Surveys occurred on 24-26 April, 17-19 July and 2-4 October 2018 at station st1 (120 m local water depth), located in Lower Lake Zurich, ~3 km from its deepest location (Fig. S1). During these surveys, we collected profile data with microstructure probes (Table S1), while a mooring system collected high-resolution current and temperature time-series within the thermocline (Table S2). We used hydrographic and biological data collected at the deepest point of the lake (st0, Fig. S1) as part of a long-term monitoring program to facilitate data interpretation over a broader seasonal context. An additional microstructure survey was conducted on 24-25 September 2019 to assess the role of nighttime convection for the seasonal breakdown of the P. rubescens layer.

Microstructure turbulence measurements
Microstructure profiles were collected with three different loose-tethered microstructure profilers deployed at st1 (2018) and st0 (2019). These included a SCAMP (Self-Contained Autonomous Microstructure Profiler, Precision Measurement Engineering) deployed in April 2018, a RSI VMP-500 (Vertical Microstructure Profiler, Rockland Scientific International) deployed in July and October 2018, and a RSI microCTD deployed in September 2019 (Table S1).
The three instruments were equipped with two FP07 thermistors that measure temperature fluctuations at a resolution of 10 −5 • C and with a response time of ≈ 10 ms (Sommer et al. 2013). The FP07 thermistors sampled at a frequency of 100 Hz (SCAMP) and 512 Hz (VMP and microCTD). Additionally, lower-resolution high-accuracy temperatureconductivity-depth (CTD) and fluorescence sensors were mounted on all three microstructure profilers. Chlorophylla (chl-a) concentrations were derived from the latter. Lake water salinity and density ( ) were calculated from the in situ temperature and conductivity measurements using lakespecific equations (Wüest et al. 1996). The buoyancy of the instruments was regulated to achieve a free-falling profiling speed of 10 to 20 cm s −1 (SCAMP and VMP) and 30 to 40 cm s −1 (microCTD). SCAMP and VMP were operated in free-falling mode, while the microCTD was operated alternatively in free-falling and free-rising mode, between the surface and 30 m depth.
The dissipation rate of thermal variance was calculated as where T = 1.4 × 10 −7 m 2 s −1 is the molecular diffusivity of heat. The Batchelor wavenumber ( k B ) was calculated from fitting the power spectra from the vertical temperature gradient data to the Kraichnan spectrum (Kraichnan 1968) using the maximum likelihood method, which explicitly considers the noise curve of the thermistors (Ruddick et al. 2000). The rate of turbulent kinetic energy (TKE) dissipation ( ) was determined from the Batchelor wavenumber, . These quantities were calculated over partially overlapping 1 m segments spaced at 25 cm. Bad fits were discarded according to the misfit criteria proposed in Ruddick et al. (2000). The estimation of also requires consideration thermistor spatial resolution vis a vis the sensor's time-response. Although this effect was corrected with a single-pole function (Nash et al. 1999), the results can be highly uncertain for high values, when temperature fluctuations are expected at large wavenumbers. For a single pole correction and frequencies higher than the inverse of the time response f = 1∕ [Hz], the spectral correction exceeds ~ 60%. Taking into account the speed of descent for different instruments and the sensor time-response, f corresponds to wavenumbers of ∼670 cpm and ∼280 cpm, for SCAMP/VMP and microCTD probes, respectively. This would equate to ≈ 6 × 10 −6 W kg −1 and ∼2 ×10 −7 W kg −1 . We interpret these values to represent safe upper limits for the detection of dissipation rates with the SCAMP/VMP and microCTD sensors. Cable strumming can cause perturbations of the descending profiler. These can result in extreme overestimation of dissipation rates during VMP sampling. Spikes interpreted to be cable strumming were removed. Similarly, due to limitations in microCTD performance during ascent when operated from a small boat under moderate winds, lower parts of profiles were often contaminated. These values were also removed. Water column stability was characterised with the buoyancy frequency, N, defined as N 2 = −g∕ ( ∕ z) , where g is gravitational acceleration. The vertical turbulent heat diffusivity was calculated using the Osborn and Cox (1972) model given by: This equation expresses the fact that turbulence stirs the large-scale temperature gradient and thereby transfers variance to small scales where it is dissipated at a rate . This boosts the action of molecular diffusion and increases the rate of mixing.
The vertical scale of turbulent overturns was estimated using the Thorpe length-scale, L T = � ⟨d � 2 ⟩ , where d ′ are the Thorpe displacements, i.e. the distance that a fluid parcel needs to be displaced in order to obtain a monotonic, stable density profile (Thorpe 1977;Lorke and Wüest 2002). Other relevant length-scales for turbulence are the Ozmidov ( L O = ( N −3 ) 1∕2 ) a n d t h e Ko l m o g o r ov s c a l e s ( L K = ( 3 −1 ) 1∕4 ), at which inertial forces enter into equilibrium with buoyancy and viscosity, respectively. The ratios of these length-scales, i.e., buoyancy Reynolds number ( Re b = (L O ∕L K ) 4∕3 = −1 N −2 ) and the turbulent Froude number ( Fr T = (L O ∕L T ) 2∕3 ), are used to characterise the balance of forces and development of turbulence, and are related to the mixing efficiency (Ivey and Imberger 1991).

High-resolution mooring
A high-resolution (HR) mooring (Sepúlveda Steiner et al. 2021; see Fig. S1) was deployed within the thermocline (~12 m depth) at station st1 to characterise velocity and temperature fluctuations within the thin layer and its surroundings (Table S2). The mooring consisted of two fins. Five RBR temperature loggers sampling at 2 Hz were installed at different depths along the upper fin. The second fin held a Nortek Aquadopp HR profiler (2 MHz) installed facing upwards for measuring current velocities at 1 Hz. The Aquadopp profiled 2.06 m divided into 104 bins of 0.02 m size along beams in a 25 • orientation relative to vertical that covered the approximate vertical range of the temperature measurements. The HR-mooring was deployed for periods of several hours during several days in April, July and October 2018 close to st1 (Fig. S1). Due to operational issues such as vibrations along the mooring, some deployments generated only poor data. For this research, we analysed the July and October datasets that contained data of sufficient quality (Table S2). Time-mean ENU (eastward-northward-upward or u-v-w) velocities were extracted by low-pass filtering at a cutoff period of 5 mins, which corresponds to about two times the local buoyancy period ( 2 ∕N ), and by rotating the measured along-beam current velocities into ENU coordinates. TKE dissipation rates ( ) were estimated using the structure function method directly applied to along-beam velocities as described in Wiles et al. (2006). For this calculation we used 1024 point half-overlapping segments, corresponding roughly to 17 mins.

FluoroProbe and CTD profiles
As part of a long-term monitoring program, bi-weekly profiles ( n = 25 ) of temperature, P. rubescens chl-a and other variables were measured near the deepest location of the lake at st0 (Fig. S1) between 10 January and 12 December 2018 with a TS-16-12 FluoroProbe (bbe Moldaenke GmbH, Kronshagen, Germany). Additionally, during the convectionoriented sampling in September 2019, 11 FluoroProbe casts were performed between 14:00 on the 24th and 10:40 on the 25th of September 2019.

Meteorological data
Hourly values for meteorological variables were obtained from the COSMO-1 (COnsortium of Small-scale MOdeling) weather model run by the Swiss Federal Office of Meteorology and Climatology (MeteoSwiss) covering the Swiss territory with a horizontal resolution of 1.1 km. Data was extracted for the grid-bin closest to st0. Heat fluxes through the lake surface (including shortwave and longwave radiation, as well as latent and sensible fluxes) were calculated following the derivation of Fink et al. (2014). Hourly surface water temperatures required for this calculation were interpolated from observations.

Analysis of basin-scale stability and mixing
Water column stability was quantified as the negative value of the background potential energy, i.e. the amount of mechanical energy that would be required to remove the existing stratification and fully mix the water column: where is the volume-averaged water density, A(z) is the lake area at depth z (where z is defined as positive upwards and zero at the surface), A 0 = A(0) is the surface area, z max is the maximum depth of the lake and z c is the depth of the lake's volume centroid. Lake bathymetry was obtained through the Swiss Federal Office of Topography (https:// www. swiss topo. admin. ch/).
The potential of wind forcing to overcome the restoring effect of stratification, disturb the metalimnion and produce diapycnal mixing was quantified as the lake number, L N (Imberger and Patterson 1989;Robertson and Imberger 1994): where z TH is the thermocline depth and u * = ( ∕ ) 1∕2 is the friction velocity. Here, = air C 10 W 2 10 is the wind stress, where W 10 is the wind speed 10 m above the lake surface, C 10 is the wind drag coefficient, calculated according to Wüest and Lorke (2003), and air is the air density. Bi-weekly values of Sc were interpolated into hourly values to match the resolution of wind data for this calculation. The relative importance of sub-surface convective and wind mixing was quantified at 1-h intervals with the Monin-Obukhov lengthscale, L MO = u 3 * ∕( B 0 ) , where = 0.41 is the von Karman constant and B 0 is the surface buoyancy flux. Basin-scale vertical heat diffusivity ( K BS T ) was calculated with the heatbudget method (Powell and Jassby 1974) for hypolimnetic layers, and using a shear and convective velocity scaling for the surface boundary layer. Supplementary Text 1 describes these methods in detail. Figure 1 shows the evolution of atmospheric forcing and lake temperature at the monitoring station (st0) for Lake Zurich in 2018. From January to March, the lake lost heat to the atmosphere (Fig. 1a). The water column was weakly stratified in January and the surface mixed layer (SML, defined as the

Seasonal stratification and P. rubescens thin layer
depth where temperature was 0.5 °C lower than the near-surface value) deepened progressively. This continued until full mixing by mid-March (Fig. 1b) following a 10-day period of significant heat loss and relatively strong winds. The holomixis period was short with solar radiation dominating the heat budget after mid-March. As a result, a weak shallow thermocline had already developed by early April. Thereafter, surface temperature and water column stratification progressively increased to 26.5 °C, as observed on 8 August.
During this warming period, a strong thermocline, with a temperature gradient ranging between 1.0 and 3.0 °C m −1 , occurred just below the SML and showed a deepening trend from 2 to 10 m depth. Net cooling prevailed from August onwards and SML temperatures progressively decreased. The SML deepened steadily from 10 to 20 m depth in December, but the maximum thermal stratification remained strong (2-3 °C m −1 ) until the end of September. It then declined Fig. 1 Seasonal time-series of atmospheric and hydrographic conditions, and P. rubescens chl-a concentration in Lake Zurich during 2018. a Hourly (red-blue filling) and weekly mean (black) net heat flux between the lake and the atmosphere (positive upwards) and weekly mean wind speed (dashed) derived from the meteorological COSMO model; b vertical distribution of water temperature; and c vertical distribution of P. rubescens chl-a concentration (μg L −1 ) (contours) and maximum concentration (bars). In b, c, the black-white thick dashed line represents the depth of thermocline (maximum temperature gradient) and the cyan-red tick dashed line, the surface mixed-layer (SML) depth. In c, the dotted thin lines rep-resent the daily PAR fluxes resulting in P. rubescens neutral buoyancy (0.28 mol photon m −2 d −1 ) and zero net growth (0.36 mol photon m −2 d −1 ) (Walsby et al. 2004). Above the upper panel, vertical lines represent the sampling dates and the color shading represents the different stages of the seasonal evolution of P. rubescens thin layer from growth and settling during spring (blue: 20 March to 30 May), maximum stable concentration in early-mid summer (orange: 17 June to 8 August) to decline in late-summer and autumn (green: 21 August-31 October). Turbulence-oriented surveys are highlighted with darker shading (color figure online) gradually as the SML expanded through the deeper part of the metalimnion, which was more weakly stratified. Planktothrix rubescens was observed throughout the SML as relatively low chl-a concentrations (chla < 4 μg L −1 ) during the mixing period in winter and occurred throughout the water column after holomixis (Fig. 1c). From May onwards, P. rubescens communities formed a thin layer located within the main thermocline between 10 and 17 m depth. Planktothrix rubescens regulate their vertical position according to the daily dosages of photosynthetically available radiation (PAR; Walsby et al. 2004). The records described here found P. rubescens mostly located between daily PAR isolines for which P. rubescens acquires neutral buoyancy (0.28 mol photon m −2 d −1 ) and for which the net growth becomes zero (0.036 mol photon m −2 d −1 ; thin dashed lines in Fig. 1c; Walsby et al. 2004). The layer remained stable, with peak concentrations between 30 and 40 μg L −1 , until October, when the thermocline deepened and approached the layer. From mid-October, P. rubescens concentration in the SML began to increase but the peak persisted until the end of the month, when the thermocline transgressed the layer, and the latter dispersed.
For the rest of the year, the P. rubescens population was located exclusively in the SML at moderate concentrations of ~ 10 μg L −1 .

Seasonal changes in stratification and mixing environment
Our turbulence sampling detected three stages in the seasonal development of the P. rubescens layer. The first microstructure and HR-mooring survey was carried out between 24 and 26 April 2018 (blue, Fig. 1), shortly before formation of the thin layer. The second survey occurred between 17 and 19 July (orange), when the thin layer had been established for more than 2 months and water column stability approached its maximum value. Finally, the third sampling took place between 2 and 4 October (green), when some P. rubescens chl-a had been detected in the SML but the chl-a peak remained sharp and well defined in spite of its reduced amplitude.
During the sampling in April, heating outpaced nighttime cooling and hourly wind speeds were maximal (6 m s −1 ) on the 25th, such that no microstructure profiles could be In a-c, microstructure profiles are indi-cated with black ticks, and the duration of mooring deployments is indicated with color filling. In panels c-f the depth range of the P. rubescens thin layer, defined as the range where chl-a is larger than an e-folding fraction of the maximum value (chl-a > max (chl-a)/e, is indicated by filled rectangles for July (orange) and October (green). Individual profiles for chl-a are shown in f as semi-transparent lines. Individual profiles were projected onto isothermal coordinates to remove the effect of internal wave motions (color figure online) Page 7 of 22 38 carried out on that day. Except for this day, wind remained < 4 m s −1 (Fig. 2a). Thermal stratification was still moderate ( N max ≈ 0.02 s −1 ), with a shallow thermocline located above 5 m depth (Fig. 2d, e). Chl-a measured with the SCAMP was low ( < 5 μg L −1 ) and no thin layer was observed (Fig. 2f). In July, nighttime cooling intensified but daytime heating still dominated the net heat budget. Maximum wind speeds of > 5 m s −1 were observed on the 19th (Fig. 2b). The water column was thermally stratified from the surface ( N ≈ 0.02 s −1 ), and a strong thermocline extended from 7 to 17 m depth with maximum stratification ( N max = 0.075 s −1 ) at 9 m depth (Fig. 2d, e). A ~2 m thick, vertically symmetric thin layer of chl-a (max. 30 μg L −1 ) occurred at 13 m, 4 m below the stratification maximum (Fig. 2f). In October, nighttime cooling was prominent and daytime heating had weakened. The relatively low wind speeds (2 m s −1 ), increased on 4 October (6 m s −1 , Fig. 2c). A homogeneous SML, not observed in the previous surveys, extended down to 10.5 m depth, and rested on an even sharper thermocline ( N max ≈ 0.1 s −1 at 11.5 m depth, Fig 2d, e). A weaker chl-a maximum (13 μg L −1 ) was observed at ~12 m depth coinciding with maximum N. The chl-a peak had broadened (4.5 m) to become asymmetric with a sharp upper edge (Fig. 2f).
April data gave elevated values of 2 × 10 −7 W kg −1 in the upper 8 m. Below this depth, dropped sharply to very low values ( ≤ 10 −10 W kg −1 ; Fig. 3a). In July and October, showed a similar vertical distribution, with larger values between 4 × 10 −9 and 4 × 10 −8 W kg −1 in the upper 14 m, and a second maximum in the depth range of the main The 90% confidence intervals from bootstrap sampling are represented as faded horizontal lines. The depth range of the P. rubescens thin layer, defined as the range where chl-a is larger than an e-folding fraction of the maximum value (chl-a > max (chl-a)/e, is indicated by filled rectangles for July (orange) and October (green) (color figure online) thermocline (10-13 m). From the base of the main thermocline, decreased to 10 −10 W kg −1 at 25 m depth.
The Thorpe scale, L T , describing the size of turbulent overturns, was 1-3 cm in the upper 10 m during April and gradually increased with depth to up to 10 cm at 20-25 m depth (Fig. 3b). In this deeper layer, L T was similar to the Ozmidov scale ( L O ), but L O was 10 times larger than L T in the upper 10 m. In July, L T , was ~3 cm in the upper 5 m but exhibited smaller thickness below ( L T = 1 mm to 1 cm). During this sampling, L O was also about one order of magnitude larger than L T throughout the water column. October data gave relatively large L T and L O values (10 cm to 1 m) that resembled each other in the upper 12 m. These declined sharply below 12 m to the order of ∼1 cm with L T being generally smaller than L O .
April data gave a buoyancy Reynolds number of Re b > 100 in the upper 8 m indicating energetic turbulence (Fig. 3c). At greater depths, Re b < 15 , except for a local peak at 12 m. Smaller values of Re b < 15 indicate that turbulence is buoyancy-suppressed (Ivey and Imberger 1991;Bouffard and Boegman 2013). In July, Re b exceeded 100 in the upper 5 m but declined rapidly below this depth to Re b < 15 in deeper layers. October data showed a sharp transition for Re b at 11 m depth from high values ( Re b ≈ 1000 ) in the SML to molecular values within the thermocline.
The turbulent heat diffusivity ( K T in Eq. 2) followed the same vertical distributions as those exhibited by Re b for the three sampling dates (Fig. 3d). K T approached (April) or fell below (July, October) molecular levels in the depth range where Re b ≲ 15 , including the metalimnion. April and July data gave  Table 1 show probability distributions and mean values for turbulence-related parameters derived from the microstructure profiles within the chl-a maximum during July and October 2018. The vertical extent of the layer was defined as the depth range where chl-a exceeded the e-folding fraction of the maximum (chl-a > max (chl-a)/e). This range should contain 85% of the layer biomass assuming a Gaussian distribution.
In July, energy dissipation rates ranged over three to four orders of magnitude ( from 10 −10 to 10 −6 W kg −1 ) within the depth range of the thin layer (11.75-14.25 m, Table 1; Fig. 4a). This parameter followed a lognormal distribution slightly skewed towards high values. As a consequence, the mean values for of 1.8 [1.6 − 3.4] × 10 −8 W kg −1 ( 90% confidence intervals in brackets) exceeded median values ( 0.99 × 10 −8 W kg −1 ). It is worth noting that only about 42% of the fits to the Kraichnan spectrum met the quality criteria. This indicates that the shape of the temperature gradient spectra was possibly influenced by fine structures under the strongly stratified conditions (Gregg 1977;Luketina and Imberger 2001). In October, the distribution centered on smaller values (median 0.12 × 10 −8 W kg −1 , mean 0.25 [0.21 − 0.78] × 10 −8 W kg −1 ) but also exhibited a small number of outlying peaks of large values ( > 10 −8 W kg −1 ). These suggest a sharp vertical gradient for within the plankton layer (Fig. 3a). Nearly half of the fits ( 47% ) from the October data failed to meet the goodness of fit criteria.
The Ozmidov length-scale ( L O ) exhibited a tight symmetrical distribution for July and October with respective median values of 1.3 and 0.60 cm, and respective mean values of 1.  Fig. 4b). However, L T was 0 for ~90% of the datapoints indicating that overturning motions were undetectable (the median of L T was zero in both samplings). When detectable, the Thorpe scale tended to fall below the Ozmidov scale, i.e., L T < L O . For typical overturns detected ( L T > 0 ) in July and October, L T equaled 0.23 [0.14 − 0.35] cm and 1.5 [0.7 − 2.8] cm, respectively. The ≈ 10 −8 W kg −1 gives a Kolmogorov length-scale on the order of 3 mm. This estimate approaches the size of the detected overturns in July suggesting effective damping by viscous forces.
Within the thin layer, Re b < 15 was found for 87% and 99% of the samples for July and October, respectively (Fig. 4c). July and October data gave median values of 4.8 and 0.86 and

Short-term variability
Microstructure profiles measured on a repeated basis established patterns of some short-term temporal variability in turbulence and chl-a profiles. Figure 5 shows this variability for the days of 19 July and 4 October, when winds reached relatively high speeds (~5 m s −1 , Fig. 5a, f). For the July date, the isothermals remained mostly flat during the 5 h of sampling ( Fig. 5a-e) indicating weak internal wave activity. The P. rubescens layer persisted unperturbed between the 8 and 18 ºC isotherms of the thermocline (Fig. 5e). Despite moderate wind speeds, the Monin-Obukhov length-scale ( L MO ) characterising the penetration of wind mixing below the surface remained shallow ( < 5 m) due to the stabilizing heat flux. The distribution showed little variation around the mean estimate for the profile ( Fig. 5b; blank spaces represent rejected fits for the theoretical temperature-gradient spectrum). The upper 10 m showed some variability with turbulent patches having ≈ 10 −7 − 10 −6 W kg −1 . The L T record illustrates the scarcity and small size of overturns within the main thermocline (empty bins represent L T = 0 in Fig. 5c). In the upper 10 m, estimates of L T ≈ 10 cm showed very little variation. Below 15 m depth, intermittent turbulent patches with L T ≳ 10 cm occasionally appeared. Vertical diffusivity (Fig. 5d) approached or fell below molecular levels throughout the day except in the relatively turbulent patches of the upper SML, but none the turbulent patches reached the P. rubescens layer (Fig. 5d). The profiles collected on 18 July, a day of lower wind speeds, showed even weaker variability (Fig. S3).
Conditions varied more in October when early morning profiles captured the transition from nighttime cooling to daytime heating and wind speeds increased during the morning (Fig. 5f). Over the 6 hours of sampling, the thermocline rose gradually by 2 m (Fig. 5g-j) indicating likely a basin-scale response to enhanced winds. The chla maximum appeared between the 8 and 18ºC isotherms and tracked a stable upwards movement (Fig. 5j). Estimates of turbulence remained relatively uniform between the upper part of the SML and the top of the P. rubescens layer until 9am. These indicated moderate levels of mixing with ≈ 10 −8 − 10 −7 W kg −1 , K T ≈ 10 −2 m 2 s −1 and large turbulent overturns L T > 1 m (Fig. 5g-l). With the onset of warming, L MO diminished ( < 5 m), overturns decreased ( L T < 1 m) and and K T declined, particularly in the lower part of the SML. As a result of the relatively strong winds, and K T remained large in the upper SML down to 7-8 m depth. This daytime wind-driven sub-surface mixing however did not reach the thermocline or cyanobacterial layer.

High-resolution mooring observations within the thermocline
The high-resolution measurements by moored instruments revealed high-frequency internal wave variability not captured by the profilers. The patterns detected included cmscale vertical internal wave-driven isothermal displacements and velocities on cm s −1 scales (Fig. 6). In July, the sampling range (10.84-12.9 m depth) spanned a region included in the main thermocline with a relatively uniform and pronounced temperature gradient (roughly 3.5 °C over 2.5 m; Fig. 6a, b). The overall amplitude of fluctuations in current velocities and isothermal displacements varied little during the 8-h deployment. In October, the sampling range (10.13-12.19 m depth) included part of the SML at the beginning of the deployment. Measurements detected only the upper, sharper edge of the P. rubescens layer, located between the 18ºC and 15ºC isotherms (Fig. 6e, f). For the first two hours, background currents were relatively weak with no fluctuations of the isotherms. After 10:30, current velocities increased gradually, and high-frequency fluctuations commenced. The intensification (up to a maximum of 5 cm s −1 at about 13:00) coincided with an expansion of the pycnocline such that it occupied the whole sampling range by the end of the deployment. This event likely reflects rising winds that morning, which also appeared in the microstructure profiles.
Horizontal f low velocities documented vertical structure within a narrow, meter-scale range during both deployments. The shear in the mean flow ( S 2 = ( u∕ z) 2 + ( v∕ z) 2 ) can trigger shear instabilities to drive turbulence and mixing if it overcomes the stabilizing effect of stratification. This situation occurs when the gradient Richardson number ( Ri = N 2 ∕S 2 ) falls below the critical value of 1/4 (Miles 1961). High-resolution measurements show that, despite the existence of background shear, unstable Ri were not often obtained, due to the strong stratification. Specifically, conditions in July and October met stable criteria ( Ri > 1 ) 88% and 80% of the time (Fig. 6c, g). Furthermore, most unstable values occurred in October within the weakly-stratified SML. These exerted little effect on the thermocline, where the P. rubescens layer occurred.
Energy dissipation rates derived from turbulent velocity fluctuations using structure functions ranged between 10 −10 and 10 −7 W kg −1 in July and rose to values on the order of 10 −10 − 10 −6 W kg −1 during October (Fig. 6d, h). In July, exhibited intermittent structure between calmer (from 10:30 to 11:00 and between 12:00 and 13:00) and more dissipative conditions. The periods of lower dissipation took place under conditions of reduced vertical shear when Ri < 1 events were not detected (Fig. 6c). In October, began with weak values during the first 2 h, when internal wave activity was weak and strengthened thereafter. The average was 1.4 [1.4 − 1.5] × 10 −8 W kg −1 (median: 0.91 × 10 −8 ) and 6.8 [6.3 − 7.3] × 10 −8 W kg −1 (median: 2.9 × 10 −8 ) for July and October, respectively. Table 1 lists these values as well as their respective ranges. The values in July showed very good agreement with microstructure estimates for the same depth range (Fig. 3a). Values derived for L O and Re b also resembled those derived form profile observations (Table 1). In October, high-resolution velocities gave somewhat larger values for , but, Re b remained in the range indicating a molecular regime (mean 21 [19 − 23] , median: 10).

Nighttime sampling: the role of convective mixing
The asymmetric vertical distribution of chl-a in October and enhanced turbulence during the early morning of 4 October 2018 (before solar heating outpaced heat lost) suggest that convective mixing in autumn could potentially disrupt the P. rubescens layer and cause the decline observed in the longterm data (Fig. 1). It was not possible to test this idea with the 2018 dataset, since daytime measurements do not capture layer dynamics during fully developed convective mixing.
To overcome this limitation, a nighttime survey was conducted with a microstructure profiler, a CTD and a FluoroProbe between 14:40 UTC (16:40 local) on 24 September 2019 and 10:40 on 25 September 2019 (Fig. 7). This field campaign occurred under conditions of low wind speeds (2 m s −1 ) and moderate heat loss (Fig. 7a). Observations recorded a 10 m SML with T = 18 °C and a sharp pycnocline between 10 and 15 m depth (Fig. 7b). Consistent with the previous results, data failed to provide good fits to the theoretical temperature-gradient spectrum for the thermocline. As a consequence, could not be determined in several instances (Fig. 7b). Turbulent overturns were rare (Fig. 7c), and vertical diffusivity approached or fell below molecular values (Fig. 7d). The SML, below the Monin-Obukhov depth gave initially low values for ( 10 −10 W kg −1 ) with L T < 1 m and K T approaching molecular scale. After 17:00, a mixing layer with enhanced ( 2 × 10 −8 W kg −1 ), L T ( > 1 m) and K T ( > 10 −4 m 2 s −1 ) penetrated beneath L MO . At about 21:00, the mixing reached the thermocline, which also deepened slightly overnight. The intensity and downward penetration of the convective mixing declined after 04:00 and completely ceased after sunrise when the heat flux was negative (into the lake).  Figure 8 shows effects of convective mixing on P. rubescens distribution as determined by the FluoroProbe measurements. The initial profiles in the early afternoon (14:40) are plotted in light blue in every panel for clarity, and bars representing the L T are included to illustrate the vertical extent and intensity of mixing. Planktothrix rubescens appears as a well defined peak in the afternoon and at relatively low concentration in the SML. Overnight, convection reached the top of the thin layer. The maximum P. rubescens layer concentration progressively declined and the community's concentration in the SML increased due to filament entrainment by convective plumes eroding the layer from above. These results clearly show the disruptive effect of nighttime convection on the structure of the P. rubescens metalimnetic layer in early autumn.

On the nature of suppressed thermocline mixing
Independent estimates of TKE dissipation with microstructure and high-resolution current measurements in the highly stratified seasonal thermocline of Lake Zurich revealed moderate levels of dissipation, ≈ 10 −8 W kg −1 . Dissipation exhibited a degree of intermittency spanning three orders of magnitude. The turbulent heat diffusivity, K T , exhibited persistently small values of the order of 10 −8 m 2 s −1 . Their scale relative to the molecular heat diffusivity means that turbulent fluxes across the thermocline were negligible. These estimates indicate very limited mixing efficiency ( = K T N 2 ∕ ≈ 0.002 ) wherein most of the TKE is lost to viscous dissipation and does not produce mixing. Data also revealed infrequent unstable Richardson number values ( Ri < 1∕4 ) and vertical turbulent overturns ( ≲ 10% occurrence) within detection limits. Detected overturns were characterised by L T ≤ L O ≈ 2 cm, which fell below typical values in the ocean thermocline (Gregg 1987) and even the lowest estimates from hypolimnetic observations in lakes (Etemad-Shahidi and Imberger 2001;Lorke and Wüest 2002). Vertical turbulent overturns store available potential energy, 0.5L 2 T N 2 (Dillon and Park 1987), which is transferred to the mean density profile when overturns collapse. The scarcity of overturns and their small size thus provide strong indication of very weak turbulent mixing in the thermocline.
Different sources of evidence, including observations, laboratory experiments and simulations, indicate that turbulent fluxes become negligible when the buoyancy Reynolds number ( Re b = (L O ∕L K ) 4∕3 ), which measures of the size of the inertial subrange ( L O > l > L K ), falls below O(10 − 100) (Ivey and Imberger 1991;Gibson 1999;Shih et al. 2005). Our thermocline observations fall mostly within this range and thus explain the relative weak mixing observed. Another parameter used to estimate mixing efficiency of turbulence is the turbulent Froude number, Fr T = (L O ∕L T ) 2∕3 . This parameter represents a ratio between the turbulent eddy size required to achieve an equilibrium between inertial and buoyancy forces ( L O ) and the observed vertical eddy size ( L T ). When Fr T approaches or falls below one, turbulent eddies begin to disrupt background stratification, and the system experiences its maximum potential for turbulent mixing (Mashayek et al. 2017). Because L T values fell far below L O for our dataset on the Lake Zurich thermocline, Fr T was very large. This means that the turbulent eddies were too small to disturb the background stratification. Instead, low mean values for Re b and L T (the latter sometimes close to the Kolmogorov scale) suggest strong dampening of turbulent motions by viscosity (Ivey and Imberger 1991). Such situations can produce only very weak turbulent fluxes, consistent with those appearing in the observations described here.
Overall, the results presented here comport with previous turbulence measurements in small and medium size lakes, which showed only very weak turbulent fluxes occurring in their interiors (Saggio and Imberger 2001;Etemad-Shahidi and Imberger 2001;Boegman et al. 2003;MacIntyre et al. 2009;Weck and Lorke 2017). Mixing at the basinscale depends most strongly on processes occurring at the boundaries (Goudsmit et al. 1997;Wüest et al. 2000;Mac-Intyre et al. 1999). With direct measurements using eddy covariance, Etemad-Shahidi and  and Saggio and Imberger (2001) quantified turbulent fluxes in the thermoclines of two stratified lakes. These researchers found similar situations to that described here for Lake Zurich, wherein high dissipation rates ( ≈ 10 −7 W kg −1 ) occurred only with small scale overturns and negligible vertical fluxes. The authors interpreted turbulent dissipation as driven by wave-wave interactions, which resulted in highintermittency and weak vertical fluxes rather than formation of overturning billows due to shear instabilities (Itsweire et al. 1993;Teoh et al. 1997). Turbulence properties of the Lake Zurich metalimnion also resembled those reported in association with thin plankton layers in the ocean. These layers are frequently found within sharp thermoclines, where Re b (McManus et al. 2003;Broullón et al. 2020) and vertical diffusivities ) are small ( Re b < 100 and K T ≲ 10 −6 m 2 s −1 ) suggesting weak vertical turbulent mixing.
Turbulence characterised by low Re b and only occasional small scale overturns, suggests anisotropic turbulent motions. These conditions violate assumptions underlying methodologies used to calculate both from microstructure and high-resolution current measurements (Dillon and Caldwell 1980;Smyth and Moum 2000;Bluteau et al. 2011) such that the large values must be interpreted with caution.

Page 15 of 22 38
About half of the temperature-gradient spectral fits were discarded because the observational spectra did not follow the expected shape according to the established fit criteria. This indicates that conditions were possibly not turbulent in all cases, and fine-scale structures influenced the spectra (Gregg 1977). Dissipation rates derived from high-resolution velocities also depend on a multiplicative fit constant. Under unisotropic conditions, this could deviate from the prescribed isotropic value and generate uncertainties of > 50% in the estimates (Jabbari et al. 2016). Agreement between the two different methods and the consistency with previous observations allow for a certain degree of confidence in our results however. Substantial errors associated with estimates furthermore would not compromise first order interpretations of suppressed metalimnetic mixing.

Basin-scale mixing and thin layer persistence
In addition to its well known seasonal persistence, our highresolution measurements revealed the exceptional stability of the P. rubescens thin layer at sub-daily scales (Fig. 5). This only abated when nighttime convective mixing penetrated into the layer during the September 2019 sampling. The persistence appears due to the reduced turbulent fluxes in the metalimnion. However, our sampling strategy may have undersampled the effect of turbulent mixing on the P. rubescens thin layer for two reasons. First, the forcing conditions during the four surveys may not have captured the full range of weather variability. Second, the P. rubescens community occurs throughout the lake (Salcher et al. 2011;Garneau et al. 2015) including along slope and near-shore regions. As such, the likely effects of boundary mixing (Goudsmit et al. 1997) on the community would not appear in offshore observations like the ones conducted here.
To investigate the role of unresolved mixing, we used hourly meteorological data and bi-weekly water column observations to calculate the lake number ( L N ). We also estimated a basin-scale turbulent diffusivity using the heat budget method and a surface boundary layer scaling (see Supplementary Text 1 for details) (Fig. 9). The lake number indicates the potential of wind stress to significantly disturb the metalimnion by driving upwelling and exciting internal motions. Significant disturbance occurs when L N < 10 . L N exhibited values <10 almost every day over the winter period until April (Fig. 9b). This indicates that wind forcing could easily overcome the existing weak stratification and drive turbulent mixing. From April, the number of days when L N < 10 progressively declines such that no such days occurred in June and throughout the entire summer. This indicates strong damping of interior motions and mixing by stratification. Mixing events reappeared from October onward as the deepening of the SML progressively eroded stratification. The high summer L N values indicate that internal wave dynamics in the lake generally follow the linear regime described by Horn et al. (2001) of reduced energy transfer to small scales and reduced mixing. This confirms the interpretation that our microstructure surveys accurately represent summer conditions in Lake Zurich.
Consistent with the seasonal evolution of L N , spring time measurements gave higher mean basin-scale vertical heat diffusivity ( K BS T , Fig. 9a) below the SML (~1 ×10 −5 m 2 s −1 ) than that estimated for summer and autumn periods ( 10 −6 − 10 −5 m 2 s −1 ). Mixing was suppressed at ~ 12 m depth ( ~15 m) during summer and fall. Data from these times gave K BS T values as low as 1 − 2 × 10 −6 m 2 s −1 . Despite large L N values during summer months, the metalimnetic basin-scale diffusivity exceeded molecular scales by an order of magnitude. This indicates weak but detectable mixing.
To understand how P. rubescens filaments navigate mixing to maintain their vertical position, we compared the vertical restoring velocities that the P. rubescens filaments can achieve when displaced from their depth of neutral buoyancy ( w buoy ), with the velocities required to compensate for dispersion by turbulent mixing ( w eq ) (Fig. 10). This required an equilibrium velocity that was calculated by integrating the advection-diffusion equation assuming steady state and zero net growth (Stacey et al. 2007;Steinbuck et al. 2009): where C represents the concentration profile. In this case we used the P. rubescens chl-a profiles from the bi-weekly monitoring at the beginning of each period (Fig. 1). Figure 9 shows vertical diffusivity profiles used as K BS T . The restoring velocity represents the balance between buoyancy and drag forces: where r s is the equivalent spherical radius of a filament ( ∼ 2 m; Walsby 2005), fill is the density of the filaments, and ≈ 3.97 is a form factor for cylindrical filaments. We assumed that for short term perturbations, the density of the filaments equals the density of water at the P. rubescens chla peak. However, laboratory cultures have shown that if the filaments are exposed to higher (lower) light levels for sufficiently long periods of several hours or days, they gain ballast (buoyancy) to a maximum (minimum) anomaly of fil − = +27 kg m −3 (-30 kg m −3 ). In this situation, the restoring velocities would reach maximum values according to Eq. (6) ( w max buoy ≈ 0.58 m d −1 ; Walsby 2005). Here, we will refer to the latter as a long-term response.

Fig. 10
Vertical velocity ( w eq , triangles) required for P. rubescens to outcompete turbulent diffusion, as derived from the advection-diffusion balance (Eq. 5). Upward (downward) velocities are shown as filled upward-pointing (empty downward-pointing) triangles. Stokes velocities ( w buoy ) for P. rubescens filaments having their neutral density at the P. rubescens peak are shown as black solid lines (i.e. short-term response). Maximum P. rubescens stokes velocities (i.e. long-term response) for maximum/minimum density anomalies (+27 / −30 kg m −3 ) observed for P. rubescens in lab cultures are shown as dashed lines, positive (resp. negative) anomalies are used above (resp. below) the peak. The P. rubescens chl-a profile at the beginning of each period is shown in grey. Spring, summer and autumn periods are as defined in Fig. 1 For all the periods, the short-term restoring velocities ( w buoy < 0.1 m d −1 ) fell below the equilibrium velocities (Fig. 10), suggesting that episodic mixing events disrupt the layer on a transient basis. This situation changes when considering the long-term response ( w max buoy ) indicating that the buoyancy regulation capacity can generally restore the layer after mixing events. In spring, when P. rubescens lived below 10 m depth without layer formation, equilibrium velocities of 0.1-0.5 m d −1 approached but generally fell below restoring velocities (Fig. 10a). This indicates favorable conditions for the P. rubescens motility in rising from deeper zones to form the layer. The time-scale for layer formation by buoyancy alone would equate to ≈ (|z mean |∕2)∕w max buoy ≈ 44 days. This estimate shows good agreement with observations since stratification began at the end of March and the plankton layer was first observed in early May. Vertical migration thus appears to initiate the layer.
In summer and autumn, equilibrium velocities ( w eq ≈ 0.1 m d −1 ) below the layer were several times smaller than w max buoy . This explains the stability of the layer. Equilibrium velocities increased towards the upper edge of the layer in both periods but remained w eq ≲ w max buoy during summer. In autumn w eq > w max buoy , which is consistent with entrainment of the layer by sub-surface mixing.

The role of convective mixing vs. wind mixing
This analysis corroborates direct observations and suggests that metalimnetic mixing driven by internal motions in the interior of the basin or at its boundaries, does not generally reach strengths necessary to overcome P. rubescens mobility. Penetration of sub-surface mixing can disrupt the layer from above in autumn when the SML deepens across the metalimnion. Microstructure observations also suggest that, at least at the time of sampling, direct SML wind-driven mixing does not penetrate deep enough to produce entrainment of the P. rubescens layer. Nighttime convective mixing however generates overturns large enough to penetrate depths necessary to disrupt the P. rubescens layer in autumn. This suggests that convective mixing is responsible for the layer breakdown in autumn.
To assess the generality of this conclusion, we quantified the frequency of conditions favorable for wind and convective entrainment at the base of the SML by comparing the Monin-Obukhov length ( L MO ) and the SML depth ( z SML ) (Fig. 9c). When L MO < 0 , the lake is warming (stable conditions). The reverse applies for lake cooling (unstable conditions). Wind can potentially produce entrainment in the metalimnion when |L MO | > z SML under either stable or unstable conditions. Fig. 9c shows that conditions for wind-driven entrainment occurred frequently only in late spring (April-June), when the SML was shallow. The frequency of wind-driven entrainment peaked at ∼80% in early April. Strong wind events occurred with 20-30% frequency (see figure caption). From June to December |L MO | > z SML events were rare ( 10% frequency). Over that period, unstable conditions potentially conducive to convective entrainment of the thermocline become more and more prevalent ( 50% in July to 80% in December). These results indicate that the deepening of the SML in late summer and autumn occurs due to convective mixing and not wind mixing. In summer, the P. rubescens layer is located at the core of the metalimnion, ∼5 m below the base of the SML. Entrainment cannot disturb it at this position. Once the SML reaches the top of the P. rubescens layer in autumn, however, nighttime convective mixing can entrain the layer and cause its seasonal breakdown.

Comparison with Lake Geneva
Planktothrix rubescens is a widespread species found throughout the world, particularly in perialpine lakes. For those lakes in which P. rubescens accumulates as a metalimnetic thin layer in summer, these bacteria become the dominant planktonic species for several years or decades. Lakes Zurich, Hallwil (Stöckli 2012) and Baldegg (Buergi and Stadelmann 2000) in Switzerland, and Lac de Bourget in France (Jacquet et al. 2005) show this trend. However, in Lake Geneva (Switzerland/France), P. rubescens occurs only in modest concentrations and mainly during autumn months (CIPEL 2019). Furthermore, P. rubescens populations have declined or disappeared in some lakes (i.e., Lac de Bourget, Jacquet et al. 2014), while they remain abundant in others that possess apparently similar nutrient regimes (i.e., Lake Zurich, Posch et al. 2012). To elucidate the reasons behind these differences, Jacquet et al. (2005) analysed a series of bio-physical parameters from two lakes with similar environmental conditions, including a recent history of reoligotrophication, but also having strongly divergent concentrations of P. rubescens (Lac de Bourget and Lake Geneva). These researchers interpreted higher water column stability as contributing to the abundant P. rubescens observed in Lac de Bourget.
Here, we extended this comparison between Lake Geneva and Lake Zurich using meteorological and hydrographic data collected in 2018. Wind speeds were generally weaker over Lake Zurich (2.6 m s −1 ) relative to those affecting Lake Geneva (3.8 m s −1 ) (Fig. 11a). Bulk water column stability, calculated as the one-dimensional version of the Schmidt stability ( Sc ′ , see caption of Fig. 11) exhibited similar maximum summer values but estimates for Lake Geneva showed a delay of more than one month with respect to Lake Zurich (Fig. 11b). These two parameters conferred a different pattern in lake number. Disruptive events ( L N < 10 ) continued to occur regularly over the summer period in Lake Geneva (Fig. 11c) and with relatively high frequency until the end of June. Lake Zurich experienced no disruptive events indicating only weak disturbance of the metalimnion by basin-scale, wind-driven motions.
A series of temperature profiles over the spring-autumn period indicate the different thermal and mixing dynamics for the two systems ( Fig. 11d-g). These showed a faster pace of warming and stratification build up in Lake Zurich that became particularly evident in June (Fig. 11e), when Lake Zurich exhibited a warmer epilimnion and more intense metalimnetic stratification. They also suggest that Lake Zurich experiences weaker summer metalimnetic mixing, as indicated by a stronger and thinner thermocline during late summer and autumn months (August, October). The amount of heat accumulated below 20 m over summer also indicates weaker mixing (Fig. 11d-g). This accumulation arises primarily due to turbulent mixing and fell below that observed in Lake Geneva by a factor of three (78 vs. 281 MJ m −2 ).
Metalimnetic mixing ultimately depends on the rate of mechanical energy supplied by the wind estimated as wind work, P 10 = W 10 (where and W 10 are wind stress and wind speed at 10 m above the lake surface, respectively; Wüest et al. 2000). The lakes also showed average P 10 values for the stratified period (April-October) that differed by a factor of three (55 mW m −2 vs. 171 mW m −2 ). This demonstrates a strong link between differences in wind forcing and the intensity of metalimnetic mixing.
Early onset of the metalimnion is also thought to be critical for the formation of the P. rubescens thin layer (Yankova et al. 2016). This onset took place more than one month earlier in Lake Zurich than in Lake Geneva due to the faster spring warming in the former. Analysis of a 10-year time-series on lake heat content (Fig. S4) revealed a consistent time-lag of ~56 days throughout the time-series confirming this finding. A combination of a stronger thermal inertia (because of the greater depth) and stronger wind forcing enhance turbulent (latent and sensible) heat fluxes (Fig. S5) to contribute to the observed intensity and duration of winter cooling in Fig. 11 Comparison of the stratification and mixing dynamics of Lake Zurich (red) and Lake Geneva (blue) during 2018. a Hourly and weekly wind speed, b one-dimensional Schmidt stability ( Sc ′ ), c Lake number ( L N ) and d-g temperature profiles collected in mid-April, June, August and October 2018. The changes in depth-integrated heat content between 20 and 120 m ( H z>20 m ) with respect to the April profile are shown in d-g. Dates of the temperature profiles in d-g are marked with a cross in b. The datasets used for Lake Geneva are obtained at a single sampling point referred to as SHL2 corresponding to the deepest and pelagic part of the lake. Data were obtained bi-monthly, except for the winter period, for which sampling is performed once a month (Rimet et al. 2020). The one-dimensional Schmidt stability was calculated as Sc � = − ∫ 0 m 120 m g( (z) − )(z − z c ) dz . Compared to Eq. (3), here we do not consider the hypsometric function and we stop the integration at the maximum common depth ( ∼120 m). This was done to make the values for both lakes comparable (color figure online) Lake Geneva. Lake size and the intensity of wind forcing thus appear to be factors supporting P. rubescens growth and development. P. rubescens thrives in mid-sized lakes deep enough to sufficiently isolate the metalimnion from surface and bottom enhanced turbulence, but shallow enough so that thermal inertia allows for an early onset of stratification. Mean and maximum depths of lakes Hallwil, Baldegg, Bourget and Zurich fall in intermediate ranges of 28-45 m and 48-145 m, respectively, relative to Lake Geneva which has mean and maximum depths of 155 m and 309 m, respectively.

Conclusions
This study investigated the interaction between vertical mixing and the seasonal persistence (May-October) of an annually recurring thin layer of toxic cyanobacteria (P. rubescens) in Lake Zurich using four microstructure surveys and high-resolution mooring deployments. Our measurements revealed that strong metalimnetic stratification inhibits overturning motions and vertical turbulent fluxes which provides a very stable environment for the P. rubescens thin layer to form and persist. Comparable sets of measurements available from the oceanographic literature describe similar turbulent environments for a variety of thin plankton layers which frequently appear in association with strong pycnoclines and weakly turbulent or laminar conditions. However, while oceanic features are ephemeral due to frequent disruption by episodic turbulence, weak mixing conditions can persist in Lake Zurich during summer. Our analysis showed that neither entrainment due to surface wind-driven turbulence nor metalimnetic mixing, related to the dissipation of internal lake motions reach the strengths or frequencies necessary to disrupt the layer. This explains its extraordinary seasonal persistence. Instead, the layer breakdown in autumn occurs due to convective deepening of the surface mixed-layer.
A comparative analysis of Lake Geneva, which does not host a P. rubescens layer, indicates that low wind forcing (underpinning the weak metalimnetic mixing) and an early onset of stratification (allowed by weak thermal inertia compared to larger lakes) favor the proliferation and stability of P. rubescens in the form of a thin layer. These two facts would explain why P. rubescens can form thin layers and thrive in medium-sized lakes. Persistent and thin metalimnetic chlorophyll maxima and weak interior turbulent mixing are common features in many lakes. Our results illustrate their connections and underline the relevance of convective mixing in the seasonal evolution of plankton blooms in stratified waters.