Hypolimnetic oxygen depletion in a deep oligomictic lake under climate change

Dissolved oxygen (DO) concentration is a fundamental metric to describe climate-related alterations in deep lakes. Increasing water temperatures enhance thermal stratification, leading in temperate basins to a growing isolation of deep waters. This leads to the depletion of hypolimnetic DO, which adds up to limited nutrient circulation and restricted replenishment of the trophogenic layers. With vanishing convective mixing, it is commonly believed that the only source of hypolimnetic DO replenishment will be represented by deep intrusions of cold oxygenated waters from the tributaries. In this study, we first analyse the 1993–2020 long-term observed trends of DO concentrations in the subalpine deep oligomictic Lake Maggiore (Italy/Switzerland). Then, through an algorithm calculating daily intrusion depths and mass discharges of DO for the major tributaries, we show that deep insertions are suppressed for increasing winter water temperatures and residual thermal stratification. Turbulent entrainment is proved fundamental for DO replenishment, leading to mass discharges of DO released into the deep hypolimnion up to more than two orders of magnitude larger than the original ones from the tributaries. Last, we discuss the results of simulations made through a one-dimensional coupled ecological–hydrodynamic model about the possible effects of a full turnover on DO concentrations in the deep hypolimnion. Two cases are displayed, with the turnover taking place either now or with an anoxic hypolimnion deriving from decades of isolation due to severe climate warming. Through this study, climate warming is shown to be a fundamental driver of DO in Lake Maggiore, its depletion harming both water quality and the ecosystem.


Introduction
Dissolved oxygen (DO) is a fundamental indicator of water quality (Ito and Kazuro 2015;Rogora et al. 2018), adopted as one of the key parameters to assess the ecological status of a lake according to the EU Water Framework Directive (European Commission 2000). Furthermore, alongside water temperature, hypolimnetic DO concentration is useful for tracking climate-related alterations in deep lakes. Increasing air temperatures can lead to longer stratification periods and a higher degree of thermal stability, which cause prolonged isolation of the hypolimnion from the atmosphere (Rhodes et al. 2017;Mammarella et al. 2018). Continuous stratification extending for several years, with winter turnovers involving only the upper part of the water column, ordinarily results in significant depletion of hypolimnetic DO (Missaghi et al. 2017). Deep-water anoxia may have severe consequences on lake water quality (Liu et al. 2019), promoting phosphorus release from the sediments (Zhang et al. 2015), internal nutrient loading (Komatsu et al. 2007) and accumulation of methane and sulphide (Posch et al. 2012). Deterioration of water quality in turn affects fish habitat, with variations in temporal and spatial distributions (Stefan et al. 2001;Foley et al. 2012).
Hypolimnetic DO depletion has been either observed or predicted for the future for several lakes around the world (Sahoo et al. 2012;Fukushima et al. 2017;Schwefel et al. 2016Schwefel et al. , 2018Fenocchi et al. 2020) and its monitoring is critical for managing and preserving lake ecosystems under climate change (Ito and Kazuro 2015). In case of hypolimnion isolation due to climate-change-induced meromixis, it is generally concluded that the only source of hypolimnetic DO replenishment would be represented by winter deep intrusions of colder oxygenated waters from the tributaries . Deep-water DO renewal by inflows depends on the plunging depths and the volumes of river intrusion plumes (Schwefel et al. 2016;Råmån Vinna et al. 2018). Flood events are the major triggers of oxygen-rich deep insertions (Lenzi and Marchi 2000;Fink et al. 2016), as the associated high suspended sediment concentrations, increasing plume densities, allow larger depths to be reached by density currents prior to the achievement of neutral buoyancy. However, climate warming would affect not only the convective mixing frequency of lakes, but also the discharges and temperatures of their tributaries (Råman Vinnå et al. 2018), thus impacting on the frequency, volume and intrusion depth of the underflows, and in turn on deepwater DO renewal. Fink et al. (2016) studied the effects of potential changes in the hydrological regime due to climate warming on underflows in Lake Constance, using an inflow intrusion algorithm. They showed that Lake Constance will face in the next decades a decrease in deep-water DO concentration due to both the reduction of convective mixing and the change in the seasonal discharge regime of inflowing River Rhine, overall reducing the frequency of deep-water renewal by ~ 10%. Råman Vinnå et al. (2018) examined the indirect impact of climate change on Lake Geneva and Lake Biel brought about by effects on their main tributaries, River Rhȏne and River Aare, respectively. They employed a coupled river-lake model addressing variations in river discharge regime and suspended sediment concentration, changes in lake stratification dynamics and alterations in river intrusion depths and volumes, the latter through an inflow intrusion algorithm. Results showed that river discharges are expected to decrease in summer and increase in winter, thus enhancing the warming of rivers in summer and smoothening it in winter. This would increase the volume of riverine water sinking to deep lake layers in winter, supporting the reoxygenation of the hypolimnion.
In this paper, we studied the dynamics of hypolimnetic DO renewal in the subalpine deep oligomictic Lake Maggiore (Northern Italy/Southern Switzerland). We first analysed the 1993-2020 long-term observed trends of hypolimnetic DO concentrations, highlighting the influence of increasing water temperature and water-column stability. Then, we applied an algorithm calculating the 1993-2020 daily intrusion depths and mass discharges of DO for the four main tributaries of Lake Maggiore, i.e. River Maggia, River Ticino, River Toce and River Tresa, taking into account turbulent entrainment. We will also discuss in the end the results of simulations performed with the GLM-AED2 model (General Lake Model-Aquatic EcoDynamics; Hipsey et al. 2013Hipsey et al. , 2019 of the possible effects of a full turnover on DO concentrations in the deep hypolimnion, considering two scenarios: (1) a complete-mixing event taking place now and (2) a full turnover occurring after decades of hypolimnion isolation brought about by climate warming, as predicted by the Swiss Climate Change Scenarios CH2011 (CH2011 2011Fischer et al. 2015).
The overall aim of our analysis was to provide insight into the future evolution of DO dynamics in Lake Maggiore under climate change. These evaluations may be relevant also for the other large deep lakes south of the Alps (Lake Como, Lake Garda, Lake Iseo and Lake Lugano), which share morphological, geographical and climatic characteristics with Lake Maggiore (Salmaso and Mosello 2010;Rogora et al. 2018).

Case study
Lake Maggiore is a deep fluvio-glacial lake at the foothills of the Southern Central Alps (Fig. 1). About 80% of its surface area lies in Italy, the rest belonging to Switzerland, while the watershed is shared almost equally between the two countries. Lake Maggiore drains the upstream watershed of River Ticino, its only outflow, which includes Lake Lugano, Lake Orta and Lake Varese. The lake has 33 tributaries, the main ones being River Ticino, River Toce, River Tresa and River Maggia (Fig. 1), which have been considered in this study. Within the 1993-2020 study period, these tributaries equalled together 71% of the discharge outflowing from Lake Maggiore through the Miorina Dam (Fig. 1a), which regulates lake levels for water storage purposes. The main morphometric and hydrological features of Lake Maggiore and of its major tributaries are listed in Table 1.
Lake Maggiore is classified as holo-oligomictic (Ambrosetti and Barbanti 1999), full turnovers occurring only at the end of particularly cold and windy winters. Completemixing events within the 1993-2020 study period occurred at the end of winters 1999, 2004 and 2005. Homeothermy of the water column also occurred at the end of winter 2012 ), yet such condition did not last long enough to produce a full chemical homogenisation of the water column. Lake Maggiore suffered from eutrophication in the 1960s and 1970s, reaching the limit between mesotrophy and eutrophy (Ruggiu et al. 1998). The recovery process started in the 1980s, through the establishment of wastewater treatment plants and the reduction of phosphorus in detergents (Kamenir and Morabito 2009). At present, the lake is oligo-mesotrophic, with total phosphorus (TP) concentrations at spring turnover around 12 µg L −1 and primary production dominated by diatoms (Morabito et al. 2012;Rogora et al. 2021).
Lake Maggiore and its tributaries have been routinely monitored since 1981 by the Water Research Institute of the National Research Council of Italy (CNR-IRSA) of Verbania Pallanza (formerly Institute of Ecosystem Study, CNR-ISE), within the framework of the limnological campaigns funded by the International Commission for the  Protection of Italian-Swiss Waters (CIPAIS; www. cipais. org). This regular monitoring has produced long-term datasets for the Ghiffa site of maximum depth ( Fig. 1) and the mouths of most tributaries. Together with Lake Como, Lake Garda, Lake Iseo and Lake Orta, Lake Maggiore belongs to the "Southern Alpine Lake" site (LTER_EU_IT_008) of the Italian and European LTER (Long-Term Ecological Research) networks.

