Background

Reservoirs are important sources of drinking water in most parts of the world. Like other water bodies reservoirs are impacted by climate change. This is reflected, for example, in physical effects such as increasing surface water temperature [1, 2], decreasing ice cover duration [3,4,5], changing stratification [6] or in biological effects such as changes in the phytoplankton community [7] and the increasing risk of cyanobacteria blooms [8, 9]. Climate change also amplifies processes leading to eutrophication of water bodies [10] which might reinforce global warming [11]. Most of these changes are associated with decreasing water quality [12, 13] and threaten reliable drinking water production.

Reservoirs respond differently to climate change compared to lakes because storage and outflow are actively managed [14]. The operational parameters associated with reservoir control are: the withdrawal rate or quantity, the withdrawal schedule, and the withdrawal depth [15]. The withdrawal depth directly influences storage or dissipation of heat and material, thermal stability, and thus resistance to mixing [15]. In drinking water reservoirs, adaptation of withdrawal depth is used as a tool to optimize raw water quality for drinking water production [16]. Reservoir operation strategies that include an adaptation of withdrawal depths were recently examined by a number of studies focusing on the thermal properties of the reservoir [17, 18] or the temperature regime of the downstream river [19, 20].

The impact of the withdrawal depth from a single outlet has been assessed by comparing the thermal structure of a lake with surface discharge and a similar reservoir with ground release [21], by comparing periods of hypolimnetic withdrawal with periods of withdrawal at intermediate depths in the same reservoir [22], or by computer simulations [17, 23,24,25]. In contrast to surface release, withdrawal of water from deeper layers causes an increase in internal heat energy and results in a shorter period of stratification. So far, modeling studies investigated the effects of the withdrawal depth on reservoir characteristics, like e.g., the duration of inverse stratification [17], the thermal regime [24] or the depth of the thermocline [23]. The studies all concluded that the investigated variables are affected by the withdrawal depth.

Until now, most studies [17, 24,25,26] considered withdrawal from a single constant depth layer throughout the year. However, practical and realistic reservoir management needs to consider multiple factors simultaneously and thus the depth of raw water extraction varies over time. Moreover, it is common to release water to the downstream river from a different depth than the water which will be used for drinking water production. Furthermore, only a few of the existing studies took into account the impact of interannual variations in the external forcing, by e.g., simulating over several years.

With this study we investigated whether dynamic withdrawal strategies are an effective means to mitigate the impact of climate change on drinking water reservoirs. We therefore focused on three questions:

  1. 1.

    What is the impact of the different withdrawal schemes on the thermal structure of reservoirs?

  2. 2.

    How does the impact depend on reservoir characteristics and external forces like, e.g., wind speed or withdrawal rate?

  3. 3.

    What is the possible contribution of dynamic withdrawal to operation strategies targeted at minimizing the negative impact of climate change?

Specifically, we studied the effects of four realistic withdrawal schemes applied to three German reservoirs by means of a one-dimensional hydrophysical model (GLM [27]). In contrast to most existing studies, the effect of the alternative schemes was demonstrated for a wide range of hydro-meteorological forcing (16 years).

Methods

Study area

The investigated reservoirs, Eibenstock (ES), Lichtenberg (LB) and Saidenbach (SB) are located in the low mountain range Westerzgebirge, Germany (Fig. 1) and are managed by the State Reservoir Administration of Saxony (LTV). The main purpose of all three reservoirs is to provide raw water for drinking water production, and to protect downstream areas from major floods. Due to socioeconomic factors a noticeable decrease in the average quantity of raw water withdrawal can be observed in ES and SB during the period from 1990 to 2017, whereas the raw water withdrawal in LB remained constant. Specifics on the location and morphological features of the reservoirs are provided in Table 1.

Fig. 1
figure 1

