Seasonal discharge response to temperature-driven changes in evaporation and snow processes

This study analyses how temperature-driven changes in evaporation and snow processes influence the discharge in large river basins. Using a distributed efficient hydrological model at high spatio-temporal resolution, we investigate the relative contribution of snow and evaporation. Comparing two 10-year periods (1980s and 2010s) in the Rhine allowed to determine the contribution of changes in snow, evaporation and precipitation to the discharge. Around half of the observed changes could be explained by the changes induced by snow (11%), evaporation (19%) and precipitation (18%), while 52% was driven by 5 a combination of these variables. Increased temperature scenarios show that seasonal changes in snow-dynamics could offset a fairly constant negative change in relative runoff induced by evaporation, but not during the melt season. This study shows how the combined effect of temperature-driven changes affect discharge. With many basins around the world depending on meltwater, correct understanding of these changes is vital.

showed that snow depth decreased over the majority of Europe. The Rhine river basin covers many different types of land cover, and is therefore representative for north-western Europe. Several studies have investigated the response of the basin under different climate scenarios (e.g., te Linde et al., 2010;Hurkmans et al., 2010;Pfister et al., 2004;Shabalova et al., 2003;Middelkoop et al., 2001), and the importance of melt water from snow and ice (Stahl et al., 2016). Yet, none of the studies have investigated the separate and combined response of evaporation and snow processes under increased temperature scenarios. 30 Spatially distributed modelling becomes increasingly viable, due to the increased computational power, gains in model performance when adding spatial information (Comola et al., 2015;Lobligeois et al., 2014;Ruiz-Villanueva et al., 2012), and increased availability of high resolution data (e.g., Huuskonen et al., 2013;Cornes et al., 2018;van Osnabrugge et al., 2017;C3S, 2017). However, the choice of spatial resolution can affect the model parameters (Melsen et al., 2016), and the sign of the simulated anomalies (Buitink et al., 2018). Besides, when finer spatial resolutions are used, the timestep should be 35 reduced as well, as the space and time dimensions are linked (Blöschl and Sivapalan, 1995;Melsen et al., 2016). However, simulations ran at high spatial (and temporal) resolutions usually greatly increase computational demand (for example, the study by Mastrotheodoros et al. (2020) took more than 6 × 10 5 CPU hours). This not only requires usage of high performance clusters, but also has considerable effects on the climate via increased power consumption (Loft, 2020). There is need for innovative hydrological models which can run on high spatio-temporal resolution without excessive computational demands, This study investigates the hydrological response to temperature-driven changes in evaporation and snow processes, testing our main hypotheses that both seasonal changes in snowmelt and enhanced evaporation will aggravate low flows, and that the changes will increase with temperature under realistic warming. We simulate the Rhine basin at high spatial (4 km) and temporal (1 hour) resolution using a calibrated version of the computationally efficient dS2 model, which is based on the 45 simple dynamical systems approach (Kirchner, 2009;Teuling et al., 2010). The model was run for two decades and increased warming scenarios, to show the response of the basin to changes in snow processes, evaporation and precipitation. Simulations performed at high spatial and temporal resolutions ensure that small scale variability is accounted for. By separating the temperature-driven effects on evaporation and snow processes, we can understand and quantify the relative importance and interaction of each process.

Methods
To understand how the discharge of large river basins responds to different changes in climate forcing, we used the computationally efficient distributed dS2 model (Buitink et al., 2019) to simulate the Rhine basin. The Rhine basin was selected because the climate and basin heterogeneity are representative for north-western Europe and many other basins globally. The model is based on the simple dynamical systems approach (Kirchner, 2009), and is extended with snow and routing modules. As dS2 55 requires actual evaporation data as input, we ran a soil moisture model prior to the rainfall runoff model to simulate the translation from potential evaporation (PET, calculated using the Penman-Monteith equation (Monteith, 1965)) to actual evaporation (AET). Since rootzone depth is an important yet highly uncertain parameter, we included simulations with rootzone depths ranging from 25 to 125 cm, with increments of 25 cm. Details on the two models and the calibration procedure can be found in the supplementary information (Text S1, Text S2).