Data
Monthly samples of water physical and chemical variables for Lake Maggiore at the Ghiffa site ( Fig. 1) (Wetzel and Likens 2000). Linear interpolation in space and time was adopted to transform these data into daily values at ∆z = 0.01 m vertical resolution, as needed by the inflow intrusion algorithm adopted here. These regularly interpolated values of water temperature and DO were also employed to assess their evolution within the study period in the 0-20 m epilimnion and in the 200-370 m deep hypolimnion, calculating daily volume-weighted averages through the hypsometric curve extracted from the 2008 bathymetry of Lake Maggiore (Fig. 1b) by the Italian Navy (CNR-IRSA 2021, unpublished data). Mean daily discharges of River Maggia at the Locarno Solduno station, River Ticino at the Bellinzona station and River Tresa at the Ponte Tresa Rocchetta station, just downstream of the Rocchetta Dam that regulates the levels of Lake Lugano, were obtained from the Swiss Federal Office for the Environment (UFAM). Mean daily discharges of River Toce at the Candoglia station were supplied by the Regional Environmental Protection Agency of the Piedmont region of Italy (ARPA Piemonte). These stations (Fig. 1a) were selected as reference ones for the considered tributaries as complete daily series for the 1993-2020 period were available.
For River Ticino, River Toce and River Tresa, the discharge measured at the reference stations was corrected to represent the inflows to Lake Maggiore more accurately. These adjustments were made employing the slopes of the linear regressions of the reference complete discharge series against incomplete ones available for the river mouths (Fig. 1a). To such aims: (1) the 2000-2003 discharges of River Ticino at the Riazzino station obtained from UFAM were employed to correct the values measured at Bellinzona, thus including the inflow from Creek Morobbia; (2) the 2002-2020 discharges of Creek Strona at the Gravellona station prior to the junction with River Toce given by ARPA Piemonte were employed to supplement the discharges at Candoglia, creating an incomplete series for the mouth of River Toce; and (3) the 1998-2013 discharges of Creek Margorabbia at the Germignaga station prior to the junction with River Tresa measured by CNR-IRSA were adopted to integrate the discharges at Ponte Tresa Rocchetta, originating a partial series for the mouth of River Tresa. The resulting correction coefficients of the reference discharges ranged 1.12 -1.17. For River Tresa, the adopted estimation of its discharge released into Lake Maggiore was made possible by the Creva Dam, located 3.7 km upstream of the river mouth ( Fig. 1a), not ordinarily providing flow regulation. This was verified comparing the 2011-2019 discharges released by the Creva Dam supplied by the managing hydroelectric company Enel Green Power against those at the Ponte Tresa Rocchetta station, finding negligible differences. The resulting discharge series of the main tributaries released into Lake Maggiore in the study period are displayed in Fig. 2. The largest flood events involved River Toce, the tributary with the widest and rainiest watershed, in September 1993, October 2000 and October 2020. The mean values of these discharge series are reported in Table 1.
In situ measurements of water temperature and sampling for chemical analyses, from which salinity was calculated as done for Lake Maggiore, were performed monthly near the mouths of the considered tributaries for the 1993-2020 period. Optical turbidity measurements were contextually made from 2002 onwards. All these monthly data were transformed into daily ones to feed the inflow intrusion algorithm.
Daily water temperatures of the considered main tributaries were obtained from the available monthly measurements through the air2stream model (Toffolon and Piccolroaz 2015), which estimates the temperatures of streams through a hybrid approach; i.e. it employs a physically based equation coupled to a stochastic calibration of the included parameters against water temperature observations. Daily water discharges and air temperatures are needed as physical drivers for air2stream. The mean daily air temperatures at the Verbania Pallanza station (Fig. 1a) were employed for all tributaries, data being provided by CNR-IRSA for 1993-2008 and by ARPA Piemonte thereafter. Few missing air temperature values were retrieved from the measurements at the ARPA Piemonte Cannobio station (Fig. 1a), linearly regressed over the Verbania Pallanza series. The full 8-parameter version of air2stream was used, solving the driving equation at daily resolution through the 4th-order Runge-Kutta method. Parameter calibration was performed minimising the Root Mean Square Error (RMSE) 4 Page 6 of 20 between model and observations through the Particle Swarm Optimisation (PSO) algorithm, employing 1000 particles and 2000 iterations (Islam et al. 2019). RMSE values of 1.08 °C, 0.96 °C, 1.08 °C and 1.51 °C against the monthly measurements near the river mouths were obtained for River Maggia, River Ticino, River Toce and River Tresa, respectively. The results of air2stream were compared for River Ticino with the available 1997-2020 mean daily water temperatures measured by UFAM at the Riazzino station, which lies 4.1 km upstream of the river mouth, obtaining RMSE = 0.81 °C. For reference, RMSE = 0.95 °C results between the monthly samples near the river mouth taken by CNR-IRSA driving air-2stream and the UFAM measurements at Riazzino. This shows that model errors are comparable to the uncertainty due to different measurement locations within the final reach of the tributaries. The slightly larger error obtained for River Tresa should be ascribed to its complex short 13.3 km long reach, in which two dams and the junction with Creek Margorabbia 0.6 km upstream of the river mouth are present (Fig. 1a). These alterations cause the river temperatures to elude the physical equations and the stochastic parameterisation of air2stream more than for the other tributaries.
For River Maggia, River Ticino and River Toce, daily salinity values were obtained from the monthly samples through power-law regressions of the latters against water discharges. Exponents − 1 < b < 0 were found, highlighting a definite dilution effect with discharge for these mountain rivers, in which solute content is mostly related to weathering. Such behaviour was not found for River Tresa, the salinity content of which rather echoes that of the epilimnion of Lake Lugano. An annual cycle, following the annual vertical stratification dynamics of Lake Lugano (Holzner et al. 2009), was evident throughout the considered period. Salinity in River Tresa decreases in the April-September period, in which stratification separates the epilimnion of Lake Lugano from the more solute-rich deep layers (Wüest et al. 1992;Holzner et al. 2009), and rises in the remnant part of the year, in which some mixing increasingly occurs. Daily salinities of River Tresa were, thus, estimated through linear interpolation of the monthly samples. The autocorrelation function of the interpolated salinities of River Tresa, put in the standardised variable form, confirmed a definite annual cycle, yielding a R = 0.66 autocorrelation coefficient for lag = 365 d.
Measured turbidity data for the main tributaries were converted into suspended sediment concentrations through linear regression laws obtained from occasional samples for which gravimetric assessments of suspended sediments were performed in addition to optical measurements. Daily values were then estimated through power-law regressions of the monthly values against water discharges, leading in all cases to b > 1 exponents (Syvitski et al. 2000).

