Simultaneous Prediction of Stream-Water Variables Using Multivariate Multi-Step Long Short-Term Memory Neural Network

1. Introduction

With the ever-increasing demand, surface water has become the most crucial resource among the communities [1,2,3,4]. From the scale of industrial and electricity generation to the scale of agricultural and drinking purposes, the availability of surface water needs to be abundant with the appropriate quantity and quality [5,6,7]. Among all the surface water bodies, Stream-Water (SW) is considered the most important source to provide countless benefits to human beings [8,9,10]. Conveying not only the drinking water to the communities and irrigation water for agricultural purposes, but streams can also significantly wash away the wastes, and provide habitat for wildlife and hydroelectricity. Often, it is used for several recreational purposes e.g., fishing, swimming, and boating [11,12]. The selected SW quantity variables are discharge and water level of the stream is highly influential on the overbank flooding in the surrounding area, demand of water supply, and fluvial ecology. Aquatic life is impacted significantly due to the temporal dynamics between the seasons of the stream discharge and water level [13,14,15,16,17]. Surface water quality standards are established by examining surface water quality measures such as dissolved oxygen, pH, turbidity, toxic substances, and aquatic macroinvertebrate life [18,19]. According to the New Jersey Department of Environmental Protection (NJDEP) (2012), Integrated Water Quality Report, at least one designated use is classified as Not Supporting (NS) in every sub watershed in Trenton [20]. Selected SW quality parameters in this study are temperature, Dissolved Oxygen (DO), pH, turbidity and Specific Conductance (SC). The temperature being the most important ecological factor is directly correlated to the physical, chemical, and biological properties of water [21,22,23,24]. DO is essential for aquatic life to survive, with differing oxygen concentration tolerances among species and life stages [25]. The pH and SC have a substantial impact on the other metrics of overall water quality, both constructively and adversely. According to the previous studies, the positive correlation between them and nitrate ions, ammonia, phosphorus, calcium, and magnesium, or even the detrimental influence of high pH on exotic species invasions, could induce disruptions in natural ecosystems [26,27,28]. Turbidity is the measure of relative clarity of water caused by suspended particles or dissolved whereas high values can significantly reduce the aesthetic quality of streams and influence natural migrations of species [29]. Assessment of turbidity improves the evaluation and indication of fecal contamination in water bodies such as Escherichia coli, the most common water infection [30,31].

Traditional physics-based numerical models (e.g., HEC-RAS, MIKE) involve spatial and temporal discretization for the entire computational space to compute SW variables which require high computational efforts [32,33]. Various numerical scheme (e.g., Finite Volume, Finite Element, Finite Difference) is used to solve partial differential equations i.e., Navier-Stokes equation coupled with the conservation of mass equation. The cost of spatial and temporal discretization increases exponentially with the increase in the required resolution and accuracy [34,35,36]. Input data for the physics-based river models consist of a significant amount of morphological, operational, and measured data. Data preprocessing for the physic-based models can be daunting depending on the spatial and temporal tags of the target variables. Physics-based numerical models require measurable and empirical parameters to estimate the target variables. A substantial amount of work to retrieve the parameters and constants through extensive laboratory-based experimentation and calibration is also a prerequisite which make these models computationally costly for practical implication with varying scale [37].

