Articles | Volume 24, issue 10
Research article
28 Oct 2020
Research article |  | 28 Oct 2020

Anthropogenic influence on the Rhine water temperatures

Alex Zavarsky and Lars Duester

River temperature is an important parameter for water quality and an important variable for physical, chemical and biological processes. River water is also used by production facilities as cooling agent. We introduced a new way of calculating a catchment-wide air temperature using a time-lagged and weighed average. Regressing the new air temperature vs. river water temperature, the meteorological influence and the anthropogenic heat input could be studied separately. The new method was tested at four monitoring stations (Basel, Worms, Koblenz and Cologne) along the river Rhine and lowered the root mean square error of the regression from 2.37 C (simple average) to 1.02 C. The analysis also showed that the long-term trend (1979–2018) of river water temperature was, next to the increasing air temperature, mostly influenced by decreasing nuclear power production. Short-term changes in timescales < 5 years were connected with changes in industrial production. We found significant positive correlations for the relationship.

1 Introduction

River water temperature (Tw) greatly influences the most important physical and chemical processes in rivers and is a key factor for river system health (Delpla et al.2009). Tw also confines animal habitats (Gaudard et al.2018; Isaak et al.2012; Durance and Ormerod2009), regulates the spread of invasive species (Wenger et al.2011; Hari et al.2006) and is therefore an important ecological parameter. River water is not solely important from an environmental perspective but is also of very significant interest for the economy. Especially for energy-intensive industries such as power plants, oil refineries, paper or steel mills, river water is an important resource. Its availability is a basic requirement for the location of the facilities (Förster and Lilliestam2010). As a cooling agent, given a 32 % energy efficiency, 68 % of the energy is discharged through the cooling system into the respective stream (Förster and Lilliestam2010). This leads to a significant heat load, even on large rivers such as the Rhine (IKSR2006; Lange2009). As a consequence, anthropogenic effects such as industrial heat input, river regulation or streamside land change can contribute significantly to the heat budget of a river and, furthermore, on Tw (Cai et al.2018; Gaudard et al.2018; Råman Vinnå et al.2018). The natural influences on Tw are the following: (1) meteorology, including sensible heat flux, latent heat flux and radiative heat fluxes; (2) source temperature, which describes the origin of the water, e.g., snow-fed, glacier-fed and groundwater-fed; (3) hydrology, which influences the water temperature through the amount of water and the flow velocity, together with the change in riparian vegetation; and (4) ground heat flux.

Dependent on data availability, computing power, accuracy and the questions asked, Tw can be modeled in different ways. The common options are statistical models and physically based models.

A physically based Tw model (Sinokrot and Stefan1993) usually parameterizes or estimates the meteorological and ground heat fluxes and adds anthropogenic heat input. Each modeled heat flux is then applied to the water mass and initialized with the starting and boundary conditions of source temperature and discharge. However, it is difficult to obtain a good estimation of these different terms over a larger catchment area. Hybrid models are in between physically based and statistical models. They use the physical formulation of fluxes but determine their parameters stochastically (Piccolroaz et al.2016). Hybrid models can reproduce river water temperatures better than simple statistical models (e.g., linear regression; Toffolon and Piccolroaz2015). Their approach includes more parameters and, thus, is more complex. However, a simple hybrid model with three parameters is comparable to a statistical model with the same number of parameters. Statistical models use air temperature (Ta) as a proxy for sensible, latent and radiative heat fluxes (ground heat flux can be neglected) and establish a TaTw relationship through regression. Ta is rather easily available from meteorological networks or reanalysis products. This is an established method, and depending on the complexity, linear or exponential models are used (Stefan and Preudhomme1993; Mohseni et al.1998; Koch and Grünewald2010). Generally, exponential models deliver better results with temperature extremes. However, they lack the distinct separation between contribution to Tw from anthropogenic heat input and natural influences. Using linear models, Markovic et al. (2013) show that between 81 % and 90 % of the Tw variability can be described by Ta. Furthermore, the authors showed that 9 %–19 % can be attributed to hydrological factors (e.g., discharge). The study was conducted for the Danube and Elbe basins using data from 1939 to 2008. These two rivers have comparable discharges and catchment areas to the Rhine river, which could mean these results are transferable. These, although simple, linear models are able to clearly separate the different influences on Tw. Another development is spatial statistical models. They correlate various landscape variables (e.g., elevation, orientation, hill shading, river slope, channel width, etc.) across the catchment area and aim to statistically determine their influence on Tw at a certain point. These correlations can be across any distance and do not have to satisfy flow connection or direction in the river system. As a prerequisite, a detailed knowledge about the river system and its characteristics is needed (Jackson et al.2017a, b). An improvement on spatial statistic models was recognizing rivers as a network of connected segments with a definite flow direction (Hoef et al.2006; Hoef and Peterson2010; Isaak et al.2010; Peterson and Hoef2010; Isaak et al.2014). The correlation of the variables (e.g., Ta, Tw discharge, etc.) which influence other Tw is weighed on their flow connectivity and Euclidean distance or flow distance. These models can also include time lag considerations using temporal auto correlation (Jackson et al.2018). Artificial neural networks (ANN) are a subset of the statistical models and are used when an incomplete understanding of most contributing processes is given (Hassoun1995). ANN use a sample data set to train artificial neurons the relationship between input (e.g., air temperature) and output (Tw) (Zhu et al.2018).

We used a simple linear regression model (transferable to other streams) to investigate the temperature changes in the Rhine, over 40 years, which had been influenced by 12 nuclear power plants (NPPs) along the Rhine. These NPPs had caused, for decades, the largest part of anthropogenic heat input (Lange2009). The nuclear power production increased in the 1970s and 1980s and reached a peak in the mid-1990s. After the Fukushima disaster (2011), the German government decided to exit nuclear power production, and the first NPPs were shut down. After this political decision, a distinct drop in nuclear power production was visible on top of the already decreasing production rates. By July 2019, eight NPPs remained operational in the catchment area of the Rhine that were using (partly) river water as a cooling agent. In this publication, we hypothesized that, next to environmental factors, the long-term decrease in power production, which is coupled to a decreasing use of river water as cooling agent, has a long-term (> 10 years) impact on Tw of the Rhine. Short-term economic changes, observable in the change in the gross domestic product (GDP), may influence Tw on shorter timescales (< 5 years). As several industrialized hot spots are present along the river, this impact might be spatially heterogeneous. Using the nuclear power production and GDP data, we also investigated the varying anthropogenic impact on Tw along the Rhine at four monitoring stations (Basel, Worms, Koblenz and Cologne).

2 Methods

We investigated the change in anthropogenic heat input and its spatial and temporal heterogeneity along the Rhine, combining ideas from spatial correlation models, to develop a new method for calculating a representative catchment air temperature (Tc). Tc and discharge at the measurement station Q were used in a multiple linear regression TcTw (Eq. 1). The resulting regression coefficients a1, a2 and a3 describe the magnitude of the respective influences (anthropogenic heat input and meteorological and hydrological influences).

