Assessing the performance of a suite of machine learning models for daily river water temperature prediction

Study sites and data

The measurements of daily water temperature (T_w), flow discharge (Q) and air temperature (T_a) were obtained from seven rivers in Europe and USA: (i) one river in Croatia, (ii) three rivers in Switzerland, and (iii) three rivers in USA. The main characteristics of the rivers with the period of records, and the calibration and validation periods are briefly summarized in Table 1.

The rescaled adjusted partial sums (RAPS) method

In order to detect and quantify trends and fluctuations in time series, the widely used rescaled adjusted partial sums (RAPS) method (Bonacci, Trninić & Roje-Bonacci, 2008; Bonacci & Oskoruš, 2010; Basarin et al., 2016) was employed in this study. The RAPS method can highlight trends, shifts, data clustering, irregular fluctuations and periodicities in the time series (Bonacci, Trninić & Roje-Bonacci, 2008), and is defined as: (1)${\displaystyle RAP{S}_i=\sum _{n=1}^i\frac{T_n-T_m}{S_T}}$

Descriptions of the variables are summarized in Table 2.

Machine learning models

In this study, we compared three machine learning models: feedforward artificial neural network (FFNN), GPR and DT. Detailed descriptions of these models can be found in Zhu, Nyarko & Hadzima-Nyarko (2018).

Air2stream

Air2stream is a hybrid model which combines a physically based structure with a stochastic calibration of parameters (Toffolon & Piccolroaz, 2015; source code available at https://github.com/spiccolroaz/air2stream). The model has been tested in various river systems, and generally it performs well. The equation for the air2stream model can be expressed as (the 8-parameter version): (2)${\displaystyle \frac{d{T}_w}{dt}=\frac{1}{{\theta }^{a_4}}\left[a_1+a_2{T}_a-a_3{T}_w+\theta \left(a_5+a_{6}cos\left(2\pi \left(\frac{t}{t_y}-a_7\right)\right)-a_8{T}_w\right)\right]}$

Descriptions of the variables can be found in Table 2. In this model version, T_w is estimated from T_a, Q and a sinusoidal term.

Assuming that the effect of flow discharge can be approximately retained using only a constant value and by combining the constant terms and those proportional to T_w, a 5-parameter model version can be obtained (Toffolon & Piccolroaz, 2015). In this model version, T_w is estimated from T_a and a sinusoidal term. (3)${\displaystyle \frac{d{T}_w}{dt}=a_1+a_2{T}_a-a_3{T}_w+a_{6}cos\left(2\pi \left(\frac{t}{t_y}-a_7\right)\right)}$

Disregarding the second term on the right hand of Eq. (2) and assuming that the effect of the flow discharge can be approximately retained using only a constant value, a 3-parameter model version can be obtained (Toffolon & Piccolroaz, 2015). In this model version, T_w is estimated from T_a only. (4)${\displaystyle \frac{d{T}_w}{dt}=a_1+a_2{T}_a-a_3{T}_w}$

The detailed information about the model can be found in Toffolon & Piccolroaz (2015) and Piccolroaz et al. (2016).

Performance indices

Model performances were evaluated using four indicators in this study: the coefficient of correlation (R), the Willmott index of agreement (d), the root mean squared error (RMSE), and the mean absolute error (MAE). (5)${\displaystyle R=\left[\frac{\frac{1}{N}\sum _{i=1}^N\left({\mathrm{O}}_i-{\mathrm{O}}_m\right)\left({\mathrm{P}}_i-{\mathrm{P}}_m\right)}{\sqrt{\frac{1}{N}\sum _{i=1}^N{\left({\mathrm{O}}_i-{\mathrm{O}}_m\right)}^2}\sqrt{\frac{1}{N}\sum _{i=1}^N{\left({\mathrm{P}}_i-{\mathrm{P}}_m\right)}^2}}\right]}$ (6)${\displaystyle d=1-\frac{\sum _{i=1}^N{\left({\mathrm{P}}_i-{\mathrm{O}}_i\right)}^2}{\sum _{i=1}^N{\left(\left|{\mathrm{P}}_i-{\mathrm{O}}_m\right|+\left|{\mathrm{O}}_i-{\mathrm{O}}_m\right|\right)}^2}}$ (7)${\displaystyle RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left({\mathrm{O}}_i-{\mathrm{P}}_i\right)}^2}}$ (8)${\displaystyle MAE=\frac{1}{N}\sum _{i=1}^N\left|{\mathrm{O}}_i-{\mathrm{P}}_i\right|}$

