Deep learning-based algorithms for long-term prediction of chlorophyll-a in catchment streams

doi:10.21203/rs.3.rs-1643745/v1

An accurate estimation of harmful algal blooms is imperative for surface water protection. Site-dependent models have been developed using deep-learning (DL) methods for chlorophyll-a (Chl-a) simulations. These DL models have been developed and evaluated using data from one site and hence are challenging to apply to other sites. Moreover, several factors affecting Chl-a concentration in stream waters were not considered. Herein, we developed and compared the performance of site-independent DL models for Chl-a simulations in the Nakdong River, South Korea. The DL-based models we developed can incorporate irregularly measured water quality observations, static physical features, and climate data measured at constant time steps. Our results indicate that the attention-based explainable DL model exhibited the best performance, with an R² value of 0.85 at the training site and 0.52 at the test site, and can be used to develop regional water quality models.

Deterioration of surface water quality due to harmful algal blooms is a global concern ¹. Harmful algal blooms adversely affect aquatic life, with repercussions for the fishery industry. Surface water quality maintenance costs are increasing owing to the increasing frequency of occurrence of algal blooms and their widespread effects ². This increase in algal blooms requires urgent, timely, and reliable predictions to enable effective management interventions ¹. Chlorophyll-a (Chl-a) is considered a proxy of algal concentration in surface waters ³ and is used to estimate the total phytoplankton biomass as a parametric variable to quantify the trophic state of the water body ⁴. Thus, the accurate estimation of Chl-a by simulating its concentration can provide useful information regarding surface water quality.

The concentration of Chl-a in surface waters has been modeled using theory-based numerical models, such as enhanced stream water quality models (QUAL2E) ⁵, soil and water assessment tools (SWAT) ⁶, and CE-QUAL-W2 ⁷. These models estimate Chl-a concentrations using equations that reflect the underlying Chl-a dynamics. Zhang et al. ⁸ simulated Chl-a concentrations to understand the response of a lake to the changes in the climate and hydrological conditions. Kim et al. ⁹ modeled the eutrophication of a reservoir to evaluate the importance of point sources. However, process-based models are limited by their complexity, low predictive performance, large calibration time, and inability to incorporate the dynamic input features that affect the target variable ¹⁰.

Physical attributes like stream length, slope, width, the surface area of the sub-basin, and flow within sub-basins can significantly influence stream water quality and chemistry ¹¹. These attributes vary from one sub-basin to another within the catchment. The concentration of phosphorus and nitrogen, which determine water quality, can also vary within a catchment from one site to another site ¹². These site-dependent factors can be divided into time-variant data such as flow rate and time-invariant data like stream morphology. The water quantity and quality variables contribute to Chl-a concentrations in the stream. However, observations of water quality and quantity time-series data have not been conducted concurrently and continuously. Sampling and water quality experiments are time-consuming and expensive and, therefore, cannot be conducted at a frequency as high as flow rate observations. Park et al. ³ showed that, for high-performance Chl-a prediction models, it is pertinent to incorporate hydrological, geochemical, and ecological variables that affect algal growth in stream waters, thus requiring a modeling approach that considers data with different sampling frequencies.

With the rise in computational power and data availability in recent years, data driven machine-learning (ML) and deep-learning (DL) techniques are gaining attention. ML offers an alternative approach for developing high-performance prediction models for Chl-a simulation. Furthermore, DL-based models can effectively incorporate input data with different time steps ¹³. Several studies have been conducted to develop Chl-a prediction models using ML-based algorithms^14,15. A literature review shows that majority of the data-driven models developed for Chl-a prediction in surface waters have used classical ML algorithms, such as support vector machines, random forest, and XGBoost ^15,16. Recently, more studies have used advanced algorithms such as long short-term memory (LSTM), convolution neural network (CNN), and CNN combined with LSTM (CNN-LSTM) ^16,17. However, most ML-based studies have been conducted at monthly or weekly time steps ^3,18. Studies with daily or sub-daily time steps have focused on a limited simulation period, from one week to three years ^15,19. Moreover, these studies focused only on the development of site-specific DL models. These models were evaluated using data from the same area where they were trained. Using this approach, developing a regional model requires a separate neural network (NN) for each site, which is computationally expensive ²⁰. Thus, a site-invariant model that can be generalized to different regions for long-term Chl-a simulations is absent.

In this study, we propose a DL-based approach that considers different types of input data, such as static sub-basin characteristics, continuously measured climate observations, and discontinuously measured water quality observations (Fig. 1) (see Methods for details). This was accomplished using separate blocks of NNs to process continuous, discontinuous, and static data (Fig. 1a). Continuous and static data were processed using a single NN block. We compared the performance of six DL algorithms: LSTM, CNN, temporal convolution network (TCN), CNN-LSTM, LSTM-autoencoder, and input attention-based LSTM (IA-LSTM) in this block (Fig. 1b). Discontinuous water quality data were processed using a separate NN block (Fig. 1c). The outputs from these two NN blocks were concatenated and processed to predict daily Chl-a concentration (Fig. 1d). We used data spanning 16 years (2004-2020) from four sub-basins of the Nakdong River in South Korea. The sub-basins Hwangji (HG), Bonghwa (BH), and Dosan (DS) were used for model training and the Andong (AD) site was used for model validation (Fig. S1). To ascertain the impacts of static physical features, streamflow, and the attention mechanism on water quality, the performance of the LSTM model was further analyzed in different scenarios: (1) the use of static physical features to initialize hidden and cell states in LSTM (LSTM Cond), (2) removal of static physical features from input data (LSTM NoCond), (3) removal of streamflow data (LSTM NoFlow), and (4) use of the self-attention mechanism (LSTM SA).

