Successive-Station Streamflow Prediction and Precipitation Uncertainty Analysis in the Zarrineh River Basin Using a Machine Learning Technique

Abstract

Precise forecasting of streamflow is crucial for the proper supervision of water resources. The purpose of the present investigation is to predict successive-station streamflow using the Gated Recurrent Unit (GRU) model and to quantify the impact of input information (i.e., precipitation) uncertainty on the GRU model’s prediction using the Generalized Likelihood Uncertainty Estimation (GLUE) computation. The Zarrineh River basin in Lake Urmia, Iran, was nominated as the case study due to the importance of the location and its significant contribution to the lake inflow. Four stations in the basin were considered to predict successive-station streamflow from upstream to downstream. The GRU model yielded highly accurate streamflow prediction in all stations. The future precipitation data generated under the Representative Concentration Pathway (RCP) scenarios were used to estimate the effect of precipitation input uncertainty on streamflow prediction. The p-factor (inside the uncertainty interval) and r-factor (width of the uncertainty interval) indices were used to evaluate the streamflow prediction uncertainty. GLUE predicted reliable uncertainty ranges for all the stations from 0.47 to 0.57 for the r-factor and 61.6% to 89.3% for the p-factor.

1. Introduction

Efficient management of water resources requires precise prediction of river streamflow to evaluate the effect of climate and land-use alterations, as well as increased agricultural irrigation, on regional aquatic systems [1]. Lake Urmia (LU), located in northwestern Iran, has significant socio-economic importance. However, the water level of the lake has decreased by up to 5 m over the past few decades due to the excessive use of available water and climate changes [2,3]. The changes in river discharges flowing into the lake are primarily responsible for the changes in its water level, since rivers contribute more to the inflow of the lake than groundwater and precipitation [4]. Zarrineh River (ZR) basin is the largest and most crucial sub-basin of the LU basin, providing more than 41% of the environmental flow into LU [3]. Therefore, developing a dependable model to predict the ZR streamflow is of great significance in assessing changes in LU’s water level. Nevertheless, the complex and non-linear behavior of the hydrological system’s components and insufficient data in the region hinder streamflow prediction [5,6].

In recent years, there has been widespread research on river streamflow prediction using process-driven and data-driven methods [7,8,9]. Process-driven methods are practical techniques for understanding fundamental mechanisms of hydrological phenomena, but they require a vast number of high-resolution inputs, including meteorological data, hydrological data, vegetation coverage, soil characteristics, and topographic data [10,11]. Therefore, low computational capacity and regions with unreliable and scarce input data often limit process-based model development.

Data-driven methods of machine learning techniques can efficiently convert the nonlinearity of the input–output relationship with no familiarity of the physical procedures [12,13,14,15]. There are several machine learning techniques that can be used for river flow prediction, including Artificial Neural Networks (ANNs) [16], Support Vector Machines (SVMs) [17], Random Forests (RFs) [18], Gaussian Process Regression (GPR) [19], and Long Short-Term Memory (LSTM) [20]. For instance, in [20] a Long Short-Term Memory (LSTM) is proposed based on model for streamflow forecasting in a river with multiple dams. Moreover, Xu et al. [19] proposed a hybrid model for river flow forecasting that combines Gaussian process regression with an improved differential evolution algorithm. Similarly, Liu et al. [18] evaluated the effectiveness of a Random Forest model for predicting daily streamflow. However, the high nonlinear relationship among the inputs and model output limits these data-driven models’ performance due to their simple structure [21,22].

Recently, deep learning (DL) techniques have been efficaciously utilized to address time-series forecasting problems [23,24,25]. DL approaches are capable of simulating more multidimensional purposes than non-deep neural networks by employing multiple neuron layers in a neural network structure [26]. The Gated Recurrent Unit (GRU) was proposed to simplify the structure of a Long Short-Term Memory model and solve the vanishing gradient problem in RNNs [27]. The GRU can produce better results by improving the prediction performance compared to other RNN networks by shortening the computation time. GRU networks have demonstrated a significant performance in dealing with nonlinearity and huge quantities of data with a simpler structure and higher computational speed than other variants of RNN [9,28]. The GRU has found many successful applications in the hydrology field, particularly in river streamflow prediction [29,30].