Descriptions of the variables can be found in Table 2.

Model versions

The models were compared using three different versions: (i) version 1: FFNN1, GPR1 and DT1 using only T_a as input variable, (ii) version 2: FFNN2, GPR2 and DT2 with T_a and DOY as input variables, and (iii) version 3: FFNN3, GPR3 and DT3 with T_a, Q, and DOY as input variables. For comparison, the 3-parameter, 5-parameter, and 8-parameter version air2stream models were used correspondingly.

Seasonal dynamics of T_a, T_w and Q

Detailed statistics for the variables selected (T_w, Q and T_a) can be found in Table 3. Mean annual T_w for Rhône at Sion and Dischmabach at Davos are smaller compared with T_w for other river stations since Sion lies at the bottom of a populated Alpine valley, and Davos is located in a steep glacial valley (Dischma). Mean annual T_a for Dischmabach at Davos is extremely colder than the other stations. Mentue at Yvonand, Dischmabach at Davos, Fanno and Irondequoit have smaller annual mean flow discharge (Q < 10.0 m³/s).

Figure 2 shows the temporal variations of annual averaged T_w and T_a for the studied river stations. It can be seen that the increasing rates of T_w are smaller than that of T_a in general. T_a displays an increasing trend except for Mentue at Yvonand and Fanno, and the increasing rates range from 0.0384 to 0.0562 °C/year. T_w displays an increasing trend except for Mentue at Yvonand and Irondequoit, and the increasing rates vary between 0.0052 and 0.0449 °C/year. Generally, T_w and T_a display the same varying pattern. However, for Fanno, T_w has an increasing trend with a decreasing trend of T_a, while for Irondequoit, T_w has a decreasing trend with an increasing trend of T_a, which indicates that the temporal variations of T_w may not be explained only by T_a. It is shown in Fig. 3 that the studied river stations present different characteristics for RAPS values. For Drava at Botovo, increases of T_w and T_a started in 2006, and during the period 1991–2005, a trend of decreasing temperature is indicated (Fig. 3A). For Drava at Donji Miholjac, increases of T_a also started in 2006, however, a trend of decreasing T_w is evident until 2013, which explains the smaller increasing rate (0.0052 °C/year) at this station (Fig. 3B) compared with that of Botovo (0.0254 °C/year). A trend of decreasing T_w and T_a is shown for Mentue at Yvonand (Fig. 3C), confirming the results in Fig. 3C. For Rhône at Sion, increases of T_w and T_a started in 2000 and 1992 respectively, and during the periods 1984–1999 and 1984–1991, a trend of decreasing temperature is noticed (Fig. 3D). A trend of increasing T_w and T_a is shown for Dischmabach at Davos (Fig. 3E). For Cedar, increases of T_w and T_a started in started in 2013 and 2012 respectively, and during the periods 2001–2012 and 2001–2011, a trend of decreasing temperature is presented (Fig. 3F). For Fanno and Irondequoit, the trends of RAPS values are not quite significant (Figs. 3G and 3H).

Figure 4 presents the seasonal dynamics of T_w, T_a and dimensionless flow discharge (DQ: the ratio of daily flow discharge to annually averaged flow discharge) for the eight river stations through the climatological year. Results showed that the seasonal variations of T_w and T_a are almost synchronous for Drava at Botovo and Donji Miholjac, Mentue at Yvonand, Fanno and Irondequoit, indicating the insignificant role of Q. The damped responses of T_w to variations in T_a can be found in Rhône at Sion, Dischmabach at Davos and Cedar, especially for Rhône River at Sion since it is impacted by cold water releases from high altitude hydropower reservoir, indicating the significant role of Q in these stations.

Model performances

All the four models (GPR, DT, FFNN and air2strteam) were applied and compared for the eight stations, and modelling performances were evaluated using RMSE, MAE, R and d values. The developed models were able to successfully predict T_w using T_a, Q, and DOY as inputs. Tables 4–11 present the statistical indicators for the eight stations. For the eight studied stations, the FFNN, GPR and DT models performed with relatively similar accuracy when only T_a was used as predictor, with DT model performing slightly worse than FFNN and GPR. When DOY was included as model input, modelling performances of the three FFNN, GPR and DT models dramatically improved in terms of higher R, d and lower RMSE and MAE values (Tables 4–11). The comparison between the improvement in river T_w prediction obtained using model version 3 and version 2 indicates that the inclusion of discharge provided a lower gain than including DOY for all the cases (Tables 4–11). The relative importance of Q increases for the studied station in the Rhône, Dischmabach and Cedar Rivers, which suggests that the role of Q on river T_w modelling is more relevant for high altitude rivers (cold water release from hydropower reservoirs, snow melting and etc.). However, for lowland hydrological stations impacted by hydropower reservoirs, such as the Drava River at Botovo, the influence of Q is not that remarkable.