Inflow intrusion algorithm
An algorithm was implemented to calculate the daily intrusion depths of the four considered tributaries and the mass discharges of DO released by them at those depths. The algorithm accounts for the turbulent entrainment of lake water by the inflowing plume, following the scheme present in the open-source GLM one-dimensional (1D) hydrodynamic model (Hipsey et al. 2019), based in turn on the one employed by its ancestor DYRESM (Fischer et al. 1979;Imberger and Patterson 1981). The code was written in the MATLAB® language, starting from the standalone implementation of the DYRESM inflow mixing algorithm by Pilotti et al. (2014), and is available upon request. Four main layers were considered to spatially classify the resulting intruded mass discharges of DO: epilimnion (0-20 m), metalimnion (20-80 m), upper hypolimnion (80-200 m), and deep hypolimnion (200-370 m).
Our algorithm takes into consideration the effects on plume density of the salinities of both the lake and the tributaries, as well as of the suspended sediment concentration of the tributaries. Null background suspended sediment concentration was assumed for lake water, being Lake Maggiore a deep basin, in which turbidity is mostly ascribed to phytoplankton and organic matter (Marchetto et al. 2004), which have densities comparable to that of water. A sink term was added to model sedimentation occurring as the plume descends along the lake shore, as in the inflow intrusion algorithms in Fink et al. (2016) andRåman Vinnå et al. (2018). Salinity was included despite the small values characterising Lake Maggiore and its inflowing waters (sampled S = 0.092-0.126 psu for the lake and S = 0.019-0.312 psu for the tributaries within the study period), as it nevertheless influences the intrusion depths, especially for the extremely low thermal gradients that characterise the lower layers throughout the year and most of the water column in winter.
The full technical details of the inflow intrusion algorithm are given in the Appendix.
The maximum annual intrusion depths of the considered tributaries obtained from the algorithm were tested for correlation against the minimum annual values of the Schmidt stability (Schmidt 1928) of the lake, which was computed through the Lake Analyzer v3.4.0 MAT-LAB® package (Read et al. 2011) considering lake temperature and salinity. These tests allowed determining the relationship between inflow intrusions and the residual stratification of the lake, the latter depending on winter climate (Ambrosetti and Barbanti 1999), considering that air temperature warming boosts winter residual stability, thus limiting the depth of convective mixing in the lake Rogora et al. 2018).

