A hybrid model for river water temperature as a function of air temperature and discharge

In order to test performances of the model, we have selected three rivers in Switzerland characterized by different hydrological conditions. The three cases are briefly described below and in table 1 (for further details please refer to the supplementary material). For a more extensive description of the Swiss database see, e.g., Schaefli et al (2015).

• 1.  
  The river Mentue flows at low altitude on the Swiss plateau through a sparsely inhabited area with predominant agricultural land use, unaffected by strong anthropic thermal alterations. This case is representative of the typical conditions of a small 'natural low-land' river.
• 2.  
  The station of the river Rhône at Sion lies at the bottom of a populated alpine valley. Starting from the beginning of the 20th century its hydrological regime has been altered by the construction of a large high-head hydropower storage system, and the river is now affected by strong hydro- and thermo-peaking (Meile et al, 2011). This case is considered as representative of a strongly 'regulated' river.
• 3.  
  The river Dischmabach is located at high altitude in a steep glacial valley that is uninhabited and used for mountain pastures, with a significant influence of snow melting. This case is taken as a representative example of a 'snowmelt-fed' river.

The model is based on a lumped heat budget that considers an unknown volume V of the river reach, its tributaries (implicitly considering both surface and subsurface water fluxes), and the heat exchange with the atmosphere, as depicted in figure 1. The variation of water temperature T_w in this volume is governed by the equation

$$\begin{eqnarray}&&\rho \;{c}_{p}\;V\;\displaystyle \frac{{{\rm{d}}{T}}_{w}}{{\rm{d}}{t}}=A\;H+\rho \;{c}_{p}\left(\displaystyle \sum _{i}{Q}_{i}\;{T}_{w,i}-Q\;{T}_{w}\right),\end{eqnarray} \tag{ 1 }$$

where t is time (hereafter expressed in days), ρ is water density, c_p its specific heat capacity at constant pressure, A is the surface area of the river reach, H is the net heat flux at the river-atmosphere interface, Q is the discharge at the downstream section, Q_i and T_wi are discharge and temperature of the i-th contributing water flux (tributaries and, possibly, groundwater), and V is the total volume that reacts to the heat fluxes. We notice that V is not necessarily limited to the surface water body, possibly including also the hyporheic region, so that the explicit inclusion of heat fluxes at the stream-river bed interface is not necessary.

Zoom In Zoom Out Reset image size

The first term on the right hand side of equation (1) contains the net heat flux H, which is primarily a function of short- and long-wave radiation, and latent and sensible heat fluxes. Following Caissie (2006), we assume that air temperature, T_a, can be used as a proxy for all these processes. The overall effect is included in the model in a linear form by using a Taylor series expansion (similarly to as in Piccolroaz et al (2013), see details in the supplementary material), such that

$$\begin{eqnarray}&&H={\left.H\right|}_{{\bar{T}}_{a},{\bar{T}}_{w}}+{\left.\displaystyle \frac{\partial H}{\partial {T}_{a}}\right|}_{{\bar{T}}_{a},{\bar{T}}_{w}}\\ &&\quad \times ({T}_{a}-{\bar{T}}_{a})+{\left.\displaystyle \frac{\partial H}{\partial {T}_{w}}\right|}_{{\bar{T}}_{a},{\bar{T}}_{w}}({T}_{w}-{\bar{T}}_{w}),\end{eqnarray} \tag{ 2 }$$

where ${\bar{T}}_{a}$ and ${\bar{T}}_{w}$ are the long-term averaged values (hereafter indicated by an overbar) of air and water temperatures, respectively. This relation can also be rewritten as

$$\begin{eqnarray}&&H=\rho \;{c}_{p}\left({h}_{0}+{h}_{a}{T}_{a}-{h}_{w}{T}_{w}\right),\end{eqnarray} \tag{ 3 }$$

where h₀, h_a and h_w are parameters that can be directly derived from equation (2).