Results at Botovo station are reported in Table 4. Figure 5 illustrates the scatterplot of measured and calculated values of T_w using the best models with T_a, Q, and DOY as input variables (version 3). During the validation phase, inspection of statistical metrics for different models shows that the RMSE and MAE are lowest for the three versions of the air2strteam model (Table 4), and the best accuracy was obtained using air2strteam3 (RMSE = 0.891 °C and MAE = 0.724 °C). Averaged RMSE and MAE values for the three air2strteam models are 0.986 °C and 0.797 °C respectively. Compared to the three machine learning models, air2strteam3 is more accurate, followed by GPR3 (RMSE = 1.302 °C and MAE = 1.042 °C) and FFNN3 (RMSE = 1.307 °C and MAE = 1.040 °C) with similar accuracy, and finally the DT3 is ranked in the third place with the highest RMSE (1.366 °C) and MAE (1.086 °C) values respectively.

Results at Donji Miholjac station are reported in Table 5. Figure 6 illustrates the scatterplot of measured and calculated values of T_w using the best models with T_a, Q, and DOY as input variables. Two important points must be highlighted. Firstly, using the air2strteam, there is no significant difference between air2strteam1, air2strteam2 and air2strteam3. The R and d values were slightly improved between air2strteam1 and air2strteam3 (0.3% and 0.1%). In addition, the RMSE and MAE did not change significantly between the two models: 8.9% reduction for the RMSE, and 9.27% reduction for the MAE. Secondly, by comparing the accuracy of the three machine learning models, it is clear that the FFNN, GPR and DT models worked with slight difference. In addition, DT1 possess the lowest accuracy with high RMSE (2.732 °C) and MAE (2.134 °C), slightly higher than GPR1 (RMSE = 2.671 °C, MAE = 2.090 °C), while the FFNN1 was the best model with lowest RMSE and MAE values.

The results at Mentue station are summarized in Table 6. Figure 7 illustrates the scatterplot of measured and calculated values of T_w using the best models with T_a, Q, and DOY as input variables. As can be seen from Table 6, during the validation phase, the FFNN1, GPR1 and DT1 were practically identical. The difference in the RMSE and MAE values were small and marginal. The FFNN1 model was slightly more accurate than the GPR1 and DT1, with MAE of 0.982 °C and RMSE of 1.307 °C. However, using only T_a as input variable, air2strteam1 was more accurate than the three machine learning models. The difference in the performance of these algorithms for the prediction of T_w is discussed hereafter. The RMSE and MAE values were decreased from 1.307 °C to 0.933 °C (28.62%), and from 0.982 °C to 0.679 °C (30.86%) compared to FFNN1. Compared to the GPR1, the RMSE and MAE values were decreased from 1.311 °C to 0.933 °C (28.83%), and from 0.992 °C to 0.679 °C (31.55%), respectively. Finally, air2strteam1 decreased the RMSE and the MAE of the DT1 by 31.60% and 34.59%, respectively.

Analysis of model performances at Rhône station (Table 7) showed that, for the models using only T_a as input variable, the three machine learning models were more accurate than the air2strteam model and FFNN1 (R = 0.927 and d = 0.962) and GPR1 (R = 0.926 and d = 0.957) performed the best and slightly better than the DT1 model (R = 0.923 and d = 0.956). According to the RMSE and MAE values, the air2strteam1 was the worst model with RMSE and MAE values equal to 1.044 °C and 0.812 °C, respectively. The FFNN1 gave the lowest values (RMSE = 1.044 °C and MAE = 0.812 °C) and improved the air2strteam1 by a 17.24% and 19.70% reduction in the RMSE and MAE values, respectively. Figure. 8 illustrates the scatterplot of measured and calculated values of T_w using the best models (version 3) at Rhône station.