Field observations
During the 1993-2020 period, continuous warming of mean daily air temperatures was observed at the Verbania Pallanza station (Fig. 3). The positive trend ∇T a = + 0.051 °C a −1 , significant at the α = 0.001 level according to the modified Mann-Kendall test for autocorrelated data (Hamed and Rao 1998;henceforth mMK), determined a ΔT a = + 1.44 °C within the period.
Air temperature rise caused the warming of Lake Maggiore throughout the water column within the study period (Fig. 4). In the 0-20 m layer, a continuous warming was observed at ∇T w = + 0.035 °C a −1 rate, significant at the α = 0.001 level according to mMK, resulting in an overall 1993-2020 warming of ΔT w = + 0.97 °C. The 200-370 m layer warmed at a ∇T w = + 0.008 °C a −1 rate, significant at the α = 0.05 level for mMK, with ΔT w = + 0.22 °C overall warming for 1993-2020. Within the deep hypolimnion, warming is shaped by complete-mixing events, steeper than the average warming phases occurring in 1993-1999 and from 2006 onwards.
The differential warming which occurred between the epilimnion and the hypolimnion caused a significant increase of the Schmidt stability of the water column (Fig. 5), with linear trend ∇St = + 0.11 kJ m −2 a −1 , significant at the α = 0.01 level according to mMK. This corresponds to ΔSt = + 3.09 kJ m −2 a −1 within 1993-2020.
During the study period 1993-2020, a significant negative trend in the observed 200-370 m volume-weighted concentrations of DO occurred in Lake Maggiore (Fig. 6), significant at the α = 0.05 level for mMK. An unprecedented decrease was observed after 2015, reaching the absolute minimum (4.4 mg L −1 , 37% as saturation) in January 2019, the previous minimum having been reached in December 1998 (5.4 mg L −1 , 45% as saturation), at the end of a depletion phase lasting since 1991, when the preceding full The average volume-weighted DO concentrations in the water column measured at spring maximum circulation showed stable values ranging 8-10 mg L −1 for 1993-2015, while from 2016 onwards concentrations of 7-8 mg L −1 were attained. Despite remaining a well-oxygenated lake throughout the study period, an increase of the volume of Lake Maggiore affected by DO levels < 6 mg L −1 is evident from Fig. 7. We evaluated that these DO concentrations affect in the end about 9.0 km 3 of deep waters, corresponding to almost 25% of the lake overall volume. By comparison, in 1998-1999, when DO concentration reached the previous minimum, only 15% of the lake volume had DO concentrations < 6 mg L −1 .

