Coupled daily streamflow and water temperature modelling in large river basins

Introduction Conclusions References

conflicts may arise between protecting freshwater ecosystems by enforcing ecological water temperature standards and risks to thermoelectric power production due to cooling water shortages. This has been reported for example for the River Rhine (Rutten et al., 2008), and for the Loire and Rhone rivers (Manoha et al., 2008) during the warm summers of 2003 and 2006 in Europe. In addition, significant impacts of water temperature of the River Rhine on electricity prices were found when water temperatures are above 22-23 • C (Boogert and Dupont, 2005).
For effective management of water and freshwater ecosystems, estimates of river discharge and water temperature at high temporal resolution, preferably on daily basis, are required. Both statistical (data-based) and process (physically based) models have been used to estimate river discharge and water temperature using climatic forcings. Statistical models (e.g. regression, stochastic models and neural networks) are appealing because they require only limited input variables (e.g. Mohseni et al., 1998;Muttiah et al., 1997;Chenard and Caissie, 2008;Ahmadi-Nedushan et al., 2007;Augustin et al., 2008). However, they are fitted for a specific historical period and, therefore, are limited in their application for forecasting and scenario studies, such as climate change impact assessments. In contrast, process models represent the physical processes that affect river discharge and water temperature and have been particularly useful for predictions of the effects of anthropogenic perturbations of model forcings and boundary conditions (land use change, thermal pollution, flow regulation) (e.g. Haag and Luce, 2008;Risley et al., 2010;St-Hilaire et al., 2003) and climate change (e.g. Morrison et al., 2002;Ferrari et al., 2007).
Although water temperature is generally most sensitive to the heat exchange processes at the air-water surface interface, changes in streamflow significantly affect water temperatures due to changes in thermal capacity (Edinger et al., 1968) and travel time, and dilution capacity for thermal effluents (Sinokrot and Gulliver, 2000;Moatar and Gailhard, 2006;Webb et al., 2003;. Therefore, an integrated hydrological and river temperature modelling approach is preferred which includes both heat exchange processes at the air-water surface interface and changes in thermal capacity and travel times due to streamflow changes. Although hydrological and process-based water temperature modelling approaches have been successfully applied for small-scale catchments and subbasins (e.g. Caissie et al., 2007;Haag and Luce, 2008;St-Hilaire et al., 2000), considerably less work has been done at large scales. To our knowledge, water temperature simulations on macro-hydrological scale have only been performed by van Beek et al. (2012). Limited studies have simulated both river discharge and water temperature for long (> 20-30 yr) time periods, which is required for scenario analyses and climate change impact assessments. In addition, realistic simulations of water temperature and discharge of rivers with different basin char-acteristics and anthropogenic impacts are needed to address large-scale water management issues.
In this study, we test the performance of an integrated framework with a physically (process) based hydrological and water temperature model to simulate daily river discharge and water temperature of large river basins in different hydro-climatic regions and with different anthropogenic impacts. A spatial resolution of 1/2 × 1/2 • was used as in Haddeland et al. (2011) and for which global forcing data are available . Several macro-scale hydrological models have simulated river discharge at this spatial resolution (e.g. Arnell, 1999;Oki et al., 2001;Alcamo et al., 2003), but most studies focus on monthly or annual mean estimates of river discharge.
Our modelling framework is based on the Variable Infiltration Capacity (VIC) macro-scale hydrological model (Liang et al., 1994) and the particle tracking River Basin Model (RBM) for water temperature (Yearsley, 2009) (Sect. 2.2). The modelling framework was applied globally; however, we focus on rivers situated in different hydro-climatic zones and characterized by different anthropogenic influences (Sect. 2.1). The performance of VIC-RBM was tested for a historical period 1971-2000 (Sects. 3.1 and 3.2). In addition, we tested the modelling framework for the Rhine and Columbia during warm, dry summer periods, and during the second half of the 20th century (Sect. 3.3). The sensitivity of simulated water temperatures to the boundary conditions (headwater temperature estimates) was studied (Sect. 3.4), as well as the sensitivity to streamflow (Sect. 3.5). The overall performance and major uncertainties of the hydrological and water temperature modelling approach for large-scale applications are discussed in Sect. 4.

Study basins
We focused on the Columbia, Mississippi (North America), Paraná (South America), Rhine, Meuse, Danube (Europe), Orange (Africa), Ob, Yenisey, Lena (Arctic Asia), Mekong, Yangtze, Yellow (Southeast Asia) and Murray-Darling (Australia) to test the performance of the hydrological and water temperature modelling approach. These basins are situated in different hydro-climatic zones with different anthropogenic influences and, therefore, represent different hydrological and thermal regimes ( Table 1). The Columbia, Mississippi, Rhine, Meuse, Danube, Yangtze and Yellow basin are situated in temperate climate zones and are influenced by transient runoff (mix of rainfall and springtime snowmelt), while the Murray-Darling and Orange are mainly fed by rainwater. The Paraná and Mekong are rain (monsoon) fed rivers located in tropical climate zones, while the Ob, Yenisey and Lena are Arctic rivers that are strongly affected by meltwater during spring and summer. River discharge and water Table 1. Major characteristics of selected study basins, data sources and number of monitoring stations with river discharge (Q) and water temperature (T w ) data used for validation of modelling approach. temperatures of the Columbia are heavily influenced by reservoirs, and parts of the Rhine, Meuse, Danube and Mississippi have a high level of thermal pollution due to cooling water discharges from thermoelectric power plants. The river basins vary in size, from almost 3.0 million km 2 (Mississippi and Ob) to 36 000 km 2 (Meuse). Another important criterion for selecting these study basins was the availability of monitoring stations with daily river discharge and water temper-ature records suitable for evaluating the performance of the modelling framework. Figure 1 shows linkages between the component models in the hydrological and water temperature modelling In brief, climate forcings and soil and vegetation parameters are used as input into VIC, resulting in simulated surface runoff and baseflow. The output (surface runoff and baseflow) is then provided to an offline routing model to simulate channel flows, depth, width and flow velocity on a stream reach basis. A routing model with a reservoir scheme simulates river discharge in the strongly regulated Columbia River. Climate forcings include air temperature, shortwave and longwave radiation, vapour pressure, density, pressure and wind speed disaggregated to the VIC grid cell and RBM reach level at a 3-hourly time step. In addition, daily channel flows, width, depth and flow velocity are used to force RBM. Other required inputs are an ordered stream network with defined river reaches (Yearsley, 2012), estimates of anthropogenic point heat sources and daily headwater stream temperature estimates (boundary conditions). The integrated modelling system simulates streamflow and water temperature in each of the grid cells.

Variable Infiltration Capacity (VIC) model and routing model
The Variable Infiltration Capacity (VIC) (Liang et al., 1994) is a grid-based macro-scale hydrological model that solves both the surface energy and water balance equations. The model represents subgrid variability in vegetation, elevation, and soils by partitioning each grid cell into multiple land cover (vegetation) and elevation classes. The soil column is commonly divided into three soil layers. Surface runoff and baseflow are routed along the stream network to the basin outlet with an offline routing model that uses the unit hydrograph principle within the grid cells and linearized St. Venant's equations to simulate river flow through the stream channel (Lohmann et al., 1998). For the Columbia River, which is highly affected by dams and reservoirs, we used the reservoir scheme of Haddeland et al. (2006), which is combined with the routing scheme of Lohmann et al. (1998) to obtain a more realistic representation of streamflow below the major reservoirs. The reservoir scheme runs at a daily time step, but was originally developed for analyses at coarser time scales. Hence, we calculated 10-day moving averages of daily regulated river discharge. Information about daily river depth, width and velocity is required for the water temperature simulations. The original VIC routing model (Lohmann et al., 1998) was therefore modified to calculate hydraulic characteristics based on power equations relating mean velocity, cross-sectional area and width to river discharge (Leopold and Maddock, 1953). Allen et al. (1994) obtained coefficients for these equations by fitting empirical relationships with river discharge using data from 674 stream discharge stations across the United States (Eq. 1a, b). As these stations are situated in a wide range of hydro-climatic zones, the assumption was made that these fitted relations can be applied to estimate the hydraulic characteristic of rivers in other regions and under different flow conditions as well. Flow velocity was estimated based on river discharge and cross-sectional area (Eq. 1c).
For the river reaches controlled by reservoirs, we assumed water surface elevation and, as a result, the depth (D res ) and width (W res ) to remain constant in time. In these river reaches Eq. (1c) becomes VIC and its routing model have been applied in the recent past at spatial scales ranging from 1/16 • (Elsner et al., 2010) to 1 • (Nijssen et al., 2001). The temporal resolution is flexible between hourly to daily step. The 1/2 • spatial resolution used in this study was selected as a compromise between the ability to resolve variations in river basin contributing areas and channel variations, and computational efficiency.

Stream temperature model RBM
RBM is a process-based one-dimensional stream temperature model that solves the 1D-heat advection equation using the semi-Lagrangian (mixed Eulerian-Lagrangian) approach (Yearsley, 2009). Because of the large-scale application, the advection term dominates, and dispersion was for that reason neglected. Water temperature is calculated for a specific stream segment based on the upstream water temperature and inflow into the stream segment, the dominant heat exchange at the air-water surface, and the inflow and temperature of water advected from tributaries. RBM was developed for subbasins of the Columbia River and has been applied by Yearsley (2012) to the Salmon subbasin (36 325 km 2 ) on a 1/16 • spatial resolution. In this study, modifications were made to apply RBM to larger river basins characterized by different thermal and hydrological regimes and anthropogenic impacts. To use RBM for thermally polluted river basins, modifications were made to incorporate anthropogenic heat discharges of thermoelectric power plants as advected heat sources. This results in the following 1D-heat advection equation: Fig. 1. Flowchart of the hydrological and water temperature modelling framework, presenting the links between the hydrological model (VIC), routing model, process-based water temperature model (RBM), water temperature regression model used to assess headwater temperatures (boundary conditions), and model input and output.
where ρ w = density of water (kg m −3 ); C p = specific heat capacity of water (J kg −1 • C −1 ); T w = water temperature ( • C); A x = cross-sectional area of river at distance x (m 2 ); H air-water = heat flux at air-water interface (J m −2 s −1 ); w x = stream width at distance x (m); Q trb = advected flow from tributaries or subsurface (m 3 s −1 ); T tbr = the difference between advected temperature from tributaries or subsurface, T tbr , and T w ( • C); Q effl = advected flow from heat dumps of thermoelectric power plants (m 3 s −1 ); T effl = the difference between the advected temperature from heat dumps of thermoelectric power plants, T effl , and T w ( • C); x = longitudinal distance along the axis of the river (m); t = time (s). The net exchange of thermal energy across the air-water interface (H air-water ) is determined using a one-dimensional implementation of the stream energy balance equation of Wunderlich and Gras (1967): where H air-water = net exchange of thermal energy across the air-water interface (J m −2 s −1 ); H s = shortwave solar radiation (J m −2 s −1 ); H rs = reflected shortwave solar radiation (J m −2 s −1 ); H a = longwave atmospheric radiation (J m −2 s −1 ); H ar = reflected atmospheric radiation (J m −2 s −1 ); H evap = evaporative heat flux (J m −2 s −1 ); H cond = conductive or convective heat flux (J m −2 s −1 ) (the flux resulting from temperature differences between the atmosphere and river); H back = blackbody radiation from the water surface (J m −2 s −1 ).

Estimation of the boundary conditions (headwater temperatures)
As part of the study described in Yearsley (2012), two methods for headwater temperature estimation were compared for the Salmon River (subbasin of the Columbia). One method uses daily soil temperature from VIC, and another method uses a nonlinear water temperature regression model (Mohseni et al., 1998) based on air temperature. Overall, the performance of the RBM model did not improve by using soil temperature to estimate headwater temperature. Given the widespread use of the regression model of Mohseni et al. (1998), we decided to use the latter approach for this study to estimate headwater temperature. The water temperature  Mohseni et al. (1998) describes the Scurve relationship between weekly water temperature and weekly air temperature according to where µ = lower bound of water temperature ( • C); α = upper bound of water temperature ( • C); γ = measure of the slope at inflection point (steepest slope) of the S-shaped relation ( • C −1 ); β = air temperature at inflection point ( • C); T w head = headwater temperature ( • C); T air = air temperature ( • C); tan θ = slope at inflection point (-). The four parameters of the regression model and time lag were fitted for 333 Global Environment Monitoring System (GEMS)/Water stations globally for the period 1980-2000 using least squares regression. We applied the nonlinear water temperature regression model on a daily time step by including a lag effect, as water temperature variations lag behind air temperature fluctuations at short time scale (hourly, daily basis) (Erickson and Stefan, 2000;Jeppesen and Iversen, 1987). For each station, the optimal lag parameter was estimated by calculating correlation coefficients between water temperature and smoothed air temperature (T air smooth ) for various lag parameter values (λ) using Eq. (6), and selecting the λ for which the correlation coefficient was highest.
In a next step, the fitted parameter values µ, α, γ and β were interpolated using ordinary kriging, resulting in 1/2 × 1/2 • interpolated grids. The time lag at which water temperature variations follow air temperature variations increases with stream depth (Stefan and Preudhomme, 1993) and thus with river discharge. The lag parameter was, therefore, spatially interpolated using gridded river discharge simulations produced by VIC in combination with an empirical relationship between lag parameter and river discharge (fitted for all stations). An overview of the mean and range (minimummaximum) in fitted parameters for all study basins (Table 2) shows that the fitted values of the Mohseni parameters vary between the different study basins (in particular µ). The lag parameter (λ) is generally constant within and between the different basins (between 0.09 and 0.12 for all basins).

Application of hydrological and water temperature modelling framework
To test the performance of the modelling framework, the VIC-RBM framework was applied globally for the period 1970-2001 (including a spin-up period of one year). The models were forced with daily (24 h mean) values of precipitation, minimum and maximum surface air temperature and wind speed from the global gridded 1/2 × 1/2 • meteorological dataset developed within the EU FP6 Water and Global Change (WATCH) project (Weedon et al., 2010. VIC was applied using the elevation and land cover classifications (elevation, vegetation, and soil characteristics) described in Nijssen et al. (2001), disaggregated to 1/2 × 1/2 • spatial resolution. In their study, calibration on soil characteristics was performed for selected large river basins globally (including the Mississippi, Columbia, Danube, Paraná, Yellow, Yangtze, Mekong, Yenisey, Lena and Ob). Calibrated parameters values were subsequently transferred to other basins based on climate characteristics. The global DDM30 routing network (Döll and Lehner, 2002) was used for the lateral routing of streamflow and to create an ordered river network for the RBM water temperature simulations. For the Columbia basin where river discharge and water temperature are highly impacted by reservoirs, we used information about dams from the University of New Hampshire updated according to the World Register of Dams (ICOLD, 2003), as described by Haddeland et al. (2006).
To get realistic water temperature simulations in thermally polluted river basins (Mississippi, Rhine, Meuse, Danube), estimates of thermal discharges of thermoelectric power plants are required as input into RBM. We used gridded (1/2 × 1/2 • ) estimates of global thermoelectric water consumption and water withdrawal for the 20th century (Flörke et al., 2011;Voß and Flörke, 2010) to estimate return flows from thermoelectric water diversions (Q effl ). Because gridded data for the difference in temperature between cooling water temperature and river water temperature were not available, we assumed that the temperature of return flow (T effl ) was on average 3 • C higher than the inlet river water temperature (T w ). This value was also selected based on an average estimate for the Rhine River (Icke et al., 2006) and based on standards for heat discharges in the United States, which are written under the requirements of the Clean Water Act, and limit the T effl to 3 • C for most states. In addition, overall best results of daily simulated water temperature were obtained under a T effl of 3 • C when we tested this for the thermally polluted basins Mississippi, Rhine, Meuse and Danube with values ranging from 2 to 10 • C (van Vliet et al., 2012). Using information about the dominant cooling type in each grid cell, Q effl and T effl , gridded (1/2 × 1/2 • ) datasets of thermal discharge were calculated for the period 1971-2000 and these were used as input into RBM.

Evaluation of hydrological and water temperature modelling framework
Observed daily river discharge and water temperature records for selected monitoring stations in the study basins were used to evaluate the VIC-RBM simulations. Daily mean series of river discharge were provided by the Global Runoff Data Centre (GRDC; http://www.bafg.de/GRDC/EN) for the period 1971-2000. For the Yangtze River, we tested Table 2. Mean values and range (minimum-maximum) of interpolated parameters of nonlinear water temperature regression model and time lag used for estimating daily headwater temperature in the study basins.  (Table 1). In general, the water temperature observations represent daily instantaneous (spot) measurements, taken approximately 0-1 m below the water surface around mid-day. These instantaneous water temperature measurements that are taken at the water surface were related to simulated water temperature, which are cross-sectional averages of mean daily water temperature. Although there are vertical variations in water temperature, previous studies have shown that instantaneous observations of water temperature taken near surface are generally representative of the mean water temperature as vertical and lateral mixing of water is often very strong in large rivers (Mackay and Mackay, 1975;Liu et al., 2005).

Study basin
To quantify the performance of VIC and RBM for daily river discharge and water temperature simulations, we used the root-mean-square error (RMSE) and mean bias (BIAS). In addition, the Pearson correlation coefficient (r) was calculated to quantify the linear dependence between simulations and observations values. For river discharge, normalized values of RMSE and BIAS were calculated (NRMSE and NBIAS henceforth) by dividing by the mean observed river discharge values. The equations for the selected performance coefficients are where P i = predicted value at time step I ; O i = observed value at time step I ; O = average of daily observed value; P = average of daily predicted value; n = number of data pairs to be compared. For the Columbia and Rhine basins, more detailed and longer-term daily water temperature datasets were available. This allows a validation over the simulated water temperature trends over the entire 1971-2000 period and for warm, dry summers, specifically, when critically high water temperatures and low water availability occur. We focused on the warm summers of 1992 and 1994 in the Rhine, and the summers of 1998 and 1999 in the Columbia. These summers were selected, as highest water temperature values were observed, considering the average of all water temperature records in the river basin.

Sensitivity of simulated water temperature to headwater temperature
For coarse spatial resolution, uncertainties in the estimates of the boundary conditions are expected to propagate over large distances. A sensitivity analysis was therefore performed to assess the impact of uncertainties in headwater temperature estimates on simulated water temperatures at different spatial resolutions: 1/2 • , 1/4 • and 1/8 • . We focused on the Rhine and Meuse basins in Western Europe, because these basins are the smallest study basins and have reasonable running times at 1/8 • resolution. The routing and water temperature simulations for the Rhine and Meuse on 1/4 • and 1/8 • were performed by using the river routing networks derived from HY-DRO1K (Wu et al., 2011). We compared water temperature simulations produced by using an overestimated headwater temperature of +2.0 • C with simulations based on the original gridded headwater temperature estimates (reference case) at 1/2 • , 1/4 • and 1/8 • resolutions for the period 1971-2000.

Sensitivity of simulated water temperature to river discharge
In addition to headwater temperature, we assessed the impact of uncertainties associated with the hydrological model output and changes in river discharge on the simulated daily water temperature. We compared simulated water temper-ature for the reference case with simulated water temperature under a change in streamflow of −25 %, −50 %, +25 % and +50 %. Simulations with RBM were performed for the period 1970-2000 (including one year spin-up) assuming a constant decrease and increase in both daily simulated runoff and baseflow from VIC of −25 %, −50 %, +25 % and +50 % compared to the reference conditions.

Performance of daily river discharge simulations
The spatial patterns of simulated mean annual river discharge of the study basins (Fig. 2) generally show a close correspondence with the mean observed river discharge (small circles). For some downstream stations in the Orange and Murray-Darling basins, VIC overestimated river discharge. Part of this overestimation can be explained by anthropogenic water withdrawals (e.g. for agriculture, energy, manufacturing and domestic water use) which are relatively high in these basins. This results in lower observed river discharge values compared to the simulated values (which do not include anthropogenic water extractions). This overestimation is also reflected by relatively high values in NBIAS (> 2) and high values in the NRMSE (> 3) for both the Murray-Darling and the Orange river basins (Table 3). For the Ob, Yenisey, Lena, Mekong and Yangtze, a slight underestimation was found resulting in small negative values of NBIAS. However, values of NRMSE were generally low and r was relatively high (r > 0.75 for most of these basins). The use of the reservoir scheme resulted in a distinct improvement (significantly smaller bias; p < 0.05 using paired t-test) in the simulated river discharge of the highly regulated Columbia River (Fig. 3). This is reflected by a lower value of mean NBIAS and NRMSE (+0.3 and 1.4, respectively) compared to the simulation without the reservoir scheme (+0.5 and 2.0, respectively). Although the onset of the discharge peak in spring is somewhat too early, e.g. at The Dalles (Fig. 3), the hydrologic regime is represented more realistically when the reservoir scheme is included. Although the hydrologic regimes of some other rivers, like the Mekong and Ob, are also slightly impacted by reservoirs, we obtained a quite realistic representation of daily river  discharge with VIC without the use of a reservoir scheme for these rivers (mean NBIAS = −0.1; r = 0.91 for Mekong; NBIAS = −0.1; r = 0.76 for Ob; Fig. 3). Daily variability in river discharge was slightly underestimated for some upstream stations in the River Rhine (Fig. 3), but in general a realistic representation was found (mean NBIAS = +0.1; r = 0.76). This indicates that the hydrological model is suitable for simulating daily discharge in river basins situated in different climate zones (temperate, tropical and arctic) and with different anthropogenic impacts.

Performance of daily water temperature simulations
The spatial patterns of simulated mean annual water temperatures within the study basins averaged over the period 1980-2000 (Fig. 4) show pronounced increases in water temperature from the upstream to the downstream parts of most river basins (except for the Lena, Ob and Yenisey and Paraná River that flow to high latitude). In general, the simulated values of the grid cells correspond closely with the observed mean annual water temperatures for the different stations (circles) along the streams. For the Columbia River a significant improvement (p = 0.03) was found by using corrected geometry-streamflow relations (Eq. 2; Sect. 2.1.1) for the grid cells where reservoirs are located. Without these corrected relationships, the onset of the rising and falling limb in the simulated thermal regime is too early in the season (see Fig. 5, station Grand Coulee), because the calculated depth and width are underestimated, resulting in an underestimation of the thermal capacity of the stream segment. Furthermore, flow velocity is overestimated, which can result in greater influence of uncertainties in headwater temperature estimates on simulated water temperature. The improvement in model performance was also reflected by lower values of RMSE (mean value of 2.8 • C versus 3.5 • C) and higher values of r (0.88 versus 0.77) for a model run with corrected relations (Eq. 2) compared to the run based on uniform geometry-stream flow relations (Eq. 1).
The implementation of point sources of heat effluents also resulted in a significant (p < 0.05) improvement in model performance for thermally polluted rivers like the Rhine, Meuse, Danube and Mississippi. Without implementation of heat effluents, the simulated water temperatures are underestimated (negative bias) compared to observed water temperature, as these reflect the "naturalized" water temperature (see Fig. 5, Rhine, Koblenz). The improvement was reflected by decreases in negative BIAS, lower RMSE and slightly higher values of r for these river basins.
For some of the tropical and arctic basins, like the Mekong and Lena, only a limited number of water temperature measurements was available to test the performance of RBM on a daily time step. However, the simulated water temperature series generally fell between the observations and the variability in water temperature throughout the year was well simulated, as shown for the Mekong (Pakse) (Fig. 5). This was also found for the other eight water temperature monitoring stations along the Mekong, although slightly overestimations occurred for the most upstream stations (mean BIAS = +1.5 • C). For the Lena, which is strongly affected by meltwater, the annual cycles in water temperature were simulated realistically during the snowmelt period. However, the steepness of the falling limb during August-October was on 4313 Fig. 4. Spatial patterns of simulated (grid cells) and observed (circles) mean annual water temperatures for study basins. average too high and the decrease started too early in the season. This might be explained by an underestimation of the discharge peak for the Lena during summer (reflected by negative NBIAS; Table 3), and associated underestimation of the thermal capacity. Due to ice and meltwater inflow, water temperatures in spring were slightly overestimated for some years, but overall the timing and magnitude of the rise in water temperature of the Lena during summer were simulated realistically for most of the years during the evaluation period.
The scatter plots and histograms of the simulated versus observed daily water temperature (Fig. 6) show that simulated water temperature values match the observed values reasonably well for most of the stations. For some stations the correlation coefficients are high (r > 0.80; Yellow, Murray-Darling) or very high (r > 0.90; Snake (Columbia), Missouri, Arkansas (Mississippi), Rhine, Meuse). For the Mekong and Orange, the correlations were somewhat lower (r = 0.72 and r = 0.78), although the distributions of daily simulated and observed values correspond closely. In addition, the seasonal signal in water temperature for both rivers is weaker, resulting in a lower signal-to-noise ratio and thus a lower correlation coefficient.

Long-term water temperature 1971-2000 and performance for warm, dry summers
For the evaluation of the simulated water temperature for the Columbia and Rhine basins during the entire 1971-2000 period, we focused on the stations Anatone and Lobith for which long-term daily observed water temperature series were available. Annual mean (± one standard deviation), and annual maximum values in simulated and observed water temperature are shown for both stations for the period 1971-2000 (Fig. 7, left). The annual simulated water temperatures match closely with the observations for the entire period. The simulated daily water temperature and river discharge during the warm summers of 1998 and 1999 in the Columbia River and summers of 1992 and 1994 in the Rhine River also showed an overall realistic performance of the modelling approach. For the Rhine, the variability in river discharge and water temperature was slightly overestimated (Fig. 7  right). In addition, the simulated water temperature values were slightly underestimated for some downstream stations in the Rhine during the period with the highest water temperatures (mean BIAS of −1.2 • C). This was also found for some stations in the Columbia and Snake rivers (overall mean BIAS for all stations of −0.8 • C). River discharge was also simulated realistically (mean NBIAS for both summers of +0.05 for the Columbia and +0.18 for the Rhine), resulting in an overall a realistic performance of the hydrological and water temperature modelling framework for warm summers.

Sensitivity of water temperature to headwater temperature estimates
The impact of a positive bias in headwater temperature of +2.0 • C generally shows the largest impact in the upstream parts of the Rhine and Meuse rivers and declines in the downstream direction more rapidly for the finer-resolution simulations (1/4 • ) compared to the coarse-resolution simulations (1/2 • ) (Fig. 8a, b). In particular for the Meuse, which is the smallest basin, the impact of biases in headwater temperature estimates differs for the 1/2 • , 1/4 • and 1/8 • simulations. Assuming a positive bias in headwater temperature of +2.0 • C, the impact on simulated water temperature 175 km downstream in the Meuse is on average +1.0 • C (51 %) for the 1/2 • Fig. 6. Scatter plots and histograms of daily simulated river temperature versus daily observed river temperature for selected stations in the study basins for 1980-2000 period. Histograms on the vertical axis are for simulated values and histograms on horizontal axis for the observed water temperature.
simulations compared to +0.4 • C (20 %) for the 1/8 • simulations. For the Rhine, the impact at 175 km downstream is larger; +1.4 • C (71 %) and +1.2 • C (59 %) for the 1/2 • and 1/8 • resolution simulations, respectively. This is due to the higher flow velocities of the Rhine compared to the Meuse, which results in shorter travel times from headwaters to the downstream site. In addition, higher water depths for the Rhine compared to the Meuse result in slower response rates to atmospheric conditions and, consequently, larger propagation of uncertainties in head water temperatures along a longitudinal section.
Higher-resolution simulations also resulted in an overall higher quality of the daily water temperature simulations, although the differences are small. The mean RMSE for the stations in the Rhine decreased from 2.3 • C (at 1/2 • spatial resolution), to 2.1 • C (1/4 • ) and 2.0 • C (1/8 • ). The correla-tion coefficients between the observed and simulated daily values were also higher for most stations along the Rhine for the 1/8 • compared to 1/2 • resolution simulations, although the differences were very small (mean r of 0.944 (at 1/2 • spatial resolution), 0.952 (1/4 • ) and 0.953 (1/8 • ); Fig. 8c).

Sensitivity of water temperature to river discharge simulations
Results of the sensitivity analyses showed moderate impacts of changes in river discharge on mean annual water temperature (average value for all basins of +0.2 • C for −50 % change in river runoff and −0.1 • C for +50 % change; Fig. 9). However, pronounced impacts in the low and high water temperature range were found. A decrease in river runoff of 25 % and 50 % results in significantly ( Fig. 8. Impacts of spatial resolution on propagation of uncertainties in headwater temperature estimates on simulated water temperature along the river course for the Rhine and Meuse, and correlation coefficients between daily simulated and observed water temperature for stations in the Rhine (from upstream station Diepoldsau to most downstream station Lobith) at 1/2 • , 1/4 • and 1/8 • spatial resolution.
impact of 50 % runoff decrease is −0.4 • C) and significantly (p < 0.01) higher maximum water temperature during summer (average impact of 50 % decrease is +1.2 • C). This is due to a smaller thermal capacity and travel time which increases the sensitivity to atmospheric warming and cooling. The impacts of changes in river runoff on minimum wa-ter temperatures are very limited for basins at high northern latitude because minimum water temperature values remain around freezing point. An increase in streamflow has an opposite impact on water temperature. A +50 % river runoff increase results in an average impact for all basins of +0.3 • C in minimum water temperature and −0.6 • C in maximum water Fig. 9. Impacts of changes in river runoff of −50 % and +50 % on mean, minimum and maximum annual water temperature.
temperature, which is also significant (p < 0.01). Probability distributions of daily water temperature for the reference case and under a change in river runoff of −50 %, −25 %, +25 % and +50 % also show highest impact of changes in streamflow in the low and high water temperature range (Fig. 10).
In particular, decreases in streamflow result in substantial increases in water temperature in the high range. For the Danube station, a decrease in river runoff results in higher (rather than lower) minimum water temperatures, and strong increases in high water temperatures were found. This was also found for several other stations in thermally polluted basins, and this could be explained by a reduced dilution capacity for thermal effluents under decreasing runoff.

Discussion and conclusion
We used a physically based modelling framework with the VIC macro-scale hydrological model and process-based water temperature model RBM. The modelling framework was modified to include impacts of reservoirs and heat effluents of thermoelectric power plants and was tested for large river basins in different hydro-climatic zones and with different anthropogenic impacts. Based on our analysis, we conclude that the coupled hydrological and water temperature modelling framework is suitable to simulate daily river discharge (median normalized BIAS = 0.3; normalized RMSE = 1.2; r = 0.76) and water temperatures (median BIAS = −0.3 • C; RMSE = 2.8 • C; r = 0.91) realistically on daily time step over long (> 20 yr) periods and on large spatial scales. A similar performance was found during critical periods (warm, dry summers), which indicates that the modelling approach has potential www.hydrol-earth-syst-sci.net/16/4303/2012/ Hydrol. Earth Syst. Sci., 16, 4303-4321, 2012 for risk assessments and studying climate change and other anthropogenic impacts on daily river discharge and water temperature in large river basins. In addition, the modelling framework shows possibilities for incorporating other water quality parameters. Yearsley (2012) compared the performance of the VIC-RBM modelling framework applied to the Salmon (subbasin Columbia) with other previous catchmentscale water temperature modelling studies, and concluded that the modelling framework performs as well or better than statistical water temperature models and within the range of site-specific applications of process-based models. Van Beek et al. (2012) simulated water temperature on a global scale (without calibration) with mean absolute errors in daily simulations ranging from 1.6 • to 7.6 • C, which are comparable or slightly higher than those obtained in our study.
As the focus of our study is on global river basins, local conditions such as effects of topography, vegetation and groundwater recharge, which can significantly influence river discharge and water temperature in small streams (e.g. Brown, 1969;Sridhar et al., 2004;Cristea and Burges, 2010), were disregarded. Although this contributes to uncertainties in river discharge and water temperature simulations, impacts of processes like groundwater advection, energy exchange between river bed and water interface, dispersion of heat and local conditions (topography and vegetation) in main rivers are relatively small at this large scale (Liu et al., 2005;Caissie, 2006).
For river flow, major sources of uncertainties are from meteorological forcing data and the parameterization of the soil and land cover (vegetation) characteristics. Uncertainties in simulated river flow then affect simulated water temperatures, especially during warm, dry conditions. Results of the sensitivity analyses showed significant impacts of river discharge (thermal capacity) on water temperature in the low and especially high water temperature range (Sect. 3.5; Figs. 9 and 10). These results correspond with previous physically-based and statistical water temperature modelling studies that have found a pronounced impact of river discharge on especially high temperatures Sinokrot and Gulliver, 2000;Bartholow, 1991).
For water temperature, we also found a relatively high sensitivity of simulated water temperatures to the boundary conditions (headwater temperatures) on a 1/2 • spatial resolution (Sect. 3.4). This highlights the importance of realistic estimates of headwater temperature for large-scale applications and coarse spatial resolutions. The effects of headwater temperature are larger in the upstream parts of the basins, although the magnitude of impact also increases with higher flow velocity due to the shorter travel time from headwater to the downstream site (Yearsley, 2012). Both the scale and time of travel from the headwaters determine the propagation and impact of incorrect values of the boundary conditions on the simulated water temperatures downstream. Increasing the spatial resolution would probably improve the quality of the water temperature simulations, by decreasing the impact of biases in headwater temperature estimates on the downstream reaches. However, only relatively small improvements in model performance were found for the Rhine and Meuse on 1/4 • and 1/8 • compared to 1/2 • while storage of input and output data and running times drastically increased.
We conclude that the integrated physically based VIC-RBM modelling framework is suitable to simulate daily river discharge and water temperature in large basins realistically. The modelling approach has potential for decision support (for example for water quality planning on a large scale) and potential to perform risk analyses and study climate change impacts for large river basins on a continental and global scale.