Data-informed predictive models provide an efficient alternative approach to forecast and monitor both the SW flow and quality parameters where only observed data can be used for the prediction instead of using many environmental factors required by the physics-based models. They offer reduced computational effort while simplifying complicated system and predict the outcomes using the observational only without any physics-based equations [38,39,40,41]. In recent days, Deep Learning (DL), an advanced sub-field of artificial intelligence, has become a popular choice in predictive modelling in the field of water resource management [42,43]. However, the traditional deep neural network algorithms (e.g., Multilayer Perceptron (MLP) do not have the ability to learn sequential data because they cannot store previous information, resulting in a constrained prediction capability for long-term time series, e.g., temporal distribution of the water table depth [44,45]. The MLP algorithms needs complex procedures in the data pre-processing stage to obtain good performance in predicting the target variables [46,47,48]. While the comprehensive data pre-processing can bolster the ability of a MLP model to learn the observed data, subjective user intervention is still necessary, e.g., selecting the number of reconstructed components [49]. In addition, the pre-processing require substantial amount of time as many reconstructed components needs to be calculated [50].

The Long Short-Term Memory (LSTM) is a special type of neural network which stores extended sequential data in the hidden memory cell for further processing [51,52]. LSTM performs well in processing long term sequential data, utilizing its sophisticated network structure specifically designed to carry the temporal linkage of the time series data. Water quality and quantity data have not been widely investigated in previous work employing LSTM. Therefore, compared with aforementioned MLP model, the proposed LSTM model only requires a simple data pre-processing method [53]. LSTM neural network is recurrent in nature, where the connections between units form a directed cycle allowing data to flow both forwards and backwards within the network. Therefore, the model is capable of preserving the past information and use them further for future prediction. LSTM model have already been used as a very advanced model in the field of DL, e.g., speech recognition, natural language processing, automatic image captioning and machine translation [44,54,55]. However, only a few studies have applied Recurrent Neural Networks (RNNs) or LSTMs to forecast multivariate time series data in the field of water resource [56,57,58]. The objective of this research is to untangle the pattern of the temporal distribution and linkage among the aforementioned SW variables and perform predictive analysis on the using the previous observed data. To accomplish the goal, a comprehensive Exploratory Data Analysis (EDA) is conducted to investigate the temporal dynamics of the SW variables and LSTM prediction is performed to predict the future values based on past records. Following sections of the paper demonstrate the study location, data source and collection, EDA, LSTM prediction, performance evaluation and possible future directions.

2. Data and Methods

2.1. Study Location and Workflow

The monitoring station used in this research is located along the Central Delaware River, in Trenton City, Mercer County, NJ [59]. It's positioned 450 feet upstream from Trenton's Calhoun Street Bridge, 0.5 miles upstream from Assunpink Creek, 0.9 miles north of Morrisville, PA. The Hydrologic Unit number for this station is 02040105 based on the USGS water resources and it is located at 40°13'18" N, 74°46'41" W coordinates referenced to the North American Datum of 1983 with the 6,780 mi² of drainage area (Figure 1).

Entire workflow of the EDA and LSTM prediction tasks are divided into three distinction stages. In the first step, data collection from the USGS web portal, an exploratory analysis of the SW variables and feature engineering to transform the data for training/testing the LSTM algorithm is conducted. Variables used in the analysis are listed in the Table 1. Activities in the first step are categorized as the transformer. After investigating the dataset and performing data transformation on the variables, LSTM neural network is trained using the data prepared in the first step to perform predictive analysis. LSTM neural networks regression model is assessed using several error matrices. These activities are categorized as estimator. The LSTM algorithm is tuned and optimized by altering the hyperparameters to reduce the errors in the prediction and achieve satisfactory performance. The third step namely evaluator, the model is deployed to predict the recession rate for a new set of target variables. Model performance is further improved through iterative incorporation and validation of the input variables.

The LSTM workflow of predicting the SW performance indicator is illustrated in the Figure 2. In the first step of the workflow in he time series of the SW variables are retrieved from the USGS National Water Information System: Web Interface [59]. The range of the time series data for all the variables was different due to the various recorded duration. Mean values of the SW variables are used in this study. The range of the data used in this research is from 02/25/2006 to 03/08/2022 with observed data of seven years. Historically, for the years 1898 to 1906, peak discharges were measured at Lambertville, NJ, 14.3 miles upstream of the Calhoun Street bridge. The maximum discharge was recorded on 8/20/1955 with the amount of 329,000 ft³/s and the minimum discharge was 1,180 ft³/s, on 10/31/1963. Extreme flooding occurred on 10/11/1903 when the water level reached an elevation of 28.5 ft above NGVD of 1929 which resulted in a discharge amount of 295,000 ft³/s [59].

2.2. Multivariate Exploratory Data Analysis (EDA)

In the second activity in the Figure 2, a detailed EDA is performed to perceive the attributes and characteristics of the multivariate dataset. EDA is the pivotal step to perform preliminary investigations on data in order to obtain satisfactory performance of the LSTM model. Internal temporal distribution of all the SW variables is explored through multiple visual techniques and numerical indices. EDA consists of process of performing initial investigation of the SW variables to understand the hidden pattern of the distribution of the variables. EDA is further divided into several activities. They are denoted as the descriptive statistics, outliers/extreme values detection, and normality check. Descriptive statistics provides a great approach to delineate the distribution of the values of the SW variables with the number of data points, mean, standard deviation, percentiles, interquartile range, and range. Full multivariate descriptive statistics is shown in the Table 2. To demonstrate the normality in the variables, histograms with density plot is used as a visual representation and Pearson’s Coefficient of Skewness (PCS) is use as a numerical indicator of skewness. Numerical imputation is performed to make the dataset consistent by filling missing values with the Nearest Neighbors (NN) of the datapoints [60].

Water quantity variables e.g., discharge and water level and water quality variable, turbidity show higher non-normality compared to the other SW variables. Numerical measure of non-normality/skewness, PCS values of discharge, water level and turbidity are also higher than the PCS values of other SW variables which indicates relatively less normality.

Figure 3 illustrates the linear connection between two SW variables. Low values of the linear coefficients delineate the overall non-linearity among several variables is high. The direction of the linear relationships is found to be both positive and negative. A few variables are approximately linearly correlated e.g., discharge and water level, temperature and DO, SC, and pH.

2.3. Feature Engineering (FE)

FE is performed after a successful preliminary investigation on the dataset using EDA. Without a successful FE, LSTM method may not yield to a satisfactory performance with minimum error. An adequate optimization through the iterative gradient descent cannot be reached without a successful scrutiny of the dataset. Therefore, a comprehensive feature engineering is carried out to transform the variables most suitable for the learning algorithm of LSTM [61,62]. FE in this research involves imputation, data transformation, data standardization and splitting the dataset into training, testing, and validation sets. Imputation in preformed to fill the null values so that the entire dataset becomes consistent. Null values or observation were found in every series due to sensor malfunctions. These cells in the dataset are imputed with the values of the NN of the blank cell. However, due to the reduction of the size of the dataset as a result of the exclusion of the observations, it is implemented in this study. After a successful imputation with median values of the variables, distribution of the variable series is checked visually and numerically to confirm the normality. PCS is used as an indicator of the normality of the variables. As the distribution of the values of discharge and water level is highly skewed to the left showing a significant non-normally, the neural network regression algorithms without appropriate data transformation does not contribute to satisfactory outcomes with good optimization [63,64,65]. Data transformation is performed to decrease the non-normality of discharge and water level. In this research, logarithmic transformation is implemented to transform the distribution of the features more to the normal distribution to reach the maximum performance of the LSTM model. In the Figure 4 distributions of the transformed and observed discharge and water level series can be seen. PCS values increases for all the transformed series compared to the observed series showing the increase in the normality.

The LSTM recurrent neural network model uses the gradient descent technique where the feature value affects the step size of the technique. Smooth progress towards minima in gradient descent requires the update of the steps at the same rate for all the feature values. To establish the LSTM model's training and testing dataset, all of the values are normalized utilizing Equation 1.

$$X_{norm}=\frac{X-X_{min}}{X_{max}-X_{min}}$$

(1)

X denotes the variable of interest and subscript norm, max and min represent the normalized variable, maximum and minimum value of the values of the variable. The entire normalized variable series is split into a training set and a testing set that is used to test/evaluate the model.

2.4. Long Short-term Memory (LSTM) Recurrent Neural

LSTM has become a very popular algorithm to deal with time series data in the DL forecasting where variables are dependent on the previous information along the series [67,68]. LSTM can capture the long-term dependencies and linkage among the variables [69]. High computational effort and time is needed by recurrent backpropagation to learn to store long-term information because of the decaying error backflow [70]. Hence, the concept of the vanishing gradient problem in recognizing long-term dependency of RNN was introduced [71]. LSTM feedback connections is the principal component of processing and recalling long-term information and a unique feature which makes it different from the traditional feedforward neural network.

Both long-term memory (c[t-1]) and short-term memory (h[t-1]) are processed in a typical LSTM algorithm through the utilization of multiple gates to filter the information. For an unchanged flow of gradients, forget and update gates update the memory cell state [72,73]. Three gates i.e., input gate i_g (pink), forgot gate f_g (red) and output gate o_g (violet) and cell state (green) handle the information flow by writing, deleting, preserving past information, and reading respectively (Figure 5). Hence, LSTM is capable in memorizing information at different lead times making it suitable in time series prediction within a certain interval [74]. In forget gate, long-term information enters and passes through a filtration where unnecessary information is discarded. The forget gate filter out unnecessary data by using the sigmoid activation function where the range of the function is 0 and 1 for open and close status, respectively. Input gate filter and quantify the significance of a new data coming as input to the cell. Like the forget fate, input gate filters out information by using binary activation functions and controls the flow of both long-term and short-term information. The output gates regulate the value of the upcoming hidden state which is a function of the information on previous inputs.

In this research, a neural network with one LSTM hidden unit accompanied by a dense layer connecting the LSTM target output at the last time-step (t-1) to a single output neuron with non-linear activation function. The LSTM model was trained using the DL library, Keras in Python, the ReLU activation function, and the coefficient of determination ($R^2$), and Nash Sutcliffe model efficiency coefficient (E) loss function. To predict the SW variable of a time-step in future e.g., daily/weekly, values of the variables at previous time-steps are used. Hyperparameters are tuned to maximize the performance of the LSTM model through iterative trial and error approach. In this study, Keras, python library which offers a space search for machine learning algorithms is used to find the best combination of the hyperparameters [75,76]. Hyperparameter of the LSTM algorithm considered in this study are the size of epoch and batch and number of neurons.

2.5. Model Evaluation and improvement

In the model evaluation step in the fifth activity in Figure 2, the performance of the LSTM model is evaluated using top three standard error matrices e.g., $R^2$, and the E. Error matrices provide numeric values as the model performance indicator by comparing the observed and predicted values. The $R^2$ value is used to evaluate the LSTM model in showing the model performance improvement. The highest $R^2$ score corresponds to the best predictive accuracy. In addition to the $R^2$ and E are used to illustrate the model response due to the variation in the lead time. The better the model fits the data, the closer the $R^2$ value is to 1. The third evaluation metric refers to the E and is one of the most widely used metrics for evaluating a hydrologic model's performance. $E$ can be classified as one of the scaled forecasts that compare the predicted error to the observed error [77,78]. Having a positive value for E, shows that the prediction and computed recession rate value is better than simply selecting the average observed value. In order to ensure that training is effectively fit, it is vital to choose a model's hyperparameters and consider carefully the controlling parameters that impact the learning rate of the model. By optimizing the size of the epoch, batch, and neurons in the stochastic process of the LSTM neural network, the LSTM algorithm's performance in predicting SW parameters is further improved.

3. Results and Discussion

LSTM neural network are used to predict the multivariate SW variables based on the previous time series data. Several lead time durations e.g., 6 hours, 12, hours, 1 day, 3 days, 1 week, 2 weeks, 3 weeks and 1 month are used to forecast the future values of variables simultaneously. Predicted values are compared to the observed dataset to quantify the error matrices. $R^2$, and E are used to estimate the error from the predicted SW variables. Model performance is improved by increase in the number of epochs. Error matrices are obtained through multiple models runs to demonstrate the linkage between the model performance and the lead times. Hyperparameters of the LSTM model are adjusted to optimize the model performance considering a set of batch size, epoch size and total number of neurons.

3.1. Predicted and Observed SW variables

Th output from the LSTM algorithm is compared to the observed values of the SW variables through visual representation in the Figure 6. Both the observed and predicted values of the SW variables are plotted for the entire time series against the number of observations. The overall distribution of the predicted values of SW variables is approximately identical to the observed data providing a satisfactory performance of the LSTM algorithm. The error metrices recorded for all variables and full time series shows LSTM performed well in case of both train/test set. After the LSTM model is trained, the entire observed dataset is divided into training and testing sets with the proportion of 80% and 20%. Training dataset is used to train the model and testing dataset is used to evaluate the model performance. In the Figure 6, orange portion of the plot illustrate training portion of the dataset whereas the green portion shows the testing portion. Dashed lines with blue color show the observed data.

3.1. Model Evaluation Matrices and Improvement

The performance of the LSTM neural network is evaluated using three error matrices e.g., Root Mean Square Error (RMSE), the $R^2$ and E. The performance of the model also evaluated and improved through increasing the number of iterations i.e., epoch in the neural network. The value of error matrices is observed with the increase in the number of epochs in the Figure 7. The model performance increases significantly from the very beginning of the iteration for both the train and test scenarios. The trend of change in the increase in the $R^2$ values reach a near-steady state after 20 epochs. Local decease in the performance i.e., increase in the RMSE value can be seen after 20 epochs. Further The performance variation with the change in the lead time regarding RMSE is represented in the Figure 7.

Error matrices e.g., $R^2$and E are documented for several lead times. Lead times are pivotal parameters of LSTM algorithm towards model performance. Lead time values are 1 day, 2 days, 3 days, 4 days, 5 days, 6 days and one week. The values of $R^2$and E decreases with the increase in lead time, showing the degradation in the model performance with increase in the lead times (Figure 8). Therefore, selection of the lead times should be based on the model performance and necessity.

Accuracies in the LSTM prediction for all SW parameters are presented in the Figure 9 with the help of coefficient of determination, $R^2$. Observed and predicted values of SW variables from the LSTM prediction are plotted to determine the $R^2$ value. The range of the $R^2$ value for all SW parameters 0.552 to 0.953 delineates satisfactory performance from LSTM model prediction overall. The best prediction with minimum error is found for the DO prediction with the $R^2$ value of 0.953.

3.3. Hyperparameters Optimization

To achieve the best possible model configuration for the predictions, hyperparameter optimization of the LSTM algorithm is of utmost importance. The optimized parameters that are only effective for the SW parameters are taken into account in this study because the neural network tuning procedure is stochastic. At first step, the epoch size is optimized keeping the batch of 4 and a single neuron. A set of increasing number of epochs (50, 75, 100, 125 and 150) is chosen to observe the response in LSTM model performance. In the same fashion, a set of batch sizes (1, 3 and 5) and neurons (1-5) are selected with a constant epoch size of 1000 to see the performance improvement. Epoch size of 2000 is selected to further the optimization process with the batch size and the number of neurons as it provides the highest $R^2$value. A batch size of 1 and 5 neurons provides further highest $R^2$values which can be seen in the Figure 10. Therefore, an epoch size of 2000, batch size of 3 and number of neurons of 5 is considered as the best hyperparameter combination for LSTM model prediction.

4. Conclusions

Multivariate prediction of the SW variables under both the water quantity and quality categories at a point location using the observed data can be highly beneficial to the water managers and decision makers to perceive the future flooding, irrigation works and fluvial ecology and aquatic life. Unlike the physics-based numerical models where additional terrain, meteorological data and human interventions are pre-requisite, the proposed approach relies only on the previously recorded data of the variables. The LSTM framework to predict the SW variables can be highly beneficial for the nearby community where the short-term prediction of the dynamics of the SW variables daily/weekly/monthly in future play a critical role. In the prediction of water quantity i.e., discharge and water level can substantially aid to prepare for the flood inundation, irrigation work and water supply demand. Prior knowledge of water quality of the SW can be highly beneficial to manage the aquatic life. As the model proposed uses only previous observed data of the variables, additional burden of input data and human intervention is not required for prediction work. Several approaches through physics-based numerical modelling techniques are proven inefficient in terms of real-time forecasting and computational efficiency. On the contrary, the application of the data-informed predictive models is highly efficacious in predicting various SW variables without taking complicated differential equations and assumptions into consideration. LSTM algorithm is capable of preserving both the short- and long-term pattern of the time series to forecast. Traditional physics-based numerical modelling tool requires assumptions, other correlated variables, and expensive calibration of the parameters.

This study contributes to a reproducible template to investigate the uniqueness of the temporal dynamics of SW variables through extensive EDA. Hidden pattern of the distribution of SW variables over seven years of data is discovered various up-to-date data exploration tools which is a mandatory requirement for the satisfactory training of LSTM algorithm. After a successful training step, LSTM is tuned and optimized through an explicit iterative performance record which can further be transferred to forecast SW variables in the identical geographical location. The performance of the LSTM algorithm in predicting the river discharge illustrates the algorithm is highly suitable to the discharge time series. Several error matrices show promising performance with minimum error. The proposed LSTM configuration is proved to offer satisfactory performance for the SW variables with the lead times up to one week. However, increasing the lead time increases the error in prediction limiting the performance of the LSTM model. Physics-based models are also incompetent in real-time prediction where the proposed LSTM can easily be coupled with the sensor and cloud to predict the SW variables in real time. Computational time may increase exponentially with the increase of the size of the dataset. Principle parameters obtained after the training process with minimum error are number of neurons, batch, and epoch size. The parameters optimized to obtain the best the LSTM configuration after training the model can be transferable in the similar climatic and geographic regions. For instance, if the distribution of the values of the SW variables are identical e.g., the difference among the PCS values being negligible, the parameters of trained LSTM model can be transferred and used for predictive analysis in a different location. However, we should not use our LSTM model in an area where the distribution of the features values through time is dissimilar. Future research should be conducted to incorporate high performance computing and cloud-based operations to obtain smart predictive tool to utilize the revolution in data storage capability and computational efficiency. LSTM models with different configurations should also be applied in different geographical and climatic locations to investigate the transferability of the model.