Exploratory data analysis

The climate and hydrological time-series data, static physical features, and irregularly observed water quality data are shown in Figs. S2‒S4. These plots depict the seasonal and spatial patterns of observed variables during the simulation period at the study site. Static variables like channel depth and width increased from the HG to the AD sites (Fig. S4), indicating that the channel’s size increased downstream owing to the erosion caused by flushing ²¹. However, the channel slope decreased from the HG to the AD site, indicating that the lower course of the Nakdong River in the study area has a gentle slope (Fig. S4).

The pairwise Pearson correlation coefficients ²² between the input and output variables and among the input variables are shown in Figs. S5‒S8 as heatmaps. We found that the correlation between streamflow and total phosphorus concentration was strong at the BH and DS sites, whereas no such linear correlation was observed at the HG and AD sites. The concentration of Chl-a positively correlated with dissolved oxygen, nitric nitrogen, total dissolved nitrogen, and total nitrogen at the HG and AD sites (Fig. S5, Fig. S8). However, these correlations were negative at the BH and DS sites (Figs. S6‒S7). The positive correlations between water quality variables and streamflow indicated similarities in the hydro-chemical processes in the BH and DS sub-basins. Positive correlations were also observed between Chl-a and temperature, solar radiation, and daylight duration.

Hyperparameter optimization

The convergence plots of the hyperparameter optimization are shown in Fig. S9. While the x-axis indicates the number of iterations during Bayesian optimization, the y-axis shows the minimum validation loss obtained for each optimization iteration. The best convergence was obtained with the IA-LSTM model, whereas minimal convergence was obtained with the TCN model (Fig. S9a). Although the convergence patterns of the CNN, CNN-LSTM, and LSTM-autoencoder were different, these models converged at the same values. The LSTM, CNN-LSTM, and LSTM-autoencoder models exhibited a significant decrease in the validation loss at the start of the optimization; however, they showed no improvement with further optimization steps, which was depicted by the flat convergence plots after the first ten iterations. The flat lines suggest that the validation performance did not improve with further optimization iterations. The TCN, CNN, and IA-LSTM models exhibited a reduction in validation MSE in small steps. IA-LSTM showed a reduction in the validation MSE from 0.013 to 0.012 mg/L during the first five iterations and then from 0.012 to 0.009 mg/L between 5 to 30 iterations. Fig. S9b, which includes the convergence plots of the LSTM-based models for different scenarios, showed that all the LSTM-based models converged to a better MSE than the LSTM model without streamflow data. The simple LSTM model exhibited the best results and converged during the first five epochs. The other models showed a small but steady convergence with smaller improvements.

The training patterns of all six models during hyperparameter optimization are shown in Fig. S10. The loss curves for each iteration of hyperparameter optimization were plotted in light gray, whereas the loss curves of the best three models were highlighted in blue, green, and yellow for each model. These curves indicate a loss-reduction pattern during the training process. The training patterns of the LSTM-based models for the five scenarios are shown in Fig. S11. The training and validation loss curves for LSTM NoFlow are the shortest indicating that the model’s training was mostly aborted due to the ‘early stopping’ mechanism, owing to a lack of improvement in validation loss after the first ten training epochs.

We calculated hyperparameter importance using functional analysis of variance (fANOVA)²³ (Fig. S12) and partial dependence (PD) plots ²⁴ (Fig. S13-S18). The PD plots are graphical tools used to visually explain the importance of each hyperparameter ²⁵. In contrast, fANOVA quantifies the variance in the objective function explained by each hyperparameter in the form of a numerical value ²⁶. In PD plots, the steepness of the curves indicates the sensitivity of a particular hyperparameter ²⁴. The parameters with a gentle or horizontal slope had the least effect on the model prediction. The PD plots suggested that the parameter ‘lookback’ and learning rate were sensitive for LSTM, CNN-LSTM, LSTM-autoencoder, and IA-LSTM models (Fig. S13-S18). The CNN and TCN models were the most sensitive to activations in the dense layer (Fig. S14-S15). In general, Figs. S12 and S13-S18 demonstrated the agreement between the results of the PD plots and the fANOVA method for the quantification of hyperparameter importance. Hyperparameters with higher relative importance in the fANOVA had steeper slopes in the PD plots. The PD plots also indicated the optimized values of the hyperparameters at the convergence point of the optimization, which is indicated by the vertical red lines in the PD plots. For example, in Fig S14, the PD plots for the hyperparameter ‘lookback’ indicated that convergence was achieved when the lookback value was close to 5. The hyperparameter ‘lookback’ determined the number of historical time steps used by the model to make predictions ²⁷.

Chlorophyll-a simulation using deep-learning models

The simulation results of the six models for the prediction of Chl-a concentration at the validation site are shown in Fig. 2. The performance of each model for each training site is shown in Figs. S19‒S21. We observed that the DL models underestimated the peaks more frequently at the validation site than at the training site. The predicted Chl-a concentration in Fig. 2 was underestimated for peaks exceeding 20 mg/L. In contrast, the predicted Chl-a peaks followed the observed Chl-a at the training sites during peak events. At the HG site, the LSTM and CNN-LSTM models underestimated peaks higher than 10 mg/L, whereas the CNN and TCN models predicted these peaks accurately (Fig. S19). At BH, all the models predicted peaks in Chl-a during the entire simulation period, except for two peaks in 2009 and 2010 (Fig. S20). The ability of the DL models to simulate Chl-a was improved at the DS site. All models correctly simulated Chl-a during all peak events (Fig. S21).