(1) T w = a 1 + a 2 T c + a 3 Q .

Using an improved calculation method for Tc, which includes catchment-wide averaging with river size weighing and a time lag, the regression should deliver a better estimate for a1, a2 and a3.

The model was run on a Tw time series, from 1979 to 2018, measured at four Rhine stations, namely Basel in Switzerland (CH) and Worms, Koblenz and Cologne in Germany (DE). From 1979 to 2018, several changes in anthropogenic heat input to the Rhine catchment area occurred, making it an interesting data set to study. Webb et al. (2003) and Markovic et al. (2013) have shown that Q is inversely related to Tw and an important factor in the TcTw relationship. Additionally, it may function as a measure of how fast the water mass responds to changes in Tw. Ground heat flux, ground water influx and heat generation due to friction were not included in this model because of the comparably small influence (Sinokrot and Stefan1993, for the Mississippi; Caissie2006, as a review article).

Using the multiple regression (Eq. 1), we especially investigated the change in a1 over time, which we call, in this study, the Rhine base temperature (RBT). This temperature represents Tw, without the influence of meteorology (Ta) and discharge (Q). RBT was defined as an indicator for industrial heat input and the use of Rhine water as cooling agent, in case both are mostly independent of Ta and Q.

2.1 Water temperature and discharge

We used a data set of daily averaged Tw and Q from 1979 to 2018 provided by WSA (2019), BfG (2019), LfU (2019) and BAFU (2019). The original data sets had a 10 min sample frequency and were averaged to a daily output. Table 1 lists the respective stations along the Rhine, stream kilometers, data availability, the important tributaries upstream and the data provider.

BAFU (2019)LfU (2019)BfG (2019)WSA (2019)

Table 1Monitoring stations used in this study from Switzerland (Basel) to the lower Rhine region (Cologne, Germany). The location is shown as a Rhine kilometer (KM), and the time period, the important upstream tributary and the data source are listed.

Download Print Version | Download XLSX

Tw was measured by platinum resistivity sensors (Pt100). The accuracy of these sensors is commonly ±0.5C, but the precision, which describes the ability to detect temperature changes, is 0.05 C. As we focused on the change of Tw over time and did not compare the absolute temperature, the accuracy was not essential, and the precision of the sensors was sufficient for this study. Measurement uncertainties (e.g., depth and location of the sensor) did not influence the calculation with respect to the aim of this study, as long as the measured Tw was a linearly dependent proxy for the average river temperature. Q was provided as daily averages in m3 s−1 by the source in Table 1 and usually calculated from a river stage nearby.

The original data sets were provided by the state and federally operated monitoring stations, which usually run backup measurement systems. They verified the data, and we additionally screened the data set for suspicious features. Missing data points for up to 1 week were linearly interpolated. Longer or recurring data outages were not given.

2.2 Air temperature

Ta is retrieved from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis Model ERA5. It provides an hourly time resolution of the 2 m Ta on a 14 by 14o grid. The data set is available from 1979 to 2018. We took the hourly Ta output and calculated a daily mean for each grid point between 1979 and 2018 to fit the time resolution of Tw.

2.3 Nuclear power plants

The annual electrical power production (EPP) by NPPs is available from the International Atomic Energy Agency (IAEA) Power Reactor Information System (PRIS; IAEA2019). A total of 12 NPPs (1986–1988) were online in the Rhine catchment area, and eight remained operational by July 2019. All shutdowns were undertaken in Germany. In this study, separate reactor blocks of the same plant NPPs were combined.

The heat input (HI) by NPPs to the Rhine was calculated for each monitoring station using the conversion factor c and the yearly EPP; see Eq. (2). NPPs with an exclusive river water cooling system have a conversion factor of 3, which is based on the power efficiency of electricity generation (Lange2009). Other factors are estimated depending on the cooling system and personal communication with technicians from NPPs. If no conversion factor was available, a constant HI was assumed (Lange2009).

(2) HI [ GW ] = c EPP [ GWh ] 365 24 [ h ] .

The NPPs, their conversion factor and, if applicable, the constant HI are shown in Table 2. The time series of upstream HI by NPPs for each monitoring station is shown in Fig. 1.

Figure 1Heat input by upstream nuclear power plants (NPPs) from 1969 to 2018 at each monitoring station.


Table 2NPPs included in this study. The conversion factor describes the conversion from electrical power production (EPP) to heat input (HI). If cooling towers are installed, a constant heat input was used for the calculation based on Lange (2009). Note: NA – not available.

Download Print Version | Download XLSX

Calculated temperature change

We calculated the expected change ΔTw based on a change in HI (Δ HI) from NPPs using the average discharge Q, the heat capacity of the water cp and the water density ρ; see Eq. (3).

(3) Δ T w = Δ HI c p Q ρ .

This approach follows the idea that the contribution of NPPs significantly alters the Tw and only influences the RBT fraction.

2.4 Gross domestic product (GDP)

The GDP for the adjacent German federal states is obtained from Arbeitskreis VGdL (2019a, 2007b). Due to changes in the calculation method of the GDP before and after the German reunification (1990), two separate data sets were used. For this study, only the GDP change in the secondary sector (construction and production) was taken into account.

The RBT, if compared to the GDP, was filtered using a 10th order Butterworth bandpass filter. The sampling rate of the GDP was 1 yr−1. We used 1.1 yr−1 as higher and 0.05 yr−1 as lower cutoff frequencies for RBT. This means that signals with a periodicity larger than 20 years and lower than 0.9 years were excluded from the calculations and display. The reasoning for this was to make the RBT data comparable to the yearly data of the GDP change. The low-frequency cutoff was canceling long-term trends, as the GDP change was only related to the previous year. The high-frequency cutoff was used to dampen fast, alternating RBT signals in comparison to the slow sampled GDP data.

2.5 Rescaled adjusted partial sums

Rescaled adjusted partial sums (RAPSs) were used to visualize trends in time series which may not be clearly visible in the unprocessed data set. Equation (4) shows the calculation of the RAPS index X, using time series Y.

(4) X k = i = 1 i = k Y i - Y σ Y .

Y is the average over the total time series, σ is the standard deviation of the whole time series and Yi is the ith data point in Y.

A change in the slope of the RAPS index only indicates a change in the slope of the original time series. A negative RAPS slope does not indicate a negative slope in the original time series. Garbrecht and Fernandez (1994) and Basarin et al. (2016) used this method to investigate trends in hydrological time series.

2.6 Catchment area

