Seasonality modulates wind-driven mixing pathways in a large lake

Turbulent mixing controls the vertical transfer of heat, gases and nutrients in stratified water bodies, shaping their response to environmental forcing. Nevertheless, due to technical limitations, the redistribution of wind-derived energy fuelling turbulence within stratified lakes has only been mapped over short (sub-annual) timescales. Here we present a year-round observational record of energy fluxes in the large Lake Geneva. Contrary to the standing view, we show that the benthic layers are the main locus for turbulent mixing only during winter. Instead, most turbulent mixing occurs in the water-column interior during the stratified summer season, when the co-occurrence of thermal stability and lighter winds weakens near-sediment currents. Since stratified conditions are becoming more prevalent –possibly reducing turbulent fluxes in deep benthic environments–, these results contribute to the ongoing efforts to anticipate the effects of climate change on freshwater quality and ecosystem services in large lakes. During the summer season, turbulent mixing in Lake Geneva is strongest in the interior water-column because stratification limits the reach of wind-driven mixing, suggest meteorological, hydrodynamic and turbulence measurements over a full year.

L ake thermal regimes are being seriously altered worldwide due to ongoing climatic change 1,2 , with implications for the health of these ecosystems 3,4 . By controlling the vertical exchange of heat, greenhouse gases, oxygen, and nutrients [5][6][7] , turbulent mixing in stratified lakes regulates the strength of deep stratification 8 and the thermal and ecological response to external forcing [9][10][11] . A comprehensive characterisation of the processes conducive to mixing in lakes is imperative to refine predictions of the fate of lake ecosystems under changing environmental conditions. Although mixing in the stratified interior of lakes relies almost exclusively on wind stress as an energy source 12 , energy pathways towards mixing are often tortuous, mediated by nonlinear interactions 13,14 and dependent on several morphological and environmental factors [15][16][17] . This complexity stems from the fact that, in the bulk of stratified interior waters, wind-driven turbulence occurs only indirectly via excitation of internal hydrodynamic motions (basin-scale seiches, internal waves, gyres, etc. 18 ). Internal motions decay at the lake boundaries by friction with the bottom 19,20 or degeneration through interaction with slopes 21,22 , or within the water-column due to shear instabilities [23][24][25] , all generating turbulent dissipation and diapycnal mixing but with distinct ecosystemic consequences 26,27 .
The detailed empirical quantification of energy dissipation and mixing pathways is hampered by the technical and operational challenges of measuring turbulent fluxes in the field over long periods. The standing knowledge builds upon a few integrative field studies aiming at resolving all relevant components of the turbulent kinetic energy (TKE) budget [28][29][30][31] , interspersed with a larger number of process-oriented campaigns targeting one or a few specific processes 25,[32][33][34][35] . Integrative studies have historically represented major breakthroughs and facilitated the parameterisation of internal lake mixing in one-dimensional computational models [36][37][38] , which have been instrumental for predicting the long-term response of lakes to climatic changes [39][40][41] . Despite their relevance, integrative budgets have so far had limited temporal coverage. Therefore, little is known about the influence of seasonal variability of stratification and wind on the energy pathways and the intensity of turbulent fluxes.
We used a newly developed 100 m 2 floating research infrastructure (LéXPLORE 42 ) to carry out an integrative year-round TKE budget in the deep, large, oligomictic Lake Geneva (Switzerland-France, surface area 582 km 2 , mean and maximum depths 154 and 309 m, respectively, Supplementary Fig. 1). We quantified the temporal variability of the rates of energy input, dissipation and mixing, in the water-column and bottom boundary layer (BBL), between April 2019 and April 2020. The dataset 43 includes continuous records of meteorological variables, full water-column moored measurements of temperature and current velocity (with acoustic Doppler current profilers, ADCP, moored near the surface), and near-bed temperature, current velocities, and turbulence (measured with a high-resolution current profiler). The continuous observations were interspersed with 447 regular (approximately weekly) microstructure turbulence profiles. This unique dataset is used here to unravel the mixing pathways across seasons depending on varying stratification and wind conditions. Our findings substantiate that energy pathways are modulated by seasonal stratification, with the benthic boundary layers being highly turbulent in winter but isolated from wind-driven mixing during the stratified season.