Tributary intrusions
The intrusion depths of the considered tributaries follow an annual cycle (Fig. 8), the largest depths being attained in the second half of winter, when riverine waters reach their lowest temperatures (Fenocchi et al. 2017). Deep intrusions of the main tributaries have often occurred in such period up to 2014, reaching the bottom for the two main inflows River Ticino and River Toce in winters 2000,2004,2006 and 2012 (Fig. 8b and c). A strong decrease in the intrusion depths is evident since 2015, the 200-370 m layer having been reached only for River Tresa in winter 2018 (Fig. 8d). The only relevant exception to the annual cycle corresponds to the October 2000 severest flood of River Toce (Fig. 8c), which determined a peak daily discharge Q = 2719 m 3 s −1 , reaching an intrusion depth of 330 m due to high suspended sediment concentration, the intrusion depth before and after the event being ~ 20 m. Compared to the other tributaries, River Tresa (Fig. 8d) intrudes deeper in winter, reaching the deep hypolimnion in most of the years up to 2014, and shallower in spring and summer, often intruding at the surface. The deeper winter plunges are due to the steeper intrusion angle than the other tributaries (see Appendix), while the shallower spring and summer insertions are due to its water coming from the epilimnion of Lake Lugano, which in the stratified season is warmer than that of Lake Maggiore (Fenocchi et al. 2017). Despite the exceptionally warm 2007 winter leading to the shallowest intrusion depths within 1993-2020, the last 6 years represent the longest span of the study period in which the plunges of all considered tributaries did not reach the bottom.
Despite the discharge series of the considered tributaries (Fig. 2) generally showing two peaks in spring and autumn due to the climate of the area (Saidi et al. 2013), the actual discharges that reach the intrusion depth display an annual cycle which resembles the one of the intrusion depths (Fig. 8), with maximum values in the second half of winter, due to the prevalence of turbulent entrainment effects. These cause the ratio between daily intruded and inflowing discharges to exceed 10 for all tributaries in winter except for the warmest 2007 and 2020 years, going beyond 10 2 for bottom insertions. The same ratios result for the mass discharges of DO which reach the intrusion depth, being more than two orders of magnitude larger than the inflowing ones from the tributaries when the bottom is reached. Figure 9 shows the total volumes of water and masses of DO intruded every year by the considered tributaries into the four main layers into which we divided the water column coming from the tributaries themselves and the lake layers above, i.e. excluding the entrained water from the same layer. This allows evaluating the amount of alien water and DO being brought annually into the layer of interest by the intrusions of the main tributaries. From Fig. 9, the absence of deep-water renewal and reoxygenation from river plunges after 2014 is self-evident. Comparing the volumes of the considered layers of Lake Maggiore obtained from the hypsometric curve (Table 1) with the intruded volumes from above (Fig. 9a), it can be noticed that: (1) the surface 0-20 m layer was fully renewed in 15 years out of the 28 of the study period, the renewal ratio never dropping below 60% in the remnant years; (2)  For all the considered tributaries, the daily intrusion depths, plotted against daily Schmidt stability (Fig. 10), display a clear inverse correlation. Apart from the October 2000 major flood of River Toce, the deepest intrusions are obtained for the lowest stabilities in winter, whereas the shallowest insertions result under maximum stratification conditions. Averaging these two variables over each day of the year (data for year 2000 were removed for River Toce) reveals for each tributary a hysteretic pattern (Fig. 10). This is due to the swifter warming and cooling rates of rivers compared to the lake. In autumn, tributaries gradually start to intrude deeper due to air temperature cooling, with summer stratification slowly collapsing, while in spring the warming of air temperature determines a quick increase of river temperatures, causing a prompt reduction of intrusion depths.
Another inverse correlation is obtained plotting the maximum annual intrusion depth of each tributary against the minimum annual stability of the lake water column (Fig. 11), ultimately revealing that deepest intrusions are obtained  for years in which the winter residual stability is the lowest. Power-law regressions well describe the link between these two variables, with Pearson correlation coefficients ranging R = − 0.870 to − 0.799 for River Maggia, River Ticino and River Toce (p value < 10 -6 ). For River Tresa, a smaller R = − 0.548 (p value < 10 -2 ) results for the powerlaw regression, being hampered by large intrusion depths being reached also during warmer winters, due to the steeper intrusion compared to the other tributaries (see Appendix).

Field observations
The analysis of long-term data for 1993-2020 showed that Lake Maggiore is experiencing increasing water temperatures (Fig. 4), at different rates between the epilimnion and the hypolimnion, in response to air temperature warming (Fig. 3). The deepest 200-370 m layer is warming at a much slower rate than the surface 0-20 m one, showcasing a "sawtooth" pattern with complete-mixing events almost resetting temporary increases (Livingstone 1997). These differential warming rates are comparable with those of the other deep perialpine lakes (Riffler et al. 2015;Rogora et al. 2018) and are causing increasing stability of the water column (Fig. 5), enhancing the resistance to mixing. Being spring turnovers, both the complete and the incomplete ones, fundamental for DO distribution through the water column, Lake Maggiore is showing decreasing deep-water DO concentrations, with 200-370 m volume-weighted values dropping consistently below 6.0 mg L −1 from 2016 onwards (Fig. 6). Even if DO concentrations below 6.0 mg L −1 were occasionally registered in such layer also in 1977, 1991, and 1998, just prior to complete-mixing events restoring higher values (Calderoni et al. 1997;Rogora et al. 2018), this is the first time that such DO concentrations endured throughout several years, attaining values that are the lowest of the whole recorded data series starting in the 1960s (Calderoni et al. 1997;Rogora et al. 2018). This has been caused in the first instance by the estimated maximum annual convective mixing depths published in the CIPAIS reports (CNR-IRSA 2021) being < 100 m since 2014, leading to the consumption of the DO brought by the relevant partial turnover of winter 2012. The benefit of such event vanished from autumn 2016, when a strong depletion was recorded (Figs. 6 and 7).