The catchment area was calculated using the HydroSHEDS database (Lehner et al.2008). The 1125 by 1125 gridded data set provides the information, at each grid point, according to which cell the water of a grid cell is drained. By selecting a starting location, e.g., Koblenz at 50.350 N and 7.602 E, it was possible to iteratively identify all grid points draining into this location. These grid points represent the catchment area of this location (in the example from Fig. 2 at Koblenz). By counting the iteration steps, the distance a water drop travels to reach the monitoring station of Koblenz was determined. This was done for each of the four stations. Additionally, the accumulation number (ACC) was obtained from the data set. It defines how many cells in total were draining into a particular cell, and it is a measure for the size of a river. Finally, a grid, which defines the catchment area, the ACC and the hydrological distance, was established, spanning the whole catchment area. Figure 2 shows the catchment area, the hydrological distance and the calculated flow time to the Koblenz monitoring station.

Figure 2Catchment area of the Koblenz monitoring station as an example. The colors show the hydrological distance between the monitoring station and each grid point of the catchment area. The second y axis shows the time it takes in days, in our setup, to flow from a certain grid point to the monitoring station, based on the hydrological distance. The flow speed is 0.4 m s−1 and was defined to be constant in space and time. All monitoring stations are marked by a cross. The other markers show the location of the NPPs.

Figure 3The catchment area of the Koblenz monitoring station. (a) Number of grid points (accumulation number – ACC) which are upstream of each specific grid point. (b) ACC w, distance and ACC-weighed grid cells.

2.7 Multiple regression

A multiple linear regression was used to separate the anthropogenic heat input a1 and meteorological a2 and hydrological a3 contributions to the river water temperature. Tw was regressed with Tc and river discharge Q. Their regression coefficients a2 (Tc slope) and a3 (Q slope) represent the magnitude of the respective influences. The offset a1 (RBT) combines all other influences which were not related to a change in Tc or Q. We hypothesized that the RBT is directly linked to heat input by power plants, which are NPPs in this study, and other industrial facilities.

Instead of taking Ta directly at the monitoring station, we improved Eq. (1) by a time-dependent and weighed average of Ta (x,y,t) over the whole catchment; see Eq. (5). There were spatial coordinates (x,y) in the catchment area, and the subscript 0 marks the location of the monitoring station.

(5) T w x 0 , y 0 , t 0 = a 1 + a 2 T c x 0 , y 0 , t 0 + Δ t x , y + a 3 Q x 0 , y 0 , t 0 .

The new representative catchment temperature was called Tcx0,y0,t0. It was a weighed average of the whole catchment area (x,y), which is defined by the measurement point (x0,y0). The difference between the measurement time t0 and the reading of Ta is called time lag Δt(x,y) and depends on the hydrological distance between the measurement point and the reading. Δt(x,y) is negative and points to a moment in time before the measurement.

2.7.1 Time lag

A change in Tw is slower than a change in Ta. The time lag Δt describes this lag and is commonly used in water temperature models.

A reason for the occurrence of Δt is that the water mass mixing capability, heat capacity and surface area cause strong thermal inertia. Changing Tw through new meteorological conditions and heat fluxes takes time. Therefore, linear and exponential models include either a fixed Δt for Ta (Eq. 6) or an average of Ta, including a time span from before (Eq. (7); Stefan and Preudhomme1993; Webb and Nobilis1995, 1997; Haag and Luce2008; Benyaha et al.2008).


A second reason for the mismatch is advection. Rivers, in this case the Rhine, exhibit current velocities which enable their water to cover significant distances on timescales larger than days. Therefore, it is necessary to take the change in Ta, in space and time, during advection into account. This is especially important for daily averaged Tw (Erickson and Stefan2000). Pohle et al. (2019) average 8 d of hydroclimatic variables over the whole catchment area; see Eq. (8). However, this approach does not include the characteristics of flow path and flow speed.

(8) T c x 0 , y 0 , t 0 = x = 0 , y = 0 , t = t 0 x = n , y = m , t = t 0 - 8 T a x , y , t .

We combined the general idea of a time lag and averaging Ta over the whole catchment area (see Eqs. 6, 7 and 8), but in this study, each grid point was linked to a specific time lag Δt(x,y) which is dependent on a fixed flow speed v and the hydrological distance s(x,y) to the measurement point (see Fig. 2). The distance was obtained from the discharge map (Sect. 2.6) and calculated with v, as described by Eq. (9). The new Δt(x,y) represents the mismatch by advection but not specifically the mismatch through thermal inertia. The thermal inertia would be independent of s(x,y) and a constant added to Δt. However, we are of the opinion that a sufficient part of the thermal inertia time lag was included in our representation of Δt(x,y).

(9) Δ t x , y = - s x , y v .

2.7.2 Weighing coefficients

Tobler (1970) proposed that close spatial and temporal conditions tend to be more highly correlated than those further away. This led to the introduction of the weighing factor w. A linear decreasing weighing factor from 1 to 0 was used; 1 marks the grid point closest (smallest Δt) to the monitoring station and 0 the point farthest away (largest Δt). As the sizes of the catchment areas were different for the four monitoring stations, four weighing coefficient tables were calculated. As an example, Table 3 shows the weighing coefficient for Koblenz.

Table 3Weighing factors for the distance and the resulting Δt for the monitoring station of Koblenz. Δt is calculated from distance and flow speed; see Eq. (9). The weighing coefficient is linearly correlated to Δt.

Download Print Version | Download XLSX

A catchment-wide hydrological flow model, estimating the flow speed at every grid point for every hydrological scenario, was not used. It was not available yet for every grid point of the catchments, and the focus of this study was to create a simple setup that is also transferable to other river catchments. Therefore, using a constant flow speed of 0.4 m s−1 was decided on. This flow speed was determined by calculating the root mean square error (RMSE; whole data set) of the ACC t model, with a step-wise reduction in the flow speed from 1.4 to 0.3 m s−1. The lowest RMSE was obtained at Koblenz at 0.4 m s−1. The RMSE and Nash–Sutcliffe coefficient (NSC) coefficients at all flow speeds and all stations are shown in the Supplement. For Basel and Worms, slower flow speeds lowered the RMSE further. We did not include this as it would create unreasonable low flow speeds. The reason for the low flow speeds, with lowest RMSE, might be the thermal inertia. As thermal inertia is not explicitly represented in Δt (Eq. 9), a smaller flow speed could compensate for that, especially in smaller catchment areas.

The weighing coefficient w is combined with ACC. ACC is used as a second coefficient which overweighs grid points with large accumulation and, therefore, large water masses. This ensures a balance between the large number of low ACC grid points, which carry less water, and the influence of Ta on large water masses. Figure 3b shows ACC w over the whole catchment area of Koblenz.

Figure 4ACC bins (x axis) vs. the relative contribution (y axis). The grid points are binned by their ACC value. The red bars show the relative contribution (largest contribution normalized to one) by the number of grid points in this bin only. The white bars show the distribution using the number of grid points in this bin and weighing ACC w.



Combining Δt, ACC w weighing and the gridded temperature reanalysis data of Sect. 2.2, we proposed a new 3D x,y,t averaging of Ta, which is shown in Eq. (10).

