Anthropogenic influence on the Rhine water temperatures

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 timelagged 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.


Introduction
River water temperature (T w ) greatly influences the most important physical and chemical processes in rivers and is a key factor for river system health (Delpla et al., 2009). T w also confines animal habitats (Gaudard et al., 2018;Isaak et al., 2012;Durance and Ormerod, 2009), 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 Lilliestam, 2010). 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 Lilliestam, 2010). This leads to a significant heat load, even on large rivers such as the Rhine (IKSR, 2006;Lange, 2009). 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 T w (Cai et al., 2018;Gaudard et al., 2018;Råman Vinnå et al., 2018). The natural influences on T w 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 groundwaterfed; (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, T w can be modeled in different ways. The common options are statistical models and physically based models.
A physically based T w model (Sinokrot and Stefan, 1993) 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 Piccolroaz, 2015). Their approach includes more parameters and, thus, is more complex.
Published by Copernicus Publications on behalf of the European Geosciences Union.

5028
A. Zavarsky and L. Duester: Rhine water temperatures 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 (T a ) as a proxy for sensible, latent and radiative heat fluxes (ground heat flux can be neglected) and establish a T a → T w relationship through regression. T a 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 Preudhomme, 1993;Mohseni et al., 1998;Koch and Grünewald, 2010). Generally, exponential models deliver better results with temperature extremes. However, they lack the distinct separation between contribution to T w from anthropogenic heat input and natural influences. Using linear models, Markovic et al. (2013) show that between 81 % and 90 % of the T w variability can be described by T a . 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 T w . 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 T w 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 Peterson, 2010;Isaak et al., 2010;Peterson and Hoef, 2010;Isaak et al., 2014). The correlation of the variables (e.g., T a , T w discharge, etc.) which influence other T w 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 (Hassoun, 1995). ANN use a sample data set to train artificial neurons the relationship between input (e.g., air temperature) and output (T w ) (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 (Lange, 2009). 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 T w of the Rhine. Short-term economic changes, observable in the change in the gross domestic product (GDP), may influence T w 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 T w along the Rhine at four monitoring stations (Basel, Worms, Koblenz and Cologne).

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 (T c ). T c and discharge at the measurement station Q were used in a multiple linear regression T c → T w (Eq. 1). The resulting regression coefficients a 1 , a 2 and a 3 describe the magnitude of the respective influences (anthropogenic heat input and meteorological and hydrological influences).
T w = a 1 + a 2 · T c + a 3 · Q. (1) Using an improved calculation method for T c , which includes catchment-wide averaging with river size weighing and a time lag, the regression should deliver a better estimate for a 1 , a 2 and a 3 . The model was run on a T w 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 T w and an important factor in the T c → T w relationship. Additionally, it may function as a measure of how fast the water mass responds to changes in T w . 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 Stefan, 1993, for the Mississippi;Caissie, 2006, as a review article).
Using the multiple regression (Eq. 1), we especially investigated the change in a 1 over time, which we call, in this study, the Rhine base temperature (RBT). This temperature represents T w , without the influence of meteorology (T a ) 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 T a and Q.

Water temperature and discharge
We used a data set of daily averaged T w 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.
T w was measured by platinum resistivity sensors (Pt100). The accuracy of these sensors is commonly ±0.5 • C, but the precision, which describes the ability to detect temperature changes, is 0.05 • C. As we focused on the change of T w 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 T w was a linearly dependent proxy for the average river temperature. Q was provided as daily averages in m 3 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.