Location, shape of surface, and hypsographic curves (total volume plotted against local water level) of the investigated reservoirs: 1. Eibenstock Reservoir, 2. Saidenbach Reservoir, and 3. Lichtenberg Reservoir. Map data were taken from the R package mapdata and GADM (https://gadm.org). Reservoir contours based on data from the Saxon State Office for Environment, Agriculture and Geology (LFULG)

Full size image
Table 1 Main morphological features, withdrawal depths, and used abbreviations (Abv.) of the three investigated reservoirs
Full size table

All three reservoirs are equipped with withdrawal structures that allow water to be discharged from different depths. The reservoirs ES and LB have withdrawal structures with extraction outlets at fixed depths. ES has 6 different outlets, the lowest 1.4 and the highest 43.3 m above the reservoir bottom. LB has 5 different outlets, the lowest 2 and the highest 25 m above the reservoir bottom. SB has one extraction outlet 2 m above the reservoir bottom and one flexible raw water intake that can be used to extract water in a range from 17 to 30 m above the reservoir bottom. This intake structure was constructed in 1989 and started operation in 1990. LB and SB release water to their downstream river through the bottom outlet; whereas in ES it is possible and common practice to release water to the downstream river from one or more of 6 outlets, located at the same heights as the raw water extraction horizons.

Model input

For the model simulations we used daily meteorological data, which have been measured on-site by the Reservoir Administration (LTV). Missing data were estimated by inverse distance weighting (power one) from observations of nearby weather stations operated by the German weather service (DWD) and obtained from their climate data center [28].

Daily hydrological data were provided by LTV and, wherever necessary, the inflow was adapted so that computed reservoir volume matched the observed one. Inflow temperature was estimated using the empirical method of Horn et al. [29], where stream temperature \(T_w\) for each day of the year d is approximated by a harmonic function with an autoregressive component,

$$\begin{aligned} T_w (d)&= b_1 + b_2 \cdot f_{\text {sin}}(d) + b_3 \cdot \left( f_{\text {sin}}(d) \right) ^2 + f_{ar}(d),\nonumber \\ f_{\text {sin}}(d)&= \sin \left( \frac{2 \pi (d + b_6)}{365.25}\right) , \nonumber \\ f_{ar}(d)&= \frac{1}{m}\sum _{j=0}^{m}T_a(d-j) \cdot \left( b_4 + b_5 \cdot f_{\text {sin}}(d) \right) , \end{aligned}$$
(1)

where \(b_1\) to \(b_6\) are calibration parameters, \(T_a(d)\) is the air temperature at a given day of the year d, and \(f_{ar}\) is an autoregressive term considering the air temperature of the m previous days. We calibrated the parameters separately for ES, LB, and SB using a standard Nelder–Mead optimizer [30]. The best-fit models had root mean squared errors of 1.57 °C, 1.60 °C, and 1.13 °C and used the air temperature of the last 4, 8, and 2 days (m) for ES, LB, and SB, respectively. After parameter fitting, we approximated inflow temperatures for every day.

Temperature profiles in the reservoirs are measured by the LTV in regular intervals. For ES 25 years, for LB 24 years and for SB 39 years of data were available. From this data monthly average temperatures for 3 m and 25 m below the water surface were calculated. Years with fewer than 11 monthly values were excluded and for the remaining years annual average temperatures were calculated. From these monthly and annual averages linear temperature trends were estimated and tested for significance (Mann–Kendall test [31], p < 0.05).

From the temperature profiles the Schmidt stability index was calculated using the R package rLakeAnalyzer [32] to estimate onset and end of summer stratification for every year. Based on the approach of Engelhardt and Kirillin [33] summer stratification was defined as days where Schmidt stability exceeded a threshold of 50 J/m\(^2\). In order to exclude inverse stratification we also added a temperature threshold of 5 °C for total average temperature. Due to a measurement interval of 2–4 weeks, we estimated the date of onset as the midpoint between the last observation with stability below the threshold and the first observation above. Similarly, we defined the end of stratification as the midpoint between the time of last observation above the threshold and the next observed point.

Daily observations of ice coverage were available for ES, LB, SB starting from 1997, 1976, and 1975, respectively. Information on withdrawal quantity and elevation was available on a daily basis for LB, SB, and ES starting from 1993, 1985, and 2000.

Model setup

We used the open source, General Lake Model (GLM) version 2.4.1 [27]. This one-dimensional model is widely used and has been successfully applied in previous studies [e.g., 34,35,36]. In the Multi-lake Comparison Project (MLCP [37]), GLM was applied to a variety of different lakes and reservoirs using a common calibration and assessment procedure. In the majority of test cases minimal parameter adjustment was sufficient to reliably reproduce observed temperature profiles and stratification patterns [37].

For modeling purpose we used the time period from 2000-01-01 to 2016-12-31 for all three reservoirs, as for this period all necessary data were available. We used cloud_mode = 3 [38] and albedo_mode = 3 [39] because compared with the other modes they gave the lowest root mean square error. Bed warming and wind fetch modules included in GLM were disabled because of insufficient data availability. Like in the MLCP five model performance measures were calculated: root mean square error (RMSE), Nash–Sutcliffe efficiency (NSE), Pearson correlation coefficient (r), percentage relative error (PRE), and normalized mean absolute error (NMAE) [cf. 40].

For calibration the first 8 years of data (2000-01-01 until 2007-12-31) were used for all three reservoirs. We used the same values for mixing parameter (e.g., coef_mix_conv, coef_wind_stir, coef_mix_shear) as the MLCP and, similar to the approach described therein, calibrated four parameters (Table 2). First the two parameters wind_factor and Kw, and then strmbd_slope and strm_hf_angle were calibrated manually. For both pairs of parameters, first adequate initial guesses for the parameters were chosen from a priori knowledge. Then, parameter sets containing 200 combinations of the parameters in the range of ±30% of the initial guesses were created using Latin hypercube sampling and the model performance for all combinations was checked. The quality parameters mentioned above were calculated for all model runs and the set with lowest RMSE was chosen to start a second iteration with another 200 random parameter sets of ±10% the best parameter values from the previous run. From this second run the parameter set with the lowest RMSE was chosen as the final set. For model validation the last 9 years of data (2008-01-01 until 2016-12-31) were used and the five goodness-of-fit indicators mentioned above were calculated.

Table 2 Description and units of the four parameters used for model calibration
Full size table

Management strategies

For each reservoir, we simulated and compared four strategies of management. The latter solely differ with regard to the depths of the outlet structures from which drinking water is withdrawn or water is released to the downstream river, respectively. These include the strategies as they are currently realized, in each of the reservoirs. The four strategies are:

  • Strategy 1 low: raw drinking water is withdrawn from the lowest possible outlet structure at all times; downstream river is fed from bottom outlet throughout the year.

  • Strategy 2 high: raw drinking water is withdrawn from the lowest possible outlet structure at all times; downstream river is fed from the highest possible outlet structure throughout the year.

  • Strategy 3 dyn1: raw drinking water is withdrawn from dynamically adapted depths; downstream river is fed from bottom outlet throughout the year. This strategy is currently effective for reservoirs LB and SB.

  • Strategy 4 dyn2: both raw drinking water and downstream water are withdrawn from dynamically adapted depths. Reservoir ES is currently managed according to this strategy.

Strategy low is the classical approach that is used in many reservoirs. Strategy high preserves cold water for drinking water production by withdrawing the downstream water from warmer layers. Strategy dyn1 is the approach realized in most reservoirs in Saxony, including LB and SB, where raw water extraction depth is varied, e.g., in order to evade plumes of algae or high turbidity after storm events. In all three realizations of dyn1 we used the withdrawal depth that was used in practice. Strategy dyn2 is additionally varying the depth of water released to the downstream river. A growing number of reservoirs including ES apply this strategy because it saves cold hypolimion water for drinking water production and decreases thermal pollution in the downstream river. For LB and SB, we used the adaptive withdrawal module included in GLM [see [20], outlet_type=5]. The module was configured such that water is withdrawn from a depth where the temperature matches the inflow temperature. Along with the target temperature, at a daily resolution, the model is provided with upper and lower boundaries for withdrawal depths. If the target temperature is outside these boundaries the closest temperature within is chosen. Raw water withdrawal was set to the depth that was used in practice (same as in dyn1).

The model was run once for each strategy and reservoir. The first year of model output was generally ignored so as to eliminate the effect of the initial state. For every model run we calculated a set of characteristic quantities, called features, for comparison of the different scenarios (Table 3). The R package rLakeAnalyzer [32] was used to calculate the internal heat energy and Schmidt stability. The onset and end of summer stratification were calculated in the same way as for the observed data (See “Model input”).

Table 3 Abbreviation and explanation of quantitative hydrophysical features used for statistical analysis (doy: day of year)
Full size table

Statistical evaluation

In order to test if the effect of the withdrawal strategy varies with meteorological (e.g., wind speed), and hydrological forcing (e.g., withdrawal rate), a multiple linear regression approach was used. We calculated candidate predictor variables (Table 4) that we used as predictors. A principal component analysis (PCA) was used to determine a set of predictor variables that explained most of the variation in the observed features (Table 3). Forward and backward stepwise predictor selection was then used to identify significant terms. The Bayesian Information Criterion [41] was used to identify the optimal models in terms of parsimony and explanatory power. The stepwise predictor selection was started with a model setup that contained the predictor variables selected from the PCA, two additional (categorical) variables to account for the reservoir and the management strategy, and interaction terms between all predictors and the two categorical variables (Eq. 3). For a single reservoir a linear model for characteristic feature y that considers interactions with the applied strategy can be written as:

$$\begin{aligned} y = \beta _{0} + \sum _{i=1}^{n}\beta _{i} x_{i} + \varvec{\beta }_{\varvec{s}} \varvec{s} + \sum _{i=1}^{n} \varvec{\beta }_{\varvec{i:s}} \varvec{s} x_{i} . \end{aligned}$$
(2)

Here, \(\beta _0\) is the intercept, \(\beta _i\) are the linear coefficients for the n predictor variables \(x_i\), \(\varvec{\beta _s}\) is a vector of linear coefficients corresponding to the management strategy, \(\varvec{s}\) is a unit vector accounting for the management strategy, and \(\varvec{\beta _{i:s}}\) is a vector of linear coefficients accounting for interaction between predictor variable \(x_i\) and management strategy. The full model, additionally accounting for the different reservoirs, is then written as:

$$\begin{aligned} y &= \beta _0 + \sum _{i=1}^{n}\beta _i x_i + \varvec{\beta _{s}} \varvec{s} + \sum _{i=1}^{n} \varvec{\beta _{i:s}} \varvec{s} x_i + \varvec{\beta _{r}} \varvec{r}\\ &\quad + \sum _{i=1}^{n} \varvec{\beta _{i:r}} \varvec{r} x_i + \varvec{r}^\intercal \varvec{B_{r:s}} \varvec{s}, \end{aligned}$$
(3)

where \(\varvec{\beta _r}\) is a vector of linear coefficients corresponding to the reservoirs, \(\varvec{r}\) is a unit vector accounting for the reservoir, \(\varvec{\beta _{i:r}}\) is a vector of linear coefficients accounting for interaction between predictor variable \(x_i\) and reservoir, \(\varvec{r}^\intercal\) is the transposed unit vector of reservoirs, and \(\varvec{B_{r:s}}\) is a matrix of linear coefficients for the interaction between reservoir and management strategy.

The remaining interaction terms in the minimal adequate linear models can help to understand which of the predictor variables are important for the effect of the applied management strategies on the observed features and how the effects of the predictor variables are influenced by the reservoir. Using the fitted linear coefficients of the minimal adequate models these effects can be quantified.

Table 4 Candidate predictor variables with units and abbreviation used in the statistical models
Full size table

Results

Observed trends

Fig. 2
figure 2

Annual mean temperatures (1), annual temperature trends for every month (2), and onset and end of summer stagnation for every year (3) for the three investigated reservoirs. The temperature trends are both for 3 m and 25 m below the surface. In (2) black stars denote significant trends, tested with the Mann–Kendall test. In (3) the whiskers indicate the maximum uncertainty for onset and end of stratification based on the observation dates

Full size image

The observed annual and monthly temperatures, summer stratification, and ice-off showed coherent trends. The annual temperature trends in 3 m depth were 0.03 K/a (ES) and 0.05 K/a (LB and SB). The corresponding annual trends in 25 m depth were 0.03 K/a (LB) and -0.02 K/a (ES and SB) (Fig. 2.1). The increase in 3 m depth was always positive, but the effect varied between months (Fig. 2.2). In all three reservoirs the largest trends were observed between April and July, the period when summer stratification usually begins. A negative temperature trend in 25 m depth can be seen in ES and SB from July to November, peaking around the time when summer stratification usually ends (see “Discussion”). All three reservoirs showed a slight positive trend from December to April (Fig. 2.2). In all reservoirs, the onset of summer stagnation shifted to the beginning of the year by about 0.4 d/a (Fig. 2.3). Likewise, the end of summer stagnation shifted to the end of the year, but the change was stronger in ES and SB. The ice-off in LB and SB occurred earlier by about 0.5 d/a. In ES no trend in the day of ice-off was observed (non-significant slope of linear model, p < 0.05).

Validation of simulated temperature profiles

Fig. 3
figure 3

Water temperatures for the whole simulation period showing observed (top) and simulated (bottom) temperatures for Lichtenberg Reservoir (LB). Modeled results are based on the currently used strategy of water withdrawal (see Section Management strategies)

Full size image

GLM simulates the observed temperatures with sufficient accuracy (for LB see Fig. 3, for the other reservoirs see Additional file 1: Figures S1 and S2). Model performance indicators for both, calibration and validation phase (Table 5), are within the range reported by the MLCP [37]. The negative PRE shows a slight underestimation of simulated temperatures in all three reservoirs. The superior performance of ES and LB compared to SB can be explained by the quality of the used meteorological data. For ES and LB we used direct measurements, whereas for SB only interpolated data were available. Scatter plots of simulated against observed temperatures are shown in Fig. 4 and the best-fit values of the calibrated parameters are provided in Table 6. The calibration of strmbd_slope and strm_hf_angle only slightly improved the model performance.

Fig. 4
figure 4

Comparison between observed and simulated temperatures over all depths for the three reservoirs. Each point represents an observation–simulation pair with unique depth and date

Full size image
Table 5 Goodness-of-fit measures for the three reservoirs Eibenstock (ES), Lichtenberg (LB), and Saidenbach (SB) for the periods of model calibration (calib.) and validation (valid.), 8 years each
Full size table

Except for one year in SB the model correctly predicted the occurrence or absence of ice in all three reservoirs. The first day of ice cover was predicted with a mean absolute error of 4.8, 5.7, and 5.4 days and the last day of ice cover was predicted with a mean absolute error of 11.6, 10.9, and 8.8 days for reservoirs ES, LB, and SB, respectively.

Table 6 Best-fit values of the calibrated model parameters for light extinction coefficient (Kw), wind factor (wind_factor), stream bed slope (strmbd_slope), and stream half-angle (strm_hf_angle) for all three reservoirs
Full size table

The calibrated value of wind_factor of ES is close to 1, whereas for LB and SB it is around 0.5 (Table 6). This reflects the fact that for ES observed wind data, measured directly at the reservoir, were used whereas for LB and SB the data was interpolated from nearby DWD stations. The closest used DWD stations were located at a distance of about 10 km (SB), and 20 km (LB), but have altitude differences of approximately 130 m (SB), and 200 m (LB). The reservoirs are located in valleys of the low mountain range Erzgebirge, where the local wind fields can differ substantially. This can explain the small scaling factors for LB and SB. For all three reservoirs, the calibrated values of Kw seemed plausible, when compared with values estimated from observed Secci disk depths SD, using a relationship of \(Kw = 1.7/SD.\) [43].

Modeled impact of management

Fig. 5
figure 5

Annual average values of temperature in 25 m depth (left) and end of summer stratification (right) for all three reservoirs and all four withdrawal strategies. For description of the reservoirs see Table 1, for description of the strategies see Section Management strategies

Full size image

A change of the withdrawal strategy affected all six hydrophysical features, but the effect size varied between features, reservoirs, and years. In some cases, the difference between two years was larger than the difference between two strategies for a single year. Also, the ranking of the strategies in terms of the effect size differed between years and reservoirs. Figure 5 illustrates this for hypolimnion temperatures and the end of summer stagnation. Generally, the internal heat energy per unit surface area, the temperature in 25 m depth, and the end of summer stagnation showed the strongest response to an alteration of the withdrawal strategy. By contrast, the temperature in 3 m depth, the onset of summer stagnation, and the ice-off were less sensitive (Fig. 6).

Fig. 6
figure 6

Internal heat energy per unit surface area (1), temperature in 3 m depth (2), temperature in 25 m depth (3), onset of summer stagnation (4), end of summer stagnation (5), and day of ice-off (6) for all four scenarios and all three reservoirs. The individual boxplots represent annual averages of 16 years. The scenarios annotated with an asterisk indicate the reference scenario for the respective reservoirs. For description of the reservoirs see Table 1, for description of the strategies see Section Management strategies

Full size image

On average, the withdrawal of downstream water from the uppermost horizon (strategy high) resulted in the lowest internal heat energy in all three reservoirs (Fig. 6.1). For LB and SB the withdrawal of downstream water from the bottom outlet (strategy low) maximized the internal heat energy, while for ES the adaptive drinking water withdrawal (dyn1) resulted in the highest internal heat energy. We observed similar responses in the temperature in 25 m depth as in the internal heat energy per unit surface area. The maximum average difference between the strategies was around 1 K in all three reservoirs (Fig. 6.3), but there were years with maximum temperature differences between strategies of up tp 2 K in ES and up to 1.5 K in LB and SB (Fig. 5). On a monthly scale, we found the largest monthly differences of up to 6 K during the summer stratification (see Additional file 1: Figures S3 to S5). Regarding the end of summer stagnation the largest difference between the scenarios in SB was 16 days, in LB 9 days, and in ES 3 days (Fig. 6.5).

For the temperature in 3 m depth the largest annual average effect of changing the withdrawal strategies in all three reservoirs was about 0.2 K. At a monthly scale we found larger differences of up to 2 K, that were caused by time shifts at the beginning or end of the stratification period (see Additional file 1: Figures S3 to S5). For the onset of summer stratification we found the largest effects between the strategies in LB (3 days) and the smallest in ES and SB (each 1 day). The day of ice-off proved to be almost insensitive to a change of the withdrawal strategy and had a maximum difference of 1 day in all three reservoirs (Fig. 6.6).

External and internal forcing

From the candidate predictor variables (Table 4), we selected the following six as predictor variables by means of the PCA: average summer air temperature, average winter air temperature, average wind speed, average reservoir volume, total downstream water withdrawal, and total raw water withdrawal. The minimal adequate linear models selected by minimizing BIC had coefficients of determination (R\(^2\)) of 0.574 for the begin of summer stratification to 0.92 for the internal heat energy per unit surface area (Table 7). Compared to the full model the BIC improved by more than 100 units while the R\(^2\) only slightly decreased (maximum decrease 0.067) for all six characteristic features. The stepwise model selection excluded the variable representing the applied withdrawal strategy and all interactions with it for the characteristic features “onset of stratification” and “day of ice-off”, these were also the variables where we saw the smallest scenario effect (Fig. 6) and where the statistical models had the lowest explanatory power.

Table 7 Parameters of linear regression (minimal adequate model according to Bayesian information criterion [41]) between modeled reservoir features and external predictor variables
Full size table

The interactions with the management strategy that remained in the models were only with hydrological predictor variables, namely the withdrawal rate of downstream water, the withdrawal rate of raw water, and the reservoir volume. The meteorological predictor variables that remained in the minimal adequate models showed only interactions with the reservoir variable.

Interpreting the parameters of the linear models can give insight on the average effect of changing withdrawal strategy on the observed features. For example, changing the withdrawal strategy in ES from dyn2 (the currently applied strategy) to high would on average increase the temperature in 3 m depth by \(0.0142 - (-0.154) = 0.1682\) K (Table 7). In all three reservoirs the temperature in 3 m depth decreased by 0.936 K per 1 m/s average annual wind speed. In terms with interactions, the default (i.e., intercept) parameters represent reservoir ES or strategy dyn1. So, the average end of ice cover in reservoir ES is 20.8 days earlier in the year per Kelvin of average winter temperature, while in reservoirs LB and SB it is \(-20.8 + 9.79 = -11.01\) (LB) and \(-20.8 + 9.41 = -11.39\) (SB) days earlier per Kelvin.

Discussion

Observed trends in hydrophysical features

The observed summer surface temperature trends in all three reservoirs were larger than the global average (0.034 K/a) reported by O’Reilly et al. [2]. This is in concordance with other studies [2, 44] according to which lakes in mountainous areas and with partial ice cover warm faster. For SB we calculated similar values as reported in previous studies by Jäschke et al. [45] (in April; Jäschke et al.: 0.09 K/a here: 0.07 K/a).

The negative temperature trend in 25 m, observed in ES and SB, can be attributed to an extension of the stratification period. Compared to LB, the end of summer stratification shifted more towards the end of the year in ES and SB. This time shift could partly be caused by the decrease in raw water withdrawal that can be seen in ES and SB but not in LB, where raw water extraction was rather constant throughout the considered period of time. Previous studies found stronger impacts of climate change on the stratification in deeper lakes [46]. So, an additional reason for the stronger shift of stratification in ES and SB, could be their larger depth and volume.

The observed shifts towards the beginning of the year in the day of ice-off in ES and SB are within the lower range of values reported for other German reservoirs by Wilmitzer et al. [5] and similar to values reported for lakes in Poland by Skowron [4]. The fact that we did not observe an effect for ES could be due to the shorter time series analyzed (21 a) compared to the other reservoirs (43 a). Additionally, ES is located at the highest altitude of the three reservoirs (Table 1), that is in agreement with Wilmitzer et al. [5] who reported smaller shifts in day of ice-off at higher altitudes.

Modeled impact of management

By changing the withdrawal strategy the internal heat energy of the reservoirs can actively be influenced. We saw a larger impact of changing the withdrawal strategy for temperatures in deep layers (25m depth) than close to the surface (3 m depth). The effect on the end of summer stratification was smaller in ES than in LB and SB, we attribute this to the larger volume of ES.

If we compare the difference between the strategies low and high with other studies that investigated the impact of changing withdrawal depth [17, 23, 25] the results are in concordance. Withdrawing water from deeper outlets increases the internal heat energy, decreases the duration of summer stratification, and increases hypolimnic temperatures. But the predicted management impact in our study was smaller than in previous studies. Similar to the results reported by Mi et al. [17], we found that withdrawing water from higher outlets (strategy high) increased inverse stratification and caused a prolongation of the ice cover duration. But, in both cases we found that the latter effect was small (maximum effect of 1 day). We attribute the smaller effect sizes to the fact that, in our case, the raw water was withdrawn from a second outlet and only the withdrawal depth of the downstream water had been changed.

For the management strategies where the withdrawal depth varied over time (scenarios dyn1 and dyn2), we found similar patterns as reported by Weber et al. [20]: in all reservoirs the dyn2 strategy caused a drop in the hypolimnion temperature and thus maximized the storage of cold deep water. In ES, similar to the results found by Weber et al. [20] the thermocline moved slightly upwards and the Schmidt stability increased (compared to low; see Additional file 1: Figure S6). This is in contrast to LB and SB where the thermocline moved downwards and thermal stability decreased. In all three reservoirs the thermocline decreased around March and increased at the end of autumn, but the effect in March was stronger in LB and SB. We did not fully understand this behavior, but suspect that the difference could be caused by the bathymetry of the reservoirs and the depth of the outlets.

The reservoirs showed different response to the withdrawal strategies for the internal heat energy. In ES strategy dyn1 resulted in the highest internal heat energy, while in LB and SB low gave the highest internal heat energy. We suspect that this difference was caused by the withdrawal during winter, which in dyn1 was from higher outlets, resulting in the withdrawal of colder water, while during summer the withdrawal depth was about the same (from low elevations). In ES the reservoir manager chose to withdraw raw water from higher elevations during winter more often than the reservoir managers in LB and SB, this can explain why the effect is only observed in ES.

Another difference between the three reservoirs can be seen in the end of summer stagnation. Changing from strategy low to high increased the end of summer stagnation for LB and SB but decreased it for ES, while changing from strategy dyn1 to high decreased it in ES and SB but increased it in LB. These differences are related to the size of the hypolimnion volume and thus the reasons for the end of summer stagnation in the reservoirs. In ES the autumn turnover is mainly caused by convective cooling while in LB and SB the turnover is caused by shrinking or even depletion of the hypolimnion due to withdrawal. In ES strategy high decreased internal heat energy allowing for faster cooling, while in LB and SB it took longer to (completely) deplete the hypolimnion (see Additional file 1: Figure S7). In SB, with strategy dyn1 the effect of increased internal heat energy seemed to outweigh the one of decreased hypolimnion volume.

Statistical model

From the minimal adequate linear models we saw that meteorological forcing was in general important for the inter annual differences in the characteristic features, but only the hydrological forcing was important for the effectiveness of the applied withdrawal strategy. Also, the response of the three reservoirs to meteorological forcing was different, e.g., the time shift of end of summer stratification along with increased wind speed was larger in LB and SB than in ES. Or the warming effect of summer air temperature on temperatures in 3 m depth was smaller in LB and SB than in ES (Table 7). Previous studies [46, 47] found similar differences in the response of reservoirs and lakes to meteorological forcing and attributed them to the different volumes and surface areas of the reservoirs.

The modeled effects of the applied withdrawal strategy on the beginning of stratification and the end of ice cover were so small, that the minimal adequate linear models omitted these terms. These were also the two characteristic features where the linear models explanatory power (R\(^2\)) was smallest. The begin of summer stratification is partly controlled by the timing and temperature of the inflow [48]. A possible reason for the relatively weak model performance could be that we were not able to find data for such an explanatory variable for our linear model.

With increasing summer air temperature the internal heat energy decreased. A possible reason for this could be increased stratification strength in warmer years. We also saw larger values of internal heat energy with increasing rate of downstream and raw water withdrawal, whereas for the downstream withdrawal the slope of the relationship depended on the applied strategy. In strategies where (during summer) colder water was withdrawn from the bottom outlet (dyn1 and low) the warming effect was larger, whereas in cases where (during summer) warmer water was withdrawn (dyn2 and high) the warming effect was lower. Also, the effect was largest in reservoir LB, which has the smallest volume. Similarly, there was a positive relationship between the withdrawal rate of downstream and raw water and the temperature in 25 m depth, where the warming rate for the downstream withdrawal was also depending on the applied strategy. In all three reservoirs the effect of raw water withdrawal was larger compared to the effect of downstream water withdrawal.

In all three reservoirs we saw that at larger downstream withdrawal rates the end of summer stratification shifted towards the beginning of the year. But, the effect of increasing withdrawal rate depended on the reservoir and strategy. In SB increasing raw water withdrawal always shifted the end of stratification towards the beginning of the year, whereas in LB larger raw water withdrawal rates were always associated with a delayed end of stratification. In ES we observed shorter summer stratification for high and low and longer summer stratification for dyn1 and dyn2. The different response could be due to the different morphological properties of the reservoirs and highlight the non-linear reaction of stratification to changes in the thermal structure of lakes and reservoirs [46].

Using linear models, we were able to better understand the impacts of withdrawal depth and rate on the thermal structure of the three investigated reservoirs and to a certain extend to quantify the impact of the withdrawal strategy. This worked well for the integrated variable internal heat energy (R\(^2\) = 0.92) and still reasonably well for the variables end of stratification (R\(^2\) = 0.84) and temperature in 25 m depth (R\(^2\) = 0.63), which have finer spatial or temporal resolution.

Climate change and implications for management

We evaluated the suitability of mitigating the impact of climate change in drinking water reservoirs by analyzing the differences between the withdrawal strategies (Fig. 6) as well as the parameters of the linear models (Table 7) and comparing them to the measured impact of climate change (Fig. 2). According to our results, the rise in surface temperature (3 m depth) which is attributed to climate change can hardly be compensated by an adaption of the withdrawal strategy. Comparison of the median values (Fig. 6.2) suggests an effect of changing withdrawal strategy of about 0.2 K, the linear model suggest an maximum difference of about 0.17 K, that is equivalent to the average warming within 4 years (Fig. 2.1). Similarly, our results suggest that the trends in the beginning of summer stratification and the end of ice cover can not be compensated by adapting the withdrawal strategy.

For the temperature in 25 m depth and the end of stratification we found that, in LB and SB, the maximum difference between the strategies is in the same order of magnitude than the average shift observed in the last 30 years. So for these two features mitigation is, in theory, possible, but reservoir managers also need to consider water quality parameters when selecting the withdrawal depth. In SB switching to strategy dyn2 could decrease water temperature in 25 m depth by 0.5 K, but it would shift the end of summer stratification about 3 days towards the end of the year. In LB switching to strategy dyn2 could slightly decrease water temperature in 3 m depth by 0.1 K and shift the onset of stratification towards the end of the year by 2 days.

Our findings suggest, that other operational parameters like the withdrawal rate or the volume in the reservoir also affected the spatio-temporal patterns of water temperature. The linear models (Table 7), indicate that increasing the volume can be an effective measure in order to reduce water temperature in 25 m depth in all reservoirs. Also, the effectiveness of the applied withdrawal strategy depends on the withdrawal rates of raw water and downstream water.

Despite the pronounced effects of global warming on the hydrophysical structure and hence water quality of lakes and reservoirs [12], the response of water temperatures and stratification patterns to management appeared to be moderate. We expect that the predicted effects of changing withdrawal strategy on physical features also translate into moderate changes of water quality parameters [49, 50]. In particular, the altered stability and duration of stratification might trigger changes in phytoplankton composition and dynamics [51]. Especially decreased stability, stratification duration, and lower temperatures can weaken the otherwise competitive advantage of cyanobacteria in dimictic systems [52, 53]. Likewise, hypolimnetic oxygen might be stabilized at higher levels as the stability and duration of stratification is reduced [49, 54]. Nevertheless, it appears unlikely that climate-driven trends in reservoir water quality can fully be stopped or reversed just by an adjustment of withdrawal depth. Water quality models can help to quantify the extend to which adjusting the withdrawal strategy can mitigate water quality issues and should be addressed in future studies.

Conclusion

By means of hydrophysical modeling, we were able to show that an adaption of the withdrawal strategy influences the thermal regime and stratification of a reservoir. We could demonstrate that the impact of the withdrawal strategy depends on the reservoir properties. Specifically, the impact on the end of summer stratification was larger in the reservoirs with smaller volume. Using statistical models we identified the most important predictor variables impacting thermal structure and stratification in our reservoirs, which are: air temperature, wind speed, withdrawal rates, and volume in the reservoir. Our simulations did not indicate that the effect of withdrawal strategy on the characteristic features is modified by meteorological forcing. But, the effect of changing withdrawal strategy on the temperature in 25 m depth depended on withdrawal rates of water released to the downstream river. The effect on the end of summer stagnation depended on raw water withdrawal rate, and the volume in the reservoir.

There is no optimum strategy for all reservoirs that simultaneously decreases water temperature, length of summer stratification, and increases ice cover duration. Reservoir managers hence, need to consider their management objectives, prioritize them, and adapt the strategy accordingly. For example, in LB a switch to strategy dyn2 could decrease temperature in 25 m depth, but it would increase the end of summer stagnation in return (Fig. 6). We did not consider biological or chemical implications, but literature suggests that changes in stratification and temperature can affect these two. Increased stratification can lead to anoxic conditions in the hypolimnion [50], and changes in temperature and stratification can lead to shifts in the phytoplankton community [51, 52]. Reservoir managers thus also need to consider possible consequences on water quality when adapting withdrawal strategies. Further investigations could consider this by coupling a water quality model to GLM.

Compared to previous studies the impact of a change in the withdrawal depth is lower, but our approach is closer to practical management and realistic conditions. Nevertheless, for the hypolimnion temperature and the end of summer stratification the possible impact of an adaption of the withdrawal strategy is comparable to the impact of climate change observed over the last 30 years. This implies that in some reservoirs adapting the used withdrawal strategy could help to mitigate the impacts of climate change. It is important to note that changing the withdrawal strategy can only reduce the impact of climate change within a limited range. It is thus unlikely that an adaptation of the withdrawal strategy can fully compensate for the negative effects of climate change. This is especially the case under the conditions of accelerated global warming [55].