(10) T c x 0 , y 0 , t 0 = 1 n m 1 w Δ t x , y ACC x , y x = 1 , y = 1 x = n , y = m w Δ t x , y ACC x , y T a x , y , t 0 + Δ t x , y .

Tcx0,y0,t0 was calculated by weighed (ACC w) averaging Tax,y,t+Δtx,y over all grid points of the catchment area (x=1,ny=1,m), which was set by the measurement point (x0,y0). The time lag Δt was an estimate for the time it takes a water droplet from a specific grid point (x,y) in the catchment area to reach the measurement location.

Based on Eq. (10), the daily Tc was calculated for each monitoring station. This temperature represents the meteorological influence all water droplets have experienced on their way to the monitoring station and is subsequently used in the multiple linear regression.

2.7.4Tc calculation methods

Additionally, we used these four calculations methods, namely (1) wt, (2) avg t, (3) avg and (4) point, to compare the results of the linear regression to the calculation proposed in Eq. (10).

(1) In the following, only w weight (Eq. 11) and time lag are shown:

(11) T c x 0 , y 0 , t 0 = 1 n m 1 w Δ t x , y x = 1 , y = 1 x = n , y = m w Δ t x , y T a x , y , t 0 + Δ t x , y .

(2) Here, no weight and only time lag (Eq. 12) are shown, as follows:

(12) T c x 0 , y 0 , t 0 = 1 n m x = 1 , y = 1 x = n , y = m T a x , y , t 0 + Δ t x , y .

(3) A mean Ta(x0,y0,t0) over the whole catchment area at the time t0 of the measurement was calculated; see Eq. (13). Δt was not used in the following:

(13) w x , y = 1 T c x 0 , y 0 , t 0 = 1 n m x = 1 , y = 1 x = n , y = m T a x , y , t 0 .

(4) Ta(x0,y0,t0) at the location x0,y0 and time t0 of the measurement (see Eq. 14) is shown, as follows:

(14) T c x 0 , y 0 , t 0 = T a x 0 , y 0 , t 0 .
3 Results

3.1 Water temperature time series

To investigate the long-term changes, we fitted a time-dependent linear function to the time series of Tw and Ta (catchment average) of all four monitoring stations (Basel, Worms, Koblenz and Cologne). The same was also done when all four monitoring stations had an overlapping data set (1985–2018); see Table 1. Figure 5a presents the yearly averaged Tw and the linear fits for the two time periods. The average Ta of the catchment area is also shown. In Fig. 5b, the RAPS index of Ta and Tw is shown. The fit coefficients and the rate of warming per year are displayed in Table 4. The calculated Ta increased in the catchment area of all monitoring stations, and the respective slopes are shown in columns four and five of Table 4.

Figure 5(a) Yearly averages of water temperatures at the four monitoring stations (black line). The red dashed line is a fit to the available data set. The red dotted line is a fit to the overlapping time period (1985–2018). The blue line is the yearly average air temperature of the catchment area. (b) RAPS Tw (black) and Ta (blue) indexes. The triangle markers divide the RAPS index into sections, based on a slope change in the RAPS index. Each section also represent a trend change in the original Ta and Tw time series.


Figure 5 and Table 4 show that the change in Tw was found to be heterogeneous along the Rhine. The slope at Basel is approx. 6 times higher (0.049 C yr−1) than the one at Cologne (0.0084 C yr−1), when comparing only the overlapping data set. However, during the same period, Ta (0.05 C yr−1 at Basel; 0.05 C yr−1 at Cologne) displays similar behavior at these two stations, which is an indication of a similar meteorological influence. The Tw warming rate from 1985 to 2018 for Worms and Koblenz is in between those from Cologne and Basel. These two stations show similar Ta warming rates compared to Basel and Cologne. Generally, the Ta warming rates have a difference of less than 5 %. Arora et al. (2016) showed a mean Tw warming rate in northern and northeastern German rivers of 0.03 C yr−1 (1985–2000) and 0.09 C yr−1 (2000–2010). Regarding our time period (1985–2010), these values are plausible. Basarin et al. (2016) found a maximum increase in Tw at the Danube at Bogojevo (1950–2012) of 0.05 C yr−1, which matches the maximum increase at Basel. Ta increased by 0.02 C yr−1 between 1985 and 2010 in the study by Arora et al. (2016). We found a steeper slope at all stations. The reason could be the hiatus of global warming (Hartmann et al.2014), which is a flattening of the Ta increase between 1998 and 2012. This period is fully included in the Arora et al. (2016) study and our data set, but we investigated further until 2018, when the warming of Ta already increased again (Hu and Fedorov2017). Michel et al. (2020) investigated Tw at 52 river gauges in Switzerland, representing most of the Rhine catchment area at Basel. The authors reported an average Tw increase at the 52 stations of 0.037 C yr−1 (1998–2018) and 0.033 C yr−1 (1979–2018). Ta increased by 0.039 C yr−1 (1998–2018) and 0.046 C yr−1 (1979–2018). Comparing this to our results at Basel, the Ta warming rates are similar. The difference might originate from the use of meteorological stations nearby river gauges only (Michel et al.2020) instead of a reanalysis product. The difference in Tw warming (approx. 0.021 C yr−1) could be interpreted as showing that most warming might occur in the broader vicinity before the Basel monitoring station.

The R2 values make the differences between the four monitoring stations visible. Basel exhibits the largest R2 values, and these are consistently high for Ta and Tw. This is in contrast to the station at Cologne, where R2 of Tw was low and statistically not significant. The slope of Ta at Cologne is lower than at the other stations but still statistically significant. The Pearson correlation coefficients between Ta and Tw were the lowest at Cologne and the largest in Basel. For Ta, the RAPS index of all the monitoring stations showed four concurrent sections (start–1987, 1987–2000, 2000–2014 and 2014–end). Their borders are marked by the blue triangles in Fig. 5. The section between 2000 and 2014 could be a consequence of the hiatus of global warming between 1998 and 2012 (Hartmann et al.2014). Each section represents slope changes in the RAPS index and indicates trend changes in the original time series. The Tw RAPS index for Basel displayed the same pattern in the sections as the Ta index. All other stations showed other RAPS Tw to RAPS Ta patterns. This means that the Ta and Tw trends of the original time series were different at these stations. Ta cannot fully describe the trends in Tw.

We hypothesized that different meteorological conditions were not the reason for this difference. Meteorological differences should also have been visible in the Ta warming rate of the four stations, which was not the case. The RAPS analysis for Ta and Tw only coincided within the Basel data set.

Table 4Slope of linear fits and Pearson correlation coefficients to the daily temperature data at the four monitoring stations. The data set used is described in the column header. Next to the slope values are the R2 values, which are statistically significant if R2>0.19.

Download Print Version | Download XLSX

3.2 Regression