Tributary intrusions
The results presented here have shown that under warm winters with strong winter residual stratification, deep intrusions are also suppressed. This occurs as tepid winters lead to both warmer water temperatures of the tributaries and to the entrainment of warmer lake waters in the surface layers, determining shallower intrusions. On the contrary, cold winters, in which very low to null residual stratification endures in the end, promote deep intrusions. Therefore, deep intrusions (Fig. 8) and strong convective mixing act in parallel in Lake Maggiore, significant hypolimnetic DO replenishment by the former events (Fig. 9b) having been calculated for years in which either full turnovers (1999,2005,2006) or relevant partial-mixing events (2000,2004,2012) events occurred. After 2014, river intrusions in the deep hypolimnion ceased (Figs. 8 and 9) together with extended convective mixing, leading to the unprecedented low DO concentrations in the deep hypolimnion. A further proof of the strong connection between deep intrusions and extended convective mixing is given by the positive Pearson linear correlation coefficients between the maximum annual intrusion depths of each tributary and the estimated maximum annual convective mixing depths (CNR-IRSA 2021), ranging R = 0.784-0.837 for River Maggia, River Ticino and River Toce (p value < 10 -6 ) and equal to R = 0.606 for River Tresa (p value < 10 -3 ).
This synergy found for Lake Maggiore must be attributed to the suspended sediment concentration of tributaries attaining low values except for catastrophic floods, such as the October 2000 one of River Toce. Consequently, the contribution of flood suspended sediments to plunging is ordinarily overcome by water stratification. As a matter of fact, a maximum suspended sediment concentration C s = 1801 mg L −1 is predicted for River Toce for the peak daily discharge Q = 2719 m 3 s −1 attained in 2000. Riverine suspended sediment concentration appears, thus, much less relevant for deep intrusions than found for River Rhine flowing into Lake Constance (Fink et al. 2016) and River Rhône flowing into Lake Geneva (Råman Vinnå et al. 2018). The lower suspended sediment concentration of the tributaries of Lake Maggiore compared to those of Lake Constance and Lake Geneva must be attributed to the steeper slopes and shorter paths of the formers, so that larger sediment diameters are involved, being conveyed as bedload.
Increases in air temperature would act also on the discharge regimes of tributaries. In particular, the CH2011 (2011) scenarios predict that during winter less precipitation will be in the form of snow and that freshets will happen earlier. This would lead to an increase in spring discharges and consequently in the suspended sediment concentration of tributaries in a period of low stratification. However, we doubt that this would boost deep intrusions in the future, due to both rising residual stabilities with climate warming and the need for very large discharges to have significant suspended sediment concentrations in the tributaries of Lake Maggiore. Deep intrusion events occurring in seasons other than winter could instead be promoted by the increase of extreme flood events projected with climate change in the Alpine region (Glur et al. 2013). However, the deep-water DO replenishment obtained through impulse flood events would be at best of the same order of magnitude of those which were attained up to the past years in winter with deep plunging of ordinary discharges lasting weeks.

Comparison between the 2000 and 2020 flood events
To shed further light on the relevance of extreme floods for hypolimnetic DO replenishment in Lake Maggiore within the context of climate warming and ongoing deep-water DO depletion, we carried out a further analysis, comparing the October 2000 and October 2020 flood events of River Toce in terms of intrusion depths and intruded DO mass (Fig. 12). The comparison of the diagrams of the intrusion depths of River Toce against the lake Schmidt stability (Fig. 12a) first highlights the differences between the 2000 and 2020 years, largely ascribable to 20-year-long climate warming. The maximum winter intrusion depth dropped from 261 in 2000 to 72 m in 2020. The much smaller anomaly relative to the mean hysteretic pattern (Fig. 10c) of the October 2020 event compared to the October 2000 one also emerges from Fig. 12a.
The values of water discharge and suspended sediment concentration for River Toce were Q = 2719 m 3 s −1 and C s = 1801.4 mg L −1 for the peak day of the 2000 event, i.e. 15 October 2000, and Q = 1940 m 3 s −1 and C s = 910.3 mg L −1 for the peak day of the 2020 event, i.e. 3 October 2020. Schmidt stabilities were St = 22.8 kJ m −2 on 15 October 2000 and St = 29.8 kJ m −2 on 3 October 2020, resulting from the warmer epilimnion and more developed stratification in 2020 (Fig. 12b), in turn ascribed to a lesser extent to the two-week-earlier advent of flood, but mostly to the 20-yearlong climate warming. The temperatures of River Toce were T w = 7.73 °C on 15 October 2000 and T w = 10.50 °C on 3 October 2020. In this case as well, the earlier occurrence of the 2020 event and the inverse relationship between discharge and water temperature which generally holds for mountain rivers would play a part in explaining the 2.77 °C difference, yet a much larger contribution is to be ascribed to elapsed climate warming. Effects of recent extended hypolimnetic DO depletion in Lake Maggiore are last strongly evident in Fig. 12c.
We first ran the river intrusion algorithm for the peak day of the October 2000 flood event employing the lake temperature profile (Fig. 12b) and the inflow temperature of River Toce of the peak day of the October 2020 event. Our calculations show that under such conditions the 2000 flood would have reached an intrusion depth of 152 m instead of 330 m, reducing the intruded mass of DO from 23.9·10 3 to 9.3·10 3 t d −1 , with a 61% decrement. This proves that with climate change, future deep intrusions due to severe floods would surely be hampered by warming river and lake waters and increasing stratification of the latters.
Next, we ran the algorithm again for the peak day of the October 2000 event, employing the DO vertical profile of Lake Maggiore of the October 2020 flood (Fig. 12c) and leaving everything else unchanged. We found out that as an effect of the reduced entrained DO mass, the mass of DO intruded by the 2000 event during the peak day would have decreased from 23.9·10 3 to 17.7·10 3 t d −1 , with a 26% decrement. This highlights that a further future expansion of the hypoxic layer would severely limit the deep-water DO replenishment by the tributaries, turbulent entrainment having been proved fundamental for the process.
Last, we ran the algorithm for the peak day of the October 2000 flood combining the substitutions of the previous two simulations, thus reproducing what would have happened if the October 2000 flood discharge and suspended sediment concentration were released under the October 2020 environmental conditions. In this case, in addition to the reduction of the intrusion depth from 330 to 152 m, the intruded mass of DO would have dropped from 23.9·10 3 to 8.8·10 3 t d −1 , with a 63% decrement. The effect of the exchange of the lake and river temperatures, thus, prevails over that of the DO profile substitution, due to the intrusion depth and turbulent entrainment limitation.