Results
Seasonal stratification. The sampling covered a typical annual cycle of stratification and destratification conditions 44 (Fig. 1). During the spring period (April-July 2019, green bar) the net surface heat flux warmed the lake (Fig. 1a). Water-column temperatures (Fig. 1b) and the Schmidt stability (Sc, Fig. 1a), representing the strength of density stratification, increased over this period, as a broad seasonal thermocline formed, extending from the base of a shallow surface mixed layer (<10 m) to~30 m depth (Fig. 1b). During summer (July -October 2019, yellow shading), heating progressively declined, water temperatures and Sc were maximal by early August and decreased thereafter. During autumn (October 2019-January 2020, brown bar), cooling was maximal, water temperatures dropped and the surface mixed layer deepened steadily from <10 m to~70 m at the end of the period, eroding the main seasonal thermocline. Following this trend, Sc decreased, reaching minimum values at the onset of the winter period (January-April 2020, blue bar). During the first month of this period, the cooling rate progressively reduced. By mid-February, the lake experienced the minimum water-column temperature (7-8 ∘ C) and the deepest mixed layer (90 m). Thereafter, warming was re-established and a weak shallow stratification formed.
Wind energy input. The supply of mechanical energy by wind, quantified as the downward flux of wind mechanical energy (P 10 ) at 10 m above the lake surface, showed marked seasonal variations. Over the annual cycle, P 10 averaged to 172 ± 243 mW m −2 (±standard deviation; Fig. 2a, Table 1), corresponding to a mean wind speed of W 10 = 2.9 ± 2.5 m s −1 , and was minimal in summer (respectively, 89 ± 228 mW m −2 and 2.6 ± 2.0 m s −1 ), and maximal in winter (resp., 261 ± 670 mW m −2 and 3.2 ± 2.9 m s −1 ) (Figs 1c, 2b,c). The actual transfer of wind energy into internal lake motions was quantified with the rate of wind work on surface currents (RW) determined from the correlation of wind stress and current velocity corresponding to the surface expression of deep-water motions (see Methods). The mean annual RW was 0.61 mW m −2 , corresponding to an overall energy transfer efficiency of γ W = RW/P 10 = 0.36%. An alternative estimate of γ W = 0.43 ± 0.05%, corresponding to 0.72 mW m −2 , was computed with a linear regression of RW against P 10 (Fig. 3). For low P 10 ≲ 500 mW m −2 , negative RW values (i.e. wind working against lake currents) were as likely as positive, while positive RW were only clearly dominant for strong winds (Fig. 3). This suggests an important role of inertia in modulating the energy uptake by the lake, since only strong and persistent winds seem able to displace lake waters along the wind direction 45 . Local effects, contingent to the alignment of the wind and currents at the sampling location, explain the large scatter observed in this figure. Taking this variability into account, we assume γ W = 0.36 − 0.43% as a plausible range of the fraction of mechanical energy that is injected into the lake's internal motions. This small fraction is counterbalanced by the bulk of the downward wind energy flux (P 10 ) being dissipated in the atmosphere and lake surface boundary layers and temporally stored in the surface wave field, which dissipates mostly in near-shore areas 12 .
Internal motions. Despite the large seasonal variations in the wind energy flux, which was three times larger in winter than in summer, the seasonally-averaged, depth-integrated kinetic energy (KE) stored in internal motions measured with ADCP's varied weakly around a mean value of 87 J m −2 (Fig. 1d, Table 1). Indeed, contrary to the energy input, mean KE was even slightly smaller in winter (79 J m −2 ) and autumn (81 J m −2 ), compared to spring (98 J m −2 ) and summer (91 J m −2 ). Stormy wind events in spring (e.g. May 2019) and winter (e.g. February 2020) produced short-lived events of high KE that were rapidly damped, not affecting the mean KE significantly. The available potential energy (PE) contained in baroclinic internal motions, quantified using temperature measurements recorded with two moored systems ( Supplementary Fig. 2), was slightly larger than KE (the average PE:KE ratio was ≈1.26). With this ratio, we estimated that~196 J m −2 of available mechanical energy (ME = KE + PE) were stored on average as internal motions (Fig. 2a).
Concerning the spatio-temporal structure of internal motions, the measured currents were mostly zonal, following the coastline orientation at the sampling location ( Supplementary Fig. 1). Although the temporal variability of current velocities was concentrated around periods of 2-8 days year-round, both near the surface (4 m) and bottom (100 m), the overall intensity of near-bed velocities exhibited a strong seasonal modulation with lower values between June and October ( Supplementary Fig. 3). This modulation was caused by a shift in the vertical structure of currents (Fig. 1e), from a two-layer structure with a node at the center of the thermocline (~15 m) during summer to vertically homogeneous profiles in winter (Fig. 1e). The summer current structure is consistent with the first baroclinic mode (M 1 ), which dominated KE during this season, while the barotropic mode M 0 dominated KE in winter (84%) and autumn (62%) (Fig. 1d). Summer observations are compatible with coastal-trapped internal Kelvin waves 46,47 . These co-existed with a net depthintegrated westward transport during spring and summer ( Supplementary Fig. 1), explained by the contribution of the barotropic mode, which represented 33% and 56% of KE, Fig. 1 Stratification, wind speed and currents. Time-series of (a) weekly-averaged net heat flux across the surface (Q net , shading) and water-column Schmidt stability (Sc, dots represent daily averaged microstructure/CTD casts and the solid-line was computed by smoothing with a 30-days Gaussian window), (b) lake temperature smoothed with a 30-days Gaussian window, (c) wind-speed at 10 m above the lake surface (W 10 ) (weekly average in black and 10 min values in grey), (d) depth-integrated kinetic energy (KE) split in barotropic (M 0 ) and three baroclinic (M 1 -M 3 ) normal vertical modes, and (e) zonal velocity reconstructed by projecting the measured velocities to normal modes (colour scale, u NM ) and temperature (grey contours). In panel (b) the black and white markers at the bottom represent microstructure and CTD casts, respectively, and the green dashed line indicates the mixed layer depth (0.5 ∘ C temperature difference with respect to 0.5 m depth). In panel (c) dots represent weekly mean W 10 centred at days of microstructure profiling. The sampling period was divided in four 3-months periods defined as spring (green, 15 April-15 July 2019), summer (yellow, 15 July-14 October 2019), autumn (brown, 14 October 2019-13 January 2020), and winter (blue, 13 January-13 April 2020). Mean seasonal values of W 10 , P 10 and KE are shown in panels c. and d., respectively. In panel (d) the relative contribution (%) of the lowest mode (M 0 ) to KE is indicated in brackets and the mean KE values for each period are represented by horizontal black lines.
respectively. This suggests the presence of a basin-wide cyclonic gyre as a response of dominant northeasterly winds 48 . During autumn and winter periods, the net depth-integrated transport was more erratic, due to stronger and more variable wind forcing ( Supplementary Fig. 1).
Water-column interior dissipation. The fate of the injected wind mechanical energy was assessed by quantifying the rates of energy dissipation and mixing in the water-column and the BBL. The TKE dissipation rates (ε) in the water-column determined from microstructure profiles were in the range of ε int % Oð10 À10 À 10 À8 Þ W kg −1 (mean: 4.4 [4.2 − 4.7] × 10 −9 W kg −1 [99% confidence interval], median: 7.0 × 10 −10 W kg −1 ). In spring, summer and autumn, the average ε profile decreased with depth across the thermocline from~10 −8 W kg −1 in subsurface layers to~10 −10 W kg −1 at~70 m depth ( Supplementary Fig. 4). In winter, ε profiles were almost homogeneous, with some enhancement in the upper and lower-most layers. The depthintegrated TKE dissipation in the water-column interior (between the surface and bottom boundary layers) was between E int ¼ 0:04 À 0:72 mW m −2 (5 and 95% percentiles) for individual profiles and averaged to 0.27 ± 0.24 mW m −2 (± SD) during the sampling period, with little seasonal variations ( Supplementary  Fig. 5, Table 1).
The microstructure profiles were collected on 51 individual days along the annual cycle and may not be representative of the system under all conditions. We attempted to overcome this limitation by deriving a scaling relationship between E int and KE at the time of profiling using a type-II log-log linear regression 49 . This correlation suggested a significant dependence of dissipation on KE, E int / KE 0:64 ± 0:04 (Fig. 4a). This empirical relation provided hourly E int estimates based on ADCP-derived KE ( Supplementary Fig. 5b), which averaged at 0.27 ± 0.22 mW m −2 over the annual cycle, very close to the microstructure observations (Table 1). Coherently with KE, measured and extrapolated E ADCP int seasonal-mean estimates showed minor variability, mostly within <15% of the annual mean ( Fig. 2a-d, Table 1).
Bottom boundary layer dissipation. In the BBL, at 1 m above lake bed (mab), TKE dissipation rates obtained from the Aquadopp measurements varied within the range ε 1mab % Oð10 À10 À 10 À6 Þ W kg −1 , with higher values during the most energetic episodes in February 2020, when near-bed velocities reached~10 cm s −1 (Supplementary Fig. 6). The average ε 1mab during the Aquadopp deployment period (August 2019-August 2020, mean: 2.2 [1.4 − 3.6] × 10 −8 W kg −1 , median: 4.1 × 10 −9 W kg −1 ) was about one order of magnitude larger than ε int . The TKE dissipation rate integrated over the Ekman BBL and averaged during the deployment time overlapping with the study period (August 2019-April 2020) was E BBL ¼ 0:16 ± 3:38 mW m −2 . The large standard deviation reflects the strong dissipation events occurring in winter (Supplementary Fig. 5c). As a consequence, seasonally averaged E BBL increased by almost one order of magnitude from a minimum in summer (0.05 mW m −2 ) to a maximum in winter (0.42 mW m −2 ) ( Table 1).
Since near-bed observations were not available in spring 2019, an indirect estimate was computed from near-bed velocities measured with the surface-moored ADCP (see Methods). The derived mean annual BBL dissipation was E ADCP BBL ¼ 0:20 ± 0:77 mW m −2 , In panels a-c., E lam BBL is included in E BBL , and RW has been slightly adjusted to match the sum of energy sinks. In (d) grey bars indicate 99% confidence intervals. e The relative contribution of the BBL to total dissipation and mixing as a function of stratification (Schmidt stability, Sc) and wind energy flux (P 10 ). Weekly averages are used for all the quantities in (e). Pearson correlation coefficients (r) of the BBL contribution to dissipation with Sc and P 10 are reported along with the corresponding p-value.   (Fig. 2a-d). These estimates were in agreement with the Aquadopp measurements within roughly a factor of 2 (Table 1), which is coherent with the error of~40% estimated for the drag coefficient used in this calculation (C D = 0.0032 ± 0.0013). In addition to turbulent dissipation, KE in the BBL is also lost by laminar friction within the thin viscous sublayer 28,50 . This term averaged to E lam BBL ¼ 0:11 ± 0:38 mW m −2 (Fig. 2d).
Mixing. The fraction of TKE spent in vertical mixing, that is, in producing a net downward buoyancy flux (or an upward mass flux) across the stable stratification, was evaluated using the flux Richardson number, R f = b/(ε + b), where b = K T N 2 is the turbulent buoyancy flux, K T is the turbulent diffusivity 51 , and N is the buoyancy frequency (a measure of density stratification). Mean R f , computed from microstructure profiles in the water-column (Fig. 4b) and near-bed measurements (Fig. 5, Supplementary  Fig. 7), were R f ≈ 0.14 and~0.12, respectively, both close to the canonical value of R f = 0.16 − 0.20 52,53 , suggesting a relatively tight coupling between b and ε (insets in Figs. 4b, 5b). The relatively high value of R f in the BBL is at odds with previous studies suggesting that R f can be much lower in this layer, where strong turbulence normally operates in a weakly stratified environment 27,54,55 . The near-bed measurements showed that, despite episodes of strong dissipation, the BBL was moderately and persistently stratified (mean N 2 = 0.91 × 10 −5 s −2 , Fig. 5a). These findings endorse the idea that secondary circulations continuously resupply buoyancy to sloping regions 19,56-58 , maintaining stratification and allowing mixing to continue. Taking these results into account, we assumed a time-constant mixing efficiency (Γ = R f /(1 − R f ) ≈0.17, and 0.13 in the water-column interior and BBL, respectively) and estimated the time-averaged, depth-integrated buoyancy fluxes as B ¼ ΓE. Annually, water-column and bottom turbulence produced a buoyancy flux of B int ¼ 0:046 ± 0:037 mW m −2 and B BBL ¼ 0:026 ± 0:100 mW m −2 , respectively.