We fitted the multiple regression model (Eq. 5) using Tc and Q to Tw of each monitoring station for the available data set. Afterwards, we recalculated Tw,modelled using the regression coefficients a1, a2 and a3. From the comparison between the Tw,modelled and measured Tw, the RMSE and the NSC for each monitoring station was derived (see Table 5). To support the introduction of weighing coefficients ACC w and Δt, we compared five different calculations of Tc from Sect. 2.

Table 5 shows the RMSE and NSC values for all correlations. The lowest (RMSE) and highest (NSC) values are displayed in bold in Table 5. The lowest RMSE was found to be 1.02 C for ACC w+Δt (first row) at the Koblenz station. At this location, the largest NSC, of 0.97, also appeared. The flow speed was optimized for lowest RMSE at the Koblenz station (see Sect. 2.7.2). It was evident that the three methods, including a Δt, have a lower RMSE (below 2.01 C; lowest 1.02 C) than the two methods without a Δt (above 2.37 C; largest 2.97 C). The same trend held true for NSC where the Δt methods were above 0.90 and the other two were below 0.86. We think that the use of a catchment-wide Δt improved the quality of the multiple regression analysis and delivered a significant improvement in the TaTw-based modeling. Interestingly, combining ACC and the weighing factor w provided the best estimation for all stations, except for Basel. The content of Fig. 4 could explain this result. Without ACC weighing, small water masses (small ACC) may be overrepresented in the contribution to Tc. Large ACC grid points represent large water masses (rivers and lakes), and their influence on Ta may be otherwise underestimated. At Basel the fraction of lowACC grid points was relatively small compared to the other stations, as Basel is closest to the water sources and has the smallest catchment area. Therefore, the ACC weighing might have provided weaker results.

As ACC w+Δt provided the smallest RMSE, this calculation method was used for all further calculations of Tc.

In the Supplement, we provide a calculation of the regression coefficients for the year 2001 only. These coefficients were then taken as a basis for calculating Tw for each year from 2000 to 2018. The RMSE and NSC data were consistent in magnitude with the long-term regressions of this section. The RMSE at Koblenz ranged from 0.75 to 1.22 C. A lower RMSE was caused by the shorter regression period. This supports the stability and validity of the regression model.

Table 5RSME (in degrees Celsius) and NSC for all Tc calculation methods. Different Tc calculation methods and the regressions are applied over the total data set. The RMSE and the NSC are calculated between Tw and Tw,modelled. The first column contains the calculation method. The best results for each monitoring station and each calculation method are indicated in bold.

Download Print Version | Download XLSX

3.3 Rhine base temperature

The RBT was used to explain differences in the Tw warming rates of Table 4. We regressed a 2-year segment of the Tw time series and set a step size of 1 month in order to create a RBT time series over the full data set. The regression of a 2-year segment should also compensate for extreme events occurring during 1 year. These could be extremely low discharge or extreme water temperatures in which industrial and power production had to react. As the absolute RBT cannot be meaningfully interpreted, only the changes in RBT over time are shown in Fig. 6. We subtracted the last data point of each time series from the rest of the data and showed the change in RBT, a 4-year running mean and Δ RBT (Eq. 3) vs. time. The HI by NPPs is shown as a dotted blue line relating to the y axis on the right (Fig. 6).

Figure 6Rhine base temperature (RBT) from four monitoring stations (black solid line). The red dashed line is the RBT 4-year running mean. The magenta line with the “+” markers shows the ΔRBT relative to the last year. The blue dotted line is the upstream heat input (HI) by NPPs; see Sect. 2.3.


3.3.1 Long-term trends

In this study, long-term trends were visible on timescales of decades. The HI by NPPs, the 4-year running mean RBT and ΔRBT followed a similar trend in this analysis (see Fig. 6). After the maximum heat discharge from NPPs between 1996 and 1998, the HI and the RBT of Worms, Koblenz and Cologne declined. The RBT started its decline 1–2 years before 1995, which might have been triggered by the recession in 1993 and a sharp drop in the German trade balance. At Basel, the RBT and the HI remained comparably constant. Additionally, we calculated ΔTw based on the change in HI, using Eq. (3), at every station and compared it to the ΔRBT from the regression model (see Table 6). The period for each monitoring station starts at the maximum HI by NPPs for the respective station and ends in the year 2017.

Table 6Change in RBT (column three) in the time period given in column two. The start of the period indicates the maximum HI of NPPs at the respective monitoring station. The calculated ΔTw (column four) and the change in HI by nuclear power plants (column five) are also provided. The calculations were done using Eq. (3). Note: GW – gigawatt.

Download Print Version | Download XLSX

At Basel, both simulated and calculated RBT changes were negligible due to the lack of change in HI. At all other stations, the change in HI was reflected in the change in RBT. The maximum difference between simulation and calculation was found to be 0.34 C. Before 1995, Worms, Koblenz and Cologne showed an approximately 1 C offset between ΔRBT and ΔTw (see Fig. 6). This occurred during a time when the NPPs HI remained relatively stable, but the GDP increased by 30 % between 1985 and 1995 (Worldbank2020). The change in nuclear power production over a time period of 30 years or more can explain the changes in and the heterogenous warming rates of Tw along the Rhine. NPPs may also impact Tw at much shorter timescales but do not change their power output accordingly.

3.3.2 Short-term trend

Short-term changes (<5 years) in RBT (Fig. 6) are not likely to be influenced by the overall HI from NPPs, as these adopt production at longer timescales. More important are local industrial conditions, which could also include fossil fuel power plants. However, not all influences on the coefficient a1 and, subsequently, on the RBT originate from industrial production. Various potential influences are unknown and not within the scope of this publication.

For Basel, it was not possible to satisfyingly explain the short-term variations. The Rhine and its tributaries upstream are flowing through subalpine lakes and, in relation to the downstream part, are not strongly industrialized. Lakes have a complicated heat budget (Råman Vinnå et al.2018), which was not focused on in this analysis.

For all other stations, we hypothesized that local production facilities and their HI into the Rhine are responsible for the short-term changes, by comparing the RBT time series to economic data. Figure 7 shows the comparison of RBT (black line; 1 year running mean) to the changes in the GDP (blue line). A discontinuity in the GDP in 1991 is visible due to the German reunification, which is when the calculation method of the GDP changed. Therefore, the data were plotted as separate lines.

Figure 7The change in RBT (black solid line) at three monitoring stations (Cologne, Koblenz and Worms). The blue dashed line is the GDP change in the adjacent federal states. To explain trends during two time periods, the red dashed line, which is the turnover of BASF, and the red dotted line, which is the production rate of the oil refineries, were added. The triangles mark the years 2000 (burst of the dot-com bubble) and 2008 (mortgage crisis).