Taylor diagrams for the six models for all sites in the study area are shown in Fig. 3. The plots revealed that the variability in the observed Chl-a was highest in AD, with a standard deviation of 6.4 mg/L; for the remaining sites, it varied between 4.4 and 4.9 mg/L. The magnitude of the variability was measured as the radial distance from the origin of the plot. A radial dashed line is shown from the reference point to aid in the interpretation of variability. The ‘reference’ point indicates the observed Chl-a concentration. Thus, the simulations from all the models at all sites have lower variability when compared to the observed values because they all occur below the dashed line. The top-right arc in the Taylor plots (Fig. 3) represents the correlation coefficients. Models closer to the x-axis have a higher correlation than those farther from the x-axis. Thus, the best-performing models in terms of correlation were CNN for HG and BH, whereas IA-LSTM had the best performance for DS and AD sites. The simple LSTM model had the lowest correlation coefficient for all the sites. The concentric gray lines from the ‘reference’ point indicated the root-mean-squared difference. These lines signified that the CNN model had the lowest root-mean-squared error (RMSE) for HG and BH, whereas LSTM had the highest RMSE.

The RMSE and coefficient of determination (R²) values for all models at the three training and validation sites are listed in Table 1. While IA-LSTM showed the best performance among all models in terms of R² at the DS and AD sites, the CNN model showed the best performance at the BH and HG sites with R² values of 0.87 and 0.83. The highest R² value (0.52) was obtained at the validation site by using the IA-LSTM model.

Table 1

Root mean square error (mg/L) and coefficient of determination for Chl-a prediction at the four sites in the study area. Hwangji, Bonghwa, and Dosan were used for training, while Andong was used for validation.
Model Name	Hwangji		Bonghwa		Dosan		Andong
Model Name	RMSE	R²	RMSE	R²	RMSE	R²	RMSE	R²
LSTM	6.49	0.24	6.33	0.49	3.24	0.73	4.86	0.42
CNN	3.09	0.87	3.71	0.83	2.51	0.84	4.75	0.49
CNN-LSTM	5.18	0.52	5.07	0.66	2.87	0.79	4.81	0.48
LSTM Autoencoder	4.13	0.70	4.71	0.74	2.82	0.81	4.75	0.45
TCN	3.66	0.74	5.44	0.64	2.77	0.81	4.91	0.41
IA - LSTM	3.36	0.53	3.13	0.78	1.76	0.85	4.49	0.52
LSTM based scenarios
LSTM	6.49	0.24	6.33	0.49	3.24	0.73	4.86	0.42
LSTM SA	5.12	0.36	5.65	0.59	3.08	0.75	4.71	0.44
LSTM NoFlow	7.42	0.08	8.26	0.28	4.55	0.56	5.70	0.21
LSTM NoCond	6.43	0.25	6.49	0.48	3.35	0.72	5.07	0.36
LSTM Cond	5.65	0.46	6.11	0.56	3.14	0.76	4.65	0.50

Scenario analysis of LSTM model

The simulation results of the LSTM model under different scenarios for the AD site are shown in Fig. 4. These scenarios include (1) the use of static physical features to initialize hidden and cell states in LSTM (LSTM Cond), (2) the removal of static physical features from input data (LSTM NoCond), (3) the removal of streamflow data (LSTM NoFlow), and (4) the use of the self-attention mechanism (LSTM SA). Static physical features include the slope area of the sub-basins and the volume, depth, slope, and area of the channel. The simulation plots for the HG, BH, and AD sites are shown in Figs. S22‒S24. These results indicated that “LSTM NoFlow” was unable to simulate peaks in the observed Chl-a. Moreover, the simulated Chl-a failed to follow the patterns of Chl-a during low-flow events. “LSTM Cond,” which used the static physical features to initialize hidden and cell states, showed a better correlation between observed and predicted Chl-a. Figure 5 shows the Taylor plots of the LSTM-based model for different scenarios. These plots indicated that “LSTM SA” and “LSTM Cond” showed the best performances when compared to the other scenarios. Overall, the best results were obtained at the DS site.

A comparison of the performance metrics for the Chl-a prediction under different scenarios is presented in Table 1. The table shows that the highest R² and lowest RMSE values were obtained for the DS site. The average R² value at the DS site was 0.7, which was the highest compared to the HG, BH, and AD sites. Similarly, the mean RMSE at the DS site was 3.47 mg/L. Overall, the best performance in terms of R² was obtained using the LSTMCond model at the DS site with 0.76 R². On the contrary, the best performance in terms of RMSE was obtained using LSTM SA at the DS site. At the validation site, the best performance was shown by LSTMCond, with an R² value of 0.5 and an RMSE value of 4.65 mg/L.

Explaining black box ML models

The attention weights of the ‘self-attention’ mechanism in the LSTM-based model are depicted in Fig. 6. Attention weights explain the prediction of DL models, which are usually considered black boxes ²⁸. These attention weights were plotted for the four sites in the study area along with lookback steps. The plots in Fig. 6 show higher attention weights for the input data from the previous day compared to the input data with lookback steps greater than one. This is evident from the higher attention weights at the lookback step of 1. The intensity of the attention weights gradually decreases with an increase in the number of lookback steps. Thus, we can infer which historical input data are more critical for the LSTM predictions ²⁹. The peaks in the simulated Chl-a were also reflected in the attention weights at all four sites. The attention maps indicated that the LSTM model focused on the input data from the previous three days to simulate Chl-a. However, the attention weights extended beyond the lookback step of three days during periods of higher Chl-a concentrations, indicating that at higher Chl-a concentrations, the LSTM requires more than three days of input data to simulate Chl-a concentration.