Comparison with other large deep subalpine lakes
From the observations of long-term data for the study period 1993-2020 and from the calculations we performed, we can state that future climate warming would cause deep intrusions to provide Lake Maggiore with a fading deep-water oxygenation, going along with the vanishing oxygenation due to convective mixing (Fenocchi et al. 2020). Both factors are due to the same causes, i.e. water warming and the enforcement of stratified conditions in winter.
These results are coherent with what Fink et al. (2016) found for Lake Constance, for which it was found that, due to climate change, also affecting the discharge regime of River Rhine, future intrusions would reach more frequently the mid-depths of Lake Constance but less frequently its deep layers, enforcing the decline of DO renewal in the deep Fig. 12 Daily intrusion depths of River Toce against daily Schmidt stability for 2000 and 2020 (a the plain end is the beginning of the year, the arrow its conclusion); vertical water temperature (a) and DO (a) profiles of Lake Maggiore for the peak days of the October 2000 and October 2020 floods of River Toce hypolimnion due to decreasing convective mixing. Råman Vinnå et al. (2018) found that for Lake Geneva deep river intrusions would shift from summer towards winter due to a climate-change-induced variation in the hydrological regime of River Rhȏne, balancing out the effect of the reduction of convective mixing due to air temperature warming. The difference of the Lake Geneva case is due to deep intrusions currently occurring there mostly in the summer, with high suspended sediment concentrations brought about by the combination of freshet and rainfall discharges, while much lower discharges occur in winter, due to the 2127 m a.s.l. mean elevation of the catchment of River Rhȏne (Råman Vinnå et al. 2018). In such a basin, air temperature warming would cause a strong decrease in the winter snow amount and the anticipation of thawing, effectively triggering deep intrusions under the least stratified conditions of the lake. This effect would be much less evident in the Lake Maggiore catchment, whose average elevation is 1270 m a.s.l. (Saidi et al. 2013), so that air temperature warming is expected to lead to a milder increase of winter discharges. This same consideration would hold for the other deep lakes south of the Alps, due to the similar features and climatic conditions of their watersheds (Salmaso and Mosello 2010).