Air temperature
T a 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 T a on a 1 4 • by 1 4 o grid. The data set is available from 1979 to 2018. We took the hourly T a output and calculated a daily mean for each grid point between 1979 and 2018 to fit the time resolution of T w .

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; IAEA, 2019). A total of 12 NPPs (1986)(1987)(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 exclu-sive river water cooling system have a conversion factor of 3, which is based on the power efficiency of electricity generation (Lange, 2009). 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 (Lange, 2009).
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.

Calculated temperature change
We calculated the expected change T w based on a change in HI ( HI) from NPPs using the average discharge Q, the heat capacity of the water c p and the water density ρ; see Eq. (3).
This approach follows the idea that the contribution of NPPs significantly alters the T w and only influences the RBT fraction.

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 lowfrequency cutoff was canceling long-term trends, as the GDP change was only related to the previous year. The highfrequency cutoff was used to dampen fast, alternating RBT signals in comparison to the slow sampled GDP data.

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 .   Y is the average over the total time series, σ is the standard deviation of the whole time series and Y i 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.

Catchment area
The catchment area was calculated using the HydroSHEDS database (Lehner et al., 2008). The 1 125 • by 1 125 • 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. Addition-ally, 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.

Multiple regression
A multiple linear regression was used to separate the anthropogenic heat input a 1 and meteorological a 2 and hydrological a 3 contributions to the river water temperature. T w was regressed with T c and river discharge Q. Their regression coefficients a 2 (T c slope) and a 3 (Q slope) represent the magnitude of the respective influences. The offset a 1 (RBT) combines all other influences which were not related to a change in T c 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 T a directly at the monitoring station, we improved Eq. (1) by a time-dependent and weighed average of T a (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.
The new representative catchment temperature was called T c (x 0 , y 0 , t 0 ). It was a weighed average of the whole catchment area (x, y), which is defined by the measurement point (x 0 , y 0 ). The difference between the measurement time t 0 and the reading of T a 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.

Time lag
A change in T w is slower than a change in T a . 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 T w through new meteorological conditions and heat fluxes takes time. Therefore, linear and exponential models include either a fixed t for T a (Eq. 6) or an average of T a , including a time span from before (Eq. (7); Stefan and Preudhomme, 1993;Nobilis, 1995, 1997;Haag and Luce, 2008;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 T a , in space and time, during advection into account. This is especially important for daily averaged T w (Erickson and Stefan, 2000). 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.
We combined the general idea of a time lag and averaging T a 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).

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. 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  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 T a on large water masses. Figure 3b shows ACC ·w over the whole catchment area of Koblenz.

T c
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 T a , which is shown in Eq. (10).
· T a (x, y, t 0 + t (x, y)) . (10) T c (x 0 , y 0 , t 0 ) was calculated by weighed (ACC ·w) averaging T a (x, y, t + t (x, y)) over all grid points of the catchment area (x = 1, . . . ny = 1, . . . m), which was set by the measurement point (x 0 ,y 0 ). 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 T c 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.

T c calculation methods
Additionally, we used these four calculations methods, namely (1) w + t, (2) avg + t, (3) avg and (4) point, to compare the results of the linear regression to the calculation proposed in Eq. (10).
(2) Here, no weight and only time lag (Eq. 12) are shown, as follows: T a (x, y, t 0 + t (x, y)) . (12) (3) A mean T a (x 0 , y 0 , t 0 ) over the whole catchment area at the time t 0 of the measurement was calculated; see Eq. (13).
t was not used in the following: (4) T a (x 0 , y 0 , t 0 ) at the location x 0 , y 0 and time t 0 of the measurement (see Eq. 14) is shown, as follows: 3 Results

Water temperature time series
To investigate the long-term changes, we fitted a timedependent linear function to the time series of T w and T a (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 T w and the linear fits for the two time periods. The average T a of the catchment area is also shown. In Fig. 5b, the RAPS index of T a and T w is shown. The fit coefficients and the rate of warming per year are displayed in Table 4. The calculated T a 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 and Table 4 show that the change in T w 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, T a (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 T w warming  rate from 1985 to 2018 for Worms and Koblenz is in between those from Cologne and Basel. These two stations show similar T a warming rates compared to Basel and Cologne. Generally, the T a warming rates have a difference of less than 5 %. Arora et al. (2016) showed a mean T w warming rate in northern and northeastern German rivers of 0.03 • C yr −1 (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and 0.09 • C yr −1 (2000-2010). Regarding our time period , these values are plausible. Basarin et al. (2016) found a maximum increase in T w at the Danube at Bogojevo (1950-2012) of 0.05 • C yr −1 , which matches the maximum increase at Basel. T a 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 T a 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 T a already increased again (Hu and Fedorov, 2017). Michel et al. (2020) investigated T w at 52 river gauges in Switzerland, representing most of the Rhine catchment area at Basel. The authors reported an average T w increase at the 52 stations of 0.037 • C yr −1 (1998-2018) and 0.033 • C yr −1 . T a increased by 0.039 • C yr −1 (1998-2018) and 0.046 • C yr −1 . Comparing this to our results at Basel, the T a 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 T w 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 R 2 values make the differences between the four monitoring stations visible. Basel exhibits the largest R 2 values, and these are consistently high for T a and T w . This is in contrast to the station at Cologne, where R 2 of T w was low and statistically not significant. The slope of T a at Cologne is lower than at the other stations but still statistically significant. The Pearson correlation coefficients between T a and T w were the lowest at Cologne and the largest in Basel. For T a , 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 T w RAPS index for Basel displayed the same pattern in the sections as the T a index. All other stations showed other RAPS T w to RAPS T a patterns. This means that the T a and T w trends of the original time series were different at these stations. T a cannot fully describe the trends in T w .
We hypothesized that different meteorological conditions were not the reason for this difference. Meteorological differences should also have been visible in the T a warming rate of the four stations, which was not the case. The RAPS analysis for T a and T w only coincided within the Basel data set.

Regression
We fitted the multiple regression model (Eq. 5) using T c and Q to T w of each monitoring station for the available data set. Afterwards, we recalculated T w,modelled using the regression coefficients a 1 , a 2 and a 3 . From the comparison between the T w,modelled and measured T w , 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 T c 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 T a → T w -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 T c . Large ACC grid points represent large water masses (rivers and lakes), and their influence on T a 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 T c .
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 T w 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.

Rhine base temperature
The RBT was used to explain differences in the T w warming rates of Table 4. We regressed a 2-year segment of the T w 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 Table 4. Slope 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 R 2 values, which are statistically significant if R 2 > 0.19.

Station
Slope T w whole Correlated T w ↔ T a Slope T w 1985-2018 Correlated T w ↔ T a Slope T a whole Slope T a 1985-2018 data set ( • C yr −1 ) whole data set  Table 5. RSME (in degrees Celsius) and NSC for all T c calculation methods. Different T c calculation methods and the regressions are applied over the total data set. The RMSE and the NSC are calculated between T w and T w,modelled . The first column contains the calculation method. The best results for each monitoring station and each calculation method are indicated in bold. 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).

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 T w 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.
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 T w (see Fig. 6). This occurred during a time when the NPPs HI remained relatively stable, but the GDP increased by 30 % between 198530 % between and 199530 % between (Worldbank, 2020. 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 T w along the Rhine. NPPs may also impact T w at much shorter timescales but do not change their power output accordingly.

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 a 1 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.
For Worms (Fig. 7, bottom panel), we added the change in the turnover of BASF (red dashed line; BASF AG, 1989). 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 (MWV, 2003). 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).

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 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.

Conclusions
We introduced a new catchment-wide air temperature T c , which decreased the RMSE (Table 5) in a T c → T w regression. T c is a weighed (ACC ·w) average of all T a 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 T a a water droplet experienced on its way to a certain point (in this study it is a monitoring station), and it delivered better linear T c → T w estimates. This improvement in the T c → T w relationship may support the analysis of processes in the heat budget of rivers. Usually, T a 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 T w 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 T a data are often globally available. Morrill et al. (2005) showed a linear T a → T w 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 T c , we regressed four T w time series (Basel, Worms, Koblenz and Cologne) along the Rhine. The offset in the this regression a 1 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 longterm 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 T w of the Rhine could increase by 0.6-1.4 • C. This trend is plausible, according to the historical data analyzed, if the T a 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 T w 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 T a , 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 https://doi.org/10.5281/zenodo.4095396 (Zavarsky and Düster, 2020). 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 dif-ferent calculation methods for the catchment wide air temperature as described in this paper.
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.