The second term on the right hand side of equation (1) represents the difference between the heat flux leaving the control volume, $\rho {c}_{p}{{QT}}_{w},$ and the sum of all the entering heat fluxes. These are almost impossible to determine precisely, and are inherently case-specific. Noting that $Q={\displaystyle \sum }_{i}{Q}_{i},$ we introduce a fictitious mean inflow temperature ${\widehat{T}}_{w}={Q}^{-1}\;{\displaystyle \sum }_{i}{Q}_{i}\;{T}_{{wi}}.$ Such a reference temperature is likely to change in the different seasons, so we express it as a first approximation in the form of sinusoidal annual variation

$$\begin{eqnarray}&&{\widehat{T}}_{w}={\widehat{T}}_{w0}+{\widehat{T}}_{w1}\;\mathrm{cos}\left[2\pi \left(\displaystyle \frac{t}{{t}_{y}}-\varphi \right)\right],\end{eqnarray} \tag{ 4 }$$

where the reference value ${\widehat{T}}_{w0},$ the amplitude of the variation ${\widehat{T}}_{w1}$ and the phase $\varphi \in [0,1]$ are left unspecified, with t_y being the duration of a year expressed in the same units as t.

Thus, we can rewrite equation (1) as

$$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}{T}}_{w}}{{\rm{d}}{t}} & = & \displaystyle \frac{A}{V}\left({h}_{0}+{h}_{a}{T}_{a}-{h}_{w}{T}_{w}\right)\\ & & +\displaystyle \frac{Q}{V}\left({\widehat{T}}_{w}-{T}_{w}\right),\end{eqnarray} \tag{ 5 }$$

which contains two characteristic quantities. The inverse of the first ratio, V/A, is related to the flow depth D (and hence depends on the discharge), but does not coincide with it. In fact, a portion of the saturated sediments should be included, especially for shallow flows with high transparency where the incoming shortwave radiation can directly warm the river bed (Evans et al, 1998, Caissie 2006). The second ratio in the equation, Q/V, represents the inverse of the residence time in the river reach, and varies with the discharge as well. In order to allow for an explicit dependence of these two ratios on easily measurable quantities, we introduce the dimensionless reference volume, δ, and we assume a simple power law as a function of the discharge:

$$\begin{eqnarray}&&\delta =\displaystyle \frac{V}{\bar{V}}={\theta }^{1/m},\end{eqnarray} \tag{ 6 }$$

where $\theta =Q/\bar{Q}$ is the dimensionless discharge and m is related with the exponent of the rating curve (e.g., $m=5/3$ according to the Manning's normal flow relationship in a wide rectangular channel). Introducing these relations into equation (5), we obtain the final form of the air2stream model

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{T}}_{w}}{{\rm{d}}{t}}=\displaystyle \frac{1}{\delta }\left\{\phantom{{{\frac{S}{A}}^{{\frac{S}{H}}^{a}}}^{H}}\right.\phantom{\rule{-1.8em}{0ex}}\left\{{a}_{1}+{a}_{2}{T}_{a}-{a}_{3}{T}_{w}\right.\\ &&\quad +\left.\theta \left[{a}_{5}+{a}_{6}\mathrm{cos}\left(2\pi \left(\displaystyle \frac{t}{{t}_{y}}-{a}_{7}\right)\right)-{a}_{8}{T}_{w}\right]\right\},\\ &&\quad \delta ={\theta }^{{a}_{4}},\end{eqnarray} \tag{ 7 }$$

where the eight parameters a₁–a₈ of the model (7) directly follow from equations (3), (4) and (6). It is important to stress that the values of these parameters are estimated through calibration, so that all the geometrical characteristics of the river reach (length, volume, area, etc) are not required to be explicitly specified, nor the role of specific heat inputs (e.g., internal friction). This feature of the model also implies that local effects, such as that of shading, are expected to be automatically captured by the calibration of the parameters. However, we note that there are cases where air temperature and discharge alone may not be sufficient for an accurate estimate of the net heat flux. A possible example is given by rivers affected by intermittent industrial point releases generally characterized by modest water fluxes at relatively high temperature, whose effect on stream thermal dynamics may be difficult to reproduce especially when these releases are irregular.

Equation (7) resembles the model air2water developed by Piccolroaz et al (2013) for lakes and reservoirs. However, two main differences exist: (i) rivers are typically not stratified, so δ is not a function of T_w, but rather of Q; and (ii) the effect of the contributing water fluxes, both natural and regulated (e.g., affected by hydropower release), is included as it can be substantial for some types of rivers and under specific circumstances. With particular reference to the second point, here we propose a few alternative versions of the 8-parameter formulation (7).

The first variation comes from disregarding the second term on the right hand side of equation (1). The resulting 4-parameter model reads

$$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}{T}}_{w}}{{\rm{d}}{t}} & = & \displaystyle \frac{1}{\delta }\left({a}_{1}+{a}_{2}{T}_{a}-{a}_{3}{T}_{w}\right),\\ & & \delta ={\theta }^{{a}_{4}},\end{eqnarray} \tag{ 8 }$$