The attention weights indicating the importance of ammonia nitrogen from the IA-LSTM model for the four sites are shown in Fig. 7. These weights were generated for each input parameter used in the model. Thus, similar weights can be plotted for each input parameter processed using the IA-LSTM model. Figs. S23-S25 show attention maps for daylight hours, dissolved oxygen, and dissolved phosphorus, providing information on the importance of each input parameter for the prediction of Chl-a. These weights are plotted along with lookback steps. A higher weight value for the lookback step indicates the historical information that is more relevant for the prediction of Chl-a.

Model size and performance

A comparison of the number of epochs to train each model and the total number of trainable parameters is presented in Table S2. Each DL model was trained until overfitting was achieved. This was indicated by the deterioration of the validation performance. Thus, the number of training epochs for the model was different. We also monitored the size of each DL model, as indicated by the number of training parameters. A higher value for the training epochs suggested that the model required more training time. However, a smaller number of epochs indicated that the model can be trained quickly. We observed that the introduction of the attention mechanism led to faster convergence of the LSTM, thereby reducing the number of training epochs. This was indicated by the smaller number of epochs for the LSTM SA, which took 338 epochs, compared to the LSTM model, which took 491 epochs. We also observed that the LSTM-NoFlow model required 109 training epochs, which was the minimum among all models, showing that the early stopping mechanism stopped training during the early stages after 109 epochs. This is also evident from the learning curves (Fig. S11).

Exploratory data analysis

We observed strong seasonality in water quality and hydrological characteristics of the Nakdong River (Fig. S2-S4). Most precipitation occurs during the summer, from April to August, while the winter is mostly dry. The mean evapotranspiration during summer is 3.2 mm, while during winter, the average precipitation is 1.5. Similarly, the average daily streamflow during the summer was 150 m³ at the HG, BH, DS, and AD sites. Among the water quality parameters, the average value of DO during summer was 12.5, whereas, during winter, it was 7.5. This seasonality of the Nakdong River was also reported in previous studies ^{30 31}. The positive correlation between Chl-a and temperature, solar radiation, and daylight hours can be attributed to the role of light in the photosynthesis phenomenon ³². The exploratory data analysis of the time-series input data and static physical features indicates that significant differences exist in the dynamics of the flow and water quality constituents at the four sites. This is indicated by the average travel time of 0.2 at the AD site as compared to 5,7,9 at the HG, BH, and DS sites, respectively. These differences can be attributed to the location of the sites in the study area. The HG site is located upstream with no incoming tributaries or inflows. The BH and DS sites are more closely related to each other than the HG and AD sites because they have both incoming and outgoing streams. The presence of a reservoir and regulated outflow at the AD dam site resulted in changes in flow dynamics. Chl-a concentration at the dam is strongly related to the inflow and outflow from the dam, which is regulated by dam operation ³⁰. Thus, the simulation of Chl-a at these sites requires modeling of the underlying drivers that affect Chl-a concentration ²⁹.

Hyperparameter optimization

We observed a continuous reduction in training loss with an increase in the number of training epochs (Figs. S10-S11). However, this trend was not always observed with the validation loss. For the LSTM Cond and LSTM-SA models, the validation loss began to increase after 300 and 250 epochs, respectively (Fig. S11). The validation loss of the LSTM Cond increased from 0.01 to 0.025 mg/L. Similarly, the validation loss of the LSTM SA increased from 0.01 to 0.02 mg/L. This phenomenon in which the model starts memorizing instead of generalizing is called overfitting ²⁷. As a result, the training performance of the model increased at the cost of validation performance ³³. The training and validation loss curves for the LSTM NoFlow were the shortest (Fig. S11), which was indicated by 110 training epochs for this model. After 110 epochs, the training of the LSTM-NoFlow model was aborted due to the ‘early stopping’ mechanism owing to the lack of improvement in validation loss after the first ten training epochs. This indicates that streamflow data are influential in the prediction of Chl-a.

Chlorophyll-a simulation using deep-learning models

The performance of the models improved as we moved downstream from HG to DS. However, the performance deteriorated for AD, which was located farther downstream (Fig. 3). This change in performance can be attributed to the NN training method used in this study. In our training method, the three sub-basins were used sequentially. Therefore, the amount of data observed by the model increases as we move from the HG to the DS site. This increase in the data available to the DL model causes an improvement in its performance. This is indicated by the R² value of 0.73 from LSTM for the DS site as compared to 0.24 and 0.49 for HG and BH sites, respectively (Table 1). The performance of data-driven models has been reported to improve with an increase in the input data ³³. The R² was 0.52 at AD, lower than other sites because it was used as a validation site. This can be attributed to the different geological features of the AD site owing to the presence of a dam at this site ³⁴. It has already been reported that the validation performance of machine-learning models deteriorates when the validation data differ from the training data ³⁵. When comparing the models across the sites, the best results were obtained at the BH and DS sites. Table 1 shows that the R² of the simple LSTM model was the lowest for all sites compared to the other models. The IA-LSTM model performed the best at the AD site. This shows that the IA-LSTM model has better generalization ability than the other models.

The deterioration in the performance of DL models at validation sites compared to training sites can be due to the out-of-distribution (OOD) of validation data ³⁶. The OOD scenario refers to situations in which the distribution of the training and validation data vary significantly³⁵. DL models have limited capabilities when the validation data distribution differs from that of training data ³⁵. The inability of ML models to generalize OOD data is one of the main challenges and questions in active research ³⁷. This is because ML models assume that the training and validation data are identically and independently drawn from the same distribution ³³. This assumption states that both the training and validation data are independent and drawn from the same probability distribution. The differences in physical features and flow dynamics between the three training sites and the validation site can result in a violation of the identically and independently drawn assumption. Other studies have shown that physics-informed hybrid DL models can improve the prediction of streamflow on OOD data ³⁸.