Streamflow prediction models face significant uncertainties due to insufficient data information and the complexity of the hydrological system. Uncertainties arise from non-optimal model features that are difficult to detect, systematic errors or measurement errors in initial data, and the calculation system due to simplification and assumption [11,31]. Several studies have evaluated and estimated the uncertainties in streamflow prediction [31,32,33,34,35] concluded that input uncertainty is an essential factor affecting the correctness of the streamflow estimation system. Since rainfall is the most essential input of the rainfall–runoff computation, its low spatial and temporal resolution or errors in the evaluation of precipitation data lead to significant uncertainty in the streamflow prediction [31,36].

Several studies have quantified uncertainty in the streamflow prediction models using various approaches [37,38,39]. Among these methods, the Generalized Likelihood Uncertainty Estimation (GLUE) is a cutting-edge method used to estimate uncertainty in prediction models [37]. Furthermore, the GLUE method is one of the most widely used approaches used to analyze uncertainty due to its simple concept, low vulnerability to model discontinuity, and easy implementation [11,31,33,40]. This method uses the Monte Carlo (MC) approach coupled with Bayesian estimation to determine the “behavioral’’ simulations based on the threshold value of the likelihood score.

This paper aims to predict the monthly streamflow of the Zarrineh River at successive stations, from upstream to downstream, by employing the GRU network. To achieve this goal, five model structures with different input variables and various time-lags were considered for each station. The selected input variables include precipitation, temperature, and streamflow with zero to four-month lag time. In addition, the GLUE method was used to quantify precipitation uncertainty in model prediction. However, instead of using MC simulation to produce random precipitation data series, the precipitation data were obtained from General Circulation Models (GCMs) under different Representative Concentration Pathways (RCPs) to avoid the stochastic errors caused by random data generation.

2. Study Area

Lake Urmia, which has a total area of 5750 km², is the largest lake in Iran and accounts for 7% of the country’s surface water [4]. Zarrineh River (ZR) basin is the largest and most crucial sub-basin of the LU basin, providing more than 41% (i.e., 1271 MCM) of the environmental flow into LU [3]. ZR is situated to the southeast of Lake Urmia and covers an area of about 12,025 km² with a length of around 300 km, as depicted in Figure 1 [41]. However, the lake surface area has drastically decreased to one-tenth, to 500 km, with the volume of half a billion cubic meters due to the unconventional use of available water and climate changes [2,3]. The Boukan Dam is the largest and most significant dam operating in the ZR basin with a live storage capacity of 650 MCM, storing water for drinking, agricultural, and industrial uses [42]. The average annual precipitation over the basin for the last four decades was 352 mm, which classifies the region as semi-arid with a Mediterranean climate.

3. Data Collection