60
All simulations are performed at a resolution of 4×4 km and at an hourly time step. The input data were obtained from the ERA5 reanalysis dataset (C3S, 2017). This dataset is globally available on a 0.25 × 0.25 • resolution and at an hourly timestep from 1979 to present. ERA5 data was interpolated to the model grid using bilinear interpolation. We selected two periods with equal length based on the maximum distance between available decades of ERA5 data: 1980-1989 and 2009-2018, referred to as 1980s and 2010s, respectively. As mentioned earlier, we simulate the Rhine basin as the climate and basin heterogeneity are 65 representative for north-western Europe.

Results
A first comparison of average temperature, precipitation and potential evaporation reveals considerable differences between the two periods ( Fig. 1). Over the entire Rhine basin, yearly average temperature has increased with more than 1 • C, from 8.1 • C to 9.3 • C between 1980s and 2010s. Largest differences are found in the eastern Alps, where average temperature 70 has risen by 1.5 • C (Fig. 1a, d). Average precipitation is lower in the 2010s over the majority of the Rhine basin, with the yearly average precipitation sums decreasing from 1146 mm to 1066 mm (Fig. 1b, e). Spatial differences in precipitation are, however, less homogeneous over the basin than the changes in temperature and potential evaporation. As a result of the increased temperatures, average potential evaporation also substantially increased from 607 mm to 678 mm from the 1980s to the 2010s, with the largest increases occurring in the northern parts of the basin (Fig. 1c, f).

75
A thorough validation is required in order to ensure that models simulate the correct sign and magnitude of the trends . Therefore, we validated dS2 on multiple levels: discharge of the total catchment, and snow and evaporation dynamics at the local scale. Discharge validation (Fig. 2a) shows that dS2 simulates the discharge with high KGE values. Panel b shows how the average discharge differs between the two periods, with lower discharges during late summer/autumn (in line with our main hypothesis) and higher flows during late winter in the 2010s. Discharge during the 2010s does not show as high 80 discharge values in June, and shows lower discharge values occurring later in the year. Kling-Gupta efficiencies for each period and several stations within the basin can be found in the supplementary information (Table S1).
Additionally, local validation is performed with point observations from the Rietholzbach research catchment in Switzerland (Seneviratne et al., 2012), by comparing simulated actual evaporation with observed evaporation from a lysimeter, and by comparing simulated snow storage with observed snow height measurements in Fig. 3. Due to data availability limitations, we had 85 to resort to our calibration period. Since dS2 was only calibrated on discharge, this still can be interpreted as validation. Both variables are correctly represented, and show similar variability as the observations, even at hourly timescale. The simulated evaporation generally shows a smoother signal than the observations. However, the simulated evaporation is based on relatively coarse ERA5 data, which could cause the lack of small scale variability. Snow storage shows a very similar pattern. It has to be noted that snow height observations cannot be directly converted into snow water equivalent, due to e.g. compaction. Yet, dS2 Given that dS2 is not calibrated on these variables, and the difference in spatial scale of the input data, this shows that dS2 is able to correctly simulate evaporation and snow processes.
In our first experiment, we aim to understand how the individual forcing variables affect the hydrological cycle. We swapped each forcing variable of the 1980s period with their timeseries of the 2010s period. The dots in Fig. 2d represent results for a 95 single timestep from separate simulations as the difference with respect to the 1980s. By summing the differences of the three swapped simulations, we make an estimate for the 2010s case. The ratio ∆P+∆Tevap+∆Tsnow ∆total is interpreted as the explained fraction, and is set to zero when they have different signs.
Investigating the difference of the "forcing-swapped" runs gives insight on how each variable affects the discharge (see Fig.   2e). Changing only the temperature affecting snow processes shows discharge differences mostly in the first half of the year.  The fraction explained in Fig. 2f-g gives an indication on the amount of interaction between the three forcing variables.
Values close to 1 indicate that there is little interaction, as the sum of the differences is able to explain all changes. The explained fraction is lowest during spring and late summer. During these periods, the storage conditions of the basin largely control the discharge response, either through snow storage or water available to generate runoff. In spring, changes in the available snow storage are the result of interactions between temperature and precipitation. In summer, discharge is controlled by water  Using the dS2 model, separate simulations of temperature effects on evaporation, snow and their combined effect allow us to understand which variable is causing the main changes. These time series are presented in Fig. 4a, including the 1980s run as reference. Three periods are highlighted, which represent typical discharge regimes: high winter discharge, the spring meltwater peak, and late summer low flow. For each of these periods, the change in discharge shows a roughly linear relation 120 with temperature increase, confirming our hypothesis. Surprisingly, for both the maximum and minimum discharges (panels b and d), the modified snow run shows behaviour opposite to both the modified evaporation run and the combined run. In these cases, change induced by evaporation is offset by the change induced by increased snowmelt or decreased snowfall, ensuring that the combined change is reduced with respect to the evaporation induced change alone. However, in the late spring discharge peak (panel c), both evaporation and snowmelt enhance the change, as they both show a reduction in mean flow during this 125 period. As a result, the mean flow of the combined run shows an even larger reduction in mean discharge, where even the peak from the 1980s has been largely diminished (see panel a). Substantial influence of rooting depth on the evaporation simulation is visible, yet the trend direction with increasing temperatures remains equal. Shallower rooting depth values induce more soil moisture stress since less water is available, leading to higher discharges.
To understand the cause of these changes, the change in generated runoff per model pixel is shown in Fig. 5a. This figure   130 shows that the majority of the basin produces less runoff for all three periods. Only the southern regions of the basin show a different response. In the winter period, these regions produced more runoff, resulting from the increased snowmelt. In the  other periods, only a few pixels produced more runoff. These pixels correspond to the glaciers in the Alps, which produced more meltwater resulting from the increased temperatures.
interplay of these variables. In Fig. 5b, the fraction of the basin that is dominated by one of these three options is plotted against the relative change in mean discharge for each period. As expected, the majority of the basin is controlled by the change induced by a change in evaporation (84-97%). As a result, the mean discharge is reduced by ±18%. Contrasting, a limited fraction of the basin (1-8%) is dominated by the change induced by snow, yet still has a considerable effect on the mean discharge: varying between -10% and 4%. Pixels in the basin where a combination of changed induced by snow and evaporation are required, 140 take only a very small fraction of the basin (<1%). Generally, these regions are at the transition between snow dominated and evaporation dominated regions. Overall, the change induced by T snow -despite the small contributing area-substantially affects the discharge. More details on the response for each temperature increase can be found in the supplement (Text S3).