and corresponds to assuming that the mean temperature of the contributing water fluxes is approximately equal to that of the river (i.e., ${\widehat{T}}_{w}\simeq {T}_{w}$). Another case in which the model (8) is valid is when these fluxes are absent or their effect is not relevant, so that ${\widehat{T}}_{w}$ is close to the temperature of the river at the upstream end of the considered reach. In general, this model can be successfully applied if the longitudinal gradient of temperature is small.

We will see in the following analysis that the effect of δ is minor in the three examined cases, which implies that the thermal response of the river does not strongly depend on the flow depth at a daily timescale. Therefore, we identify two new versions of the model by setting $\delta =1$ (i.e., ${a}_{4}=0$): a 7-parameter model derived from (7), and a 3-parameter model derived from (8). The latter model is the simplest that can be used to simulate river temperature variations.

A further version of the model can be obtained by simplifying the 7-parameter version imposing also $\theta =1,$ i.e. assuming that the effect of the discharge can be approximately retained using only a constant value. In this way, by combining the constant terms and those proportional to T_w, we obtain a 5-parameter version

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{T}}_{w}}{{\rm{d}}{t}}={a}_{1}+{a}_{2}{T}_{a}-{a}_{3}{T}_{w}\\ &&+{a}_{6}\;\mathrm{cos}\left[2\pi \left(\displaystyle \frac{t}{{t}_{y}}-{a}_{7}\right)\right],\end{eqnarray} \tag{ 9 }$$