For Worms (Fig. 7, bottom panel), we added the change in the turnover of BASF (red dashed line; BASF AG1989). BASF is a major chemical company, and one of its largest production facilities, with an estimated HI of 500 MW to 1 GW, is located 12 km upstream (Rhine kilometer 431) from the Worms station. It was investigated whether the production and HI changes in this factory were also visible. In 1985, although the change in GDP did not indicate a large RBT change, a RBT decrease was visible. This was indicated by a turnover decrease in 1985 and 1986. After the German reunification in 1990, a negative GDP change (recession) was evident. This was followed by a BASF turnover decline and a decrease in RBT. After that, the RBT followed the zigzag movements of the GDP, and so did the BASF turnover (only shown until 2000). In particular, the economic events, such as the burst of the dot-com bubble (early 2000s) and the mortgage crisis (2008), were visible in the RBT and in the GDP when a decrease in both parameters followed. The two events are marked by triangles in Fig. 7.

Before 1990, the RBT at Koblenz did not follow the GDP trend, and it showed a rather anticyclic behavior which can not be explained yet. After 1991, the RBT followed the general trend of the GDP but did not seem to be strongly influenced by the short recession after the German reunification. Again, economic events such as the burst of the dot-com bubble (early 2000s) and the mortgage crisis (2008) displayed an influence on the RBT.

The RBT at Cologne did not seem to be strongly influenced by the recession connected to the German reunification, but after 1999, the RBT follows the zigzag trends of the GDP.

For all monitoring stations, a red dashed line was added between 1995 and 1999. This dashed line indicates the production rate of German oil refineries (MWV2003). From 1995 to 1999, German refineries ran at full capacity (100 %); usually the capacity levels did not exceed 90 %. The increase in production was clearly visible in the RBT at Cologne, where a large oil refinery is located 19 km upstream at Rhine kilometer 671 (Rhineland refinery). The RBT at Worms and Koblenz could be influenced by the output of a refinery next to Karlsruhe at Rhine kilometer 367 (Mineraloelraffinerie Oberrhein).

Figure 8Detrended and filtered RBT signal (black solid) and the GDP change (blue dashed) at Cologne, Koblenz and Worms.


3.3.3 Correlation

We correlated the GDP change and the filtered RBT signal. It was noticeable that a 480 d shift to the past was needed to obtain matching trends. This means that a change in RBT or anthropogenic HI appeared about 480 d earlier than in the GDP calculation. The shift could be caused by the following two reasons: (1) using the GDP difference of 2 consecutive years has a significance at a unspecific point of time within these 2 years. (2) The GDP is lagging behind the real economic situation which is, in this case, the industrial production. Yamarone (2012) claimed that GDP was a coincident economic indicator similar to industrial production. However, Yamarone (2012) used quarterly GDP calculations, and in this study annual data were used. The quarterly data set may react faster to changes. A second thought was that Yamarone (2012) compared the calculations of industrial production, which is an economic index, to GDP (another economic index). In this study, real-time data from industrial HI into the river was processed. This shift has not been done for Fig. 6) because a shift of 1.5 years on a 40-year timescale is negligible.

Table 7Spearman's rank correlation coefficients between RBT and GDP change for ACC w+Δt. The last column shows the significance.

Download Print Version | Download XLSX

Table 7 shows the Spearman's rank correlation coefficients of Worms, Koblenz and Cologne for the ACC w+Δt calculation method, which resulted in the lowest RMSE at Koblenz.

All correlations were found to be positive and statistically significant (p<0.05). The correlation at Koblenz was highest. Figure 8 shows the filtered RBT signal vs. the GDP change at the three monitoring stations. The RBT time series was detrended and filtered. Most of the time, the variations in the RBT (filtered and shifted) were coincident with the GDP change. The RBT peak between 1995 and 1998 was not very well represented by the GDP change, which has already been discussed earlier in context of Fig. 7.

4 Conclusions

We introduced a new catchment-wide air temperature Tc, which decreased the RMSE (Table 5) in a TcTw regression. Tc is a weighed (ACC w) average of all Ta across the catchment area, including the use of Δt for each grid point according to the hydrological distance and flow speed. In the approach, this time lag was used as an indicator for the point in time at which a water droplet was at a certain grid cell in the catchment area. As a result, one can obtain a better estimate of which Ta a water droplet experienced on its way to a certain point (in this study it is a monitoring station), and it delivered better linear TcTw estimates. This improvement in the TcTw relationship may support the analysis of processes in the heat budget of rivers. Usually, Ta data is readily available and can easily be combined with Q data for multiple linear regression analysis. Still, a sufficiently long (decadal) time series of Tw was required. Nevertheless, a linear relationship was found to be simpler than a full, physically based model which requires all meteorological fluxes as input quantities.

In the proof of concept, we focused on the Rhine catchment area, but in principle, the model can be applied to any river system around the globe if the respective long-term data are available. However, catchment area data and reanalysis Ta data are often globally available. Morrill et al. (2005) showed a linear TaTw relationship for 43 rivers with various catchment areas in the subtropics. This potentially indicates that the proposed model and procedure can be applied elsewhere. However, this still has to be verified. Future calculations may be coupled with catchment-wide hydrological models to improve the accuracy of the time lag. The time lag used in this study was based on trial and error in search of the lowest RMSE. A detailed catchment-wide hydrological flow model would be especially beneficial for setting an upper limit for the time lag and constraining its validity. It would also be interesting to estimate the importance of the advection time lag vs. the thermal inertia time lag.

With Tc, we regressed four Tw time series (Basel, Worms, Koblenz and Cologne) along the Rhine. The offset in the this regression a1 was called RBT, and its change over time was found to be an indicator for anthropogenic HI. The RBT positively correlated to long-term economic changes, such as the decrease in nuclear power production, and to short-term economic events. We showed that changes in production rates (oil refineries or chemical industry) and a change in GDP may influence the RBT and, therefore, the Rhine's water temperature. Additionally, the Spearman's rank correlation coefficient between RBT and GDP is positive and significant, providing another indication for the relation. This case study might provide a tool for a better understanding of the long-term consequences of industrial water use, and it might be used as a verification tool for reported HI. Germany has a rigorous reporting system on cooling water use. However, other countries could check if industrial HI is in accordance with legislative guidelines, without depending on official reports. Whether the ongoing COVID-19 (2020) pandemic and its impact on the economy is also visible, using the offered procedures, may be explored after the crisis.

Hardenbicker et al. (2016) estimated, using a physically based model (QSIM), that between the reference period of 1961–1990 and the near future 2021–2050, the mean annual Tw of the Rhine could increase by 0.6–1.4 C. This trend is plausible, according to the historical data analyzed, if the Ta increase remains constant. However, they used a constant anthropogenic HI by, for example, power plants and production industries, and different warming rates along the Rhine can result from changes in anthropogenic HI. Next to the global air temperature increase, the industrial use of river water will have major impacts on the Rhine water temperature in future.