The underestimation of Chl-a peaks (Fig. 2) using ML methods has been reported in previous studies ^39,40. Yajima and Derot ⁴⁰ explained the underestimation of Chl-a from the random forest model for periods of high peaks of Chl-a. A similar explanation is provided for this case. The Chl-a peaks at the validation site are higher (up to 50 mg/L) than Chl-a peaks at training sites (up to 30 mg/L). Figure 3 indicates that the CNN model showed the lowest RMSE at the HG and BH sites. Wunsch et al. ⁴¹, in a groundwater level forecasting study, also found higher prediction performance in terms of RMSE from CNN compared to LSTM for smaller datasets. For the DS and AD sites, the difference between the models in terms of RMSE was not significant.

Our results indicated that the performance of LSTM decreases when static physical features are not used as conditional variables. This is indicated by an R² value of 0.36 from LSTM NoCond, as compared to 0.42 of R² from LSTM (Table 1). The improved results from ‘LSTM Cond’ show that using static physical features of sub-basins to initialize hidden and cell state provides important contextual information to LSTM. The information regarding the physical features of a sub-basin improves the generalization ability of LSTM by increasing its prediction performance at the validation site. Several attempts have been made to incorporate static features into LSTM. Lees et al. ²⁹ used an entity-aware LSTM (EA-LSTM) to enable LSTM to learn hydrological patterns based on static catchment attributes. However, their results showed a decrease in model performance using EA-LSTM compared with simple LSTM. However, our results show that ‘LSTM Cond’ does not suffer from deterioration in performance upon incorporating static physical features. This shows that initializing the hidden and cell states of the LSTM with physical features instead of random values improves the learning ability. The performance of the LSTM Cond was better than that of the simple LSTM in terms of the RMSE and R² at all sites (Fig. 5).

Explaining back box ML models

We investigated the impact of the difference between the sites on attention maps. The patterns of attention weights from the BH and DS sites were more closely related than those from the AD and HG sites. At the BH and DS sites, the ammonia nitrogen concentration corresponding to a lookback step of one had a significant impact on the Chl-a simulation. At the AD site, Chl-a was more strongly affected by ammonia nitrogen during the peaks. Nutrients such as ammonia nitrogen in surface water increase phytoplankton biomass and productivity ⁴². Malek et al. ⁴³ also found that NNs were highly sensitive to ammonia nitrogen for Chl-a simulations. We also observed that the ammonia nitrogen from the previous three days played the most important role in simulating Chl-a at all sites.

The input attention weights for daylight hours, saturation concentration of dissolved oxygen, and total dissolved phosphorus are shown in Figs. S25‒S27. A comparison of the attention weights from the IA-LSTM model revealed the relative importance of the input variables ⁴⁴. Therefore, a comparison of Fig. 7 with Fig. S25‒S27 shows that ammonia nitrogen and the saturation concentration of dissolved oxygen play a more important role than other water quality variables in the simulation of Chl-a, which is evident from the higher attention weights for ammonia nitrogen (Fig. 7) and saturation concentration of DO compared to the lower weights for total dissolved phosphorus (Fig. S27). Several studies have highlighted that the concentration of phytoplankton increases with increasing oxygen levels at depth^45,46. This was due to the increase in oxygen produced by phytoplankton photosynthesis ⁴⁶. Similar results were obtained from process-based simulations of Chl-a in river waters. ⁴⁷. Thus, the incorporation of the attention mechanism in DL models makes them interpretable and explainable. The attention-based DL models showed that input data from the previous three days are the most important for the prediction of the daily Chl-a concentration in the Nakdong River.

Model size and performance

A higher number of training parameters indicates that the model is large, and training the model requires more computational resources ³³. The number of parameters of the CNN model was 28594, which makes it the smallest model when considering the size. However, the TCN model, with 190097 parameters, is the largest, owing to the stacking of the convolution layers within the TCN. (Fig. S28). The convolution layers in CNN and TCN have an advantage over LSTMs because convolution-based NNs can be easily parallelized ²⁷. This makes training faster for larger datasets using graphical processing units ²⁷. Therefore, a comparison of the training time taken by each model shows that the CNN model required the minimum training time among all models (Table S1). The reduction in training epochs for the LSTM SA model compared to LSTM indicates that the purpose of the attention mechanism is twofold. This helped improve the performance of the NNs and reduce the learning time.

Study area and data description

This study was conducted in the upper course of the Nakdong River, South Korea (Fig. S1). The Nakdong River is located in southeastern Korea, and is the longest river in Korea (530 km). The catchment at the study site consists of 12 sub-basins, with the Andong (AD) dam as its outlet. Several water quality observation points were located, starting from Hwangji (HG), upstream of the AD. In this study, data from the HG, Bonghwa (BH), Dosan (DS), and AD sites were used. Our input data consisted of three types of time-series data measured between 2005 and 2019. This included continuous, static, and discontinuous water quality data. Continuous data consisted of environmental data, including precipitation, minimum and maximum air temperatures, evapotranspiration, solar radiation, wind speed, daylight hours, streamflow in each sub-basin, and water temperature in the stream. We also used travel time obtained from a study that focused on simulating Chl-a using (SWAT) in the upper Nakdong River ³⁴. These data were measured using daily time steps. We considered static sub-basin characteristics such as the channel volume, width, length, depth, and slope. In addition, we considered the slopes and areas of the sub-basins. We also considered water quality observations, including ammonia nitrogen, total nitrogen, suspended solids, dissolved oxygen, biological oxygen demand, total phosphorus, and nitric nitrogen, for each of the four sub-basins. Water quality observations were performed continuously with variable time steps. Thus, these time series consist of irregularly measured observations from 2005 to 2019. Our target data were Chl-a concentrations, which were measured irregularly. The number of Chl-a observations for the HG, BH, DS, and AD sites were 175, 181, 382, and 733, respectively. Detailed statistical analyses of the input and target data are provided in Table S3‒S6 in the supplementary information.