Effects of full turnovers on deep-water DO
Having proved that for Lake Maggiore the contribution of riverine intrusions to deep-water oxygenation would extinguish together with convective mixing, we present here the results (Fig. 13) of two simulations of the possible evolution of DO concentration in the 200-370 m deep hypolimnion of Lake Maggiore in the 2020-2085 period from Fenocchi et al. (2020). In these simulations, performed with the 1D coupled ecological-hydrodynamic model GLM-AED2 (Hipsey et al. 2013(Hipsey et al. , 2019, inflows were not considered. The two simulations consider the CH2011 A2_upper worst climate change scenario for Southern Switzerland (CH2011 2011Fischer et al. 2015), which projects a rise of mean annual air temperature in 2085 from 2020 of + 4.1 °C, and differ one another by the random series of input meteorological variables, generated through the Vector-Autoregressive Weather Generator (VG; Schlabing et al. 2014). The hydrodynamic module GLM and the ecological module AED2 were calibrated for Lake Maggiore in Fenocchi et al. (2017) and Fenocchi et al. (2019), respectively. Details on the Lake Maggiore model are available in the mentioned papers. In the simulations displayed herein, the 2006-2015 series of input loads to Lake Maggiore were cycled, representing a credible constant load case (Fenocchi et al. 2020).
In In the second simulation, a single full turnover occurs at the end of winter 2070, representing what could happen if complete mixing happened with an anoxic deep hypolimnion (Fenocchi et al. 2020), brought about by the suppression of deep convective mixing due to severe climate warming (Livingstone 2003;. Such event would bring DO concentrations in the deep hypolimnion back to ~ 4.5 mg L −1 , which would yet drop to < 2 mg L −1 in less than 10 years. The same ~ 4.5 mg L −1 concentration would be attained throughout the water column with the full turnover, ordinary ~ 10 mg L −1 concentrations being restored after ~ 70 days in the 0-20 m epilimnion. Within such period, severe consequences would the expected for the ecosystem (Foley et al. 2012;Perga et al. 2015;Missaghi et al. 2017). In Lake Lugano, which has an anoxic deep hypolimnion, the sudden drop in surface DO concentration after the full turnover occurred in 2006 caused substantial fish death (Holzner et al. 2009;Rogora et al. 2018). Nutrients, notably phosphorus, which have been accumulating in the hypolimnion (Rogora et al. 2021) would also become available in the surface layers, triggering phytoplankton blooms and critical effects on water quality, lasting till the disposal of the excess load.

Conclusions
The analysis of 1993-2020 long-term data of Lake Maggiore has revealed that the lake is experiencing increasing water temperatures and an intensifying stability of the water column. Due to the high relevance of spring turnovers for DO distribution across the water column, Lake Maggiore has been recently affected by a steady decrease of DO concentrations in the deep waters. An in-depth analysis on the relevance of tributary insertions for hypolimnetic DO replenishment has revealed that in the last years maximum annual intrusion depths have decreased together with convective mixing ones, strengthening hypolimnion separation. Therefore, they would not likely be able to support hypolimnetic oxygenation with future climate change, not even during exceptional floods with high suspended sediment concentrations.
Climate warming and its effects on the thermal stability of the water column and on the intrusion depths of the tributaries hence drive the DO levels in the lake. Soon, the overall ecological state of the lake will be strongly driven by climate change and management plans will be needed to sustain water quality and avoid cancelling the past efforts to contrast eutrophication.

Data availability statement
The datasets analysed and generated during the current study and the MATLAB® code of the inflow intrusion algorithm are available from the corresponding author on reasonable request.

Appendix: Details on the inflow intrusion algorithm
The Chen and Millero (1986) equation of state is employed to evaluate water density within the inflow intrusion algorithm, including the effects of temperature, salinity and pressure. The density of the mixture of water and suspended sediments ρ m is computed from the water density ρ w as: where C s is the suspended sediment concentration and ρ s is the density of the suspended sediments, assumed as ρ s = 2650 kg m −3 due to the siliceous watershed (Morabito et al. 2012).
Saturation conditions were assumed for the DO concentrations of the tributaries due to the turbulence of lotic waters (Fenocchi et al. 2019). Saturation concentrations for the tributaries were determined through the Mortimer (1981) formula, as function of the inflow temperature and of the standard barometric pressure at the reference lake elevation 194.0 m a.s.l. Considering the generic tributary i with temperature T i,0 in [°C], the saturation concentration of DO at the lake entrance is computed as: in which the last term inside the square brackets expresses the natural logarithm of the standard barometric pressure in [atm] at an altitude of 194.0 m a.s.l.
Following the original DYRESM schematisation, the entrance of the tributary into the lake is modelled as a triangular channel continuing below the surface, α being half the bottom angle of the triangular section and φ being the slope angle of the shore in the entrance region (a schematic is available in Hipsey et al. 2019). These angles were estimated for the considered tributaries from the lake bathymetry, being α = 86.97° and φ = 5.49° for River Maggia, α = 86.69° and φ = 3.59° for River Ticino, α = 84.58° and φ = 3.74° for River Toce, α = 52.54° and φ = 12.78° for River Tresa. The steeper slope angle of River Tresa is due to its insertion in the deepest central basin of Lake Maggiore (Fig. 1b). A drag coefficient C D = 0.016 was adopted to represent bottom friction for the four throughflows (Hipsey et al. 2019).
The Richardson number Ri and the entrainment rate E are first calculated for the generic tributary i as function of the drag coefficient and of the insertion geometry (Hipsey et al. 2019): The density of the tributary and of the lake surface are then compared. If the inflow is lighter, it enters on top of the water column and entrainment is not computed. All plume variables, i.e. discharge Q, temperature T, salinity S, suspended sediment concentration C s and dissolved oxygen concentration DO, are hence the same as those of the inflow. Intrusion of the tributary otherwise occurs, its initial thickness being calculated as: where Q i,0 is the initial discharge of the inflow and g' i,0 is its reduced gravity: in which g = 9.806 m s −2 is the gravity acceleration at the latitude of Lake Maggiore and ∆ i,0 is the non-dimensional density anomaly of the tributary with respect to the lake surface, ρ i,0 being the density of the inflow and ρ l,0 the density of the lake surface. The intrusion process is iterated along its spatial development, potentially down to the lake maximum depth of 370 m, employing a ∆z = 0.01 m vertical resolution. All plume variables are initialised at the inflow values. The thickness of the inflowing plume at generic step j is computed from that at previous step j-1 as: where Δs = Δz∕sin i is the distance along the lake shore covered by the intruding plume at each step. A more accurate algorithm, considering the slope variation of the shore along the presumed intrusion path of the tributary, and hence a variable entrainment rate, is under development and will be the object of a future paper. The plume discharge at generic step j is given by: so that the increase in discharge from previous step j-1 is: The other plume variables T, S, C s , DO are updated at each step through mass balances between the plume and the entrained lake water. For C s , sedimentation between the previous step and the current one is accounted if the plume velocity at the previous step U i,j-1 < U c = 0.46 m s −1 , the latter being the critical settling velocity for medium silt given in Mulder et al. (1998) and employed in Fink et al. (2016) andRåman Vinnå et al. (2018). The plume velocity at the previous step is calculated as (Hipsey et al. 2019): The decrease in suspended sediment concentration due to sedimentation between the previous step and the current one is computed as (Syvitski and Lewis 1992): in which r = 4.7 day −1 is the sediment removal rate for medium silt in Mulder et al. (1998), employed in Råman Vinnå et al. (2018, and Δt i,j−1 = Δs∕U i,j−1 is the time taken by the plume to cover the sloping distance between the two steps. The suspended sediment concentration of the plume which enters the mass balance at the current step is therefore: A corrected density * i,j−1 is then calculated from Eq. (1) using the corrected C s * i,j−1 . Considering the generic plume variable k, its value at generic step j is obtained balancing the value at the previous step k i,j-1 ( C s * i,j−1 is used for the suspended sediment concentration), weighted by the corresponding mass flux Q i,j−1 * i,j−1 , with the value for the entrained lake water at the present step k l,j, weighted by the mass flux increment from the previous step ΔQ i,j l,j : The density at the current step ρ i,j is then computed using the calculated T i,j , S i,j and C si,j . The resulting DO concentration is checked never to overcome the absolute saturation concentration, obtained from the Mortimer (1981) formula considering the actual pressure at depth z j in [m] (Wetzel 2001): The iteration of the process requires the continuity between consecutive steps of Q i and of the internal Froude number (10) U i,j−1 = Q i,j−1 h i,j−1 2 tan i .
(11) ΔC s i,j−1 = rC s i,j−1 e −rΔt i,j−1 Δt i,j−1 , Fr i = U i ∕ √ g � i H , with g � i = gΔ i = g i − l ∕ l being the reduced gravity between the plume and the lake and H the total lake depth. For this to occur, the thickness of the throughflow to be used as initial thickness at the next step j + 1 is corrected to (Pilotti et al. 2014): where Δ i,j and Δ * i,j are the non-dimensional density anomalies of the plume with respect to lake water at depths z j and z j+1 , respectively.
The intrusion loop of the algorithm ends when either the bottom of the lake is reached, in which case the inflow intrudes below the lake water column, or when the lake density at the next step ρ l,j+1 is larger than that of the plume ρ i,j , in which case the latter enters at the depth of the current step.