The difference in the Tw warming rate between Basel and the other monitoring stations in the time series data can be explained by the change in nuclear power production and the influence of general industrial production. For the Rhine, a decreasing (except for Basel) RBT, which indicates a decreasing HI, was found. Other river catchment areas with growing energy-intensive industries might be impacted by much larger warming rates than those caused by the general increase in Ta, which means they will experiences all the consequences on the physical, chemical and biological processes.

Data availability

The data used in this study are available from Zenodo at (Zavarsky and Düster2020). These data include the daily averaged water temperatures and discharge at the Rhine monitoring stations, the yearly heat input from nuclear power plants, and the time series of four different calculation methods for the catchment wide air temperature as described in this paper.


The supplement related to this article is available online at:

Author contributions

AZ and LD did the study concept and design. AZ collected the data, developed the model, analyzed the data and wrote the paper in consultation with LD. LD provided critical feedback and helped shape the research.

Competing interests

The authors declare that they have no conflict of interest.


We acknowledge all the data providers for their help and the German Federal Ministry of Transport and Digital Infrastructure for founding the Federal Institute of Hydrology and, therefore, making this work possible.

Review statement

This paper was edited by Alberto Guadagnini and reviewed by Ina Pohle and one anonymous referee.


Arbeitskreis VGdL: Bruttoinlandsprodukt, Bruttowertschöpfung in den Ländern der Bundesrepublik Deutschland 1991 bis 2018 (Reihe 1 Band 1), Berechnungsstand: August 2018/Februar 2019, available at:, last access: 1 March 2019a. a

Arbeitskreis VGdL: Arbeitstabellen: Rückrechnungsergebnisse für das frühere Bundesgebiet nach Bundesländern 1970 bis 1991 Berechnungsstand: August 2006, available at: (last access: 1 March 2019), 2007b. a

Arora, R., Tockner, K., and Venohr, M.: Changing river temperatures in northern Germany: trends and drivers of change, Hydrol. Process., 30, 3084–3096,, 2016. a, b, c

BAFU: Water temperature and discharge Basel, Bundesamt fuer Umwelt BAFU Abteilung Hydrologie, available at:, last access: 1 March 2019. a, b

Basarin, B., Lukić, T., Pavić, D., and Wilby, R. L.: Trends and multi-annual variability of water temperatures in the river Danube, Serbia, Hydrol. Process., 30, 3315–3329,, 2016. a, b

BASF AG: Geschaeftsbericht, BASF AG Oeffentlichkeitsarbeit und Marktkommunikation, Ludwigshafen, Deutschland, available at: (last access: 1 March 2019), 1989. a

Benyaha, L., St Hilaire, A., Ouarda, T., Bobee, B., and Dumas, J.: Comparison of non-parametric and parametric water temperature models on the Nivelle River France, Hydrolog. Sci. J., 53, 640–655,, 2008. a

BfG: Water temperature and discharge Koblenz, Bundesanstalt fuer Gewaesserkunde, Koblenz, available at:, last access: 1 March 2019. a, b

Cai, H., Piccolroaz, S., Huang, J., Liu, Z., Liu, F., and Toffolon, M.: Quantifying the impact of the Three Gorges Dam on the thermal dynamics of the Yangtze River, Environ. Res. Lett., 13, 054016,, 2018. a

Caissie, D.: The thermal regime of rivers: a review, Freshwater Biol., 51, 1389–1406,, 2006. a

Delpla, I., Jung, A.-V., Baures, E., Clement, M., and Thomas, O.: Impacts of climate change on surface water quality in relation to drinking water production, Environ. Int., 35, 1225–1233,, 2009. a

Durance, I. and Ormerod, S. J.: Trends in water quality and discharge confound long-term warming effects on river macroinvertebrates, Freshwater Biol., 54, 388–405,, 2009. a

Erickson, T. R. and Stefan, H. G.: Linear Air/Water Temperature Correlations for Streams during Open Water Periods, J. Hydrol. Eng., 5, 317–321,, 2000. a

Förster, H. and Lilliestam, J.: Modeling thermoelectric power generation in view of climate change, Reg. Environ. Change, 10, 327–338,, 2010. a, b

Garbrecht, J. and Fernandez, G. P.: Visualization of trends and fluctuations in climatic records, J. Am. Water Resour. Assoc., 30, 297–306,, 1994. a

Gaudard, A., Weber, C., Alexander, T. J., Hunziker, S., and Schmid, M.: Impacts of using lakes and rivers for extraction and disposal of heat, WIREs Water, 5, e1295,, 2018. a, b

Haag, I. and Luce, A.: The integrated water balance and water temperature model LARSIM-WT, Hydrol. Process., 22, 1046–1056,, 2008. a

Hardenbicker, P., Viergutz, C., Becker, A., Kirchesch, V., Nilson, E., and Fischer, H.: Water temperature increases in the river Rhine in response to climate change, Reg. Environ. Change, 17, 299–308,, 2016. a

Hari, R. E., Livingstone, D. M., Siber, R., Burkhardt-Holm, P., and Guettinger, H.: Consequences of climatic change for water temperature and brown trout populations in Alpine rivers and streams, Glob. Change Biol., 12, 10–26,, 2006. a

Hartmann, D. L., Tank, A. M. G. K., Rusticucci, M., Alexander, L. V., Brönnimann, S., Charabi, Y. A. R., Dentener, F. J., Dlugokencky, E. J., Easterling, D. R., Kaplan, A., Soden, B. J., Thorne, P. W., Wild, M., and Zhai, P.: Observations: Atmosphere and Surface, in: Climate Change 2013 – The Physical Science Basis, edited by: IPCC, 159–254, Cambridge University Press, Cambridge, United Kingdom,, 2014. a, b

Hassoun, M.: Fundamentals of artificial neural networks, MIT Press, Cambridge, Massachusetts, 1995. a

Hoef, J. M. V. and Peterson, E. E.: A Moving Average Approach for Spatial Statistical Models of Stream Networks, J. Am. Stat. Assoc., 105, 6–18,, 2010. a

Hoef, J. M. V., Peterson, E., and Theobald, D.: Spatial statistical models that use flow and stream distance, Environ. Ecol. Stat., 13, 449–464,, 2006. a

Hu, S. and Fedorov, A. V.: The extreme El Niño of 2015–2016 and the end of global warming hiatus, Geo. Res. Lett., 44, 3816–3824,, 2017. a

IAEA: Power Reactor Information System (PRIS), Web, available at:, last access: 1 March 2019. a

IKSR: Vergleich der Waermeeinleitungen 1989 und 2004 entlang des Rheins, IKSR-Bericht, available at: (last access: 1 March 2019), 2006. a

Isaak, D. J., Luce, C. H., Rieman, B. E., Nagel, D. E., Peterson, E. E., Horan, D. L., Parkes, S., and Chandler, G. L.: Effects of climate change and wildfire on stream temperatures and salmonid thermal habitat in a mountain river network, Ecol. Appl., 20, 1350–1371,, 2010. a