A summary of the model performances at Dischmabach station is provided in Table 8 and Fig. 9. It can be observed that GPR1 slightly outperformed FFNN1 and DT1 according to all measures. A RMSE of 0.887 °C was observed in estimated T_w using GPR1, while the respective values for FFNN1 and DT1 were 0.896 °C and 0.972 °C, respectively. GPR1 produced the highest R and d between the measured and calculated T_w (R = 0.951 and d = 0.974), while the air2strteam1 yielded R of 0.945 and d of 0.972, less than the values obtained using the GRP1. To show the importance of including the DOY and Q as inputs to the models, the corresponding performances are illustrated and compared. Comparing the FFNN2, GPR2 and DT2 models with the DOY as input in addition to the T_a, it can be concluded that: (i) the three machine learning models provided the same accuracy with only marginal difference and (ii) the air2strteam2 yielded less accuracy compared to the FFNN2, GPR2 and DT2 models, with high RMSE and MAE, and low R and d values.

Results obtained using the applied models at Cedar station are illustrated in Table 9. Figure 10 illustrates the scatterplot of measured and calculated values of T_w using the best models. Overall, using only the T_a as input, the maximum R and d values of 0.967 and 0.981 respectively during the validation phase were obtained using the air2strteam1 model. In addition, the RMSE and MAE values of the FFNN1, GRP1 and DT1 models obtained were 1.211 °C and 0.933 °C, 1.248 °C and 0.978 °C, and 1.243 °C and 0.955 °C, which were respectively 13.54% and 10.075%, 16.10% and 14.21%, and 15.76% and 12.14%, greater than the values obtained using the air2strteam1 model. Thus, it was demonstrated that the air2strteam1 model was more accurate than the machine learning models. With the inclusion of the DOY as input variable combined with the T_a, the model performances varied accordingly and were significantly improved. The best accuracy was obtained with the FFNN2 model that had a significant decrease of the RMSE and MAE compared to the FFNN1 by 33.94% and 33.23%, respectively. Comparing the overall accuracy of the models using the version 2, it is clear from Table 9, that the four models provided the same accuracy with very marginal difference. Finally, as shown in Table 9, R and d values for the FFNN3, GPR3, DT3 and air2strteam3 reached up to 0.983, 0.985, 0.989, and 0.983, respectively, which were relatively higher compared to the values obtained using the versions 1 and 2.

Results at the Fanno station are reported at Fig. 11 and Table 10. The RMSE and MAE of the four models for this station were relatively low. It was indicated that the T_w calculated by air2strteam1 agreed well with the measured value, with R and d values of 0.984 and 0.989 during the validation phase. Compared to the machine learning models, air2strteam1 was the best model. Moreover, from the modeling accuracy reported in Table 10, the FFNN1, GPR1 and DT1 were less accurate than the air2strteam1. Analogously for the version 2 of models, air2strteam2 was the best model and yielded the highest R and d values (0.991 and 0.992), the lowest RMSE (0.961 °C) and MAE (0.794 °C), among the four developed models. On the other hand, except for the R values, the performances of the air2strteam2 was quite similar to FFNN2, with R, d, RMSE and MAE values being 0.985, 0.993, 0.964 °C and 0.768 °C, respectively. The DT2 model yielded the poorest accuracy with lowest R (0.982) and d (0.986) and highest RMSE (1.265 °C) and MAE (1.021 °C) values.

Table 11 presents the performance metrics for the three machine learning models and air2strteam for Irondequoit River. Figure 12 illustrates the scatterplot of measured and calculated values of T_w using the best models. Superiority of the air2strteam for predicting T_w is evident for all three versions of models. During the validation phase, the superiority of the air2strteam1 is indicated by higher R (0.984) versus 0.966, 0.967, and 0.965 for FFNN1, GPR1 and DT1, respectively. This superior performance is confirmed by lower RMSE and MAE values (1.160 versus 1.867 °C, 1.926 °C, and 1.976 °C), and (0.890 versus 1.459 °C, 1.430 °C, and 1.387 °C) for FFNN1, GPR1 and DT1, respectively. As is shown in Table 11, air2strteam2 performs best, as indicated by higher R and d, and lower RMSE and MAE values. It is important to highlight the significant improvement of the models by the inclusion of the DOY as input variable. By adding the DOY to the model inputs, the R value was increased from 0.989 to 0.993 for air2strteam2, from 0.965 to 0.985 for DT1, from 0.967 to 0.986 for GPR2, and 0.966 to 0.986 for FFNN2, respectively (Table 11). The RMSE was decreased by 35.99%, 34.58%, 34.51% and 19.05%, when using FFNN2, GPR2, DT2 and air2strteam2, respectively. Similarly, the MAE was decreased by 33.23%, 29.51%, 28.78% and 21.12%, when using FFNN2, GPR2, DT2 and air2strteam2, respectively.