Modelling daily water temperature from air temperature for the Missouri River

Linear regression model

For linear regression models, where one dependent variable exists, water temperature is modelled with linear functions. Linear regression models provide a first order estimation of the sensitivity of water temperature on air temperature (Webb & Nobilis, 1997; Kothandaraman, 1971). This simple model is shown in Eq. (1): (1)${\displaystyle T_w\left(t\right)=A+B{T}_a\left(t\right)+\epsilon \left(t\right)}$

where T_w(t) is water temperature for a given day, T_a(t) is air temperature for the same day as water temperature, A and B are regression parameters, and ε(t) is an error term.

Non-linear regression model

Whether water and air temperature relationship is linear is a key point for linear regression models, and this assumption has been questioned. For example, Mohseni, Stefan & Erickson (1998) found a non-linear relationship between water and air temperature, and they developed a non-linear regression model for water temperature prediction based on a logistic S-shaped function. The non-linear regression model proposed by Mohseni, Stefan & Erickson (1998) is shown in Eq. (2): (2)${\displaystyle T_w\left(t\right)=\mu +\frac{\alpha -\mu }{1+e^{\gamma \left(\beta -T_a\left(t\right)\right)}}}$

where α is a coefficient which estimates the highest water temperature, μ is a coefficient which estimates the minimum water temperature, β represents the air temperature at the inflexion point, γ represents the steepest slope of the logistic S-shaped function.

Stochastic model

Stochastic models are normally based on a stochastic time-series analysis function which can predict water temperatures using air temperatures (Caissie, El-Jabi & St-Hilaire, 1998). In stochastic models, water temperature is generally modeled as a function of time consisting of two entirely different components, a short-term residual component and a long-term component which is periodic in time. There are several forms of stochastic models. One form is used by Hadzima-Nyarko, Rabi & Šperac (2014): (3)${\displaystyle T_w\left(t\right)=a+bsin\left[\frac{2\pi }{365}\left(t+t_0\right)\right]+{\beta }_1{T}_a\left(t\right)+{\beta }_2{T}_a\left(t-1\right)+{\beta }_3{T}_a\left(t-2\right)}$

where a, b, t ₀, β₁, β ₂, β ₃ are regression coefficients.

Artificial neural networks

Artificial Neural Networks are often used to model complex nonlinear functions from observations, which is their biggest advantage. The most popular ANN is the multi-layer backpropagation network, with a basic structure generally consisting of three distinctive layers: the input layer, one or more hidden layers, and the output layer. In such a network, data travels forward through the layers. Data is introduced to the model through the input layer, then it is processed in the hidden layers, and finally output through the output layer. Every layer is made up of nodes or neurons which are linked to other neurons in the proceeding layer. With the exception of the neurons in the input layer, all the neurons in the other layers are made up of several components: an activation function, an offset or bias and weights. Common activation functions are logsig and tansig functions, which are sigmoid non-linear functions, and purelin is a linear function.

The training of the network is performed by comparing the desired or measured outputs with the network or model outputs. The difference between these values is then used to adjust the neuron weights. This process is repeated for a sufficiently large amount of training cycles until the desired network output accuracy is achieved.

The stability of ANN model strongly depends on the number of neurons in the hidden layers. A review on approaches to determine the suitable number of hidden neurons within a neural network was performed by Gnana Sheela & Deepa (2013), where they pointed out that the random selection of the number of hidden neurons may cause errors for the constructed models by either overfitting or underfitting. In this paper, the optimal number of neurons for the ANN model of each station was determined by performing 10-fold cross validation using the training (or calibration) dataset of each station. For a given station, the optimal number of neurons was determined by gradually increasing the number of neurons from 1 to 10 in the single hidden layer of the ANN model and calculating the mean cross validation error for each given ANN structure. The results of this procedure are displayed in Fig. 3. For the studied three stations, the optimal number of hidden neurons are respectively 3, 5 and 6.

Gaussian process regression