Discussion
In summary, the correlation between wind and near-surface water velocities indicates that annually γ W = 0.36−0.43% of the downward wind energy flux (P 10 ≈ 172 mW m −2 ) feeds the internal motions of Lake Geneva (RW = 0.61 − 0.72 mW m −2 ) (Fig. 2a), which store~196 J m −2 of mechanical energy, giving a mean energy residence time of τ = ME/RW ≈ 3.5 days. The sum of the energy sinks (~0.65 mW m −2 , equivalent to 0.38% of P 10 ) roughly balances the energy input. The sinks are split between dissipation (E int þ E BBL þ E BBL;lam ¼ 0:58 mW m −2 , represent-ing~89% of RW) and mixing (B int þ B BBL ¼ 0:072,~11% of RW). The energy sinks are almost equally shared between the water-column interior and the bottom boundary layer, with turbulent dissipation accounting for 0.27 mW m −2 (42% of RW) and 0.20 mW m −2 (30% of RW) in the respective regions, and laminar dissipation in the viscous BBL contributing with 0.11 mW m −2 (17% of RW). Because laminar dissipation (E lam BBL ) does not produce mixing, the contribution of water-column turbulence to mixing (B int % 0:046 mW m −2 ) is larger than BBL turbulence (B BBL % 0:026 mW m −2 ).
Unique to our dataset is the one-year temporal coverage, which revealed the key role of seasonal variations in both the energy input rates and the energy pathways down to dissipation and mixing (Fig. 2b,c). This result contrasts previous energy budget analyses circumscribed to short periods, mostly during the stratified seasons. Due to seasonal changes in P 10 , the wind energy transferred to internal motions was three times larger in winter (1.00 mW m −2 ) compared to summer (0.32 mW m −2 ), even though the mean energy transfer efficiency seemed to remain rather stable over the annual cycle at γ W ≈ 0.38% (Fig. 2). The water-column interior dissipation pathway reacted weakly to the changes in forcing, with a basal loss rate of~0.27 mW m −2 throughout the seasons (coherent with the annual average), while the bottom sink damped rapidly the excess energy supplied in winter. Therefore, the water-column interior sinks dominated (~85%) the energy loss in summer, whereas the bottom pathway dominated (~68%) in winter (Fig. 2b, c). Despite the variations in the energy flux, internal motions maintained a fairly constant energy level year-round and, as a consequence, the energy residence time was substantially longer during summer (6.6 days) compared to winter (2.2 days), indicating a distinct functioning of the system across the seasons. In summer, the energy residence time is close to that expected for a quasi-free oscillating system with a maximum depth of 309 m (~8 days) 28 , suggesting that stratification allows to transiently store wind energy away from the dissipative bottom environment in form of baroclinic motions. Conversely, the much shorter energy residence time in winter indicates that strong wind pulses dominate the lake's energetics in this period, through transferring large amounts of energy, which get rapidly dissipated at the BBL. Therefore, instead of the slow continuous relaxation of the stored wind energy described for summer, most winter turbulent fluxes occur via a faster pathway, as a series of short, bottom-locked highly dissipative events. Interestingly, the constancy of the mean mechanical energy content in spite of seasonally variable energy fluxes resembles the properties of self-organized criticality described for a number of complex phenomena and recently attributed to geophysical flows 59 . According to this concept, turbulence re-organises the flow in such a way that a statistically stable mean state is maintained notwithstanding the variable forcing.
The large contribution of the lake water-column interior to dissipation and mixing on an annual basis and particularly during the summer season challenges the standing view that the bottom boundary layer is the dominant energy sink below the surface layer of stratified lakes 28,60 , and confirms previous limited evidence derived from observations in large lakes 25,33 . In addition, our results provide the first evidence that the contribution of the BBL to dissipation and mixing is strongly modulated by the intensity of wind forcing (r = 0.47) and, more important, watercolumn stability (r = − 0.81; Fig. 2e). The role of wind forcing is explained by strong winds transiently energising the BBL, where dissipation has a strong dependence on the flow velocities (E BBL / u 3 ). The chief role of stratification emerges because it regulates the transition from baroclinic motions during the stratified season to depth-uniform currents in winter, as previously described for other lakes 33,61,62 . Baroclinic motions retain most of the kinetic energy above the thermocline so that, in deep lakes with a shallow thermocline, near-bed velocities are comparatively reduced during periods of stratification.
The above findings suggest that the energetics of the BBL and their contribution to total dissipation and mixing in such lakes could be affected by changes in stratification and wind conditions, making them particularly sensitive to climate change. Like other perialpine and temperate lakes worldwide 1,63,64 , Lake Geneva has experienced an overall warming and increased prevalence of stratified conditions over the past decades. To assess the potential impact of these changes in the distribution of energy dissipation and mixing, we quantified the annual occurrence of stratification values above a threshold (Sc > 15 kJ m −2 ) for which the contribution of BBL dissipation becomes small (≲ 20%, Fig. 2e) in a historical (> 1980) time-series of density profiles collected monthly to bi-weekly in Lake Geneva (Fig. 6). This analysis revealed a significant increase of the number of days of highly stratified conditions at a rate of 8.6 ± 3.1 days per decade. Simulations of a Representative Concentration Pathway RCP6.0 scenario using a 1-dimensional lake model 65 confirmed that this trend will continue over the 21st century (Fig. 6). Interpreted in the light of our results, these projections suggest that turbulence in the deep bottom boundary layers of Lake Geneva (and possibly other deep temperate lakes) could become progressively weaker over the next century, slowing down the exchanges between the sediments and the water-column.
Although our annual energy budget provides an unprecedented temporal description, its significance at the basin scale may be limited, since the measurements were collected at a single location over the shelf break. The local depth was~30% shallower than the lake mean depth and the depth-integrated quantities (e.g. KE, E int ) could be underestimated. This is unlikely the case for the stratified periods when KE and dissipation rates decayed rapidly with depth, but may be a concern in winter. In addition, since Lake Geneva lateral dimensions exceed the internal Rossby radius of deformation, internal motions are influenced by Earth's rotation 46 , which could impact the distribution of BBL energy between pelagic and near-shore zones. To assess this, we examined a 6-month record of near-bed currents at a central location (Supplementary Fig. 8). The data revealed the presence of anticyclonic near-inertial Poincaré waves, not observed in the shelf break location, and showed that the BBL was overall less energetic with dissipation rates lower by a factor of 5, too low to account for a relevant fraction of the energy input, even in winter. These results confirm that the BBL at lake slopes must represent a major energy sink during the energetic periods, and that the shelf break location of our sampling point allowed to document this contribution. Further, the representativeness of the location is strongly supported by the tight coherence between the estimates of energy input and total dissipation for the different seasons (Fig. 2d), suggesting that, despite the obvious limitations, our energy budget captures the main seasonal variations of winddriven turbulent fluxes in the stratified interior of Lake Geneva. Together with the interior mixing processes considered here, turbulent mixing in the surface boundary layer -directly driven by the input of kinetic and potential energy from the atmosphere-plays a major in controlling stratification and turbulent fluxes in the upper water -column and the exchange with the atmosphere. Those processes are tightly connected to episodic atmospheric forcing events and their full range of variability is challenging to capture using observations. The availability of novel research infrastructures, such as the LéXPLORE platform 42 , and the development of autonomous profiling instruments are set to offer new opportunities for improving our understanding of these processes in the near future.