In this study, four hydrological stations were considered for the streamflow prediction of successive stations in the Zarrineh River basin: Safakhaneh (station #1), Boukan reservoir (station #2), Qezkorpi (station #3), and Nezamabad (station #4) (see Figure 1B). Available measured monthly streamflow and reservoir outflow data for 1974–2014 were obtained from (www.yekom.com, accessed on 1 February 2023) and (www.dams.wrm.ir, accessed on 1 February 2023), respectively. In addition, the meteorological data introduced into the GRU network, including the precipitation, maximum temperature, and minimum temperature dataset for 1974–2014, were collected from Iran Meteorological Organization. A total of 17 well-spread meteorological stations were considered in the basin (Figure 1B) for reliable streamflow predictions.

The projected precipitation data from GCMs under different RCPs were used in the GLUE method to estimate precipitation input uncertainty. Future precipitation data from 2025 to 2060 was collected from all the available models and RCP scenarios in the region. In total, 93 precipitation datasets under RCP6.0, RCP2.6, RCP8.5, and RCP4.5 scenarios were obtained from the Climate Change, Agriculture and Food Security (CCAFS) data portal (www.ccafs-climate.org, accessed on 1 February 2023). The RCPs represent greenhouse gas (GHG) concentration trajectories used to understand the climate change in future, and vary from very low (RCP2.6) to very high (RCP8.5) future concentrations [43]. Because the GCM projected data contain systematic errors in their rough 3D resolution, they cannot directly be applied in climate models [44]. Hence, raw climate model outputs require bias correction to improve the fit of the projected data to the observations. Mengistu in [45] compared the raw regional climate model (RCMs) and bias-corrected RCMs against observed climate data. Bias-corrected RCMs performed better in reproducing rainfall, minimum temperature, and maximum temperature than raw RCMs, which demonstrated obvious biases in estimating climate data.

4. Model Description

4.1. Gated Recurrent Unit (GRU) Cell Structure

GRU is an advanced variant of a RNN developed to deal with the vanishing gradient problem [27]. Several studies showed that RNNs have higher performance compared to feedforward networks (FFNs) because they predict better and more stable streamflows [46,47]. Compared to other RNN networks, the GRU has a faster training process for multistep-ahead prediction without affecting its prediction performance. Thus, it is a commonly used deep learning technique, which has been utilized in many hydrological investigations, particularly streamflow forecasting.

The typical GRU cell structure is demonstrated in Figure 2. It has a memory ($h_t$), a candidate hidden layer ($h_t^{\prime }$), and two controlling gates: the reset gate ($r_t$) and the update gate ($z_t$). The memory of the current t and the previous time steps t − 1 is calculated using the reset and update gates. The update gate controls how much state information h_t−1 (h^′_t−1) is transferred to the up-to-date time step from the earlier one. More state data from the prior time step is produced by the greater number of update gates. The reset gate is applied to determine the degree to which the information from the previous state is forgotten. The lesser the reset gate, the more state information is forgotten. The update equations in the GRU cell structure are computed as per Equations (1)–(4):

$$r_t=\sigma \left(W_r{x}_t+U_r{h}_{t-1}\right)$$

(1)

$$z_t=\sigma \left(W_z{x}_t+U_z{h}_{t-1}\right)$$

(2)

$$h_t^{\prime }=\tan h(W_h{x}_t+r_t{U}_h{h}_{t-1})$$

(3)

$$h_t=\left(1-z_t\right)h_t^{\prime }+z_t{h}_{t-1}$$

(4)

where W and U are the networks’ weights matrices. The sigmoid function ($\sigma$) and the $\tan h$ function limits the output range from 0 to 1 and −1 to 1, respectively.

4.2. GLUE Theory

The GLUE is a statistical technique for uncertainty quantification of forecasting computation [37]. The GLUE method uses different variables to make numerous simulations in a model in order to describe the behavioral/non-behavioral models. The generalized likelihood function is used to identify the behavioral simulations. A higher likelihood value represents a better correlation between observed and simulated values. The behavioral models are used to quantify the model uncertainty after discarding the non-behavioral simulations. The term “behavioral” signifies the accepted models based on the available data and knowledge.

5. Methodology

5.1. GRU Model Development

Figure 3 demonstrates the GRU modeling steps undertaken in this study. The GRU network was used to predict successive-station streamflow in the Zarrineh River basin. Therefore, monthly streamflow, precipitation, reservoir outflow, and maximum and minimum temperature of 1979 to 2000 (427 data points), 2000 to 2003 (50 data points), and 2003 to 2014 (127 data points) were considered as training, validation, and testing sample data, respectively. The validation dataset was used to find the optimum model factors and avoid overfitting, and the testing dataset comprised the unseen data in the training procedure to evaluate the calculative algorithm efficiency.

The process of finding essential input variables with the most influence on the model’s output requires a trial-and-error procedure since there are no unified methods to determine them. Therefore, five model structures with different input variables and up to a four-month lag time were considered to find the structure with the best performance (see

Table 1. The model structures used in the GRU network.

Name	Model Structure
S1	$Q^t=f(Q_{us}^t,P^t,T_{max}^t,T_{min}^t,T_{avg}^t)$
S2	$Q^t=f\left(Q_{us}^{t-1},Q_{us}^t,Q^{t-1},P^{t-1},P^t,T_{max}^{t-1},T_{max}^t,T_{min}^{t-1},T_{min}^t,T_{avg}^{t-1},T_{avg}^t\right)$
S3	$Q^t=f(Q_{us}^{t-2},Q_{us}^{t-1},Q_{us}^t,Q^{t-2},Q^{t-1},P^{t-2},P^{t-1},P^t,T_{max}^{t-2},T_{max}^{t-1},T_{max}^t$ $T_{min}^{t-2},T_{min}^{t-1},T_{min}^t,T_{avg}^{t-2},T_{avg}^{t-1},T_{avg}^t)$
S4	$Q^t=f(Q_{us}^{t-3},Q_{us}^{t-2},Q_{us}^{t-1},Q_{us}^t,Q^{t-3},Q^{t-2},Q^{t-1},P^{t-3},P^{t-2},P^{t-1},P^t,T_{max}^{t-3}$ $T_{max}^{t-2},T_{max}^{t-1},T_{max}^t,T_{min}^{t-3},T_{min}^{t-2},T_{min}^{t-1},T_{min}^t,T_{avg}^{t-3},T_{avg}^{t-2},T_{avg}^{t-1},T_{avg}^t)$
S5	$Q^t=f(Q_{us}^{t-4},Q_{us}^{t-3},Q_{us}^{t-2},Q_{us}^{t-1},Q_{us}^t,Q^{t-4},Q^{t-3},Q^{t-2},Q^{t-1},P^{t-4}, P^{t-3},$ $P^{t-2},P^{t-1},P^t,T_{max}^{t-4},T_{max}^{t-3},T_{max}^{t-2},T_{max}^{t-1},T_{max}^t,T_{min}^{t-4},T_{min}^{t-3},T_{min}^{t-2},$ $T_{min}^{t-1},T_{min}^t,T_{avg}^{t-4},T_{avg}^{t-3},T_{avg}^{t-2},T_{avg}^{t-1},T_{avg}^t)$

). In Table 1, f represents the GRU networks; $Q^t$ and $Q_{us}^t$ represent the streamflow of the current and upstream station at month t, respectively; $P^t$ is the precipitation at month t; $T_{max}^t,T_{min}^t, \mathrm{and} T_{avg}^t$ are the maximum, minimum, and average temperature at month t, correspondingly; t − 1, t − 2, t − 3, and t − 4 illustrate one- to four-month lag times in the model structure. Structure S₁ and S₂ consider all the features, except the current station’s streamflow with zero- and one-month lag times, respectively. However, S₃, S₄, and S₅ structures contain all the input variables with two- to four-month lag times. A total of five model structures were used to simulate all hydrometric stations, except for the first station (Safakhaneh station), which does not include the streamflow from the upstream station.

Tuning the hyper-parameters in the GRU network is an essential step in achieving accurate prediction results [48]. However, there is no specific method available to select and optimize these parameters; therefore, the trial-and-error technique is used to discover the hyper-parameters with the best model performance on the validation dataset [49]. Therefore, a large number of experiments was performed by considering a wide range for each parameter.

Considering that the stochastic gradient descent optimization algorithm is applied to train the DL networks, a loss function is defined to repeatedly estimate the current model state. Then, the network’s weights are updated to increase the model performance on the subsequent evaluation. The present study used mean squared error (MSE) as the loss function MSE (Equation (11)):

$$MSE={\displaystyle \sum }_{t=1}^n{(Q_o-Q_s)}^2$$

(5)

where Q_o and Q_s are the observed and estimated streamflow at time t, correspondingly.

5.2. Data Normalization

Normalizing raw data is an important pre-processing step in training ML approaches. Mapping all the attribute data to the same scale avoids numerical difficulties of the model and enhances the speed and accuracy of the modeling. Zhu in [5] suggested normalizing data into the range of [0, 1] for ML techniques, specifically ANN networks. The following equation (Equation (6)) was applied in the present research for the data normalization:

$$X_{norm}=\frac{X_i-X_{max}}{Xmi{n}_{max}}$$

(6)

where X_i and X_norm denote the raw and normalized data, correspondingly. X_max and X_min represent the maximum and minimum of raw dataset, respectively.

5.3. Model Evaluation Criteria

The accuracy and reliability of streamflow prediction were evaluated using four statistical measures. The Nash Sutcliffe coefficient (NSE) (Equation (7)) is a reliable and widely used criterion for assessing the hydrological models’ performance, and indicates the ratio of the modeled data variance to the observed data variance. The range of NSE is [−∞, 1], with values closer to 1 indicating better performance [46]. The coefficient of determination (R²) (Equation (8)), which has a range of [0, 1], represents the linear relation between the observed and predicted data. The prediction model shows more reliable results if the value of R² is closer to 1. The root mean square error (RMSE) (Equation (9)) evaluates the magnitude of the difference between the observed and predicted values. The closer the value of RMSE to 0, the higher the accuracy of the prediction.

$$NSE=1-\frac{{{\displaystyle \sum }}_{t=1}{(Q_m-Q_o)}^2}{{{\displaystyle \sum }}_{t=1}{(Q_o-\overline{Q_o})}^2}$$

(7)

$$R^2=1-\frac{\sum {(Q_m-Q_o)}^2}{\sum Q_m^2}$$

(8)

$$RMSE=\sqrt{\frac{\sum {(Q_m-Q_o)}^2}{n}}$$

(9)

where n is the number of data points. $Q_m$ and $Q_o$ are the predicted and observed values, respectively. $\overline{Q_o}$ is the average value of the observations.

5.4. Bias Correction Method

The bias modification method enhances the reliability of climate model simulations by adjusting projected precipitation and temperature data to the observations [50]. Thus, simulated raw climate data are corrected based on the alterations in the mean and variability among the climate model outputs and observed data in a reference period. The general procedure of the bias correction approach is illustrated in Figure 4.

The general form of the bias modification method uses observations to correct the mean and temporal variability of the climate prediction technique outputs. This bias correction is performed by the following equation:

$$T_{BC}\left(t\right)=\overline{O_{REF}}+\frac{{\sigma }_{o,REF}}{{\sigma }_{T,REF}}\left(T_{RAW}\left(t\right)-\overline{T_{REF}}\right)$$

(10)

where T_BC is the bias-corrected GCM output, T_RAW is the raw GCM output for the historical or future period, T_REF is the GCM output from the historical reference period, and ${\sigma }_{T,REF}$ and ${\sigma }_{o,REF}$ are the standard deviation of GCM output and the standard deviation of reference observations from the reference period, respectively.

5.5. Quantification of Input Data Uncertainty Using GLUE

In the present paper, the general concept of the GLUE method used is shown in Figure 5. The first step in quantifying input data uncertainty is generating random sets of data. While previous studies have used MC simulation for this purpose [11,31,40], this study used projected precipitation data from GCMs under different RCPs for the period of 2025–2060. In total, 93 precipitation datasets were acquired from GCMs. Then, the likelihood value ($L(P|Q)$) was obtained after applying each dataset to the GRU network. The widespread likelihood quantity is defined as NSE equation (Equation (11)) [31,33]:

$$L(P|Q)=1-\frac{{{\displaystyle \sum }}_{t=1}{(Q_o-Q_S)}^2}{{{\displaystyle \sum }}_{t=1}{(Q_o-\overline{Q})}^2}$$

(11)

where $Q_o$ is the observed streamflow, $Q_S$ is the simulated streamflow, and $\overline{Q}$ is the average of observed streamflow dataset.

The comparison of the chosen threshold value $\alpha =80\%$ and likelihood value $L(P|Q)$ specifies the behavioral ($L\left(P|Q\right)\ge 80\%$) and non-behavioral ($L\left(P|Q\right)<80\%$) datasets. Then, the non-behavioral datasets are discarded and the behavioral ones determine the uncertainty interval using the greater (UL, Equation (12)) and lower boundaries (LL, Equation (13)) equations. Furthermore, the streamflow of the upper and lower limits is obtained using Equations (14) and (15).

$$UL=\frac{1+\alpha }{2}\times 100\%$$

(12)

$$LL=\frac{1-\alpha }{2}\times 100\%$$

(13)

$$Q_{UL}=Q_{max }+\frac{P_{UL}-P_{max}}{P_{min}-P_{max}} \left(Q_{P_{min}}-Q_{p_{max}}\right)$$

(14)

$$Q_{LL}=Q_{min }+\frac{P_{LL}-P_{max}}{P_{min}-P_{max}} \left(Q_{P_{min}}-Q_{p_{max}}\right)$$

(15)

where $Q_{UL}$ and $Q_{LL}$ are the upper and lower limits of the predicted streamflow, correspondingly; $Q_{min }$ and $Q_{max }$ are the minimal and maximal amounts of streamflow, respectively; $P_{min}$ and $P_{max}$ are the precipitation data consistent with $Q_{min }$ and $Q_{max }$, respectively; and, $P_{LL}$ and $P_{UL}$ are the precipitation data associated to the minor and higher boundaries likelihood values, correspondingly.

The $p\text{-}factor$ and $r\text{-}factor$ are applied to quantify the strength of the simulation and evaluate the predicted streamflow uncertainty. The $p\text{-}factor$ is the percentage of observed data in the uncertainty interval (95PPU) (Equation (16)). The $r\text{-}factor$ reflects the average width of the 95 PPU band (Equation (17)). Theoretically, the prediction is a perfect fit with the observed data if $p\text{-}factor$ and $r\text{-}factor$ are 1 and 0, respectively. A $p\text{-}factor$ larger than 50% depicts low uncertainty and a low value of $r\text{-}factor$ shows lower uncertainty in the model prediction.

$$p\text{-}factor=\frac{{{\displaystyle \sum }}_{t=1}^{n}l\left(Q_o\left(t\right)\right)}{n}$$

(16)

with $l\left(Q_o\left(t\right)\right)=\left\{\begin{array}{c}1 if Q_{LL}<Q_o<Q_{UL}\\ 0 otherwise\end{array}\right.$

$$r\text{-}factor=\frac{\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n\left(Q_{UL}-Q_{LL}\right)}{{\sigma }_o}$$

(17)

where $Q_o\left(t\right)$ is the experiential streamflow at time t and ${\sigma }_o$ is the standard deviation of the declared streamflow.

6. Results

6.1. Evaluation of GRU Networks

The hyper-parameters that require tuning include the optimizer, activation function, learning rate, number of epochs, and batch size. The epoch is a process of sending the entire dataset into the network only once to complete an iterative calculation. Each epoch contains large amounts of data; thus, they are split into small batches. The epoch and batch sizes are set to 64 and 1000, respectively. Although a significant number of epochs is selected, the callback is applied to stop the training process if the validation period performance starts to decrease. The optimization is limited in the stochastic gradient descent algorithm by using a similar learning rate for each feature. Furthermore, the Adam algorithm automatically adapts the learning rate by using the applied gradient for the variable. However, the algorithm may not locate the optima by using a small learning rate for each variable. Nevertheless, the Root Mean Squared Propagation (RMSprop), an extension of previous algorithms, uses the decaying moving average of partial gradients to focus on the most recently seen partial gradients and forget early gradients [51,52]. The RMSprop optimizer with a learning rate of 0.001 was selected. Moreover, the activation function was set to “Tanh”. Note that the GRU model with various structures uses the same hyper-parameters.

Numerous experiments were performed on the chosen range for hidden layers, neurons, and drop-out values for each structure of the GRU network in each hydrometric station. The hidden layers, neurons, and drop-out values varied between 1–5, 5–500, and 0.3–0.7, respectively. For instance, the outlet station of the Zarrineh River basin (Nezamabad station) has two hidden layers with 100 and 120 neurons in each layer, and a drop-out value of 0.4 in the S₃ model structure. The GRU network might reach sub-optimal solutions using a random start point. Therefore, ten identical runs were performed for each structure, and the final model was selected based on the replication with the best performance in the testing period.

Table 2. Performance of GRU-based streamflow forecasting models for five station structures in the Zarrineh River basin with varying monthly lag time.

		Training Phase			Validation Phase			Testing Phase
Station	Structure	NSE	R2	RMSE	NSE	R2	RMSE	NSE	R2	RMSE
Safakhaneh (#1)	S1	0.34	0.35	13.9	0.49	0.34	15.1	0.46	0.29	12.3
	S2	0.53	0.66	10.2	0.52	0.63	11.1	0.54	0.57	8.1
	S3	0.74	0.8	7.5	0.73	0.69	5.8	0.69	0.71	6.6
	S4	0.75	0.86	5.8	0.75	0.78	6.7	0.8	0.8	5.3
	S5	0.73	0.81	7.7	0.79	0.75	6.6	0.75	0.78	5.7
Boukan dam (#2)	S1	−7.9	0.79	23.4	−10.8	0.82	27.4	−12.6	0.85	20.7
	S2	0.79	0.89	35.1	0.83	0.88	31.5	0.75	0.92	28.1
	S3	0.8	0.81	35.5	0.73	0.76	27.7	0.78	0.81	25.9
	S4	0.84	0.84	31.3	0.78	0.88	28.5	0.81	0.83	24.2
	S5	0.88	0.89	26.8	0.88	0.89	23.5	0.85	0.86	20.7
Qezkorpi (#3)	S1	0.94	0.96	15.2	0.86	0.84	18.6	0.95	0.99	12.7
	S2	0.93	0.96	19.2	0.93	0.95	16.5	0.94	0.99	13.4
	S3	0.95	0.96	15.9	0.92	0.94	11.4	0.96	0.99	10.2
	S4	0.96	0.96	15.1	0.97	0.98	7.6	0.98	0.99	8.2
	S5	0.94	0.95	18.7	0.91	0.94	14.1	0.94	0.99	12.6
Nezamabad (#4)	S1	0.72	0.72	42.3	0.66	0.72	34.7	0.71	0.77	27.7
	S2	0.81	0.85	34.8	0.84	0.87	26.3	0.79	0.82	23.7
	S3	0.95	0.95	18.1	0.85	0.89	17.6	0.87	0.88	18.3
	S4	0.94	0.94	18.7	0.89	0.93	22.4	0.85	0.88	19.8
	S5	0.84	0.87	31.7	0.84	0.87	24.4	0.82	0.83	21.3

lists the results of the GRU network with five model structures for each hydrometric station. The best model structure was determined based on the statistical criteria of NSE, R², and RMSE in the validation and testing phases to obtain high and comparable performance and avoid model overfitting. The S₁ structure with no lag time shows the poorest performance among the other models in that station. However, introducing antecedent streamflow of the station and a one-month lag time of other input parameters in the S₂ model increases the model performance significantly compared to the S₁ structure. The S₂, S₃, S₄, and S₅ model structures have the same input variables with a one- to four-month lag time. All the available climate data with various lag times were considered in the model structures to obtain the best combination of these inputs and their period. In addition, lag times were chosen in order to analyze how temporal variations in inputs affect the results.

In the first station of the Zarrineh River basin, i.e., Safakhaneh, the models’ performance is enhanced with the increase in the lag time, except for the S₅ model, which showed lower results than the S4 model. This indicates that the model’s performance declines when complicating the model with excessive inputs. Overall, the S4 structure shows the best performance among the other models, with NSE, R², and RMSE of 0.75, 0.78, and 5.7, respectively, in the testing phase. While the streamflow of an upstream station is not considered in the model structures of this station, the downstream stations of the Safakhaneh benefit from the upstream streamflow. The monthly inflow to the Boukan dam was predicted using five structures, in which S1 presents inferior performance compared to the other models. However, applying the station’s streamflow and various lag times substantially improves the statistical criteria of the model. The S₅ structure with all the input variables and a four-month lag time shows the best output results, with NSE, R², and RMSE of 0.85, 0.86, and 20.7, respectively.

Considering that the Qezkorpi station is located downstream of the Boukan dam, the measured monthly outflow of the dam is used as the upstream outflow in the GRU network. All the structures depict high performance with comparable results, but with a slight improvement in the S₄ model composed of all the input variables and a three-month lag time. The evaluation criteria for the S₄ model are 0.98, 0.99, and 8.2 for NSE, R², and RMSE, respectively, demonstrating the most accurate model. The most critical station in the Zarrineh River basin is the outlet station, i.e., Nezamabad, which yields the outflow to Lake Urmia. The model illustrates satisfactory output results in all the structures with the highest model performance in the S₃ model, with NSE, R², and RMSE of 0.87, 0.88, and 18.3, respectively. Thus, the GRU network shows a significant capability to predict the successive-station monthly streamflow of the basin, particularly at the outlet station contributing to the Lake Urmia inflow.

The observed and predicted hydrograph and the scatter plot of the structure with the best performance for each hydrometric station in the training, validation, and testing phases are shown in Figure 6. The hydrographs show that the model accurately predicted the streamflow fluctuations in all the stations. In addition, the scatter plots illustrate that the streamflow is predicted with high R² in each station. Although the model shows some inconsistencies at high flows at the Safakhaneh station and Boukan dam, they performed reasonably for the low- and medium-range flows. The GRU model generally performed significantly for all the flows at the Qezkorpi and Nezamabad stations. Various climate data, the land use, and the location of stations are responsible for the inconsistency in the results for the same model structure in different stations. The results demonstrate that the model performed better for downstream stations compared to the upstream stations considering that the calibrated river flow reaches the downstream stations.

6.2. Uncertainty

The projected precipitation data from 93 GCMs were used to determine the uncertainty in the input data. The projected datasets were obtained from 2025 to 2060 and applied to the best GRU model of each station to predict ensemble streamflow. The 95 PPU plots derived from 93 precipitation datasets for each station are presented in Figure 7. The likelihood value $L\left(P|Q\right)$ was obtained by using the confidence level of $\alpha =80\%$. The results indicate that more than 92%, 96%, 98%, and 91% of the precipitation datasets satisfied the $L\left(P|Q\right)$ at the Safakhaneh, Boukan dam, Qezkorpi, and Nezamabad stations, respectively. These precipitation datasets are called behavioral and were retained to estimate the uncertainty of the input data in the GRU network (see

Table 3. Input uncertainty estimation of the GRU model for the Zarrineh River basin.

Station Names	$\mathit{L}\left(\mathit{P}|\mathit{Q}\right)$	p-Factor (%)	r-Factor
Safakhaneh	86	78.5	0.53
Boukan dam	89	89.3	0.57
Qezkorpi	91	86.6	0.52
Nezam Abad	85	61.6	0.47

). The calculated p-factor represents the number of observed streamflows falling inside the 95 PPU. The uncertainty is lower if the p-factor and r-factor are closer to 1 and zero, respectively [53]. The p-factor of all the stations is greater than 50%, showing low uncertainty in the retained datasets. The stations have a similar r-factors, but the p-factors of Boukan dam and Qezkorpi stations are higher, indicating lower uncertainty than the other stations.

7. Conclusions

The environmental studies in water related subjects are increasing in developing countries due to its vitality [54,55,56]. In Iran as a developing country, the alternation in the Zarrineh River streamflow is primarily responsible for the changes in the water level of Lake Urmia. This study uses a reliable machine learning method, i.e., GRU, to predict the successive-station monthly streamflow of the Zarrineh River basin. Through five model structures defined for each station, the structure with the most accurate results was obtained based on the three statistical criteria. The input variables in the model structures include streamflow of the current and upstream station, precipitation, and maximum, minimum, and average temperature, with a lag time of up to four months, excluding the Safakhaneh station which had no upstream streamflow. The GRU network presented significant performance in predicting streamflow, particularly at the basin’s outlet station, Nezamabad station. Furthermore, the GLUE method was applied to assess the effect of precipitation uncertainty in streamflow prediction. Therefore, ensemble streamflows were obtained by applying 93 GCM projected precipitation datasets to the GRU network. The results indicate the capability of this method to include the precipitation input uncertainty in the streamflow prediction. Most of the precipitation datasets satisfied the likelihood value considering the selected high confidence level. Furthermore, the p-factor and r-factor were used to estimate the input uncertainty by comparing the observed streamflow with the ensemble predicted streamflow. The p-factor of all the stations is greater than 50% and the r-factor is around 0.5, showing low uncertainty in the retained datasets.