Model architecture

We developed an NN architecture consisting of separate blocks to process continuous weather data, discontinuous water quality observations, and static physical features (Fig. 1a). The extracted features were further processed and refined to predict Chl-a concentrations. The complete steps of Chl-a prediction consisted of the following steps: (1) The continuous weather data, static sub-basin, and channel features were processed using the left NN block (Fig. 1b). We used six different algorithms in this block and compared their performance. All six algorithms were used for the time-series prediction problems ⁴⁸. (2) The second block comprised a dense layer to process the discontinuously observed water quality variables ²⁷ (Fig. 1c). Both networks used input data with daily time steps. (3) The third NN block consists of a concatenation layer and two dense layers. We used a concatenation layer to concatenate the features extracted from continuous and discontinuous input data. (4) These concatenated features are fed to a “Dense1 layer” whose purpose was to reduce the dimensions of the extracted features. (5) The features extracted from the Dense1 layer were further processed by an output layer that generated Chl-a concentrations in the stream (Fig. 1d). The output layer also consisted of a dense layer with only one unit.

The six DL architectures are LSTM, CNN, TCN, CNN-LSTM, LSTM-based autoencoder, and IA-LSTM. These architectures have shown promising results in hydrological and water quality modeling ⁴⁸. The LSTM architecture is a recurrent NN that can extract temporal features from time-series data ⁴⁹. The CNN architecture consists of convolution operations to extract local or spatial features from input data ²⁷. In this study, we used 1-dimensional CNN. A CNN-LSTM architecture combines CNN and LSTM to first extract features from the input data. LSTM is then used to extract temporal features ¹⁶. A TCN consists of hierarchical layers of a CNN to extract features in different hierarchies ⁵⁰. The LSTM-based autoencoder consists of two stacked layers of LSTM used to compress the input data and then reconstruct the data from its compressed form. The second layer, called the decoder, was used for prediction ⁵¹. In IA-LSTM, an attention mechanism first selects important input features that are further processed by the LSTM layer ⁴⁴. The details of these architectures are provided in Supplementary Information Text S1 and Table S6.

Model training

One of the major objectives of this study was to build a DL model that was trained and evaluated at different locations. Therefore, we used data from three upstream sub-basins (HG, BH, and DS) in consecutive order to train the NNs and data from the AD site for validation (Fig. 8). The input data were shuffled during the training process to avoid bias in the data ³³. The performance of this NN structure was evaluated on a separate sub-basin, the AD, located downstream. Thus, we trained the NN on three sub-basins located upstream and evaluated it on a separate sub-basin located downstream. We then evaluated the generalization ability of different NN architectures. The NNs were trained using the mean squared error (MSE) as the loss function (Text S2).

The NN parameters were optimized using an ADAM optimizer ⁵². During the training of the NNs, the MSE of the training data was used as the objective function. After each training iteration, we evaluated the performance of the NN on validation data. We used an ‘early stopping’ mechanism during training ³³. The training was stopped when the validation error stopped decreasing or started increasing. We quantified the similarities and differences between the observed Chl-a and the model predictions using a Taylor diagram ⁵³. Taylor diagrams provide a concise statistical summary of how well the predicted results from a model relate to observations in terms of their correlation, root-mean-square difference, and variance ⁵³. The hyperparameters of all the six models were optimized using the Bayesian optimization method ⁵⁴. The results of hyperparameter optimization were analyzed using functional analysis of variance (fANOVA) ²³ and PD plots ²⁴. Text S3 in the supplementary information provides a detailed description of the hyperparameter optimization method adopted in this study.

Incorporation of static data

The simple LSTM-based model uses continuous input data concatenated with static data as the input. However, static input data does not contain any temporal information because each static variable consists of a single constant value. Therefore, the use of static time-series data is inefficient. Hence, we adopted a strategy to initiate the hidden and cell states of LSTM with static data (Fig. 9). This approach conditions the simulations of the LSTM on the initial values of the hidden and cell states. A similar approach was adopted in other studies to incorporate static categorical input data for predictions using LSTM ⁵⁵.

It is pertinent to note that we initiated the hidden and cell states of the LSTM before training it with the data of a particular sub-basin ⁵⁵. Thus, we initiated the cell state and hidden state of the LSTM thrice during one training epoch, owing to three subbasins usage for training purposes. We call this LSTM ‘LSTM Cond.’ To test the importance of static input data, we built an LSTM model without static data. This model is called ‘LSTM NoCond.’ In LSTM NoCond, the hidden and cell states are initialized with zero values. Similarly, we built a model without using streamflow as the input data. This model is called ‘LSTM NoFlow.’ We also applied a ‘self-attention’ mechanism to the LSTM (LSTM SA) output. This attention mechanism aggregates temporal features using dynamically generated weights, allowing the network to directly focus on significant time steps in the past, even if they are far back in the lookback window ²⁸.

Data Availability

Flow data were obtained from the Water Resource Management Information System (www.wamis.go.kr). The weather data for the study site were obtained from the Korea Meteorological Administration (www.kma.go.kr). The static physical attribute data were obtained from a digital elevation model provided by the Rural Development Administration of Korea (www.rda.go.kr). The water quality and Chl-a concentration data were obtained from the Water Environment Information System of Korea.

Acknowledgments

This work was supported by the Korea Environment Industry & Technology Institute (KEITI) through the Aquatic Ecosystem Conservation Research Program, funded by the Korea Ministry of Environment (MOE) (2020003030003).

Author Contributions

