Tributaries affect the thermal response of lakes to climate change

. Thermal responses of inland waters to climate change varies on global and regional scales. The extent of warming is determined by system-speciﬁc characteristics such as ﬂuvial input. Here we examine the impact of on-going climate change on two alpine tributaries, the Aare River and the Rhône River, and their respective downstream peri-alpine lakes: Lake Biel and Lake Geneva. We propagate regional atmospheric temperature effects into river discharge projections. These, together with anthropogenic heat sources, are in turn incorporated into simple and efﬁcient deterministic models that predict future water temperatures, river-borne suspended sediment concentration (SSC), lake stratiﬁcation and river intrusion depth/volume in the lakes. Climate-induced shifts in river discharge regimes, including seasonal ﬂow variations, act as positive and negative feed-backs in inﬂuencing river water temperature and SSC. Differences in temperature and heating regimes between rivers and lakes in turn result in large seasonal shifts in warming of downstream lakes. The extent


Introduction
The thermal and hydrodynamic responses of lakes to climate change are considerably diverse. Observed responses vary on global, regional and even local scales (O'Reilly et al., 2015). Even neighboring freshwater bodies can react differently to a given increase in air temperature. This indicates that lake-specific characteristics will determine the response to climate change (for clarity and brevity, we refer to anthropogenic climate change simply as "climate change" or "climate" from now on). Local factors which affect climate warming of lakes include, among others, morphology (Toffolon et al., 2014), irradiance absorption (Kirillin, 2010;Williamson et al., 2015), local weather conditions (Zhong et al., 2016), stratification , atmospheric brightening (Fink et al., 2014a) and ice cover (Austin and Colman, 2007).
Throughflows affect epilimnion and hypolimnion temperatures of lakes. Studies of climate impact typically do not address these sorts of subtleties in lake dynamics due to a lack of data or difficulties in predicting future temperature and discharge conditions (Fenocchi et al., 2017). Several studies of large lakes suggest that major tributaries play only a minor role in climate-induced warming and deep-water oxygen renewal (Fink et al., 2014a;Schwefel et al., 2016). Mediumand smaller-scale lakes are, however, more abundant than large lakes (Verpoorter et al., 2014) and exhibit a greater degree of sensitivity to point sources of anthropogenic thermal input which can affect temperature and stratification (Kirillin et al., 2013;Råman Vinnå et al., 2017). Medium-and L. Råman Vinnå et al.: Tributaries affect the thermal response of lakes to climate change small-sized lakes also make a more significant contribution to the temperature-dependent global greenhouse gas budget (Holgerson and Raymond, 2016). Accurate prediction of climate change impacts therefore requires a more detailed understanding of small-to medium-scale lake and tributary systems.
Climate change exerts a dual influence on alpine rivers by introducing variation to both flow and temperatures. Discharge variation takes the form of less flow in summer and more flow in winter due to warmer high-altitude snow and ice melt/runoff regimes (Addor et al., 2014;Birsan et al., 2005), which also influence river temperature (Isaak et al., 2012;Van Vliet et al., 2013). Increased air temperature may also enhance erosion rates in river basins thereby supplementing river-borne suspended sediment loads (Bennett et al., 2013). River temperature and suspended sediment content determine water density and, by extension, the depth of river plume intrusions into downstream lakes or reservoirs. The depths and volumes of river intrusion plumes determine deep-water oxygen renewal, especially for deeper lakes where climate-related warming can reduce seasonal deep convective mixing and thereby deplete deep-water oxygen (Schwefel et al., 2016). Major (deep penetrating) river intrusion events typically occur due to flooding, which flush large sediment loads into the river (Fink et al., 2016). The frequency and volume of floods in the Alps are notoriously hard to predict, although a decrease in floods has occurred in association with recent warmer summers observed in the Alps (CH2011, 2011; Glur et al., 2013).
Recent model studies have identified inland waters as risk hotspots under expected climate change scenarios (IPCC, 2014). These systems require a more detailed analysis given their role in supporting fisheries, agriculture, drinking water supply, heat management and hydropower. This paper examines the complex interactions between tributaries and lakes in response to temperature increase and other modifications expected from climate change. Our objectives were to quantify the impact of specific climate change scenarios on (i) alpine tributaries and (ii) downstream peri-alpine lakes with a focus on river-borne suspended sediment concentration (SSC), water temperature, stratification and river intrusions.
We used coupled river-lake models to build on previous research by Fink et al. (2016). These authors investigated the effect of flood frequencies on deep-water renewal under established climate change scenarios. Their work did not generate tightly constrained estimates for flooding events. Our analysis therefore provisionally assumed that flood frequency does not change in the future. In addition to these sources of natural variation, our models addressed variation in river discharge regimes (i.e., daily mean level shift) under the specified A1B climate change scenario. These in turn affect SSC and thermal regimes for rivers and their associated downstream lakes. Furthermore, here we show that local point sources of anthropogenic thermal pollution can have discharge Regional climate models (RCMs) Lake model temperature Figure 1. Schematic illustration of the one-way model chain of this study. Orange models represent modeling performed by this study, while grey models represent simulated data inputs obtained from external sources. a major impact on the response of inland waters to climate change as previously suggested by Fink et al. (2014b).

Approach
The investigation of tributary influence on lake response to climate change followed these procedural steps: i. Define river temperature and SSC models for two major alpine rivers and designate a one-dimensional lake model for a large-and medium-sized peri-alpine lake.
ii. Integrate model (i) with a river intrusion scheme: Fig. 1 shows the integration of the one-way component models.
iii. Obtain and apply estimates of future regional air temperature, tributary discharge and changes in local anthropogenic thermal emissions to both river and lake models.
iv. Identify patterns in model outputs of water temperature, SSC, lake stratification and river intrusion parameters (volume and depth).
LG is a large meso-eutrophic lake resting at 372 m elevation and covering an area of 580 km 2 . It reaches a maximum depth of 309 m and holds a volume of 89 km 3 with an average hydraulic residence time of 11.5 years. Complete seasonal deep convective mixing only occurs on average every Figure 2. Study area and predicted regional air temperature increases. Elevation above sea level (green to white color ramp), locations and number of river stations (white diamonds) and atmospheric monitoring stations (red triangles), drainage area (Aare: vertical lines; Rhône: horizontal lines) and location of Mühleberg Nuclear Power Plant (MNPP, orange circle). Area covered by regional climate models (dark red dashed-dotted line) with (a) predicted air temperature increase T in the near-future (blue, 2030-2049) and far-future (orange, 2080-2099) for medium (thick lines) and upper/lower estimates (thin lines) under the A1B emission scenario (CH2011, 2011). fifth winter but is predicted to become less frequent with ongoing climate change (Perroud and Goyette, 2010;Schwefel et al., 2016). Whereas the global average lake surface temperature has increased by ∼ 0.34 • C decade −1 between 1985and 2009(O'Reilly et al., 2015, the Rhône River supplying ∼ 75 % of LG's inflow has experienced a temperature increase of ∼ 0.21 • C decade −1 from 1978 to 2002 (Hari et al., 2006).
LB is a 74 m deep, meso-eutrophic, medium-sized lake resting at an elevation of 429 m. It covers a surface area of 39.3 km 2 and holds a volume of 1.18 km 3 with a hydraulic residence time of 58 days. Complete deep convective mixing occurs every winter and effectively replenishes the oxygen-depleted deep-water. The Aare River provides ∼ 61 % of LB's inflow and experienced a 0.34 • C decade −1 increase in temperature from 1978 to 2002 (Hari et al., 2006). Several dams/lakes trap sediment along the upstream Aare course and increase sediment settling and water temperature prior to entering LB. The Mühleberg Nuclear Power Plant (MNPP), situated ∼ 19 km upstream from LB (46 • 59 N, 7 • 16 E; Fig. 2) represents a point-source of thermal pollution. Planned for decommission in 2019, the plant emits ∼ 700 MW of heat into the Aare and substantially warms the river water (Råman . The ∼ 8 km long Zihlkanal, LB's second largest tributary, supplies ∼ 32 % of the lake inflow and connects LB to Lake Neuchâtel (Fig. 2). This tributary is neglected here since it mainly transports lake surface water, which has approximately the same temperature as LB surface water and thus without net heat effects.

River models 2.3.1 Temperature
Uncertainties concerning river morphology, heat fluxes, shadowing and atmospheric conditions such as wind speed and cloudiness (Caissie, 2006) pose a significant challenge to accurately model future river temperatures. Deterministic models typically require detailed knowledge unavailable for future climate scenarios. Regressions and stochastic models rely heavily on observed natural variability in a given time frame and typically do not include inputs representing additional or interacting physical processes. On their own, these sorts of "black box" models cannot balance trade-offs between constraints available from empirical data and the complexity offered by theoretical frameworks.
To overcome these limitations, we used the hybrid model air2stream . The model combines the simplicity of stochastic models with accurate representation of the relevant physical processes affecting temperature. Similar to the neural networks approach, the model calculates river water temperature (T w ) through a Monte-Carlo-like calibration process, which identifies optimal parameters for weighting physically dependent variables. We used the eight-parameter (a 1 to a 8 ) version of the model which incorporates air temperature (T a ) and river discharge (Q) as a function of time (t).

Suspended sediment concentration
Water density and intrusion depth of river water into downstream lakes is influenced by SSC. Intensive flow events create high levels of SSCs (Rimmer and Hartmann, 2014), as can exposure/erosion of sediment sources within the river basin through the so-called hysteresis effect, in which SSC varies for the same level of discharge (Tananaev, 2012). River  -1999; b 2000-2009 discharge regimes have been predicted to change in the future (Birsan et al., 2005), suggesting that SSCs will also change.
To simulate future SSCs, we used the supply-based rating model described in Doomen et al. (2008), which Fink et al. (2016) adapted to the River Rhine. The model consists of a base level SSC (g m −3 ) function expanded to express erosion of sediment at high discharge and sediment accumulation at low discharge. The model is expressed as where b x , c x and m are adjustable parameters in combination with the threshold discharge (Q th ), which determines whether erosion or deposition occurs within the river. The parameters d 1 and d 2 control the deposition of (or erosion from) the river sediment storage (ψ (g)).
Erosion occurs if Q exceeds Q th and the river basin contains erodible sediment (ψ > 0). Sedimentation occurs if Q is smaller than Q th . The change in ψ over time can be formulated as Parameters in Eqs.

Lake model
We used the one-dimensional model SIMSTRAT (Goudsmit et al., 2002) to assess the impact of climate change on temperature and deep-water renewal in LB and LG. The model calculates heat fluxes and vertical mixing driven by wind and the internal wave field using a k − ε turbulence closure scheme. It has been adapted to and validated for multiple lakes including Lake Zürich , LG (Perroud and Goyette, 2010;Schwefel et al., 2016), Lake Neuchâtel , Lake Constance (Fink et al., 2014b;Wahl and Peeters, 2014) and LB (Råman Vinnå et al., 2017).    Table 3. One-dimensional lake model SIMSTRAT best-fit parameters and model performance statistics reported as vertical volumeweighted averaged root mean square deviation (RMSD-V).
Parameter (unit) Lake Biel Lake Geneva The model contains seven tunable parameters, including p 1 (irradiance absorption), p 2 (sensible heat flux) and K (vertical light absorption) for heat flux adjustments from the atmosphere to the lake. Momentum and kinetic energy transfer from the wind to internal waves is tunable by C 10 (wind drag). The internal seiche energy balance can be adjusted through α (production), C Deff (loss by bottom friction) and q (vertical distribution of turbulent kinetic energy). To include the effect of seasonally varying stratification strength, we followed Schwefel et al. (2016) and varied α: α S for summer (April to September) and α W for winter (October to March), where α S > α W . Here we used the best-fit parameter setup (Table 3) already established and validated for LG and LB by Schwefel et al. (2016) and Råman .
Building upon the model developed by Råman , we introduced an extended river intrusion scheme described in Appendix A1 (including sensitivity analysis). This scheme was chosen in order to include the effect of steep bathymetry on plume entrainment. Additionally, the robustness and simplicity of the intrusion scheme limits the uncertainty associated with more complex intrusion models including multiple parameters which can be hard to predict in the future. The entrainment of lake water into plunging underflows was modeled as proposed by Akiyama and Heinz (1984) with additional sedimentation of suspended load (Mulder et al., 1998;Syvitski and Lewis, 1992). The method addresses the transition of a homogenous open channel flow to a stratified underflow where entrainment and settling of sediment depend on bottom slope angle. The model scheme consists of (i) the homogenous region, where river water extends from the surface to the lake bed; (ii) the plunging region, where the plume separates from the lake surface and (iii) the underflow region, where the plume descends downslope while entraining surrounding water until it sepa-rates from the bottom and intrudes into the lake interior (Fink et al., 2016).

Data, hydrology and climate forcing
The models described above used hourly resolved data from 1989 to 2009 as inputs. For calibration/validation of river temperature, we used flow and temperature data from the Aare monitoring station no. 2085 ( Fig. 2 Fig. 2; 46 • 20 N, 6 • 55 E) for Rhône, provided air temperature data. Due to insufficient representation of high turbidity events, we calibrated/validated the SSC model with turbidity data converted to SSC with suspended sediment samples from 2013 and 2014.
The meteorological data used for SIMSTRAT included air temperature, vapor pressure, wind speed, solar radiation and cloud cover. These data were collected from the meteorological stations Cressier (no. 6354 in Fig. 2   and Schwefel et al. (2016) provide additional information on climate data inputs to the one-dimensional model. The river intrusion scheme requires as input the slope angle traveled by the river underflow, which was obtained from a 25 m resolved digital height model (DHM25). Vertical temperature profiles, sampled at the deepest location of both lakes in January 1989, were used as initial conditions.
Van Vliet et al. (2013) suggested that river discharge and air temperature should be used while predicting future river temperatures. We incorporated recent findings of climateinduced changes in air temperature and river discharge regimes to model both future river temperature and SSC. Seasonal mean predictions for air temperature increase in western Switzerland (Fig. 2) were estimated from CH2011 (2011) for the A1B emission scenario (balanced use of renewable and fossil fuels) using results from 20 regional climate models. Flow projections were obtained from published results generated by the PREVAH (PREcipitation-Runoff-EVApotranspiration HRU Model) hydrological model (Viviroli et al., 2009) using a gridded configuration as described in Speich et al. (2015) and Kobierska et al. (2011). The model explicitly incorporates changes in glacial extent, snow melt, catchment runoff, floods and low water flows (FOEN, 2012;Bosshard et al., 2013;Speich et al., 2015). The PREVAH outcomes for the 1981-2009 period have been validated with data from 65 river gauges (Speich et al., 2015), including the two gauges upstream of LG (Rhône; no. 2009 in Fig. 2) and LB (Aare; no. 2085 in Fig. 2) used here.

Scenarios
Six different model scenarios were used to propagate climate change effects through the major tributaries and their associated downstream lakes. Model scenarios LG1 to LG3 represented LG while LB1 to LB3 represented LB (Table 4). Each scenario includes three time periods: a reference period (1990-2009), a near-future period (2030-2049) and a far-future period (2080)(2081)(2082)(2083)(2084)(2085)(2086)(2087)(2088)(2089)(2090)(2091)(2092)(2093)(2094)(2095)(2096)(2097)(2098)(2099). The 20-year intervals allowed us to resolve natural variations on seasonal and shorter time scales. We initialized the models 1 year prior to the investigated period for each time frame (1989, 2029 and 2079) in order to remove effects of initial conditions. Scenarios LG1 and LB1 excluded river inflow in order to isolate lake response to climate change from potential tributary influence. Scenarios LG2, LG3, LB2 and LB3 were used to differentiate between the effects of tributary temperature and SSC, and to provide model sensitivity estimates. The LB3 scenario excluded MNPP thermal pollution from near-future and far-future time periods but not from the reference period. The LB2 scenario included thermal pollution in modeling river water temperature. Scenarios LB2, LB3 and LG3 included SSC while LG2 did not. Low SSC values found in the Aare data resulted in negligible differences between models including and excluding SSC. Because they served primarily validation and sensitivity analysis purposes, the Aare/LB model results excluding particles and including/excluding MNPP thermal pollution (LB4 and LB5) are relegated to Appendix Fig. B1 and not discussed further. Scenarios LG3 and LB3 represent expected future developments.
The unmodified air temperature and modeled river discharge, temperature and SSC were used as inputs for the reference periods. Near-future and far-future models incorporated predicted changes in air temperature and river discharge, temperature and SSC with maximum, medium or mean, and minimum values serving as envelopes for each parameter (Figs. 2a and 5). This strategy gave nine simulations (three for scenarios LG1 and LB1, which exclude rivers; i.e., no variation in discharge nor river temperature) for each nearfuture and far-future time period. Predicted results included a total of 87 model runs. Upper, mean and lower impact estimates (described and interpreted below) were derived from the nine basic model runs.

Rivers
The seasonality of predicted river discharge (Q) from FOEN (2012) varies with respect to the reference period 1990-2009 ( Fig. 5a and b). The PREVAH model shows a future decrease in mean summer discharge (1 April to 30 September) for both the Aare (−3.7 m 3 s −1 decade −1 ; no. 2085) and Rhône Table 4. Model scenarios of climate change effects for near-future and far-future time periods, including (Inc.) and excluding (Exc.) the effects of rivers and suspended sediment. Thermal input from MNPP was also included/excluded. Most likely scenarios are shown in bold.

Lake
Exc. rivers Inc. rivers Exc. suspended sediment Inc. suspended sediment Geneva LG1 LG2 LG3 The decrease in summer will be compensated by an observed increase in winter flow (1 October to 31 March) of the Aare (+3.3 m 3 s −1 decade −1 ) and Rhône (+3.7 m 3 s −1 decade −1 ). These results confirm previous findings presented by Addor et al. (2014) and Bosshard et al. (2013). Regional air temperatures from the A1B emission scenario (∼ +0.32 • C decade −1 ; CH2011, 2011; Fig. 2a) cause a predicted increase in mean annual water temperature (T ) for both the Aare (∼ +0.10 • C decade −1 ) and the Rhône (∼ +0.08 • C decade −1 ). Both rivers experience seasonal variations in temperature increase similar to that predicted for air temperatures (Figs. 2a, 5e and f). The effect is strongest for the Aare during summer with warming of up to +2.5 • C in water temperatures for the far-future time period relative to the reference period.
Thermal pollution from the MNPP in the Aare during the reference period (blue-green line in Fig. 5e; Råman  causes approximately twice as much heating in winter relative to warming from climate change in the far-future. In summer, the relationship reverses with minor MNPP warming relative to that induced by climate change. The net effect of climate warming and MNPP decommission (i.e., removal of MNPP heat from near-future and far-future time periods) on the Aare is cooling in winter and warming in summer relative to the reference period (Fig. 5c).
Like river temperatures, SSCs depend on river discharge. Our model therefore show SSC increasing in winter and decreasing in summer due to shifts in the discharge regime ( Fig. 5g and h). The model results for the Rhône exhibit a mean seasonal increase of +14 g m −3 decade −1 in winter and a decrease of −11 g m −3 decade −1 in summer. For reasons explained above (Sect. 2.2), results for the Aare show less variation, with a seasonal increase of +0.3 g m −3 decade −1 in winter and a decrease of −0.4 g m −3 decade −1 in summer. Altered temperature and SSC caused increases and decreases in water density for both rivers in winter and summer, respectively.

Lakes
Warmer air temperatures (Fig. 2a) predicted from climate change resulted in temperature increases in both LG and LB for all scenarios (Table 5). Models showed the highest warming rates in the epilimnion, intermediate values throughout the metalimnion and the lowest rates in the hypolimnion (Table 5). We defined the epilimnion, metalimnion and hypolimnion using the water column stability method described in Råman . The predicted warming of LG varied only slightly among the three different scenarios Table 5. Change in temperature, length of the stratified period and depth of the thermocline (negative values correspond to a shallower thermocline) for each scenario listed in Table 4. Estimates given as mean of the daily difference between the reference period and the far-future time period. Temperature anomalies are volume-weighted and vertically averaged. Most likely scenarios are shown in bold.

Temperature
Stratification Thermocline Scenario  Fig. 6a-c). Predicted warming of LB depends strongly on the scenario used ( Fig. 6d-f). Similar to the predicted warming patterns for rivers (Sect. 3.1), both lakes showed seasonally varying warming patterns. Reduced warming corresponds with periods of high river discharge ( Fig. 5a and b). This cooling effect occurs primarily in winter and midsummer, and focussed in depth to the level of river intrusion (Figs. D1b, d and 7c-f). Model results showed a greater degree of fluctuations of the warming in LB than in LG. This probably results from the greater influence of the Aare on LB compared to that of the Rhône on LG, as LG has a longer hydraulic residence time. Scenario LB1, which excludes river intrusion, showed only limited seasonal variation in warming ( Fig. C1c and e). According to these results, the closure of MNPP could offset climateinduced warming of LB by ∼ 25 %.
Model results show that enhanced warming of the epilimnion relative to the hypolimnion strengthens stratification ( Fig. 7g and h). This enhances the duration of stratification (for both lakes ∼ +2 days decade −1 ; Table 5) and slightly lifts the thermocline (in LB ∼ −0.1 m decade −1 and in LG ∼ −0.05 m decade −1 ; Table 5). We used the Schmidt (1928) Equation (7) incorporates gravity (g = 9.81 m s −2 ), depth (z), lake surface area (A 0 ), horizontal cross section area (A(z)), lake density (ρ(z)), maximum depth (z max ), mean lake density (ρ m ), lake volume (V ) and volumetric mean depth (z m ) defined as  Figure 6. Temperature increase T for near-future (blue) and farfuture (orange) time periods relative to reference period temperatures, displayed as mean (columns) and standard deviation (black bars) calculated from the nine basic model runs in the near-future and far-future scenarios. Epilimnion (left pair of columns), metalimnion (middle pair) and hypolimnion (right pair) in LG (a to c) and LB (d to f). Graphs represent river intrusion excluded (a and d), river-borne SSC included (c, e and f) and excluded (b). Mühleberg Nuclear Power Plant (MNPP) heat release included in (e) and excluded in (f) from near-future and far-future time periods but retained for the reference period.
The duration of stratification was determined by counting the days when temperature differed by more than 1 • C between surface (2 m depth) and deep-water (280 m for LG; 50 m for LB) (Foley et al., 2012). The maximum water column stability expression N 2 = −(g/ρ) ρ(z)/ z (s −2 ) was used to determine the thermocline depth.  (black, 1990-2009). Anthropogenic MNPP heat input entering LB has been excluded from near-future and far-future time periods but retained for the reference period. Temperature T (a and b), increase in temperature T in epilimnion (c and d) and hypolimnion (e and f) as well as increase in stability S (g and h).
The river intrusion depth is dependent on water density (temperature and SSC are dominant; dissolved solids are negligible). The Rhône is colder (Fig. 5c and d) and carries more suspended sediment ( Fig. 5g and h) than the Aare. Reference period results showed that the Rhône intruded in LG at greater depths relative to depths of the Aare intrusion into LB (Figs. 8 and D1). Given the future change in river temperature and SSC, intrusion patterns will thus change as the densities of both the Aare and Rhône increase and decrease during respective winter and summer seasons (Sect. 3.1). This explains model results showing respective deeper and shallower intrusions during winter and summer for both rivers (Fig. D1).
Model results show that warming of the Rhône generally diminishes the amount of river water penetrating beyond 200 m depth in LG (Fig. 8a). Enlarged river flow in winter enhances SSCs and counteracts heating, thereby increasing the amount of river water intruding beyond 200 m depth (nearfuture ∼ 30 %; far-future ∼ 65 %; Fig. 8b). The difference in winter heating for the Aare and LB epilimnion (Figs. 5c, e and 7c) generally increased the amount of water penetrating into the hypolimnion (Fig. 8c). Decommission of the MNPP enhances temperature differentials between LB and the Aare, thereby increasing the amount of water reaching past 30 m in LB (near-future + ∼ 80 %; far-future + ∼ 120 %; Fig. 8d). In summary, the change in river discharge regime for the Aare and Rhône results in respective increase and decrease in winter and summer water density, resulting in a summer to winter shift of the amount of river water penetrating deeper than the metalimnion for both lakes.

Rivers
Increases in air temperature expected from climate change modify tributary runoff. Less water is predicted to be bound in snow and ice at high elevation during winter and spring/summer floods will occur earlier (CH2011, 2011FOEN, 2012). The changed river discharge regime, appearing as increased flow in winter and decreased flow in summer ( Fig. 5a and b), amplifies the increase and decrease in river temperature during respective summer and winter periods ( Fig. 5e and f). Amplification results from (i) a smaller flow volume requiring less energy to heat and (ii) lower flow velocities which extend heat exposure. The PREVAH model predicts that the future discharge of the Aare in summer will be ∼ 20 % less than summer discharge in the Rhône. These results therefore suggest that the Aare summer conditions will be more impacted by climate change than the Rhône summer conditions (Fig. 5e and f).
Model results concerning discharge-dependent responses to climate-induced warming were consistent with previous findings reported by Isaak et al. (2012) and van Vliet et al. (2013). The river temperature increases predicted by this study (+0.10 • C decade −1 for the Aare and +0.07 • C decade −1 for the Rhône) were much smaller than past observed warming rates (0.34 • C decade −1 for the Aare and 0.21 • C decade −1 for the Rhône; Hari et al., 2006). These differences may reflect contrasting reference periods with past observations conducted from 1971 to 2001 and modeled observations addressing 1990 to 2099. Past observations also incorporate effects of solar brightening during the 1980s (Fink et al., 2014a;Sanchez-Lorenzo and Wild, 2012;Wild et al., 2007), which led to additional warming of air and water.
Climate change effects aside, MNPP decommissioning in 2019 is predicted to decrease the temperature in the Aare by up to 4.5 • C at station no. 2085 (Råman . The cooling effect of this plant closure primarily affects winter conditions when climate-change-induced warming is weaker and river flow is lower (Fig. 5e). The heating of the Aare and LB by MNPP heat emissions equates to approximately 1 decade of climate-induced warming of lake surface waters (O'Reilly et al., 2015;Råman Vinnå et al., 2017). This result highlights the role of point source thermal contributions in local climate impact assessments.
The amount of suspended sediment carried by rivers depends on both discharge and the amount of erodible sediment in the watershed (Fink et al., 2016). We used a supply-based sediment rating model subjected to a changing discharge regime to examine seasonal changes in suspended sediment for both the Aare and Rhône (Fig. 5g and h). Consistent with previous findings reported by Pralong et al. (2015), we predict an increase in SSC during winter and decrease in SSC in summer. This is caused by two phenomena associated with increased river discharge: (i) amplified river bed erosion linked to increased intensity of high discharge events carrying enhanced volumes of SSC and (ii) increase in the sediment available for erosion in the river catchment due to enhanced supply at low flow velocities. Figure 4 and Table 2 show that the SSC model gives robust results for the Rhône (coefficient of determination R 2 = 0.68 from 2013 to 2014) but not for the Aare (R 2 = 0.06 from 2013 to 2014). The Aare includes several sediment-trapping reservoirs/lakes upstream of station no. 2085. Thus, peaks in SSC at station no. 2085 do not reflect watershed-scale discharge events (Fig. 4) but rather local precipitation and discharge events in the headwaters of a tributary (Saane River) to the Aare (Fig. 2). This tributary hosts few sediment traps and contributes ∼ 34 % of the downstream flow at station no. 2085. Given the limited impact of SSC on the Aare water density, models show only negligible impact on river intrusion depth and corresponding intruding volumes (Figs. 8c, B1c, e and D1c). The lower reaches of the Rhône are not dammed, thus adhering more directly to model assumptions and giving clearer results (Fig. 4). High SSC events are usually associated with extreme floods (Fink et al., 2016), which are predicted to vary in alpine lake catchments with on-going climate change (Glur et al., 2013). The lack of constraints on extreme precipitation events introduces uncertainty into future flood frequency and magnitude predictions (CH2011, 2011). Shifts in river discharge regimes also depend on the amount of water bound in snow and ice as well as on the timing of spring/summer melt. Future climate scenarios predict that ∼ 30 % of the glacier mass will remain in the Aare and the Rhône catchments by the end of the 21st century (FOEN, 2012). Glacial meltwater is thus expected to continue to supply the Aare and Rhône throughout the time frames considered in this study. We thus assumed that the flood frequency remained unchanged, while the amplitude of the floods was adjusted in the future according to river discharge regime shifts predicted by FOEN (2012).

Lakes
All model scenarios showed that increased air temperature leads to warming of both lakes, especially of the epilimnion (Table 5, Fig. 6). Piccolroaz et al. (2015) showed that an increase in lake stability and earlier onset of stratification causes warming of surface waters due to the smaller volume undergoing warming and diminished heat transfer to the hypolimnion. The lake model used here showed an increase in stratification strength and a lengthening of the stratified period in both lakes (Table 5; Fig. 7g and h). Our results thus consistently support previous findings for LG reported by Foley et al. (2012), O'Reilly et al. (2015 and Schwefel et al. (2016).
Seasonal variations in warming of both epilimnion and hypolimnion ( Fig. 7a-f) surpassed the seasonality of applied changes in air temperature (Fig. 2a). The model showed a decrease in warming during winter and midsummer, which corresponds to time periods of high river discharge from the main tributaries ( Fig. 5a and b). This cooling effect was more effective for LB than for LG (Fig. 7) and appeared in all scenarios except for LB1 and LG1 (Fig. C1), both of which exclude coupled river effects. The extended seasonal variation in climate warming is thus driven by river discharge volume and temperature trends (Figs. 5 and 7). This response applies to aquatic systems in which a difference exists in temperature and heating regimes between rivers and lakes, but does not appear to affect water bodies with uniform temperature/heating regimes. Our results thus supports the hypothesis put forward by Zhang et al. (2014), stating that climate warming of lakes might be reduced and even reversed by addition of external water.
To investigate this effect, we varied the hydraulic residence time of LB and LG, while holding all other factors constant (Fig. 9). We implemented a stepwise reduction in LG size (to 1/80 of its original volume), simultaneously reducing hypsographic area but keeping maximum depth unchanged. Similar adjustments were made to LB to obtain corresponding hydraulic residence times. This stepwise approach required 972 additional model runs. These iterations showed that river water had to be cooler than lake water in order to generate a dampening effect for climate warming ( Fig. 9a and d). Deep penetrations by large riverine volumes increase the cooling of the hypolimnion (Fig. 9b). The climate dampening effect is suppressed when river and epilimnion temperature are similar. MNPP thermal input creates such conditions in the Aare and therefore largely counteract the river cooling effect of the Aare on LB (Fig. 9c). For shorter residence times (<∼ 1000 days), rivers can exert influence if a significant temperature difference exists between river and lake waters. For longer residence times (>∼ 1000 days), tributaries cannot significantly offset climate effects in downstream water bodies.
Climate-induced warming of lakes (Schwefel et al., 2016), along with changing frequency or intensity of deep penetrating flood events (Fink et al., 2016) may curtail oxygen supply to deep lakes. Recent flood analysis has also indicated that input of river-borne organic matter increases respiration, causing a paradoxical net oxygen reduction within the intruding layer (Bouffard and Perga, 2016). Models showed respective winter increase and summer decrease in river water density relative to lake stability. This creates summer to winter seasonal shifts in deep intrusion dynamics for both lakes (Fig. D1), causing a net annual increase in the river water penetrating to deeper parts of both lakes (Fig. 8). An increase in Rhône SSC in winter represented the primary driver in LG (Figs. 8a, b and D1a, b), while the dominant factor in LB was Aare river temperature, which cooled in winter by increased discharge and removal of MNPP heat (Figs. 8c,d and D1c,d). Fink et al. (2016) also found evidence that climate change will cause diminished deep river intrusion events in summer and enhanced intrusion in winter. They predicted an annual decrease in the amount of river water reaching the deepest parts of Lake Constance. The tributaries considered here differ from the Rhine River investigated by Fink et al. (2016) primarily in terms of temperature. The Rhône catchment for example rests at a mean elevation of 2127 m and includes greater glacial coverage (11 %), whereas the Rhine  Figure 9. Variation in lake hydraulic residence times (changed lake volume) in response to modeled temperature increase ( T ) in the epilimnion (grey) and hypolimnion (black) displayed as decadal mean (solid line) and standard deviation (dotted line) for LG (a, b) and LB (c, d). River-borne SSC included (b) and excluded (a, c and d), MNPP heat input included in (c) and excluded (d) from near-future and far-future time periods but retained for the reference period. Black x's mark modeled lake residence times, while full-height black rectangles mark the lakes' presentday residence times. catchment has a mean elevation of 1771 m and only 1 % glacial coverage (www.hydrodaten.admin.ch). The closure of the MNPP and associated temperature decrease contribute to increase the volume/frequency of deep intrusions (Fig. 5). While Fink et al. (2016) focused primarily on flood frequencies, our models emphasized river discharge regimes and interacting river and lake temperature regimes. The annual increase in river penetration to depth predicted by our models suggests future increase in deep-water oxygen supply in similar tributary-lake systems. This prediction applies mainly to meromictic lakes such as LG. Analogous effects in holomictic lakes such as LB, which mix completely each winter, are less significant. Similar to findings of Fink et al. (2016), our models indicate that deep-water oxygen conditions will worsen during strongly stratified conditions due to seasonal shifts in deep river intrusions from summer to winter. Concluding, as river water density increases in winter, the volume of those intrusion events, which occurred in the reference period, will increase in the future. Likewise, high discharge events, which were previously unable to penetrate into the deep, are likely to do so in the future.

Model reliability
Predictions concerning the effect of climate change on rivers and lakes depend on (i) the choice of emission scenario, (ii) the accuracies of models linking climate to hydrology and climate to heat fluxes and (iii) natural variability of the system being investigated (Raymond Pralong et al., 2015). This section describes uncertainties and limitations of our approach.
Results of long-term forecasts (beyond 2050) depend strongly on representations of global greenhouse gas emission scenarios (FOEN, 2012). Given the uncertainties in future global climate policy, we chose a median scenario, which falls between the best (e.g., RCP3PD) and the worst case scenarios (e.g., A2) in terms of greenhouse gas emissions. A1B assumes a peak population at mid-century, balanced use of renewable and fossil fuels and rapid introduction of new technologies.
Estimates of future air temperatures and river discharge were obtained from a combination of regional climate models (RCMs; CH2011, 2011FOEN, 2012). Uncertainties associated with individual RCMs were offset by combined forecasts from multiple-model chains. Numerous studies have performed detailed evaluations of uncertainty in air temperature and river discharge under established emission scenarios (RCP3PD, A1B, A2) and accounting for globalregional climate model interactions (Addor et al., 2014;Bosshard et al., 2011Bosshard et al., , 2013CH2011, 2011. The degree of accuracy with which model input parameters represent future conditions determines the accuracy of model predictions. We therefore ran the river temperature model with varying parameters to evaluate model sensitivity (Table 1) for different yet similar datasets. The air2stream parameter a 1 showed the greatest degree of sensitivity, varying within 3 orders of magnitude. The a 1 parameter, however, does not respond to variations in river discharge or air temperature (Eq. 1), which limits its sensitivity to climatic input data. The other parameters (a 2 to a 8 ) varied only within 1 order of magnitude (Table 1). The SSC model gives better results for the Rhône (coefficient of determination R 2 = 0.68 from 2013 to 2014) than for the Aare (R 2 = 0.06 from 2013 to 2014). Dams and reservoir infrastructure upstream of station no. 2085 along the Aare dampen sediment transport events and decouple them from regional discharge events (see above; Fig. 4). Given the relatively minor effect of SSC on the Aare water density, variation in the input parameter does not influence river intrusion depths (Figs. B1e-f and D1c-d). As with other vertical, one-dimensional models, SIMSTRAT cannot account for lateral heterogeneities in lakes. This inherent weakness in model design, however, does not significantly diminish the accuracy of model pre-dictions concerning LB and LG (Råman Schwefel et al., 2016).
Of special importance for climate research in lakes is the sensitivity of models to shifts in the heat budget. Forcing parameters of importance, besides air temperature, include wind speed, solar irradiance, vapor pressure and light absorption. The sensitivity of SIMSTRAT to variable forcing has previously been established for lakes in Switzerland. Schmid and Köster (2016) demonstrated how solar brightening from 1981 to 2013 increased Lake Zürich surface warming comparable to heating by increased air temperature. Schwefel et al. (2016) revealed strengthening of the thermocline and decrease in the mean lake temperature by increased light absorption in LG, whereas a decrease in absorption had the reverse effect. As of yet, reliable predictions of wind speed, irradiance and vapor pressure under future climate conditions are not available for Switzerland (CH2011, 2011). Therefore, we use long-term (1981 to 2013) data from station no. 8100 (Fig. 2) as guidance for potential annual atmospheric forcing trends ( Fig. A6; Table 6).
The sensitivity of SIMSTRAT was tested in LG by applying these trends, individually and combined, to the reference period. The increasing trend in air temperature was included for comparison, while no trend could be identified in cloud cover which was excluded. The decreasing trend in wind speed cooled the lake, while the increasing trend in irradiance and vapor pressure heated the lake comparable to air temperature (Table 6). By combining all trends, we obtained similar warming of the LG epilimnion (∼ +0.38 • C decade −1 ) as observed over land (+0.  1983-2000Gillet and Quétin, 2006). The historical effect of increased air temperature caused ∼ 40 to ∼ 70 % of the heating in the epilimnion/metalimnion and ∼ 240 % in the hypolimnion.
Here we include predictions of future temperature and precipitation. The extrapolation of observed atmospheric trends into the future is outside the scope of the present study. Yet, we expect our lake water temperature predictions for the near-future and far-future scenarios to underestimate the total heating in shallow water and overestimate warming of deep-water. Nonetheless, the solar brightening trend observed over Switzerland from 1980 to 2000, caused by a decrease in atmospheric aerosols, will not continue into the future (Sanchez-Lorenzo and Wild, 2012), thereby reducing the uncertainty of our predictions.
In this study we assumed that glacial meltwater feeding both the Aare and Rhône in summer will not disappear within the time frames considered. Loss of glacial sources would modify the discharge regime, especially in summer, which would affect accuracy of temperature, SSC and intrusion depth estimates. However, as stated in Sect. 4.1, FOEN (2012) predicts that the Aare and Rhône catchments will retain 30 % of their glacial masses by the year 2100. Table 6. Observed trends in atmospheric forcing (1981 to 2013) at station no. 8100 (Fig. 2) per decade, and modeled temperature increase in Lake Geneva (LG) with forcing trends applied to the reference period (1990 to 2009 These predictions support assumptions concerning the Aare and Rhône discharge regimes used here. Point sources/sinks of anthropogenic heat can affect inland water bodies response to climate change, as shown by the MNPP effects described here. Other changes in catchment management, such as hydropower damming would also alter river discharge regimes and by extension, temperatures, SSCs and deep-water renewal (Fink et al., 2016). Thus, the correctness of future climate change predictions depends on adequate accounting of regional anthropogenic factors affecting physical processes in the system under investigation.

Conclusion
Aquatic processes in lakes are the result of regional forcing and the upstream catchment environment. This study investigated the impact of climate change on inland waters by propagating climatic inputs through integrated fluvial-lacustrine systems. We fed predicted future climatic data into models for two connected river and lake systems in order to evaluate downstream thermal responses and how river discharge regime shifts might affect deep-water renewal in the lakes. Climate data propagated through discharge-dependent river temperature and SSC models, which were coupled to a onedimensional lake model. We applied this approach for the two peri-alpine lakes Biel and Geneva. The models showed that climate warming of rivers is enhanced in summer and diminished in winter due to future river discharge regimes with decreased flow in summer and increased flow in winter. This climate-caused alteration of the flow regime likewise increase and decreases the river-borne suspended sediment load in winter and summer, respectively.
Both lakes showed large seasonal temperature increases that could not be solely explained by climate-related (predicted) increases in air temperature. Instead, the lakes experienced a cooling effect associated with upstream tributaries, where responses to increasing future air temperatures differed from that of the lakes. The smaller Lake Biel showed stronger response to this repressive effect of climate warming compared to the larger Lake Geneva. Predicted changes in Lake Biel strongly depend on the removal of upstream anthropogenic thermal emission into the Aare River. Local anthropogenic point sources of heat can thus rival climate change in their influence on lakes. This damping of climate warming depends on the lakes hydraulic residence times and requires adequate river/lake temperature differences. Our models indicate that tributaries can exert a system-wide influence on lakes with hydraulic residence times less than ∼ 1000 days. Lake systems with longer residence times are resistant to tributary effects but may respond on a local level.
The combination of changes in river SSC and differential lake/river temperature/warming result in a seasonal shift of deep-water penetration (by rivers) into lakes. The volume of river water penetrating to deeper parts of lakes specifically decreases in summer and increases in winter. Higher rates of deep-water renewal can in turn enhance reoxygenation of the deepest reaches of lakes, which may otherwise experience lower oxygen concentrations under climate change.
Data availability. CTD profiles are available from the Office of Water Protection and Waste Management of the Canton of Bern (GBL/AWA) at: http://www.bve.be.ch/bve/ de/index/wasser/wasser/messdaten.html (GBL/AWA, 2017). Meteorological data are available from the Swiss Federal Office of Meteorology and Climatology (MeteoSwiss) at: http://www.meteoswiss.admin.ch/home/services-and-publications/ beratung-und-service/data-portal-for-teaching-and-research.html (MeteoSwiss, 2017 Figure A1. Illustration of river intrusion model. Appendix A: River intrusion model Figure A1 summarizes the river intrusion model. The depth where the river plume separates from the surface, the socalled plunge depth (h p ), depends on the slope angle (σ ), gravity (g), coefficients (S 1 = 0.25, S 2 = 0.75; Ellison and Turner, 1959), bed friction (f t = 0.02; Akiyama and Heinz, 1984), initial flow per unit width (q 0 ) = V r /W r dependent on river discharge (V r ), river width (W r = 100 m for the Aare and W r = 120 m for the Rhône) and the relative density difference (ρ = (ρ r −ρ l (z 1 ))/ρ l (z 1 )) between the homogenous river (ρ r ) and lake (ρ l ) with z 1 = surface.
The level of initial plume entrainment is treated differently on a gentle versus a steep slope. This is controlled by the two coefficients e 1 and e 2 .
where the critical slope angle σ c = f t S 1 /S 2 distinguishes between gentle and steep slope designations. The initial height of the underflow (h d ) can then be written as where γ is the entrance mixing coefficient equal to ∼ 0 for gentle slopes and increasing to larger values for steeper slopes. Here we find that a value of γ = 0.1 provides best results. The initial underflow temperature (T u ), velocity (U u ), particle content (P u ) and volume (V u ) is consequently expressed as a function of ambient lake water temperature (T l ), river temperature (T r ) and river particle content (P r ) in the homogenous region.
Once the plume has passed through the plunge region into the underflow region, we express h d , U u , T u and V u as where x is the horizontal distance between z and z + 1 and the entrainment factor (E) is expressed as a function of the entrainment constant (β = 0.0015; Ashida and Egashira, 1975) and the Richardson number (R i ).
For P u , we include a sedimentation term as proposed by Syvitski and Lewis (1992), which depends on the removal rate (r) and t = x/U u (z).
Sedimentation occurs only if the plume velocity drops below a critical settling velocity (U c ) subject to the parameter e 3 : We set U c equal to 0.46 m s −1 and r equal to 4.7 day −1 to represent medium-sized silt following Mulder et al. (1998). The plume travels downslope as long as the underflow plume density (ρ u ) exceeds ρ l (z). Once ρ u ≤ ρ l (z), the plume raises from the slope and intrudes into the lake proper. The terms T u and V u were thus added to the lake model at this depth. Calculations excluded expressions for the settling of accumulated particles following plume intrusion, assuming that these exert only minor impacts on lake temperature and density.
The sensitivity of the river intrusion depth to entrainment of ambient water into the plume was tested by propagating a range of β (Eq. A13) values from 1 to 1 × 10 −6 through model spaces composed of temperature, discharge and SSC data from the Aare (station no. 2085). Figures A2-A4 compare modeled intrusion depths to empirical estimates based on vertical temperature and light transmission data at the centre of LB (47 • 6 N, 7 • 11 E) collected shortly after major river intrusion events. Additionally, acoustic Doppler current profiler (ADCP) measurements of river plume intrusions in LB (47 • 5 N, 7 • 12 E; 2 km from the Aare inlet) were used for a temporal sensitivity analysis of the intrusion model (Fig. A5). Comparison of the modeled intrusion depth with light transmission depth (whose minimum value represents a proxy for actual river intrusion depth) suggests that β = 0.0015 offers an adequate representation of intru-  Figure A6. Annual atmospheric forcing (grey dots) of wind speed (a), shortwave radiation (b), vapor pressure (c) and Air temperature (d) at station no. 8100 (Fig. 2). The black line shows trends (Table 6).
Appendix B Figure B1. Modeled climate impact on LB excluding river-borne SSC. Temperature increase T (a and b) displayed as means (bars) and standard deviations (black lines) in epilimnion (left bar group), metalimnion (middle bar group) and hypolimnion (right bar group); mean intruding river volume (c and d) and mean river intrusion depth (e and f). MNPP thermal input included (a, c, e) or excluded (b, d, f) in near-future (blue) and far-future (vermilion) time periods but retained in the reference period (black).  Figure C1. Modeled climate impact (river intrusion excluded) on LB (left column, scenario LB1) and LG (right column, scenario LG1) shown as daily mean (thick lines) and maximum/minimum model values (thin lines) for near-future (blue, 2030-2049) and farfuture (orange, 2080-2099) time periods relative to the reference period (black, 1990-2009). Temperature T (a and b), temperature increase ( T ) in the epilimnion (c and d) and hypolimnion (e and f) as well as increase in stability ( S; g and h).

MNPP Excluded (d) LB3
Lake Geneva Lake Biel Figure D1. Modeled climate impact on mean river intrusion depth. Reference period (black), near-future (blue) and far-future (orange) time periods for LG (a-b) and LB (c-d) with (b-d) and without (a) river-borne SSC and MNPP thermal input included (c) or excluded (d) from near-future and far-future time periods relative to the reference period.