Isaak, D. J., Wollrab, S., Horan, D., and Chandler, G.: Climate change effects on stream and river temperatures across the northwest U.S. from 1980–2009 and implications for salmonid fishes, Climatic Change, 113, 499–524,, 2012. a

Isaak, D. J., Peterson, E. E., Hoef, J. M. V., Wenger, S. J., Falke, J. A., Torgersen, C. E., Sowder, C., Steel, E. A., Fortin, M.-J., Jordan, C. E., Ruesch, A. S., Som, N., and Monestiez, P.: Applications of spatial statistical network models to stream data, Wiley Interdisciplinary Reviews: Water, 1, 277–294,, 2014. a

Jackson, F., Hannah, D. M., Fryer, R., Millar, C., and Malcolm, I.: Development of spatial regression models for predicting summer river temperatures from landscape characteristics: Implications for land and fisheries management, Hydrol. Process., 31, 1225–1238,, 2017a. a

Jackson, F. L., Fryer, R. J., Hannah, D. M., and Malcolm, I. A.: Can spatial statistical river temperature models be transferred between catchments?, Hydrol. Earth Syst. Sci., 21, 4727–4745,, 2017b. a

Jackson, F. L., Fryer, R. J., Hannah, D. M., Millar, C. P., and Malcolm, I. A.: A spatio-temporal statistical model of maximum daily river temperatures to inform the management of Scotland's Atlantic salmon rivers under climate change, Sci. Total Environ., 612, 1543–1558,, 2018. a

Koch, H. and Grünewald, U.: Regression models for daily stream temperature simulation: case studies for the river Elbe, Germany, Hydrol. Process., 24, 3826–3836,, 2010. a

Lange, J.: Waermelast Rhein, Bund fuer Umwelt und Naturschutz Deutschland, available at: (last access: 1 March 2019), 2009. a, b, c, d, e

Lehner, B., Verdin, K., and Jarvis, A.: New Global Hydrography Derived From Spaceborne Elevation Data, Eos T. AGU, 89, 93–94,, 2008. a

LfU: Water temperature and discharge Worms, Landesamt fuer Umwelt Rheinland-Pfalz, available at:, last access: 1 March 2019. a, b

Markovic, D., Scharfenberger, U., Schmutz, S., Pletterbauer, F., and Wolter, C.: Variability and alterations of water temperatures across the Elbe and Danube River Basins, Climatic Change, 119, 375–389,, 2013. a, b

Michel, A., Brauchli, T., Lehning, M., Schaefli, B., and Huwald, H.: Stream temperature and discharge evolution in Switzerland over the last 50 years: annual and seasonal behaviour, Hydrol. Earth Syst. Sci., 24, 115–142,, 2020. a, b

Mohseni, O., Stefan, H. G., and Erickson, T. R.: A nonlinear regression model for weekly stream temperatures, Water Resour. Res., 34, 2685–2692,, 1998. a

Morrill, J. C., Bales, R. C., and Conklin, M. H.: Estimating Stream Temperature from Air Temperature: Implications for Future Water Quality, J. Environ. Eng., 131, 139–146,, 2005. a

MWV: Mineraloel und Raffinerien, Mineraloelwirtschaftsverband e.V., available at: (last access: 1 March 2019), 2003. a

Peterson, E. E. and Hoef, J. M. V.: A mixed-model moving-average approach to geostatistical modeling in stream networks, Ecology, 91, 644–651,, 2010. a

Piccolroaz, S., Calamita, E., Majone, B., Gallice, A., Siviglia, A., and Toffolon, M.: Prediction of river water temperature: a comparison between a new family of hybrid models and statistical approaches, Hydrol. Process., 30, 3901–3917,, 2016. a

Pohle, I., Helliwell, R., Aube, C., Gibbs, S., Spencer, M., and Spezia, L.: Citizen science evidence from the past century shows that Scottish rivers are warming, Sci. Total Environ., 659, 53–65,, 2019.  a

Råman Vinnå, L., Wüest, A., Zappa, M., Fink, G., and Bouffard, D.: Tributaries affect the thermal response of lakes to climate change, Hydrol. Earth Syst. Sci., 22, 31–51,, 2018. a, b

Sinokrot, B. A. and Stefan, H. G.: Stream temperature dynamics: Measurements and modeling, Water Resour. Res., 29, 2299–2312,, 1993. a, b

Stefan, H. G. and Preudhomme, E. B.: Stream temperature estimation from air temperature, J. Am. Water Resour. Assoc., 29, 27–45,, 1993. a, b

Tobler, W. R.: A Computer Movie Simulating Urban Growth in the Detroit Region, Econ. Geogr., 46, 234,, 1970. a

Toffolon, M. and Piccolroaz, S.: A hybrid model for river water temperature as a function of air temperature and discharge, Env. Res. Lett., 10, 114011,, 2015. a

Webb, B. W. and Nobilis, F.: Long term water temperature trends in Austrian rivers, Hydrolog. Sci. J., 40, 83–96,, 1995. a

Webb, B. W. and Nobilis, F.: LONG‐term perspective on the nature of the air–water temperature relationship: a case study, Hydrol. Process., 11, 137–147,<137::AID-HYP405>3.0.CO;2-2, 1997. a

Webb, B. W., Clack, P. D., and Walling, D. E.: Water-air temperature relationships in a Devon river system and the role of flow, Hydrol. Process., 17, 3069–3084,, 2003. a

Wenger, S. J., Isaak, D. J., Dunham, J. B., Fausch, K. D., Luce, C. H., Neville, H. M., Rieman, B. E., Young, M. K., Nagel, D. E., Horan, D. L., and Chandler, G. L.: Role of climate and invasive species in structuring trout distributions in the interior Columbia River Basin, USA, Can. J. Fish. Aquat. Sci., 68, 988–1008,, 2011. a

Worldbank: World Bank national accounts data, and OECD National Accounts data files, available at:, last access: 1 February 2020. a

WSA: Water temperature and discharge Cologne, Wasserstraßen- und Schifffahrtsamt Duisburg-Rhein, available at, last access: 1 March 2019. a, b

Yamarone, R.: Indexes of Leading, Lagging, and Coincident Indicators, John Wiley and Sons, Inc., Hoboken, New Jersey, USA,, 2012. a, b, c

Zavarsky, A. and Düster, L.: Water temperature and discharge at montoring station along the Rhine, Zenodo,, 2020. a

Zhu, S., Heddam, S., Nyarko, E. K., Hadzima-Nyarko, M., Piccolroaz, S., and Wu, S.: Modeling daily water temperature for rivers: comparison between adaptive neuro-fuzzy inference systems and artificial neural networks models, Environ. Sci. Pollut. R., 26, 402–420,, 2018. a

Short summary
River water temperature is an important parameter for water quality and an important variable for physical, chemical and biological processes. River water is also used as a cooling agent by power plants and production facilities. We study long-term trends in river water temperature and correlate them to meteorological influences and power production or economic indices.