AA was responsible for conceptualization, code development, and writing the draft. MJ were responsible for review and editing. SS was responsible for review, editing, and supervision. KHC was involved in conceptualization, funding acquisition, supervision, review, and editing.

Competing Interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Page, B. P., Kumar, A. & Mishra, D. R. A novel cross-satellite based assessment of the spatio-temporal development of a cyanobacterial harmful algal bloom. International journal of applied earth observation and geoinformation66, 69-81 (2018).
Paerl, H. W. & Huisman, J. Blooms like it hot. Science320, 57-58 (2008).
Park, Y., Cho, K. H., Park, J., Cha, S. M. & Kim, J. H. Development of early-warning protocol for predicting chlorophyll-a concentration using machine learning models in freshwater and estuarine reservoirs, Korea. Science of the Total Environment502, 31-41 (2015).
Lee, D. H., Kim, J. H., Park, M.-H., Stenstrom, M. K. & Kang, J.-H. Automatic calibration and improvements on an instream chlorophyll a simulation in the HSPF model. Ecological Modelling415, 108835 (2020).
Brown, L. C. & Barnwell, T. O. The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: Documentation and user model. (Environmental Research Laboratory, Office of Research and Development, US …, 1987).
Neitsch, S. L., Arnold, J. G., Kiniry, J. R. & Williams, J. R. Soil and water assessment tool theoretical documentation version 2009. (Texas Water Resources Institute, 2011).
Cole, T. M. & Buchak, E. M. CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model, Version 2.0. User Manual. (ARMY ENGINEER WATERWAYS EXPERIMENT STATION VICKSBURG MS ENVIRONMENTAL LAB, 1995).
Zhang, C., Brett, M. T., Brattebo, S. K. & Welch, E. B. How Well Does the Mechanistic Water Quality Model CE‐QUAL‐W2 Represent Biogeochemical Responses to Climatic and Hydrologic Forcing? Water Resources Research54, 6609-6624 (2018).
Kim, D., Kim, Y. & Kim, B. Simulation of eutrophication in a reservoir by CE-QUAL-W2 for the evaluation of the importance of point sources and summer monsoon. Lake and Reservoir Management35, 64-76 (2019).
Du, Z., Qin, M., Zhang, F. & Liu, R. Multistep-ahead forecasting of chlorophyll a using a wavelet nonlinear autoregressive network. Knowledge-Based Systems160, 61-70 (2018).
Fatehi, I., Amiri, B. J., Alizadeh, A. & Adamowski, J. Modeling the relationship between catchment attributes and in-stream water quality. Water Resources Management29, 5055-5072 (2015).
Umwali, E. D. et al. Spatio-seasonal variation of water quality influenced by land use and land cover in Lake Muhazi. Scientific Reports11, 1-16, doi:https://doi.org/10.1038/s41598-021-96633-9 (2021).
Jeong, K. et al. Prediction of biogas production in anaerobic co-digestion of organic wastes using deep learning models. Water Research205, 117697, doi:https://doi.org/10.1016/j.watres.2021.117697 (2021).
Cho, H., Choi, U. & Park, H. Deep learning application to time-series prediction of daily chlorophyll-a concentration. WIT Trans. Ecol. Environ215, 157-163 (2018).
Shin, Y. et al. Prediction of chlorophyll-a concentrations in the Nakdong River using machine learning methods. Water12, 1822 (2020).
Barzegar, R., Aalami, M. T. & Adamowski, J. Short-term water quality variable prediction using a hybrid CNN–LSTM deep learning model. Stochastic Environmental Research and Risk Assessment, 1-19 (2020).
Lee, S. & Lee, D. Improved prediction of harmful algal blooms in four Major South Korea’s Rivers using deep learning models. International journal of environmental research and public health15, 1322 (2018).
Kim, K. B., Jung, M.-K., Tsang, Y. F. & Kwon, H.-H. Stochastic modeling of chlorophyll-a for probabilistic assessment and monitoring of algae blooms in the Lower Nakdong River, South Korea. Journal of Hazardous Materials400, 123066, doi:https://doi.org/10.1016/j.jhazmat.2020.123066 (2020).
Barzegar, R., Aalami, M. T. & Adamowski, J. Coupling a hybrid CNN-LSTM deep learning model with a Boundary Corrected Maximal Overlap Discrete Wavelet Transform for multiscale Lake water level forecasting. Journal of Hydrology598, 126196, doi:https://doi.org/10.1016/j.jhydrol.2021.126196 (2021).
Ouyang, W. et al. Continental-scale streamflow modeling of basins with reservoirs: Towards a coherent deep-learning-based strategy. Journal of Hydrology599, 126455, doi:https://doi.org/10.1016/j.jhydrol.2021.126455 (2021).
Ji, U., Jang, E.-K. & Kim, G. Numerical modeling of sedimentation control scenarios in the approach channel of the Nakdong River Estuary Barrage, South Korea. International Journal of Sediment Research31, 257-263 (2016).
Pearson, K. (ed, 1895).
Hutter, F., Hoos, H. & Leyton-Brown, K. in International conference on machine learning. 754-762 (PMLR).
Friedman, J. H. Multivariate adaptive regression splines. The annals of statistics, 1-67 (1991).
Apley, D. W. & Zhu, J. Visualizing the effects of predictor variables in black box supervised learning models. Journal of the Royal Statistical Society: Series B (Statistical Methodology)82, 1059-1086, doi:https://doi.org/10.1111/rssb.12377 (2020).
Molnar, C. Interpretable machine learning. (Lulu. com, 2020).
Chollet, F. Deep learning with Python. 1 edn, 386 (Manning Publications Co., 2018).
Bahdanau, D., Cho, K. & Bengio, Y. Neural machine translation by jointly learning to align and translate. arXiv preprint arXiv:1409.0473 (2014).
Lees, T. et al. Benchmarking data-driven rainfall–runoff models in Great Britain: a comparison of long short-term memory (LSTM)-based models with four lumped conceptual models. Hydrology and Earth System Sciences25, 5517-5534, doi:https://doi.org/10.5194/hess-25-5517-2021 (2021).
Jeong, K.-S. et al. Impact of summer rainfall on the seasonal water quality variation (chlorophyll a) in the regulated Nakdong River. KSCE Journal of Civil Engineering15, 983-994 (2011).
Cha, Y., Soon Park, S., Won Lee, H. & Stow, C. A. AB ayesian hierarchical approach to model seasonal algal variability along an upstream to downstream river gradient. Water Resources Research52, 348-357 (2016).
Kindle, K. L. Expression of a gene for a light-harvesting chlorophyll a/b-binding protein in Chlamydomonas reinhardtii: effect of light and acetate. Plant molecular biology9, 547-563 (1987).
Goodfellow, I., Bengio, Y. & Courville, A. Deep learning. (MIT press, 2016).
Lee, J. et al. Dynamic calibration of phytoplankton blooms using the modified SWAT model. Journal of Cleaner Productionin revision (2021).
Berend, D. et al. in Proceedings of the 35th IEEE/ACM International Conference on Automated Software Engineering. 1041-1052.
Gallagher, L. K. et al. Sandtank-ML: An Educational Tool at the Interface of Hydrology and Machine Learning. Water13, 3328, doi:https://doi.org/10.3390/w13233328 (2021).
Bengio, Y., Lecun, Y. & Hinton, G. Deep learning for AI. Communications of the ACM64, 58-65, doi:https://doi.org/10.1145/3448250 (2021).
Lu, D., Konapala, G., Painter, S. L., Kao, S.-C. & Gangrade, S. Streamflow Simulation in Data-Scarce Basins Using Bayesian and Physics-Informed Machine Learning Models. Journal of Hydrometeorology22, 1421-1438, doi:https://doi.org/10.1175/JHM-D-20-0082.1 (2021).
Cao, Z. et al. A machine learning approach to estimate chlorophyll-a from Landsat-8 measurements in inland lakes. Remote Sensing of Environment248, 111974 (2020).
Yajima, H. & Derot, J. Application of the Random Forest model for chlorophyll-a forecasts in fresh and brackish water bodies in Japan, using multivariate long-term databases. Journal of Hydroinformatics20, 206-220 (2018).
Wunsch, A., Liesch, T. & Broda, S. Groundwater Level Forecasting with Artificial Neural Networks: A Comparison of LSTM, CNN and NARX. Hydrology and Earth System Sciences Discussions2020, 1-23, doi:https://doi.org/10.5194/hess-25-1671-2021 (2020).
O'brien, W. et al. Control mechanisms of arctic lake ecosystems: a limnocorral experiment. Hydrobiologia240, 143-188 (1992).
Malek, S., Ahmad, S. M. S., Singh, S. K. K., Milow, P. & Salleh, A. in BMC bioinformatics. 1-11 (Springer).
Jang, J. et al. Prediction of antibiotic-resistance genes occurrence at a recreational beach with deep learning models. Water Research196, 117001 (2021).
Stefan, H. G. & Fang, X. Dissolved oxygen model for regional lake analysis. Ecological modelling71, 37-68 (1994).
Zang, C. et al. Comparison of relationships between pH, dissolved oxygen and chlorophyll a for aquaculture and non-aquaculture waters. Water, Air, & Soil Pollution219, 157-174 (2011).
He, W. et al. Effects of temperature-control curtain on algae biomass and dissolved oxygen in a large stratified reservoir: Sanbanxi Reservoir case study. Journal of environmental management248, 109250, doi:https://doi.org/10.1016/j.jenvman.2019.07.021 (2019).
Lim, B. & Zohren, S. Time-series forecasting with deep learning: a survey. Philosophical Transactions of the Royal Society A379, 20200209, doi:https://doi.org/10.1098/rsta.2020.0209 (2021).
Hochreiter, S. & Schmidhuber, J. Long short-term memory. Neural computation9, 1735-1780 (1997).
Bai, S., Kolter, J. Z. & Koltun, V. An empirical evaluation of generic convolutional and recurrent networks for sequence modeling. arXiv preprint arXiv:1803.01271 (2018).
Liu, D., Jiang, W., Mu, L. & Wang, S. Streamflow prediction using deep learning neural network: case study of Yangtze River. IEEE Access8, 90069-90086 (2020).
Kingma, D. P. & Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014).
Taylor, K. E. Summarizing multiple aspects of model performance in a single diagram. Journal of Geophysical Research: Atmospheres106, 7183-7192 (2001).
Snoek, J., Larochelle, H. & Adams, R. P. Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems25 (2012).
Lim, B., Arık, S. Ö., Loeff, N. & Pfister, T. Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting, doi:https://doi.org/10.1016/j.ijforecast.2021.03.012 (2021).

(Not answered)

SI.docx
Supplementary Information

Deep learning-based algorithms for long-term prediction of chlorophyll-a in catchment streams

Status:

Version 1

Abstract

Figures

Introduction

Results

Exploratory data analysis

Hyperparameter optimization

Chlorophyll-a simulation using deep-learning models

Scenario analysis of LSTM model

Explaining black box ML models

Model size and performance

Discussion

Exploratory data analysis

Hyperparameter optimization

Chlorophyll-a simulation using deep-learning models

Explaining back box ML models

Model size and performance

Methods

Study area and data description

Model architecture

Model training

Incorporation of static data

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1