Conclusions
With an exhaustive field campaign, we described the seasonal energy pathways towards dissipation and mixing in a large lake with unprecedented temporal coverage and resolution. The fullyear duration of the measurements unveils a previously overlooked seasonal modulation of the energy sinks and shows that the benthic environment is subjected to highly variable energetic conditions. Stratification isolates the deep benthic layers from wind energy and maintains very weak levels of turbulence in summer, but the same layers, particularly in slope regions, are responsible for damping the strong energy input during the weakly stratified winter season. Seasonal variability in near-bed turbulence has important biogeochemical and ecological implications, because near-bed turbulence controls the sediment-water exchange of heat, environmentally important gases 66,67 , metals, and nutrients 68 , influences benthic microbial communities 69 , and regulates the development and duration of anoxic conditions 70 . This is particularly relevant as global warming is driving a worldwide increase in the duration and intensity of stratification 1,71,72 , which, based on the relationship found between water-column stability and BBL dissipation and mixing, could reduce the overall energy levels in the deep boundary layers.
Notwithstanding the seasonal variability of the input, dissipation, and mixing fluxes, the fraction of wind energy transferred to internal motions (RW/P 10 ≈ 0.38%) and available for mixing (B=P 10 % 0:04%) remained remarkably constant throughout the annual cycle. Similar ratios have been reported for lakes of very contrasting morphology and forcing, ranging from the most voluminous lake on Earth, Lake Baikal 29 (RW/P 10 ≈ 0.2%, B=P 10 % 0:04%), to mid-sized lakes like Lake Windermere 30 in Britain (RW/P 10 ≈ 0.13−0.64%, B=P 10 % 0:05%) and the small Lake Alpnach 28 in Switzerland (RW/P 10 ≈ 0.7%, B=P 10 % 0:06%). These results suggest that in spite of the intricate energy pathways, the strong dependence of mixing on wind energy tightly constrains the magnitude of turbulent fluxes in different systems, thus providing a simple rule of thumb to estimate the overall mixing levels of lakes exposed to variable wind forcing and assess the ecological implications 73 . Altogether, our thorough characterization of the mechanical energy pathways broadens the field-based knowledge on how mixing operates in lakes, offering a valuable benchmark for model development and assessment, which will allow to improve the predictions of the future evolution of these ecosystems.