Discussion and conclusion
We compared two periods of 10 years to investigate the relative importance of changes in temperature, evaporation and precip-145 itation. Over these periods of 10 years, most interannual variability is averaged out, allowing us to objectively investigate the effect of different temperatures on the hydrological response. Furthermore, the choice of spatial model resolution is a balance between data availability, computational time and underlying modelling concept. Here we selected a resolution of 4×4 km, so we can use the ERA5 forcing data without downscaling methods (adding uncertainty and potential errors), have short runtimes (simulating 10 years including all IO operations takes just over 5 minutes on a normal desktop), and apply the model at it's 150 proven spatio-temporal scale (±10 km 2 at hourly timestep). Contrasting, the study by Mastrotheodoros et al. (2020) used a much finer spatial resolution, but at the cost of enormous CPU times.
Temperature, evaporation and precipitation substantially changed from the 1980s to the 2010s in the Rhine. In the 2010s, basin average temperature was more than 1 • C higher, potential evaporation was almost 70 mm higher, and precipitation decreased with 80 mm. Discharge between these two periods was significantly different for 8 out of 12 months. Each individual 155 forcing variable can partly explain these discharge differences: 11% can be explained by the changed snowfall and melt dynamics, 19% is explained by the changed evaporation, 18% by the changed precipitation, and 52% is explained by interaction of these variables. Increasing the temperatures further results in decreased lower discharge values, where the role of snow driven changes shifts from enhancing to softening the discharge reduction induced by the increased evaporation. Less than 10% of the basin is dominated by the changed snow dynamics, yet it reduces the discharge with as much as 11%, depending on the season.

160
This study focusses on the Rhine basin, yet these results can be interpreted for the many different basins around the globe depending on both rain-and snowfall. With higher temperatures, increased melt from glaciers and snow packs can offset the discharge reduction from enhanced evaporation over the majority of the year. However, the season where runoff generation is reduced due to smaller snow storages (and potentially smaller glaciers) should be identified in each basin, as this part of the year is impacted the most. Many regions rely on "water towers" for their year-round water availability (Immerzeel et al.,165 2020), where the mountainous regions cover varying fractions of the basin. In many basins, more of the discharge originates from these water towers than in the Rhine basin, amplifying our results. Here, higher temperatures would likely imply even