Gaussian Process Regression models are not new in hydrology (Grbić, Kurtagić & Slišković, 2013). GPR is a full Bayesian learning algorithm which has been used in diverse applications such as experiment design, multivariate regression, model approximation (Rasmussen & Williams, 2006; Girard et al., 2003; Quiñonero-Candela & Rasmussen, 2005). A process is referred to as a Gaussian processes if it is assumed that the joint probability distribution of model outputs is Gaussian. Compared to other machine learning methods, the advantage of GPR is that it combines various machine learning tasks, which include model training, uncertainty estimation and hyperparameter estimation. In GPR, it is assumed that the output variable measurements of y can be represented as (4)${\displaystyle y=f\left(\mathit{x}\left(k\right)\right)+\epsilon }$

where x represents measurements of input variables, ε represents noise with a Gaussian distribution and variance ${\sigma }_n^2$ and f represents the unknown nonlinear function to be modelled.

The space of functions has a prior probability which is defined as a Gaussian process having mean m(x) and covariance function cov(x, x′) (Rasmussen & Williams, 2006). For a given sample of input variables x_∗, the output variable is predicted using the predictive probability distribution p(y_∗—X, y, x_∗) with mean and variance: (5)${\displaystyle \begin{array}{c}{\hat{y}}_{\ast }=m\left({\mathit{x}}_{\ast }\right)+{\mathbf{k}}_{\ast }^{\mathrm{T}}{\left(\mathbf{K}+{\sigma }_n^2\mathbf{I}\right)}^{-1}\left(\mathit{y}-m\left({\mathit{x}}_{\ast }\right)\right)\hfill \\ {\sigma }_{y_{\ast }}^2=k_{\ast }+{\sigma }_n^2-{\mathbf{k}}_{\ast }^{\mathrm{T}}{\left(\mathbf{K}+{\sigma }_n^2\mathbf{I}\right)}^{-1}{\mathbf{k}}_{\ast }\hfill \end{array}}$

where I represents the identity matrix, k_∗ represents a vector which can be defined by [k_∗]_i= cov(x_i, x_∗), k_∗ = cov(x_∗, x_∗) and K represents a covariance matrix with elements [K]_i,j= cov(x_i, x_j). Thus, it can be seen that compared to classical regression methods where only the parameters are used to determine the model prediction or output, and the model output depends on the training dataset X, y.

For precise predictions, the available data is used to determine the parameters of the mean and covariance function. The predictive probability distribution is completely defined using these parameters which are often referred to as hyperparameters. The values of the hyperparameters can be obtained by maximizing the log-likelihood function of the training datasets (Rasmussen & Williams, 2006): (6)${\displaystyle logp\left(\mathit{y}|\mathit{X}\right)=-\frac{1}{2}{\mathit{y}}^{\mathrm{T}}{\left(\mathbf{K}+{\sigma }_n^2\mathbf{I}\right)}^{-1}\mathit{y}-\frac{1}{2}log\left(\left|\mathbf{K}+{\sigma }_n^2\mathbf{I}\right|\right)-\frac{n}{2}log\left(2\pi \right)}$

where n is the number of training datasets.

Bootstrap aggregated decision trees

In this paper, an ensemble of decision trees is also employed to model the daily river water temperature. A decision tree is a predictive modeling approach popular in machine learning. Specifically, for regression, it is referred to as regression tree. For regression trees, the tree structure is built by binary recursive portioning whereby the data is split into partitions or branches. All data in the training set are initially grouped into the same partition. Then the learning procedure assigns the datasets into the first two branches or partitions, by using every possible binary split of the inputs. The sum of the squared deviations from the mean in the two separate partitions is calculated and the split that minimizes this value is selected. This splitting procedure is then applied to each of the new partitions. The learning process continues until each node reaches a user-defined minimum node size and becomes a terminal node. An ensemble, which combines the predictions from multiple decision trees makes more accurate estimations than a single tree. Since decision trees are sensitive to the specific training data, their output has a high variance. In order to decrease the variance and increase accuracy, bootstrap aggregation or bagging is used (Breiman, 1996). Bagging creates several training datasets by using bootstrap sampling, i.e., by random sampling the original training set with replacement. The regression tree algorithm is applied to each dataset, and then the average amongst the models is used to calculate the estimations for the new data.