Methods
Meteorological forcing. Meteorological parameters (air temperature, relative humidity, shortwave radiation, wind speed, wind direction, and atmospheric pressure) were recorded every ten minutes since 8 April 2019 with a Campbell Scientific Automatic Weather Station installed on the platform, 5 m above the lake level. Turbulent heat fluxes (latent and sensible) and wind stress (τ w ) at the lake surface were derived following the Monin-Obukhov similarity theory 74,75 , whereas radiative fluxes (shortwave and longwave) were derived with formulas calibrated for Swiss lakes 76 . The multiplicative factor in the longwave absorption formula was specifically calibrated for Lake Geneva to ensure that heat fluxes were coherent with the temporal evolution of heat content in the lake.  above which the contribution of the BBL dissipation becomes small (≲0.2) according to Fig. 2e. The black dots correspond to historical observations collected monthly to bi-weekly by CIPEL 77 at the central lake station (SHL2, Supplementary Fig. 1). The blue lines (thin dotted: individual forecast, thick: ensemble mean) and shading (SD) were obtained from numerical simulations using the Simstrat 1-dimensional lake model forced with downscaled Representative Concentration Pathway RCP6.0 climate projections from global climate models 65 . Linear trends are reported for both datasets.

Water
Water-column currents. Profiles of horizontal currents were recorded with two RDI Acoustic Doppler Current Profilers (ADCP) from 16 April 2019 to 15 April 2020. Until 8 August 2019, an RDI600 and an RDI300 were installed in a mooring located 8 m below the surface, looking upward and downward, respectively. Thereafter, only the RDI300 was used, installed directly at one of the LéXPLORE moon pools 42 looking downward from the lake surface. The RDI300 recorded 300 pings every 15 min over a vertical grid of 53 cells with 2 m thickness, extending down to the lake bottom, and the RDI600 (when available) recorded 1000 pings every 10 min vertical grid of 31 cells with 0.25 m thickness, extending up to the lake surface. Shorter ensemble intervals and smaller cell sizes were used during the first deployment (until 12 June 2019), occasionally yielding poor data quality. For the analysis, the data were averaged with running windows of 2 h in time and 4 m in depth.
Water-column microstructure and CTD profiles. A total of 447 deep (down to 80-100 m depth) microstructure turbulence profiles were performed approximately weekly on 51 different days between March 2019 and April 2020. Three different Rockland Scientific International Inc. (Canada) instruments were used: two small self-contained microCTD profilers (SN310, SN158) and a VMP-500 (SN028). In addition to high-resolution microstructure temperature sensors (2× FP07 fast thermistors), all three profilers were equipped with accurate conductivitytemperature-depth (CTD) and optical sensors (turbidity, fluorescence). The buoyancy of the instruments was regulated to achieve a free-falling profiling rate of~35 cm s −1 . In addition, full-depth (down to~105 m) regular CTD profiles were carried out with a Sea&Sun CTD75M probe. In total, 99 CTD profiles were collected on a regular basis (typically weekly or bi-weekly) between February 2019 and April 2020.
Near-bed currents and temperature. A near-bed moored structure was deployed from 16 August 2019 to 12 August 2020 to record high-resolution near-bed velocities and temperature. Velocities were measured with a high-resolution pulsecoherent ADCP, Aquadopp 2 MHz (Nortek). The instrument was mounted looking downward on a horizontal mast of the mooring structure at 2 mab. The instrument was located 2 m away from the main body of the structure to minimise possible distortion of the measured currents. Every hour, a 1024-datapoint burst of along beam velocities with a vertical spacing of 2 cm was recorded at 2 Hz. Since 17 December 2020, seven RBRsolo T thermistors (accuracy 0.002 ∘ C, sampling frequency 2 Hz) were installed along the main vertical mast of the structure between 0 and 3 mab, spaced by 0.5 m.
Analysis: thermal cycle and lake stability. To characterise the seasonal thermal stratification cycle, a composite annual cycle of temperature, salinity, and potential density (ρ) was computed by merging the mooring, CTD, and microstructure profiler datasets by using a running average with a 30-days Gaussian window in the top 30 m and increasing length at lower depths, where data was not recorded continuously. The water-column Schmidt stability (Sc, i.e. strength of density stratification) was calculated from the running averaged ρ profiles as: where A = A(z) is the hypsometric curve, A 0 = 580 km 2 is the surface area, z bot is the maximum lake depth (309 m) and z v is the lake's center of volume. Since CTD profiles at the LéXPLORE location were restricted to 109 m, density profiles were completed down to z bot using deep casts collected by CIPEL monthly to bi-weekly at the deepest location 77 (SHL2, Supplementary Fig. 1).
Analysis: energy input. The downward flux of wind mechanical energy at 10 m (P 10 ) was taken as a reference of wind energy supply to the lake. This was calculated as P 10 = τ w W 10 , where τ w and W 10 are the wind stress modulus and wind speed at 10 m above the lake surface, respectively 28 . The actual energy transfer was quantified with the rate of wind work (RW), calculated as the correlation between wind stress and near-surface water velocities: RW = τ x u 0 + τ y v 0 , where τ x,y are the zonal and meridional wind stresses and u 0 , v 0 are the zonal and meridional nearsurface current velocities derived from ADCP measurements. To ensure that RW provides a measure of the energy transferred to basin-scale internal motions, we followed Shimizu et al. (2008) 20 and took u 0 and v 0 as the velocities at 0.5 m depth reconstructed by projecting the observed velocity profiles into normal vertical modes (see below).
Analysis: internal motions. The kinetic energy (KE) stored in internal motions was characterised by projecting the ADCP-measured velocities onto normal vertical modes via least-squares fitting 78 . The shape of the vertical modes was computed by resolving the Taylor-Goldstein equation for the buoyancy frequency 79 , N 2 = − gρ −1 (∂ρ/∂z). The N 2 profiles and normal modes were computed daily based on the smoothed seasonal ρ profiles. Similarly, the available potential energy of internal motions was computed from temperature mooring measurements by projecting the isothermal displacements δ onto normal modes, where δ ¼ ÀðT À TÞð∂T=∂zÞ À1 was estimated from the temperature anomalies with respect to the 30-days running averaged profile (T).
Analysis: interior dissipation and mixing. Water-column profiles of turbulent kinetic energy (TKE, ε) and thermal variance (χ) dissipation rates in the watercolumn were derived from microstructure measurements by fitting the temperature gradient spectra computed in 2 m-thick overlapping bins (vertical resolution 0.5 m) to the theoretical Kraichnan spectrum following the MLE procedure 80 , as described in more detail elsewhere 73,81 . Volume-weighted and depth-integrated areal TKE dissipation rates in interior layers (E int ), away from the influence of the surface (SBL, z SBL ) and bottom boundary layers (BBL, z BBL ), were calculated from each individual ε profile as: where V ¼ R 0 À109m A dz % 48:5 km 3 is the total lake volume above the local bottom depth (~109 m). To avoid the influence of turbulence derived from direct interactions with the atmosphere, z SBL was computed by examining the Thorpe scale of overturns 82,83 , L T . The buoyancy flux, used to quantify irreversible mixing, was calculated as b = K T N 2 , where K T = 0.5χ(∂T/∂z) −2 is the turbulent heat diffusivity evaluated with the Osborn and Cox 51 model. Note that, strictly speaking, a minus sign should be added to the b definition, as the turbulent buoyancy flux is downgradient, and hence downward for stable stratification. For convenience, we take both b and ε as positive.
Analysis: BBL dissipation and mixing. In the BBL, a vertical profile of the TKE dissipation rate with 2 cm vertical resolution was computed for each Aquadopp burst following the structure functions method 84 with the relevant corrections for the bottom boundary layer 85 . Depth-integrated dissipation over the bottom Ekman layer (height: h E = κu * f −1 , where κ = 0.41 is the Von Kármán constant, u * is the friction velocity and f is the inertial frequency), E BBL , was derived from the estimate at 1 mab (ε 1 mab , computed as the average within 20 cm of this nominal height) assuming a logarithmic layer behaviour. Following the law-of-the-wall (εðhÞ ¼ u 3 Ã =κh, where h is the height above the bed), the friction velocity was computed as u Ã ¼ ðhκε 1mab Þ 1=3 with h = 1 mab, and depth-integrated dissipation was obtained as: where δ ¼ 11νu À1 Ã is the thickness of the viscous sub-layer and ν is the water kinematic viscosity. Alternatively, estimates of E BBL were derived from a friction velocity u ADCP Ã ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi C ADCP D; 9 mab q u ADCP 9 mab computed using the near-bed velocities at 100 m (9 mab) measured with the surface moored, downward looking ADCP. The drag coefficient (C ADCP D; 9 mab ) was adjusted to match the dissipation rates at 1 mab derived from the structure functions. The best-fitting drag coefficient was C ADCP D; 9 mab ¼ 0:0032 ± 0:0013. In addition to turbulent dissipation, depth-integrated laminar dissipation within the viscous sublayer 28,50 was computed as E lam BBL ¼ 11ρu 3 Ã . To estimate the contribution of BBL turbulence to mixing, buoyancy fluxes were obtained following the same procedure described for the profile data, but using an estimate of χ computed from the near-bed high-resolution temperature time-series recorded with the RBR thermistors mounted in the bottom mooring. Briefly, temperature frequency spectra were calculated for each thermistor during each Aquadopp burst. Then, the frequency spectra were transformed into wavenumber (k) spectra (S T (k)) using the mean flow velocity during the burst at the thermistor height following Taylor's frozen flow hypothesis. By taking the theoretical spectral shape in the inertial sub-range (S I T ¼ C T χε À1=3 k À5=3 , where C T is the Obukhov-Corrsin universal constant typically assumed as C T = 0.4) 86 , χ was estimated as: where the average was performed up to the wavenumber corresponding to a frequency of 0.1 Hz, where the limitation due to the sensors time-response became apparent. BBL stratification, N 2 , was derived from an average temperature gradient computed from the mean temperature profile over each burst using a linear regression. In order to ensure that Taylor's frozen flow was met, we only retained the fits where the mean flow velocity was 5 times larger than the root-mean-square of the velocity fluctuations.

Data availability
The data needed to reproduce the figures and tables in the paper are available at the Eawag Data Institutional Collection 43 (https://doi.org/10.25678/0004YA). More data (i.e. longer time series, raw and processed data) and live visualisation are available here: https://www.datalakes-eawag.ch/. The scripts needed to reproduce the analyses are maintained in the following Github repository: https://github.com/bieitofernandez/ LakeGeneva.