where the interpretation of the parameters a₁, a₃ and a₆ changes from that of the initial model. The summary of the five versions of the air2stream model is reported in table 2 (the source code is available at https://github.com/marcotoffolon/air2stream).

Table 2.  Versions of the model and performances for the three case studies.

Version	Parametersa										RMSE Tw [°C] (NSE [−])
of the	a1	a2	a3	a4	a5	a6	a7	a8	δ	θ	Case 1		Case 2		Case 3
model											calb	valc	cal	val	cal	val
8-par	x	x	x	x	x	x	x	x	$f(\theta )$	θ	0.62 (0.99)	0.76 (0.98)	0.58 (0.92)	0.75 (0.89)	0.61 (0.96)	0.62 (0.95)
7-par	x	x	x	—	x	x	x	x	1	θ	0.62 (0.99)	0.78 (0.98)	0.61 (0.92)	0.77 (0.89)	0.64 (0.95)	0.65 (0.95)
5-par	x	x	x	—	—	x	x	—	1	—	0.65 (0.99)	0.78 (0.98)	0.88 (0.83)	1.05 (0.79)	0.70 (0.94)	0.70 (0.94)
4-par	x	x	x	x	—	—	—	—	$f(\theta )$	—	0.77 (0.98)	0.92 (0.98)	0.90 (0.81)	1.04 (0.79)	0.88 (0.91)	0.87 (0.91)
3-par	x	x	x	—	—	—	—	—	1	—	0.78 (0.98)	0.91 (0.98)	0.91 (0.81)	1.05 (0.79)	0.89 (0.90)	0.88 (0.91)
${T}_{w,{eq}}$ from 8-par (6 parameters)									—	θ	0.71 (0.99)	0.83 (0.98)	0.59 (0.92)	0.75 (0.89)	0.64 (0.95)	0.66 (0.95)
${T}_{w,{eq}}$ from 5-par (4 parameters)									—	—	0.72 (0.98)	0.83 (0.98)	0.87 (0.83)	1.05 (0.79)	0.76 (0.93)	0.78 (0.93)
${T}_{w,{eq}}$ from 4-par (2 parameters)									—	—	0.96 (0.97)	1.07 (0.97)	0.91 (0.82)	1.05 (0.79)	0.94 (0.89)	0.94 (0.89)
Logistic regression (4 parameters)									—	—	0.98 (0.97)	1.10 (0.97)	0.71 (0.89)	0.86 (0.86)	0.89 (0.90)	0.89 (0.90)

^a'x' denotes parameters that are calibrated, '-' parameters that are not included in the version of the model. ^bCalibration. ^cValidation.

The successful application of the model relies on a careful calibration. The best set of parameters is identified by defining a physically reasonable range of the parameter values and solving an optimization problem. The aim is to minimize the error between simulated and measured water temperatures, defined through the introduction of a suitable metric (objective function). The input data are a continuous series of daily averaged air temperature and river discharge (the latter only for some versions, see table 2). The output of the model is the series of daily averaged water temperature, which is compared with measured values. In the most favourable case, observed water temperature is a series of continuous daily values, but the model can be calibrated also with a reduced number of data or aggregated values (e.g., weekly or monthly averages, see Toffolon et al (2014) for lakes). Air and water temperature series should be long enough to capture interannual variability, thus in general should extend for several years. As a general rule, we used 2/3 of the available series for calibration and 1/3 for validation of the model.

The objective function used for model calibration is the root mean square error (RMSE) between simulated and observed T_w. Moreover, we calculated the Nash-Sutcliffe Efficiency (NSE) index that normalizes the variance of the errors with the variance of the measurements (Moriasi et al, 2007). We recall that NSE equal to 1 corresponds to perfect agreement, while negative values indicate that the mean would have been a better predictor than the tested model.

The identification of the best set of parameters for each river is performed through a Monte Carlo-based optimization procedure in which a large number of parameter sets is sampled and evaluated in terms of RMSE. Here we use the particle swarm optimization (PSO) algorithm, a population-based stochastic optimization technique firstly proposed by Kennedy and Eberhart (1995) (see supplementary materials for further details) with a population of 1000 particles and for 1000 iterations. For each run, equation (7), or any other version of the model, is numerically integrated using the Runge–Kutta 4th order method with a time step of one day. Computed water temperatures below 0^oC have been limited to that value.

The five versions of the model were applied to the three case studies (table 1), showing a general good agreement between model results and measured water temperature. Table 2 reports the values of RMSE, which are always smaller than 1^oC in calibration (i.e., first 2/3 of the series), with only minor increases for some cases in validation. The values of NSE are also reported, which are always larger than 0.8 and in some cases very close to 1, thus suggesting that the model performs well also in relative terms.

In order to visually appreciate seasonal fluctuations in temperature and discharge of the three case studies, we calculated the reference year (daily means) of all variables (measured air and water temperature, discharge, and modeled water temperature) by averaging across all years the values for each day of the year. Figure 2 shows the performances of the 8-parameter version (the results of the other versions and the values of the parameters are reported in the supplementary material) and demonstrates that the model effectively follows the measured annual oscillation. We note that the phase is not necessarily coincident with that of T_a because of the contribution of Q and of the additional sinusoidal term, possibly characterized by a different phase.

Zoom In Zoom Out Reset image size

Interannual variability of river water temperature is examined in figure 3 during a period of three years. For all rivers the model is able to reproduce the sequence of cold and warm years, as for example the particularly warm 2006 summer and 2006/2007 winter in case 1, the warm 1987/1988 winter in case 2, and the cold 2008 spring in case 3 associated with the occurrence of larger discharge values.

Zoom In Zoom Out Reset image size

The three examined rivers behave in a markedly different way. The response of water temperature to changes in air temperature is almost linear for the low-altitude river (case 1), thus producing large differences between annual maximum and minimum water temperature. Indeed, case 1 is characterized by a relatively low discharge, which typically corresponds to a shallow depth and hence a low thermal inertia (but larger lowland rivers may behave differently). Conversely, although at the same altitude, case 2 presents a much more limited increase of water temperature, especially in summer. This is attributable to the release of cold water during hydropower production, which has been shown to affect the thermal behavior of the rivers at different temporal scales (Zolezzi et al, 2011). Since thermal alterations travel for long distances (Toffolon et al, 2010), the impact can be felt along the whole river in alpine contexts, especially if there are not other significant tributaries. The resulting effect is a marked flattening of the natural seasonal pattern (Meier and Wüest 2004). Finally, water temperature in case 3 is generally low during the whole year because of the high altitude and the contribution of snow melting waters (Lisi et al, 2015). The influence of an upstream snowfield is clearly distinguishable in the abrupt variations of the discharge due to snow melting starting in spring and slowly decreasing during summer and autumn, with almost no flow in winter. The effect of the cold melting water is similar to that of hydropower releases: a much lower response of T_w to T_a during the warm season. We also note that a large contribution of cold groundwater inputs can produce similar effects (Tague et al, 2007).

For the selected cases, and over the specific study periods, the 8-parameter model always performs better than the other versions (table 2), but the improvement is different in the three types of rivers. In particular, the effect of discharge cannot be neglected in case 2 (strong influence of hydropower generation) and case 3 (significant contribution of snow melt) because the water flowing at the river station originates from a number of colder contributions. As a result, the model depending on air temperature only is not suitable to reproduce the proper dynamics, although the 5-parameter version yields a substantial improvement without considering Q explicitly. This is particularly true for case 3, for which the phase shift induced by the additional sinusoidal term (representing the mean inflow temperature) is able to postpone the warming of the river water with respect to that of the air in spring, mimicking the effect of cold contributions due to snow melting water. Conversely, the effect of discharge is not so relevant in case 1, and the 3-parameter version already produces satisfactory results (table 2). In this case, the river does not receive large tributaries, so its temperature is approximately spatially uniform along large distances (i.e., ${\widehat{T}}_{w}\simeq {T}_{w}$). Anyway, a significant improvement is given by the 5-parameter model, which compensates the small phase lag that is however present between air and water temperature, and by the 7- and 8-parameter models, despite Q is less variable during the year than in the other two rivers.

Interestingly, the performances of the 7-parameter and 3-parameter versions are comparable to those of the 8-parameter and 4-parameter versions, respectively, for the three examined rivers. This suggests that the role of δ, which represents the variation of the thermal inertia of the river as a function of Q, is minor. This result seems to contradict the common intuition that smaller flow depths (due to smaller Q) would result in a warmer river water temperature in summer, and, conversely, in a colder river water temperature in winter. Actually, δ influences the reaction time to external forcing, but not the value at equilibrium. This can be easily shown by considering the structure of the model, as it is discussed in the supplementary material where analytical solutions are derived in idealized cases. There, we show that $\delta /{A}_{3}$ (with ${A}_{3}={a}_{3}+\theta \;{a}_{8}$ for the 8-parameter model) controls the time scale for the adaptation of T_w to the external forcing. Referring for simplicity to the 5-parameter model (9), which does not include θ (hence ${A}_{3}={a}_{3}$), and assuming $\delta \simeq 1,$ the time scale is 0.94, 0.55 and 0.36 days for cases 1, 2 and 3, respectively. This means that water temperature adapts to the external forcing in a period that is comparable with the averaging window (1 day) assumed in this analysis. Such a consideration can be supported by rewriting equation (7) as

$$\begin{eqnarray}\displaystyle \frac{\delta }{{A}_{3}}\displaystyle \frac{{{\rm{d}}{T}}_{w}}{{\rm{d}}{t}} & = & \displaystyle \frac{{a}_{1}+\theta \;{a}_{5}}{{A}_{3}}+\displaystyle \frac{{a}_{2}}{{A}_{3}}{T}_{a}\\ & & +\displaystyle \frac{\theta \;{a}_{6}}{{A}_{3}}\;\mathrm{cos}\left(2\pi \left(\displaystyle \frac{t}{{t}_{y}}-{a}_{7}\right)\right)-{T}_{w},\end{eqnarray} \tag{ 10 }$$

from which it clearly appears that when A₃ is large compared to the temperature variations (multiplied by δ, which however is $\sim O(1)$ on average) the left hand side vanishes. A more precise condition has been derived in the supplementary material by also taking the time derivative into account by means of a perturbation method: the ratio ${A}_{2}/{A}_{3}^{2}$ (where A₂ represent the amplitude of an approximately equivalent sinusoidal forcing) has to be lower than a threshold, confirming the strong inverse dependence on A₃. When such a condition holds, the thermal response becomes substantially independent of δ and an instantaneous equilibrium water temperature ${T}_{w,{eq}}$ can be adopted as an approximation (Caissie et al, 2005).

Neglecting the time derivative in equation (10) and rearranging the parameters by dividing everything by a₃, the 8-parameter and 7-parameter versions reduce to

$$\begin{eqnarray}{T}_{w,{eq}} & = & \displaystyle \frac{1}{\left(1+\theta \;{e}_{4}\right)}\left\{\phantom{{{\frac{S}{A}}^{{\frac{S}{H}}^{a}}}^{H}}\right.\left\{\phantom{\rule{-2em}{0ex}},{e}_{1}+{e}_{2}{T}_{a}\right.\\ & & +\left.\theta \left[{e}_{5}+{e}_{6}\;\mathrm{cos}\left(2\pi \left(\displaystyle \frac{t}{{t}_{y}}-{e}_{3}\right)\right)\right]\right\},\end{eqnarray} \tag{ 11 }$$

where ${e}_{1}={a}_{1}/{a}_{3}$, ${e}_{2}={a}_{2}/{a}_{3}$, ${e}_{3}={a}_{7}$, ${e}_{4}={a}_{8}/{a}_{3}$, ${e}_{5}={a}_{5}/{a}_{3}$, and ${e}_{6}={a}_{6}/{a}_{3}.$ The approximation (11) explicitly illustrates how the effect of the discharge is felt only through θ, suggesting that, besides the exchanges with the atmosphere, it is the contribution of the inflows (with a temperature that may be different from that of the receiving river) that relevantly affects water temperature. Similarly to equation (11), the 5-parameter version reduces to

$$\begin{eqnarray}&&{T}_{w,{eq}}={e}_{1}+{e}_{2}{T}_{a}+{e}_{6}\;\mathrm{cos}\left(2\pi \left(\displaystyle \frac{t}{{t}_{y}}-{e}_{3}\right)\right),\end{eqnarray} \tag{ 12 }$$

and the 4-parameter and 3-parameter versions to

$$\begin{eqnarray}&&{T}_{w,{eq}}={e}_{1}+{e}_{2}{T}_{a}.\end{eqnarray} \tag{ 13 }$$

These equilibrium relationships are characterized by a number of parameters that is smaller than the original model because a₃ is used to rescale the other parameters. Thus, equation (11) is characterized by only 6 parameters to calibrate. The other relationships (12) and (13) have only 4 and 2 free parameters, respectively, the latter representing the equivalent of a linear regression with a lower bound at 0 ^oC. The performances of the equilibrium relationships are compared with the complete models in table 2. In the same table, we also report the RMSE and NSE of the most widely used nonlinear regression model, based on the logistic function (e.g., Arismendi et al, 2014), whose formulation is explicitly presented in the supplementary material and contains 4 parameters to calibrate. The 8- and 7-parameter versions of air2stream outperform the logistic regression for all rivers, both in the differential form and under the equilibrium assumption. The other versions have higher performances, as well, with the only exception of case 2. In the regulated river, in fact, the S-shape of the logistic function is able to capture the flattening of water temperature in summer, but does not account for the actual physics governing the thermal process. The information about Q is indeed essential, that is the reason why the simplest versions of the model are less adequate. In conclusion, our model performs better than the common statistical models and moreover includes the most relevant physical factors. A direct comparison is not possible with fully mechanistic models, which would require much more information. Finally, we note that, differently from the linear and nonlinear regression models that can reproduce only a biunivocal dependence of T_w as a function of T_a, the solutions (11) and (12) allow for a more realistic hysteretic pattern T_w-T_a, despite the complexity of the formulation is equivalent.