Rainfall-Runoff modeling using transductive long short-term memory and hyperparameter optimization

2.1 ML

In 2022, Xiao et al. ⁴⁹ employed four ML techniques at the Wuzhou station, located in the central section of the Xijiang River, to develop a highly accurate runoff forecasting model. The best model for forecasting runoff and water levels was the generalized regression neural network (GRNN). Singh et al. ⁵⁰ used diverse data-driven models to predict RR in the Gola watershed, including RF, multiple adaptive regression splines (MARS), and multiple linear regression (MLR). Several graphical plots were used in their thorough analysis to effectively present the results.

In 2023, Anaraki et al. ⁵¹ compared nine separate and hybrid ML models to simulate RR processes. These models included hybrids, including KNN enhanced with a gorilla troop optimizer, as well as more conventional models, including ANN and least squares SVM. Additionally, they presented a new feature selection technique that integrated empirical mode decomposition and principal component analysis (PCA). Shekar et al. ⁵² estimated monthly streamflow in the Murredu River basin using the SWAT model with eight AI models. They included KNN regression, SVR, and LSTM.

In 2024, four distinct ML models, including support vector machines (SVMs), gene expression programming (GEP), multilayer perceptrons, and multivariate adaptive regression splines (MARS), were used by Fuladipanah et al. ⁵³ to study RRM. They examined a dataset comprising 4765 daily records from the rainfall and hydrometric stations of the Malwathu Oya watershed, collected between July 18, 2005, and September 30, 2018. Using metrics like RMSE, mean absolute error (MAE), R2, and developed discrepancy ratio (DDR), the study assessed the models and concluded that the GEP model performed best. In Slovenia's karst Ljubljanica River catchment, Sezen and Šraj ⁵⁴ used three combined conceptual models, including a snow module with ML models for hourly RRM. They combined the CemaNeige Génie Rural à 4 paramètres Horaires (GR4H) conceptual model with a wavelet-based regression tree (WRT) and a wavelet-based extreme learning machine (WELM). The performance of these hybrid models was contrasted with that of independent conceptual and ML models. Shah et al. ⁵⁵ used ML techniques with models like the ANN multilayer perceptron (ANN-MLP), RF, and MLR to predict runoff in Quetta Valley, located in the Pishin Lora Basin. They utilized remotely sensed meteorological data from 1990 to 2022, gathered from the Modern-Era Retrospective Analysis for Research and Applications-2 (MERRA-2) satellite, and modeled runoff over daily, weekly, and monthly timescales. The inputs included rainfall, humidity, and temperature, with runoff data simulated using the soil conservation service curve number (SCS-CN) method due to the absence of actual records. Iamampai et al. ⁵⁶ devised a technique to refine RR model predictions by leveraging rainfall accumulation at various intervals alongside the soil water index (SWI). This method addressed the shortcomings of models that depend solely on rainfall by incorporating changes in soil moisture and runoff. They employed RF and ANN models to simulate daily runoff, noting that accumulated rainfall was a pivotal input. Remarkably, the RF model surpassed the ANN and conceptual models in accuracy. Using soil moisture data spanning 39 years, Kantharia et al. ⁵⁷ used an adaptive neuro-fuzzy inference system (ANFIS) to estimate daily discharge in the Damanganga basin.

In 2025, a novel framework that combines ML techniques with stochastic hydrologic modeling was proposed by Houénafa et al. ⁵⁸ They employed wavelet-enhanced models, such as wavelet-based gated recurrent units (WGRU) and wavelet-based extreme gradient boosting (WXGBoost). These models improved streamflow prediction by utilizing daily discharge variability and statistical features. A parameter regionalization technique for distributed RRM was presented by Sayama et al. ⁵⁹ The Bayes theorem and conditional probability were used to direct parameter assignment. This method decreased computational demand by connecting model parameters to regional hydrological features. To simulate monthly RR processes in watersheds in southern Thailand, Kaewthong et al. ⁶⁰ compared the performance of ML algorithms, including MLR, MLP, and SVM, to the Modèle du Génie Rural à 2 paramètres Mensuel (GR2 M) conceptual model.

The analysis of different ML approaches reveals a standard issue when these methods are used for RR. The methods depend heavily on adequate data quality and sparse data, while facing challenges with model overfitting during training and geographic generalization. Additionally, these models need substantial computing resources and expert knowledge to optimize their numerous parameters effectively.

2.2 DL

In 2023, to overcome the limitations of DL and conceptual models in hydrology, Kapoor et al. ⁴⁸ combined sophisticated DL architectures, particularly convolutional neural network (CNN) and LSTM networks, with the Modèle du Génie Rural à 4 paramètres Journalier (GR4J) RR model. The hybrid method unites conceptual model simplicity with deep learning prediction power for enhanced performance and simplified model design. Their research demonstrates that hybrid models deliver improved adaptability across different hydrological settings. The hybrid models demonstrate superior performance over both conceptual models and traditional deep neural networks (DNNs), specifically in arid catchment areas. Dang et al. ⁶¹ researched to enhance RR prediction accuracy through ML and DL model optimization methods. The researchers integrated four different optimizers into LSTM models, which included root mean square propagation (RMSprop) and Adagrad as well as Adadelta and Adam. The research tested different dropout strategies through 0%, 10%, 20% and 30% dropout rates in LSTM architectures across two Mekong Delta hydro-meteorological stations in Vietnam. The research aimed to establish the best optimizer together with the best dropout rate to achieve better model performance and reduce overfitting in hydrological forecasting.

In 2024, a novel DL architecture for hourly RR modeling was presented by Ishida et al. ² A sequential framework links LSTM networks to a one-dimensional CNN (1D-CNN) through this system. The LSTM system handles short-term data while the 1D-CNN processes extended hourly meteorological data. The system achieves better predictions through CNN-extracted feature integration. The Japanese watershed of the Ishikari River served as a testing site for this approach, which proved its capability to enhance hourly model forecasts. Li et al. ²⁸ developed a process-driven DL model through the combination of a conceptual hydrological model (EXP-HYDRO) with a recurrent neural network (RNN). The PRNN-EA-LSTM model integrates an entity-aware LSTM cell, which acts as a post-processor with EXP-HYDRO as a process driver inside an RNN cell. This specialized model structure unites the accurate predictions of DL methods with the fundamental understanding of process-based models. The model demonstrates enhanced robustness when energized by sub-process variables from EXP-HYDRO across 531 catchments. Li et al. ²⁷ developed a hybrid model that utilizes LSTM together with the transformer architecture and random search (RS) for efficient flood forecasting. The model accurately reproduces flood patterns through RR data analysis from the Yellow River Jingle watershed. By doing this, it gets around the drawbacks of data-driven models and traditional physical hydrology. The adaptability of AI for RRM in the Indian Bardha watershed was investigated by Shekar et al. ⁶² Six AI models, namely ANN, KNN, XGBoost, RF, CNN, and CNN-RNN, were applied in their study, originally published between 2003 and 2009. The calibration period spanned 2003–2007, and the validation period was 2008–2009. This approach tests the model's performance based on spatial distributions of rainfall, temperature, and discharge. Wang et al. ⁶³ introduced a new method of long-term runoff forecasting by using a transformer method and baseflow separation. Temperature Forecasting Using Transformer Networks deals with the problems of computational inefficiency and error accumulation after long sequences when using LSTM. To enhance the credibility of flood prediction, Kim et al. ⁶⁴ used an RRM (Rainfall-runoff model) plus AI to address the persistent issue of flood-related damages in Vietnam. They employed a genetic algorithm (GA) and pattern search (PS) to perform the calibration of the parameters, and they made use of the Tank model to compute the peak flow rate for the flood discharge. They used some of the evaluation criteria (WSSR and SSR) for measuring the peak flood discharge, and used a rating curve and Monte Carlo simulation as criteria to transform flood discharge into floodwater levels and establish confidence bounds. They utilized the measurements they made to improve the flood level predictions (by better informing the AI models in particular DNNs and LSTMs). Yoon et al. ⁶⁵ presented a semi-supervised learning method, called self-training, to augment data-driven models in the case of RR relationships, which often suffer from a lack of paired observations of climate data and streamflow. Essentially, they applied self-training in a semi-supervised, teacher-student model by using a teacher model trained on a small number of paired samples to generate pseudo streamflow for unpaired samples. These pseudo streamflow observations, combined with the real observations, provided the student model with data to train on. The approach they used has an annealing-capable loss function to deal with the uncertainties of the pseudo streamflow. Their results showed large improvements over older fully supervised models using an LSTM network and the CAMELS dataset and provided a feasible approach to RRM with augmented data. To improve RR predictions, Wang et al. ⁶⁶ developed a novel informer neural network in combination with the empirical wavelet transform (EWT). The new technique, considering only single events of rainfall, utilizes EWT for reducing the non-linearity and non-stationarity of runoff data. To reduce the impact of data variability, the model splits RR data into three sections using fractal theory as well. GPM precipitation data and 15 years of USGS runoff data were used to evaluate the method. It was superior to traditional LSTM-based methods. The effectiveness of ML-based RRM under changing environmental conditions was investigated by Li et al. ⁶⁷ They used models such as an LSTM network, a multi-hidden-layer back propagation neural network (MBP), and category boosting (CatBoost) for RRMs on a daily time scale for the Dongwan catchment of the Yiluo River Basin. These models were repeatedly and non-real-time tuned for non-stationary hydrological data. Their ability in holding variation and irregularities of hydrological time series was evaluated utilizing determinants including deterministic coefficient, peak flow error, and runoff depth error.

In 2025, six AI models for RR analysis, including SVR, MLR, XGBoost, LSTM, CNN, and a hybrid CNN-RNN, were assessed by Shekar et al. ⁶⁸ during training and validation periods spanning 1998 to 2006. Jiang et al. ⁶⁹ introduced a model for RRM, which utilizes a convolutional neural network architecture to analyze frequency-domain transformed hydrological time series, capturing local and global data variations. Li et al. ⁷⁰ developed a model, integrating LSTM with Fourier Kolmogorov Arnold networks and convolutional feature extraction to enhance spatial data handling and improve precipitation-runoff modeling accuracy. Zhang et al. ⁷¹ developed a hybrid explainable streamflow forecasting model using CNN-LSTM-attention, focusing on typical river source regions in the eastern Qinghai Tibet Plateau, incorporating base flow to enhance high-flow predictions. Yin et al. ⁷² introduced a multi–step regional RRM approach using a pyramidal transformer (PT) with hierarchical attention (PTHA). The model integrates dynamic and static hydrological attributes to capture spatiotemporal dependencies more effectively than conventional LSTM models. Using DL techniques, Bao et al. ⁷³ created a runoff forecasting model that improved feature extraction and captured temporal dependencies in hydrological time-series data by combining transformer architectures and LSTM networks with sparse and full attention mechanisms. Zhu et al. ⁷⁴ used a hybrid DL technique that combines transformer architectures with self-attention mechanisms and CNNs (ResNet). This design allows for precise runoff prediction in large river basins by capturing intricate spatiotemporal patterns in multi-source hydrological data. To capture nonlinear hydrological patterns with enhanced spatiotemporal feature learning, Xu et al. ⁷⁵ developed a runoff prediction model that combines wavelet feature extraction, a transformer encoder for sequence representation, and a TimeBlock decoding module. To improve feature extraction over conventional LSTM and transformer approaches, Yin et al. ⁷⁶ developed a temporal-periodic transformer (TPT) model that combines temporal and periodic attention mechanisms to capture time-dependent and annually recurring patterns in streamflow. A lightweight attention-based modern CNN (LMCNN) was proposed by Jiang et al. ⁷⁷ To effectively capture multi-channel temporal and cross-variable runoff features, it combines depthwise convolution (DWConv), pointwise and squeeze-and-excitation (SE)-enhanced convolution (PSConv), and improved SE networks (SENets).

Previous studies on RRM have often emphasized the predictive ability of LSTM models. Nonetheless, many LSTM models train globally on the entire dataset. This makes it more difficult for them to record regional differences in hydrological dynamics. In certain situations, this lowers forecast accuracy because the model may overlook minute but significant variations in rainfall and runoff. We suggest TLSTM as a solution to this gap. This transductive approach uses the proximity of incoming data points to adjust weights. By incorporating local adaptability and retaining the advantages of LSTM, the method enhances the model's ability to capture short-term fluctuations. Therefore, by addressing the limitations of traditional LSTM models in representing fine-grained hydrological responses, the study makes a significant contribution.

The lack of effective and methodical hyperparameter optimization for deep learning models in RRM represents another gap in the literature. Previous research employed grid search or trial-and-error methods. These approaches often yield suboptimal results and are computationally expensive. We employ the HMS-based DE algorithm to solve this. It works well for locating balanced parameter configurations and exploring high-dimensional search spaces. In doing so, our study enhances prediction accuracy and stability while lowering computational load. To summarize, this study makes two contributions: (i) it introduces TLSTM to overcome the limitations of global LSTM training and capture localized hydrological patterns; and (ii) it uses the HMS-based DE algorithm for hyperparameter optimization to increase rainfall–runoff prediction accuracy, efficiency, and robustness.

In this work, we assess the suggested model using the CAMELS dataset. CAMELS ⁷⁸ is a standard in RRM benchmarks due to its wide scope and the variety of catchment characteristics it encompasses. For thorough model validation, this dataset offers comprehensive meteorological and runoff data. The method places a strong emphasis on detailed and varied data to ensure comprehensive analyses tailored to the distinct features of each dataset. However, it remains flexible and applicable to various hydrological datasets.

CAMELS contains 11,981 records from 674 basins across the United States. This enormous dataset, which provides daily runoff records, catchment characteristics, and meteorological forcings for every basin from October 1, 1980, to December 31, 2014, facilitates thorough model evaluation. While regional models use data from the Maurer dataset, individual basin RRM uses daily meteorological data from Daymet. These meteorological data sources were selected by accepted guidelines for RRM. These two datasets are commonly used in both individual and regional models. A total of 32 features can be extracted from each record in the dataset. For individual basin modeling, we replicate the benchmarks by selecting the same 673 basins. Similarly, for regional modeling, we utilize the same 448 basins used in benchmark studies. To maintain a fair comparison, our methodology employs the same daily meteorological forcings and catchment characteristics as those used in the benchmarks.

For building and evaluating the model, the dataset was partitioned into a 70:30 ratio for training and testing. This ratio aligns with standard practice in hydrology and machine learning. In these fields, allocating 20–30% of data for testing is recommended to balance model learning and evaluation reliability. ⁷⁹ Empirical studies, including recent streamflow modeling with machine learning, show that 70:30 splits typically provide robust generalization while maintaining sufficient unseen data for unbiased assessment. ⁷⁹ In this setting, the test set consisted of the last 30% of temporal sequences in each catchment, separated chronologically to avoid data leakage. During training, only the input features of these test sequences without target runoff values were used to adjust the model in the transductive step. The actual runoff values were reserved exclusively for final evaluation. An 80:20 split provides more training examples. However, it can reduce the statistical diversity of the validation set and increase the variance of performance estimates. Using less than 20% test data can also inflate performance metrics due to overfitting. Therefore, considering the dataset size, hydrological variability, and model complexity, the 70:30 split provided the most reliable balance for both training and evaluation.

3.2 Data preprocessing

The pre-processing stage focuses on preparing raw hydrological and meteorological data to ensure the model receives consistent, complete, and high-quality inputs for RR prediction. First, the collected rainfall, runoff, and catchment attribute datasets are systematically inspected to detect missing or anomalous records that may result from sensor malfunctions, transmission errors, or reporting delays.

The nearest-neighbor imputation method is used to fill in the missing values. The approach substitutes missing values by obtaining the nearest valid observation value to keep time series continuity. Statistical thresholds that use z-scores and interquartile ranges (IQR) help detect unrealistic measurements, including negative rainfall values and extreme runoff spikes that do not match precipitation events. The procedure removes isolated anomalies but substitutes continuous anomalies with interpolated values. The multi-step data cleaning procedure produces a final dataset that will be dependable and suitable for precise RRM.

The next step after data cleaning involves min-max scaling to bring all input features into uniform numerical ranges between 0 and 1. The training process becomes more stable because this procedure stops high-magnitude variables like runoff from overpowering smaller-scale variables such as rainfall intensity during the learning process.

A sliding-window method is used to transform the data into time-series sequences. The sliding window method produces overlapping sequences of input and output pairs through a fixed-length window that slides across the temporal data. This enables the model to learn both long-term dependencies (from accumulated hydrological memory) and short-term patterns (from recent rainfall and runoff). As inputs to the model, the pre-processed data are subsequently organized into feature matrices that integrate spatial attributes and temporal sequences.

3.3 Hyperparameter optimization

An essential component of DL is hyperparameter optimization, which has a big influence on the efficacy and efficiency of training algorithms. Hyperparameters, such as learning rate, batch size, and number of epochs, regulate the training process of the model. The ability of the model to generalize from training data to unknown data is also influenced by these parameters. The overall performance of the model, computation load, and rate of convergence can all be significantly impacted by the selection of these parameters.

Table 1 lists the hyperparameters that were modified during our investigation; the value ranges were selected based on prior DL applications in RR. These ranges provide a framework for understanding how the DE algorithm works, enabling modifications and enhancements tailored to specific project needs. To create the initial population for the DE optimization process, random initial hyperparameter values were generated within the specified ranges (defined in Table 1) using the Random Key encoding.

3.3.1 Random key

The Random Key technique is used in this paper to optimize hyperparameters. It provides ease of integration with evolutionary algorithms, simplicity, and efficiency in continuous search spaces. It facilitates a thorough investigation of various model configurations. Additionally, it facilitates parallelization, which reduces computation costs and time. It also enables global searches, which reduces the likelihood of local minima. These properties make it highly adaptable for tuning deep–learning architectures with multiple hyperparameters. This adaptability ultimately improves overall model performance. ⁸⁰

Random Key employs a method that translates solutions into a series of T numeric vectors, each having D dimensions, collectively termed p 1 , \; p 2 , \; … , \; p T . A transformation function, known as the “random key”, connects each vector representing a potential solution to the hyperparameters of the model. For optimizing C hyperparameters (detailed in Table 1), each hyperparameter c = 1 , \; … , \; C aligns with D c positions within the vector. For hyperparameters that accept continuous values, D c is set to one, making the aggregate dimensionality D equal to the sum D = ∑ c = 1 C D c . Each vector, p i , is segmented into C sections, each holding D c positions that signify potential values for a hyperparameter. For categorical hyperparameters, the process converts the D c dimensional segment of the numeric vector into a categorical value using a mapping, M A P c . In this mapping, segment elements are ranked by value, and the top element determines the categorical choice from the M A P c vector for the c t h hyperparameter. This approach enables the direct application of evolutionary operations, such as mutation, crossover, and selection, to the numeric vectors p i . This allows for an efficient conversion of results into practical sets of categorical and continuous hyperparameter values.

An example of the Random Key approach for the hyperparameter ‘number of layers’ where D c =5 is depicted in Figure 2. Essentially, the random key comprises real numbers ordered by value, a crucial mechanism for aligning vectors with predetermined arrays of options. This structure provides a clear view of feature importance by progressively elevating important features while decreasing the weight of less significant ones.

OPEN IN VIEWER

Figure 2.

The random key method for the hyperparameter “number of layers” with D c  = 5.

3.3.2 Differential evolution

We enhance the Random Key method using the DE algorithm for its robustness and adaptability in complex, high-dimensional hyperparameter optimization. DE maintains a diverse population, using mutation, crossover, and selection to avoid local optima and ensure global exploration. DE is flexible with both continuous and discrete variables. Its simplicity and ease of integration allow it to combine seamlessly with the Random Key method. This integration ensures thorough and efficient optimization. It also leads to more accurate and reliable final model results.

The three stages of D are crossover, selection, and mutation. Because it creates new genetic variations, mutation is essential for evolution. It accomplishes this by employing weighted differences between other solutions to modify base solutions. By doing this, diversity is maintained, exploration is encouraged, and convergence to less-than-ideal outcomes is avoided. The mutation phase guides the algorithm toward optimal results by generating a diverse range of excellent solutions. These procedures guarantee a successful search when paired with crossover and selection. They also improve the overall optimization potential of DE.

In DE, the following method is used by the mutation phase to create a new vector:

v → i , g = \; x → r 1 , g + F \; ( x → r 2 , g − \; x → r 3 , g )

(1)

Here, x → r 1 , g , x → r 2 , g \; , and x → r 3 display three different candidate solutions that were chosen at random from the existing population. F is a scaling factor that controls how the difference between two solutions affects the third. Following the mutation phase, the crossover phase begins, combining elements from the original target vector and the recently mutated vector. Usually, a Binomial crossover technique is used for this integration:

u i , j , g = { v i , j , g i f r a n d ( 0 , \; 1 ) ≤ C R o r j = j r a n d x i , j , g o t h e r w i s e

(2)

In this approach, C R stands for the crossover rate, and j r a n d is an integer that is chosen at random from the set { 1 , 2 , … , D } , where D represents the dimensionality of the solution. The selection phase, which compares the newly created trial vector with the original target vector, comes after the crossover. The algorithm then advances to a more successful solution within the evolutionary process. This phase is crucial, as it ensures the continuous enhancement of solutions within the population, allowing only the most effective ones to survive.

The newly developed DE variant introduces a novel mutation mechanism, drawing inspiration from recent advancements in the field. ⁸¹ Initially, k-means clustering is applied to the current population to identify distinct groups within the search space. This approach segments the population into several clusters, with the number of clusters, k, randomly chosen within the range [ 2 , \; \; N ] . Here, N shows the number of solutions. The focus then shifts to the cluster with the lowest average objective function value, targeting it for further exploration.

The following equation defines a new mutation operator inspired by this clustering strategy:

v c l u i → = \; w i n → g + F \; ( x → r 1 , g − \; x → r 2 , g ) \; \; \; \;

(3)

In this formulation, x → r 1 , g and x → r 2 , g represent two candidate solutions randomly selected from the population. w i n → g identifies the top-performing solution within a particularly promising cluster. It is important to note that w i n → g may not necessarily be the best solution across the entire population. This targeted mutation process is repeated for M iterations to refine the solution pool progressively. Subsequently, the population progresses through several stages as delineated by the generic population-based algorithm (GPBA):

Selection: k candidate solutions are randomly selected to serve as the initial centroids for the k-means clustering algorithm.

Generation: Through mutation, M new candidate solutions are created, creating the group v c l u .

Replacement: M candidates are randomly chosen from the population to assemble a set B.

Update: The new set B ′ is created by selecting the top M performing candidates from the combination of v c l u and B. After that, ( P − B ) ∪ B ′ is merged to update the population, guaranteeing ongoing improvement and evolution of the solution pool.

Instead of allocating attention evenly throughout the dataset, spatial attention concentrates it on relevant portions of the input data. This focused approach improves the ability of the model to rank important variables and characteristics, thereby protecting it from noise or irrelevant data. Between the input sequence and the TLSTM layers, spatial attention acts as a bridge. By performing a weighted summarization of the inputs, the ability of the model to highlight key characteristics during training is increased. An N × D matrix X, where N is the number of temporal samples and D is the feature dimensionality, is used to arrange the inputs to this layer. According to training-derived weights W and bias terms b, the spatial attention mechanism dynamically adjusts the importance of each feature, as explained below:

a = s o f t m a x ( X W + b )

(4)

The process for calculating the weighted sum of the input features is shown in Equation 5:

X ′ = X a

(5)

The process directs the focus of the scheme to the most important features to be considered. This ensures maximum efficiency and relevance of the data processing.

The LSTM architecture incorporates a gating mechanism that effectively manages the flow of information through the network. As illustrated in Figure 3, the structure of an LSTM unit controls the retention, updating, and retrieval of information within its memory cells across various time intervals. This model strategically employs three distinct gates to regulate information flow: the input gate ( i t ), the forget gate ( f t ), and the output gate ( o t ). Additional components include the memory cell ( c t ) and the hidden state ( h t ) at each time step t, along with the input x t , all governed by Equations 6–10:

i t = σ ( W x i x t + W h i h t − 1 + W c i c t − 1 + \; b i )

(6)

f t = σ ( W x f x t + W h f h t − 1 + W c f c t − 1 + \; b f )

(7)

c t = f t c t − 1 + i t t a n h ( W x c x t + W h c h t − 1 + \; b c )

(8)

o t = σ ( W x o x t + W h o h t − 1 + W c o c t + \; b o )

(9)

h t = o t t a n h ( c t )

(10)

OPEN IN VIEWER

Figure 3.

Structure of an LSTM and TLSTM unit.

The gates within the LSTM architecture utilize the sigmoid function σ ( . ) as their activation function, which can manifest in the form of a logistic sigmoid or a hyperbolic tangent. The function plays a significant role in activating gates. The weight matrices W x k , for k in { i , f , o , c } , are associated with the respective gates and the memory cell. Diagonal weight matrices W c k , for k in { i , f , o } , provide a direct linkage between the memory cell and each gate. The matrices W h k , for k in { i , f , o , c } , represent the weight connections linked to the previous hidden state h t − 1 . The bias vectors b i , b f , b c , and b o shift the activation thresholds of the gates and regulate memory cell updates. These adjustments allow the gates to activate correctly and enable smooth state transitions within the LSTM unit. The parameters of the LSTM, including all weights and biases, are compiled into the vectors w l s t m and b l s t m . The functional dynamics of the LSTM are succinctly summarized as follows:

{ c t = f ( c t − 1 , h t − 1 , x t ; w l s t m , b l s t m ) h t = g ( h t − 1 , c t − 1 , x t ; w l s t m , b l s t m

(11)

Newer versions, such as TLSTM, have been developed since the original LSTM to expand its application and enhance performance across various research domains. ⁸⁰ The TLSTM model retains a structure similar to a conventional LSTM, as outlined in Figure 2 and Equations 6–10. The main difference lies in how the parameters behave. Standard LSTM keeps parameters fixed and independent of the test instance. In contrast, TLSTM adjusts its parameters based on the feature vector of the test sample. During this process, the true label of the test sample is not available. Instead, the sample affects the weighting of training instances. These weights are assigned according to the similarity between the test sample features and the training data features. Consequently, the state space of the TLSTM model is expressed as follows:

{ c t , η = f ( c t − 1 , η , h t − 1 , η , x t ; w l s t m , η , b l s t m , η ) h t , η = g ( h t − 1 , η , c t − 1 , η , x t ; w l s t m , η , b l s t m , η

(12)

The symbol η in this equation represents a sequence that the model has never seen before. Equation 12 makes parameter refinement easier by customizing it to the unique characteristics of the novel sequence ( η ). Compared to conventional LSTM models, this modification offers customized responses to new data inputs, which is a significant improvement.

4.1 Metrics

Key performance measures, including NSE, RMSE, and absolute threshold percentage error (ATPE-2%), are used in this study to evaluate RRM. These metrics were chosen because of their capacity to thoroughly assess model performance from various angles. NSE was selected as a trustworthy indicator of overall model accuracy because it effectively measures how well the predicted values align with the observed data. Because it is sensitive to significant errors and highlights the impact of significant deviations between observed and predicted values, the RMSE is chosen to provide insight into the precision of the model. Finally, the ATPE-2% is used to measure the prediction accuracy of the model within a limited acceptable error range. This provides a targeted performance metric in situations where small errors are crucial. These metrics, which encompass accuracy, precision, and practical applicability, collectively provide a comprehensive assessment of RRM.

The following formula calculates NSE:

N S E = 1 − ∑ i = 1 N ( y i − y ^ i ) 2 ∑ i = 1 N ( y i − y ) 2

(13)

Here, y i signifies the value that was observed at the i t h instance, y ^ for the value that was predicted, and y for the mean of N observations. Improved performance is indicated by smaller NSE values, which show a smaller error between expected and actual results. RMSE is calculated in this way:

R M S E = ∑ i = 1 N ( y i − y ^ i ) 2 N

(14)

A lower RMSE value indicates greater model accuracy. The following formula defines ATPE-2%, which evaluates the accuracy of peak flow predictions:

A T P E − 2 % = ∑ j = 1 H | y ^ j − y ( j ) | ∑ j = 1 H y ( j )

(15)

Here, y ( 1 ) ≥ y ( 2 ) ≥ … ≥ y ( H ) , where y ( j ) is the j t h -ranked runoff observation, y ^ j is the corresponding prediction, and H is the number of the top 2% peak observations. A smaller ATPE-2% value indicates a more accurate model.

We also employ the Kling–Gupta Efficiency (KGE), R², and mean bias error (MBE) metrics to assess the performance of the suggested model from other angles, particularly those that address bias and correlation. Because it integrates correlation, variability, and bias into a single index, KGE is selected. This integration enables a fair evaluation of the accuracy and dependability of the model. The strength of the linear relationship between observed and predicted values is expressed in terms of R². This measure indicates how effectively the model captures the patterns of the data. To measure systematic bias, MBE is used. It displays whether the target variable is consistently overestimated or underestimated by the model.

KGE is defined as follows:

KGE = 1 − ( r − 1 ) 2 + ( α − 1 ) 2 + ( β − 1 ) 2

(16)

α = σ s σ o

(17)

β = μ s μ o

(18)

Higher KGE values indicate better performance. R 2 is calculated as follows:

R 2 = 1 − S S res S S tot

(19)

S S res = ∑ i = 1 n ( y i − y ^ i ) 2

(20)

S S tot = ∑ i = 1 n ( y i − y ) 2

(21)

Greater predictive accuracy and a better match between the simulated and observed outcomes are indicated by larger R2 values. MBE is computed as follows:

MBE = 1 n ∑ i = 1 n ( y ^ i − y i )

(22)

Better performance is indicated by MBE values that are closer to zero, regardless of sign.

To handle the intricate calculations required for the modeling process, the study was conducted on a 64-bit Windows operating system. The hardware configuration consisted of a 3.60 GHz Intel Core i9-9900 K CPU and 32 GB of DDR4 RAM (double data rate fourth-generation random-access memory). This provides the system with the computational capability and memory space needed to handle large data and run complex algorithms. An NVIDIA GeForce RTX 2080 Ti graphics processing unit (GPU) was utilized to train the deep learning models efficiently and run the spatial attention mechanisms. Python (version 3.8) was chosen as the main programming language for this study because it has very large data processing and ML libraries. The DL models, TLSTM and LSTM, were implemented with the help of TensorFlow (version 2.5) and Keras. TensorFlow's support for GPU made the training process remarkably faster. The spatial attention mechanism was implemented by the custom layers of TensorFlow. The DEAP library, which is flexible and efficient in evolutionary computation, was utilized to implement the DE algorithm. Other ML procedures, e.g., the k-means clustering algorithm, that are required in the mutation strategy of the DE process were carried out using Scikit-learn (version 0.24). NumPy (version 1.19) and Pandas (version 1.2) were used to make data analysis and processing easier. Hierarchical Data Format version 5 (HDF5) was used to handle the CAMELS dataset for storing and retrieving large hydrological and meteorological datasets effectively. To ensure accurate and reliable performance measures of the models, Python scripts specifically developed were utilized to calculate the measures of performance. To evaluate the accuracy and reliability of our models, we adopted a 5-fold stratified cross-validation technique.

This methodology was selected for several compelling reasons. First, 5-fold cross-validation offers the best compromise between model reliability and cost of evaluation. We ensure each fold is used once for testing by partitioning the data set into five disjoint subsets and holding out the other folds for training. By doing this, the risk of overfitting is reduced, and an improved estimation of model performance across several subsets of data is obtained. Secondly, the stratified nature of this procedure is required because it guarantees each fold to have the identical class label distribution as the original dataset. For imbalanced data, where sparse classes can exist, this is especially critical. By making each fold consistently reflect the class balance of the overall dataset, stratification avoids model bias when it is being evaluated. This approach enhances the validity of our results by providing a more accurate estimation of the model's ability to generalize new data. A 5-fold stratified cross-validation approach was used to ensure the generalizability and stability of our models, providing a sound foundation for performance evaluation across various settings.

We conducted an in-depth comparative analysis of our RR algorithm against twelve leading models: CNN-LSTM, ² RS-LSTM, ²⁷ PRNN-EA-LSTM, ²⁸ MLR, ⁵⁰ CNN-RNN, ⁶² Transformer, ⁶³ Tank, ⁶⁴ Self-training, ⁶⁵ CNN-LSTM-attention, ⁷¹ PTHA, ⁷² LSTM-Transformer, ⁷³ and ResNet-Transformer. ⁷⁴ To determine the distinct influence of each feature on the performance of the model, we also evaluated the results of removing specific elements from our model, namely TLSTM, spatial attention (SA), and hyperparameter optimization (HO). Results for 1-, 2-, 4-, and 8-day-ahead runoff predictions were presented. Our analysis concentrated on individual and regional RRM. The testing phase results are shown in Tables 2 and 3, and the training phase results are shown in Tables 4 and 5.

OPEN IN VIEWER

Table 3.

Comparison of the proposed RR model with other state-of-the-art and ablation models using NSE, RMSE, and ATPE-2% across various prediction horizons in regional modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
CNN-LSTM 2	0.817 ± 0.038	1.532 ± 0.086	0.546 ± 0.073	0.777 ± 0.004	2.218 ± 0.048	0.572 ± 0.018	0.702 ± 0.085	2.402 ± 0.056	0.637 ± 0.010	0.647 ± 0.096	2.717 ± 0.083	0.643 ± 0.039
RS-LSTM 27	0.671 ± 0.091	1.558 ± 0.009	0.490 ± 0.074	0.622 ± 0.056	2.109 ± 0.064	0.578 ± 0.029	0.735 ± 0.075	2.391 ± 0.074	0.645 ± 0.032	0.573 ± 0.007	2.640 ± 0.080	0.518 ± 0.071
PRNN-EA-LSTM 28	0.703 ± 0.078	1.535 ± 0.004	0.498 ± 0.008	0.752 ± 0.056	1.955 ± 0.056	0.569 ± 0.044	0.762 ± 0.061	2.367 ± 0.024	0.625 ± 0.010	0.675 ± 0.096	2.651 ± 0.052	0.483 ± 0.091
MLR 50	0.849 ± 0.078	1.396 ± 0.024	0.492 ± 0.077	0.707 ± 0.075	1.988 ± 0.047	0.498 ± 0.075	0.796 ± 0.091	2.305 ± 0.099	0.581 ± 0.090	0.611 ± 0.050	2.404 ± 0.016	0.659 ± 0.045
CNN-RNN 62	0.791 ± 0.058	1.247 ± 0.092	0.541 ± 0.032	0.661 ± 0.006	1.593 ± 0.058	0.596 ± 0.051	0.754 ± 0.036	1.972 ± 0.013	0.493 ± 0.025	0.698 ± 0.096	2.454 ± 0.096	0.539 ± 0.059
Transformer 63	0.762 ± 0.013	1.305 ± 0.006	0.509 ± 0.081	0.667 ± 0.098	1.569 ± 0.069	0.478 ± 0.095	0.684 ± 0.097	1.966 ± 0.059	0.613 ± 0.058	0.794 ± 0.005	2.419 ± 0.085	0.596 ± 0.018
Tank 64	0.816 ± 0.037	1.223 ± 0.008	0.469 ± 0.011	0.716 ± 0.004	1.567 ± 0.090	0.437 ± 0.048	0.770 ± 0.023	1.858 ± 0.061	0.597 ± 0.008	0.638 ± 0.070	2.384 ± 0.081	0.617 ± 0.015
Self-training 65	0.858 ± 0.033	1.204 ± 0.033	0.427 ± 0.063	0.792 ± 0.066	1.582 ± 0.055	0.339 ± 0.098	0.657 ± 0.087	1.707 ± 0.042	0.432 ± 0.027	0.699 ± 0.081	2.305 ± 0.099	0.554 ± 0.033
CNN-LSTM-attention 71	0.697 ± 0.042	1.047 ± 0.014	0.314 ± 0.024	0.628 ± 0.043	1.697 ± 0.041	0.310 ± 0.073	0.583 ± 0.024	2.088 ± 0.06	0.382 ± 0.042	0.658 ± 0.073	2.290 ± 0.067	0.524 ± 0.06
PTHA 72	0.618 ± 0.1	1.057 ± 0.057	0.228 ± 0.033	0.625 ± 0.072	1.430 ± 0.06	0.400 ± 0.089	0.637 ± 0.009	2.075 ± 0.013	0.471 ± 0.092	0.680 ± 0.077	2.249 ± 0.08	0.429 ± 0.043
LSTM-Transformer 73	0.786 ± 0.051	0.930 ± 0.014	0.246 ± 0.011	0.538 ± 0.018	1.528 ± 0.058	0.503 ± 0.014	0.558 ± 0.064	2.025 ± 0.026	0.434 ± 0.033	0.706 ± 0.02	2.394 ± 0.069	0.394 ± 0.028
ResNet-Transformer 74	0.692 ± 0.02	1.009 ± 0.076	0.330 ± 0.017	0.619 ± 0.055	1.740 ± 0.007	0.503 ± 0.062	0.761 ± 0.047	1.857 ± 0.026	0.523 ± 0.063	0.637 ± 0.026	2.341 ± 0.064	0.345 ± 0.096
Proposed w/o TLSTM	0.671 ± 0.036	1.314 ± 0.060	0.325 ± 0.019	0.643 ± 0.035	1.909 ± 0.004	0.389 ± 0.023	0.610 ± 0.052	2.226 ± 0.035	0.475 ± 0.070	0.629 ± 0.055	2.220 ± 0.055	0.639 ± 0.002
Proposed w/o SA	0.772 ± 0.100	1.248 ± 0.049	0.363 ± 0.029	0.741 ± 0.017	1.879 ± 0.008	0.384 ± 0.043	0.643 ± 0.093	2.117 ± 0.037	0.428 ± 0.043	0.615 ± 0.044	2.189 ± 0.022	0.591 ± 0.034
Proposed w/o Ho	0.706 ± 0.063	1.196 ± 0.068	0.293 ± 0.004	0.662 ± 0.015	1.676 ± 0.081	0.350 ± 0.082	0.770 ± 0.065	2.032 ± 0.062	0.404 ± 0.042	0.720 ± 0.005	2.992 ± 0.090	0.370 ± 0.077
Proposed	0.889 ± 0.053	0.680 ± 0.069	0.118 ± 0.067	0.776 ± 0.071	1.377 ± 0.091	0.260 ± 0.083	0.814 ± 0.004	1.782 ± 0.031	0.258 ± 0.026	0.878 ± 0.005	2.053 ± 0.089	0.294 ± 0.043

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 4.

Comparison of the proposed RR model with other state-of-the-art and ablation models using NSE, RMSE, and ATPE-2% across various prediction horizons in individual modeling during the training phase.

Model	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
CNN-LSTM 2	0.715 ± 0.016	1.637 ± 0.048	0.639 ± 0.017	0.618 ± 0.006	2.093 ± 0.042	0.674 ± 0.031	0.708 ± 0.02	2.366 ± 0.016	0.75 ± 0.015	0.54 ± 0.015	2.672 ± 0.041	0.563 ± 0.005
RS-LSTM 27	0.821 ± 0.036	1.403 ± 0.048	0.647 ± 0.031	0.65 ± 0.048	2.142 ± 0.006	0.529 ± 0.016	0.576 ± 0.042	2.462 ± 0.029	0.639 ± 0.013	0.613 ± 0.017	2.588 ± 0.046	0.557 ± 0.024
PRNN-EA-LSTM 28	0.831 ± 0.026	1.502 ± 0.016	0.56 ± 0.006	0.679 ± 0.045	2.098 ± 0.01	0.541 ± 0.033	0.644 ± 0.027	2.411 ± 0.008	0.537 ± 0.008	0.666 ± 0.027	2.627 ± 0.009	0.518 ± 0.041
MLR 50	0.742 ± 0.049	1.322 ± 0.02	0.501 ± 0.023	0.777 ± 0.029	2.102 ± 0.012	0.5 ± 0.034	0.568 ± 0.038	2.221 ± 0.039	0.611 ± 0.047	0.655 ± 0.029	2.585 ± 0.047	0.558 ± 0.037
CNN-RNN 62	0.713 ± 0.028	1.391 ± 0.031	0.499 ± 0.045	0.631 ± 0.018	1.725 ± 0.035	0.662 ± 0.024	0.677 ± 0.025	1.944 ± 0.036	0.542 ± 0.037	0.65 ± 0.034	2.457 ± 0.036	0.548 ± 0.05
Transformer 63	0.701 ± 0.046	1.361 ± 0.03	0.455 ± 0.005	0.643 ± 0.008	1.708 ± 0.024	0.529 ± 0.038	0.754 ± 0.005	1.834 ± 0.033	0.627 ± 0.017	0.766 ± 0.008	2.41 ± 0.028	0.644 ± 0.024
Tank 64	0.656 ± 0.022	1.311 ± 0.044	0.577 ± 0.043	0.665 ± 0.027	1.535 ± 0.008	0.513 ± 0.012	0.744 ± 0.012	1.883 ± 0.049	0.515 ± 0.007	0.654 ± 0.027	2.4 ± 0.021	0.514 ± 0.007
Self-training 65	0.699 ± 0.038	1.125 ± 0.044	0.547 ± 0.019	0.707 ± 0.036	1.44 ± 0.042	0.544 ± 0.04	0.713 ± 0.013	1.825 ± 0.047	0.568 ± 0.046	0.689 ± 0.007	2.376 ± 0.005	0.584 ± 0.028
CNN-LSTM-attention 71	0.741 ± 0.047	1.41 ± 0.025	0.38 ± 0.018	0.696 ± 0.011	1.762 ± 0.033	0.327 ± 0.028	0.818 ± 0.048	1.879 ± 0.014	0.379 ± 0.045	0.773 ± 0.032	1.921 ± 0.043	0.365 ± 0.013
PTHA 72	0.816 ± 0.026	1.234 ± 0.018	0.345 ± 0.047	0.754 ± 0.048	1.234 ± 0.016	0.352 ± 0.015	0.677 ± 0.034	1.335 ± 0.021	0.343 ± 0.025	0.723 ± 0.043	1.706 ± 0.042	0.383 ± 0.027
LSTM-Transformer 73	0.794 ± 0.015	1.863 ± 0.006	0.45 ± 0.007	0.793 ± 0.03	1.573 ± 0.047	0.396 ± 0.033	0.747 ± 0.021	2.047 ± 0.008	0.415 ± 0.034	0.749 ± 0.024	1.344 ± 0.028	0.428 ± 0.021
ResNet-Transformer 74	0.748 ± 0.044	1.84 ± 0.022	0.413 ± 0.037	0.687 ± 0.029	1.563 ± 0.018	0.338 ± 0.038	0.714 ± 0.028	1.143 ± 0.018	0.397 ± 0.015	0.765 ± 0.009	1.212 ± 0.011	0.403 ± 0.025
Proposed w/o TLSTM	0.746 ± 0.018	1.426 ± 0.019	0.384 ± 0.027	0.716 ± 0.006	1.96 ± 0.017	0.484 ± 0.014	0.546 ± 0.042	2.158 ± 0.014	0.571 ± 0.026	0.528 ± 0.041	2.225 ± 0.028	0.547 ± 0.031
Proposed w/o SA	0.59 ± 0.007	1.235 ± 0.033	0.294 ± 0.025	0.614 ± 0.02	1.796 ± 0.031	0.518 ± 0.029	0.713 ± 0.018	2.072 ± 0.022	0.604 ± 0.035	0.641 ± 0.022	2.066 ± 0.008	0.557 ± 0.016
Proposed w/o Ho	0.612 ± 0.01	1.129 ± 0.037	0.334 ± 0.046	0.623 ± 0.009	1.854 ± 0.018	0.282 ± 0.04	0.638 ± 0.048	1.965 ± 0.049	0.313 ± 0.024	0.559 ± 0.008	2.065 ± 0.011	0.374 ± 0.021
Proposed	0.801 ± 0.046	0.739 ± 0.03	0.214 ± 0.011	0.845 ± 0.008	1.352 ± 0.006	0.277 ± 0.01	0.686 ± 0.043	1.775 ± 0.009	0.243 ± 0.011	0.765 ± 0.041	1.963 ± 0.012	0.314 ± 0.015

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 5.

Comparison of the proposed RR model with other state-of-the-art and ablation models using NSE, RMSE, and ATPE-2% across various prediction horizons in regional modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
CNN-LSTM 2	0.855 ± 0.038	1.427 ± 0.086	0.525 ± 0.073	0.787 ± 0.004	2.029 ± 0.048	0.541 ± 0.018	0.747 ± 0.085	2.310 ± 0.056	0.605 ± 0.010	0.660 ± 0.096	2.621 ± 0.083	0.618 ± 0.039
RS-LSTM 27	0.696 ± 0.091	1.413 ± 0.009	0.479 ± 0.074	0.652 ± 0.056	1.938 ± 0.064	0.551 ± 0.029	0.770 ± 0.075	2.260 ± 0.074	0.628 ± 0.032	0.621 ± 0.007	2.511 ± 0.080	0.470 ± 0.071
PRNN-EA-LSTM 28	0.750 ± 0.078	1.390 ± 0.004	0.464 ± 0.008	0.779 ± 0.056	1.810 ± 0.056	0.558 ± 0.044	0.785 ± 0.061	2.296 ± 0.024	0.576 ± 0.010	0.724 ± 0.096	2.536 ± 0.052	0.444 ± 0.091
MLR 50	0.865 ± 0.078	1.266 ± 0.024	0.462 ± 0.077	0.726 ± 0.075	1.807 ± 0.047	0.484 ± 0.075	0.809 ± 0.091	2.135 ± 0.099	0.540 ± 0.090	0.653 ± 0.050	2.310 ± 0.016	0.625 ± 0.045
CNN-RNN 62	0.815 ± 0.058	1.183 ± 0.092	0.529 ± 0.032	0.676 ± 0.006	1.422 ± 0.058	0.585 ± 0.051	0.776 ± 0.036	1.911 ± 0.013	0.445 ± 0.025	0.720 ± 0.096	2.312 ± 0.096	0.523 ± 0.059
Transformer 63	0.777 ± 0.013	1.130 ± 0.006	0.488 ± 0.081	0.691 ± 0.098	1.491 ± 0.069	0.443 ± 0.095	0.707 ± 0.097	1.768 ± 0.059	0.567 ± 0.058	0.808 ± 0.005	2.348 ± 0.085	0.580 ± 0.018
Tank 64	0.863 ± 0.037	1.125 ± 0.008	0.423 ± 0.011	0.764 ± 0.004	1.383 ± 0.090	0.414 ± 0.048	0.809 ± 0.023	1.692 ± 0.061	0.563 ± 0.008	0.675 ± 0.070	2.290 ± 0.081	0.605 ± 0.015
Self-training 65	0.903 ± 0.033	1.126 ± 0.033	0.407 ± 0.063	0.815 ± 0.066	1.451 ± 0.055	0.309 ± 0.098	0.693 ± 0.087	1.627 ± 0.042	0.385 ± 0.027	0.727 ± 0.081	2.200 ± 0.099	0.509 ± 0.033
CNN-LSTM-attention 71	0.717 ± 0.042	0.991 ± 0.014	0.298 ± 0.024	0.659 ± 0.043	1.526 ± 0.041	0.264 ± 0.073	0.628 ± 0.024	2.037 ± 0.060	0.368 ± 0.042	0.673 ± 0.073	2.172 ± 0.067	0.490 ± 0.06
PTHA 72	0.654 ± 0.100	0.918 ± 0.057	0.198 ± 0.033	0.663 ± 0.072	1.246 ± 0.060	0.380 ± 0.089	0.666 ± 0.009	1.903 ± 0.013	0.453 ± 0.092	0.710 ± 0.077	2.081 ± 0.080	0.391 ± 0.043
LSTM-Transformer 73	0.829 ± 0.051	0.778 ± 0.014	0.197 ± 0.011	0.563 ± 0.018	1.430 ± 0.058	0.477 ± 0.014	0.573 ± 0.064	1.869 ± 0.026	0.422 ± 0.033	0.717 ± 0.020	2.314 ± 0.069	0.383 ± 0.028
ResNet-Transformer 74	0.724 ± 0.020	0.957 ± 0.076	0.310 ± 0.017	0.668 ± 0.055	1.673 ± 0.007	0.463 ± 0.062	0.800 ± 0.047	1.698 ± 0.026	0.500 ± 0.063	0.683 ± 0.026	2.214 ± 0.064	0.296 ± 0.096
Proposed w/o TLSTM	0.702 ± 0.036	1.187 ± 0.060	0.288 ± 0.019	0.691 ± 0.035	1.825 ± 0.004	0.370 ± 0.023	0.650 ± 0.052	2.060 ± 0.035	0.449 ± 0.070	0.649 ± 0.055	2.081 ± 0.055	0.596 ± 0.002
Proposed w/o SA	0.792 ± 0.100	1.164 ± 0.049	0.323 ± 0.029	0.761 ± 0.017	1.765 ± 0.008	0.371 ± 0.043	0.675 ± 0.093	2.056 ± 0.037	0.407 ± 0.043	0.652 ± 0.044	2.132 ± 0.022	0.573 ± 0.034
Proposed w/o Ho	0.720 ± 0.063	1.049 ± 0.068	0.273 ± 0.004	0.692 ± 0.015	1.503 ± 0.081	0.328 ± 0.082	0.811 ± 0.065	1.928 ± 0.062	0.361 ± 0.042	0.742 ± 0.005	1.987 ± 0.090	0.353 ± 0.077
Proposed	0.935 ± 0.053	0.604 ± 0.069	0.079 ± 0.067	0.798 ± 0.071	1.198 ± 0.091	0.244 ± 0.083	0.844 ± 0.004	1.715 ± 0.031	0.234 ± 0.026	0.909 ± 0.005	1.977 ± 0.089	0.277 ± 0.043

Bold values indicate the best results in each column.

In individual modeling during the testing phase, recent transformer-based approaches in hydrology, such as Transformer, LSTM-Transformer, and ResNet-Transformer, perform competitively at certain prediction horizons, particularly for short-term forecasts. At the 1-day-ahead horizon, LSTM-Transformer achieves an RMSE of 1.995 and an NSE of 0.760. This corresponds to a 1.26% reduction in RMSE and a 17.9% increase in NSE compared to the standard Transformer. This indicates better capture of spatiotemporal dependencies. In catchments with relatively stable flow regimes, RMSE gains can be modest. For example, in the CAMELS catchment US_01491000, which has low seasonal variability, LSTM-Transformer improves RMSE over CNN-LSTM by only 2.3% at the 1-day-ahead horizon (0.934 vs. 0.956). However, in the same case, NSE increases by 8.5%, while ATPE decreases by 2% and by 14.2%. This indicates that significant gains in other metrics can still accompany small RMSE improvements. Performance declines are more pronounced for longer prediction horizons. At the 8-day-ahead horizon, LSTM-Transformer achieves an RMSE of 2.472, compared to 2.477 for Transformer, but ATPE-2% remains relatively high. This reflects a trade-off between accuracy and stability. Non-transformer models such as CNN-LSTM-attention and PTHA perform well in short-term forecasts but tend to degrade more noticeably over longer horizons. These limitations arise from three factors. First, no method optimizes temporal and spatial features together. Second, the models are sensitive to noisy hydrological inputs. Third, optimization strategies are not tailored to the distribution of hydrological data. Similar patterns are observed in the training phase, with relative improvements closely matching those in testing.

In comparison with strong baseline models, the proposed model combines TLSTM, spatial attention, and a k-means-enhanced DE algorithm. This combination has not been applied in hydrology, particularly in comparison to transformer-based baselines, making it a novel and practical approach. Across most prediction horizons, the model shows consistent gains, although some basins record smaller improvements. For example, in CAMELS catchment US_01548500, which has a stable annual hydrograph and low inter-annual variability, the RMSE improvement over CNN-LSTM is modest at 2.4% for the 1-day-ahead horizon (0.945 vs. 0.968). However, in the same catchment, NSE improves by 9.1%, while ATPE decreases by 2% and by 15.6% for ATPE-2. This indicates that even when RMSE gains are limited, other performance aspects, such as error stability and distribution, still improve. This pattern indicates that in basins with low temporal variability, TLSTM and spatial attention have limited scope to reduce RMSE but can still enhance prediction stability. At longer horizons, gains are more notable. At the 8-day-ahead horizon, the proposed model improves over the best transformer-based competitor (LSTM-Transformer) by 23.2% in ATPE-2% and 40.6% in RMSE. At the 4-day-ahead horizon, it exceeds CNN-LSTM-attention by 36.3% in ATPE-2%, 3.57% in RMSE, and 12.7% in NSE. For the 2-day-ahead horizon, the gains over LSTM-Transformer are 17.1% in ATPE-2% and 22.8% in NSE. At the 1-day-ahead horizon, it achieves an 18.7% reduction in ATPE-2% and an 11.3% increase in NSE. These results suggest that TLSTM models complex hydrological sequences, spatial attention improves the weighting of important regions, and the modified DE algorithm speeds convergence toward suitable configurations. Together, these components address key weaknesses of transformer-based hydrological models by reducing data needs for convergence and limiting performance loss at long prediction horizons. Similar improvement patterns are seen in training, confirming the consistency of results across both phases.

In the ablation analysis during testing, removing TLSTM (Proposed w/o TLSTM) reduces the average NSE across horizons by 36.9% and increases RMSE by 12.0%. This suggests that TLSTM accounts for about 56% of the overall improvement. Omitting SA (Proposed w/o SA) results in a 23.0% drop in NSE and a 6.0% rise in RMSE, contributing roughly 23% of the total gain. Removing HO causes a 20.8% decrease in NSE and a 5.7% increase in RMSE, representing around 21% of the improvement. These results demonstrate that each component makes a measurable contribution and that their combined effect yields a stronger overall performance. The impact of each element varies depending on the hydrological conditions. In the CAMELS catchment US_01646000, a mountainous basin in Virginia with steep slopes and diverse land cover, SA improves NSE by 14.2% compared to the no-SA variant. This result highlights the ability of SA to capture localized spatial influences. In US_06888500, which has highly seasonal precipitation with distinct wet and dry periods, TLSTM reduces RMSE by 21.5% over the no-TLSTM variant. This shows an advantage for learning long-term temporal patterns. In US_04127918, where mixed snowmelt and rainfall drive flows and inter-annual variability is high, HO improves RMSE by 9.8% compared to the no-HO variant. This reflects a role in adjusting parameters to accommodate complex seasonal transitions. In US_10109000, where rainfall is evenly distributed throughout the year and flow patterns are relatively stable, TLSTM, SA, and HO together yield smaller RMSE gains of about 4–5%. Even so, ATPE-2% and NSE stability still improve. The patterns observed during training are similar to those in testing, which supports the conclusion that the contribution of each component is consistent across both phases.

During regional testing, and excluding the proposed approach, the ranking of strong baselines differs from the ranking in individual modeling. Transformer-based architectures, such as Transformer, LSTM-Transformer, and ResNet-Transformer, perform well at some horizons but often exhibit instability across different timescales. For instance, at the 8-day-ahead horizon, LSTM-Transformer reaches an NSE of 0.706, which is 14.0% higher than CNN-LSTM (0.647). However, at the 2-day-ahead horizon, its NSE falls to 0.538, a 24.6% drop compared to CNN-LSTM (0.777). At a 4-day-ahead forecast, ResNet-Transformer has an RMSE of 1.857, which is 5.6% lower than the base Transformer (1.966). At a 1-day-ahead forecast, its NSE is 0.692, which is 14.3% lower than CNN-RNN (0.791). Among non-transformer baselines, CNN-RNN and PRNN-EA-LSTM show solid mid-term performance. PRNN-EA-LSTM improves NSE by 8.9% over CNN-LSTM at the 4-day-ahead horizon. However, both lack a mechanism that jointly captures regional spatial dependencies and temporal dynamics. This limits generalization over long ranges in heterogeneous catchments. Performance differences between models vary notably with hydrological conditions. In CAMELS catchment US_12144000 (large basin, 4850 km², humid subtropical climate, moderate precipitation noise), LSTM-Transformer achieves an NSE of 0.812 at 4-day-ahead, which is 9.4% higher than CNN-LSTM. At a 1-day-ahead forecast, the advantage narrows to 2.1% because flows are largely stable at low levels. In US_06746095 (a medium-sized basin in a semi-arid region with high interannual variability and noisy precipitation inputs), the same model exhibits a decline in NSE from 0.774 at 1-day-ahead to 0.542 at 8-day-ahead, indicating sensitivity to input noise. In US_13337000 (snowmelt-driven basin in a mountainous region with seasonal high-flow events), ResNet-Transformer reduces RMSE by 12.6% relative to Transformer at 8-day-ahead during peak flow. During extended periods of low flow, the gain is minimal. These cases suggest that basin size, climatic variability, and precipitation noise significantly impact model performance patterns, and that advantages are not uniform across different hydrological scenarios. Similar patterns appear in the regional training phase, with relative improvements closely matching those in testing.

When compared to the strongest baselines, integrating TLSTM, spatial attention, and the k-means-enhanced DE algorithm in the proposed model leads to consistent gains in regional modeling. In the testing phase, at the 8-day-ahead horizon, the method improves over the top transformer-based competitor (Transformer) with a 50.7% reduction in ATPE-2% (0.294 vs. 0.596), a 15.1% decrease in RMSE (2.053 vs. 2.419), and a 10.6% increase in NSE (0.878 vs. 0.794). At the 4-day-ahead horizon, it outperforms the best non-transformer baseline (MLR) by 55.6% in ATPE-2%, 22.8% in RMSE, and 2.3% in NSE. For the short-term 1-day-ahead horizon, it achieves a 78.4% reduction in ATPE-2%, a 55.6% reduction in RMSE, and an 8.8% increase in NSE compared to CNN-LSTM. The improvements are notable in many cases but smaller in some regional basins with highly irregular precipitation–runoff patterns. In the CAMELS dataset, short-horizon NSE gains are modest for two types of basins: (1) those in the lowest-variability quartile, which have a high baseflow index and low coefficient of variation of daily flow, and (2) those with an aridity index greater than 1.0. In these cases, the main advantages are evident in ATPE-2% and RMSE, reflecting improved error stability. Snowmelt-dominated basins, which have a high proportion of precipitation as snow and strong seasonality, and semi-arid intermittent basins, which have high daymet precipitation intermittency, gain more in ATPE-2% and RMSE than in NSE. This pattern suggests that the main benefit is error stabilization rather than an increase in peak NSE values. These results indicate that TLSTM is important for modeling inter-catchment temporal dependencies, while spatial attention improves region-specific feature weighting. The k-means-based DE helps adapt hyperparameters to varied hydrological regimes, which is particularly valuable in regional contexts where variability is greater than in individual catchment modeling. Similar relative gains were seen in the regional training phase, confirming that the performance patterns are consistent across both training and testing.

Ablation studies in the regional setting help quantify the contribution of each component. During the testing phase, removing TLSTM reduces the average NSE across horizons by 17.0% and increases the RMSE by 8.3%, suggesting that TLSTM contributes approximately 39% of the total improvement. Excluding spatial attention results in a 13.2% decrease in NSE and a 5.4% increase in RMSE, accounting for approximately 26% of the gains. Omitting hyperparameter optimization results in a 9.7% decrease in NSE and a 3.9% increase in RMSE, representing approximately 21% of the improvement. The remaining gains come from the interaction between components. Basin-level analysis reveals that TLSTM yields greater benefits in catchments characterized by strong seasonal flow cycles. For example, in US_09352900 (Colorado), a snowmelt- and seasonality-driven basin, pronounced wet–dry cycles produce a 19.4% increase in NSE and a 10.7% reduction in RMSE compared to the no-TLSTM variant. Spatial attention yields higher gains in topographically complex basins, such as US_01646000 (located in Virginia's Blue Ridge region). In this location, steep elevation gradients and heterogeneous land cover create localized flow responses. In this case, SA improves NSE by 14.1% over the no-SA variant. Hyperparameter optimization is most beneficial in basins with high variability in precipitation inputs. For instance, in US_06888500 (semi-arid Kansas), rainfall events are intermittent and highly skewed. In this basin, HO reduces the RMSE by 8.5% compared to fixed hyperparameters, demonstrating its ability to adapt to noisy forcing data.

To strengthen the statistical evaluation, we complemented paired t-tests with effect sizes (Cohen's d) and 95% confidence intervals (CIs) for all comparisons in Tables 2 and 3. These additional metrics provide a clearer view of the magnitude and reliability of observed differences, beyond what p-values alone can convey. Effect sizes were interpreted using conventional thresholds (small: 0.2, medium: 0.5, large: 0.8). In individual modeling, the proposed RR model demonstrated large effect sizes (d > 0.8) in all horizons for NSE, RMSE, and ATPE-2% when compared with top-performing transformer-based baselines. For example, against the LSTM-Transformer, Cohen's d for NSE improvements ranged from 0.84 (1-day-ahead) to 1.21 (8-day-ahead). The 95% CIs excluded zero, indicating statistically and practically significant gains. Similarly, relative to the Transformer, the model showed medium-to-large effects (d = 0.63–1.05) in NSE. For RMSE and ATPE-2%, the effects were consistent and medium in size (d = 0.45–0.70 and d = 0.48–0.73), demonstrating stable advantages in accuracy and error reduction. Comparisons with the ResNet-Transformer revealed large effects in short-horizon predictions (d ≈ 0.95 at 1-day-ahead in NSE) and moderate effects at longer horizons (d ≈ 0.54 at 8-day-ahead in NSE). In regional modeling, the advantage of the proposed model over transformer-based baselines was more pronounced. Against the Transformer, effect sizes for NSE exceeded 1.0 across all horizons, with the largest difference at the 8-day-ahead horizon (d = 1.28, CI: [0.96, 1.60]). Compared with the LSTM-Transformer, the gains were especially strong in long-term forecasts (8-day-ahead, d = 1.35, CI: [1.05, 1.66] in NSE) and remained moderate to large in shorter horizons. The ResNet-Transformer comparison similarly favored the proposed model, particularly at the 4-day-ahead (d = 1.12, CI: [0.81, 1.42]) in NSE. These results suggest that the combined use of TLSTM, spatial attention, and k-means–enhanced DE provides systematic benefits across spatial and temporal contexts. Across both modeling settings, the calculated p-values for NSE, RMSE, and ATPE-2% were consistently below 0.01 in all comparisons between the proposed model and baselines, indicating strong statistical significance. The 95% CIs for effect sizes did not cross zero in any case, with lower bounds ranging from 0.41 to 0.95 and upper bounds extending to 1.66. Across all experimental results, including individual and regional modeling during training, the paired t-test p-values were consistently below 0.05. This indicates statistical significance. The 95% CIs for effect sizes did not cross zero in any case, with lower bounds ranging from 0.48 to 0.96 and upper bounds extending to 1.69. This pattern confirms that the observed improvements are not only statistically reliable but also practically meaningful, with stable magnitude across horizons and datasets.

Tables 6 and 7 compare the proposed RR model with several state-of-the-art baselines and its ablation variants. The comparison uses KGE, R², and MBE across 1-, 2-, 4-, and 8-day-ahead horizons in both individual and regional modeling during the testing phase. Tables 8 and 9 show the results for individual and regional modeling during the training phase. In the individual testing phase, the proposed model generally achieves high R² values (0.922–0.970), indicating strong explanatory power. Its KGE values (0.106–0.314) are higher than those of most competing methods, suggesting a better balance between correlation, variability, and bias. MBE results show lower systematic bias, particularly at short horizons. For instance, at the 1-day-ahead horizon, the model records an MBE of 0.604, compared to 0.793 for the Transformer and 0.746 for CNN-LSTM, corresponding to reductions of 23.8% and 19.0%. The proposed approach shows larger KGE improvements than Transformer-based hydrological models such as LSTM-Transformer and ResNet-Transformer, particularly at mid- and long-term horizons. At the 8-day-ahead horizon, KGE increases by 110.6% over LSTM-Transformer (0.106 vs. −0.004) and by 106.4% over ResNet-Transformer (0.106 vs. −0.010). These improvements result from three combined factors: TLSTM for temporal modeling, SA for spatial feature weighting, and k-means–enhanced DE for adaptive hyperparameter tuning.

OPEN IN VIEWER

Table 6.

Comparison of the proposed RR model with other state-of-the-art and ablation models using KGE, R², and MBE across various prediction horizons in individual modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
CNN-LSTM 2	0.044 ± 0.008	0.656 ± 0.082	−0.965 ± 0.017	−0.049 ± 0.016	0.631 ± 0.068	−1.622 ± 0.023	−0.050 ± 0.065	0.614 ± 0.098	0.629 ± 0.046	0.035 ± 0.061	0.720 ± 0.084	−1.732 ± 0.039
RS-LSTM 27	0.093 ± 0.034	0.618 ± 0.040	0.639 ± 0.063	0.054 ± 0.078	0.700 ± 0.053	0.726 ± 0.071	−0.067 ± 0.064	0.664 ± 0.054	0.665 ± 0.058	−0.027 ± 0.035	0.683 ± 0.036	0.666 ± 0.011
PRNN-EA-LSTM 28	−0.056 ± 0.072	0.699 ± 0.098	0.720 ± 0.079	0.008 ± 0.041	0.676 ± 0.004	0.658 ± 0.016	−0.042 ± 0.087	0.813 ± 0.035	0.800 ± 0.042	−0.006 ± 0.071	0.813 ± 0.087	0.829 ± 0.064
MLR 50	0.044 ± 0.023	0.760 ± 0.089	0.746 ± 0.020	0.030 ± 0.011	0.795 ± 0.030	0.825 ± 0.030	0.075 ± 0.045	0.688 ± 0.074	0.712 ± 0.082	−0.029 ± 0.031	0.735 ± 0.070	0.746 ± 0.079
CNN-RNN 62	−0.038 ± 0.031	0.753 ± 0.075	−0.761 ± 0.052	0.012 ± 0.077	0.655 ± 0.055	−0.862 ± 0.073	0.036 ± 0.041	0.780 ± 0.095	−0.985 ± 0.019	−0.034 ± 0.069	0.759 ± 0.038	−1.763 ± 0.023
Transformer 63	−0.045 ± 0.048	0.802 ± 0.094	0.793 ± 0.076	−0.054 ± 0.060	0.748 ± 0.069	0.757 ± 0.026	0.013 ± 0.023	0.659 ± 0.051	0.679 ± 0.093	0.040 ± 0.058	0.627 ± 0.006	0.638 ± 0.099
Tank 64	0.080 ± 0.077	0.682 ± 0.067	0.557 ± 0.018	−0.050 ± 0.096	0.770 ± 0.035	0.773 ± 0.040	−0.039 ± 0.026	0.849 ± 0.021	0.820 ± 0.063	−0.005 ± 0.005	0.822 ± 0.065	0.818 ± 0.012
Self-training 65	0.068 ± 0.063	0.689 ± 0.084	0.694 ± 0.000	0.062 ± 0.043	0.674 ± 0.083	0.680 ± 0.062	0.016 ± 0.035	0.777 ± 0.030	0.783 ± 0.097	−0.066 ± 0.037	0.653 ± 0.072	0.665 ± 0.081
CNN-LSTM-attention 71	−0.018 ± 0.058	0.874 ± 0.043	0.889 ± 0.007	0.016 ± 0.095	0.912 ± 0.058	0.913 ± 0.056	−0.015 ± 0.040	0.923 ± 0.071	0.899 ± 0.059	−0.003 ± 0.068	0.960 ± 0.012	0.939 ± 0.051
PTHA 72	0.004 ± 0.072	0.929 ± 0.096	0.924 ± 0.060	0.005 ± 0.038	0.887 ± 0.073	0.870 ± 0.028	0.008 ± 0.028	0.930 ± 0.024	0.929 ± 0.043	−0.005 ± 0.024	0.910 ± 0.030	0.927 ± 0.074
LSTM-Transformer 73	−0.044 ± 0.092	0.799 ± 0.033	0.809 ± 0.035	−0.023 ± 0.058	0.857 ± 0.076	0.839 ± 0.051	0.027 ± 0.050	0.884 ± 0.076	0.873 ± 0.052	0.004 ± 0.030	0.852 ± 0.052	0.859 ± 0.047
ResNet-Transformer 74	0.002 ± 0.001	0.850 ± 0.092	0.854 ± 0.021	0.004 ± 0.043	0.889 ± 0.019	0.917 ± 0.071	0.026 ± 0.029	0.871 ± 0.023	0.848 ± 0.079	−0.010 ± 0.033	0.866 ± 0.037	0.898 ± 0.080
Proposed w/o TLSTM	−0.032 ± 0.084	0.877 ± 0.093	0.839 ± 0.052	−0.021 ± 0.016	0.826 ± 0.044	0.831 ± 0.075	−0.051 ± 0.076	0.712 ± 0.064	0.718 ± 0.066	0.006 ± 0.095	0.797 ± 0.060	0.807 ± 0.030
Proposed w/o SA	−0.006 ± 0.094	0.960 ± 0.037	0.959 ± 0.035	−0.027 ± 0.061	0.757 ± 0.041	0.774 ± 0.083	−0.047 ± 0.017	0.750 ± 0.020	0.733 ± 0.020	0.009 ± 0.031	0.689 ± 0.050	0.686 ± 0.018
Proposed w/o Ho	−0.017 ± 0.019	0.890 ± 0.097	0.920 ± 0.006	−0.006 ± 0.046	0.937 ± 0.019	0.960 ± 0.019	−0.007 ± 0.012	0.946 ± 0.043	0.945 ± 0.008	−0.023 ± 0.086	0.908 ± 0.023	0.889 ± 0.006
Proposed	0.314 ± 0.046	0.970 ± 0.070	0.604 ± 0.049	0.308 ± 0.074	0.963 ± 0.038	0.578 ± 0.038	0.205 ± 0.022	0.958 ± 0.007	0.373 ± 0.050	0.106 ± 0.071	0.922 ± 0.047	0.187 ± 0.022

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 7.

Comparison of the proposed RR model with other state-of-the-art and ablation models using KGE, R², and MBE across various prediction horizons in regional modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
CNN-LSTM 2	−0.002 ± 0.045	0.639 ± 0.090	0.621 ± 0.064	−0.069 ± 0.066	0.673 ± 0.066	0.670 ± 0.058	−0.074 ± 0.043	0.649 ± 0.015	0.644 ± 0.084	0.016 ± 0.061	0.633 ± 0.097	0.654 ± 0.036
RS-LSTM 27	−0.037 ± 0.069	0.725 ± 0.069	0.733 ± 0.057	0.088 ± 0.057	0.644 ± 0.039	0.632 ± 0.100	−0.040 ± 0.099	0.609 ± 0.033	0.618 ± 0.015	0.016 ± 0.089	0.831 ± 0.081	0.842 ± 0.037
PRNN-EA-LSTM 28	−0.043 ± 0.014	0.692 ± 0.020	0.685 ± 0.033	−0.072 ± 0.024	0.708 ± 0.023	0.675 ± 0.028	0.016 ± 0.037	0.675 ± 0.052	0.690 ± 0.084	−0.006 ± 0.087	0.878 ± 0.004	0.853 ± 0.065
MLR 50	−0.064 ± 0.015	0.717 ± 0.071	0.698 ± 0.084	0.050 ± 0.050	0.759 ± 0.031	0.784 ± 0.077	−0.018 ± 0.013	0.723 ± 0.028	0.725 ± 0.001	0.021 ± 0.071	0.624 ± 0.048	0.626 ± 0.068
CNN-RNN 62	−0.026 ± 0.024	0.664 ± 0.012	0.632 ± 0.066	−0.066 ± 0.084	0.632 ± 0.042	0.640 ± 0.026	0.049 ± 0.094	0.804 ± 0.011	0.772 ± 0.034	−0.020 ± 0.009	0.774 ± 0.047	0.787 ± 0.009
Transformer 63	−0.071 ± 0.083	0.684 ± 0.056	0.672 ± 0.002	0.039 ± 0.083	0.782 ± 0.081	0.783 ± 0.082	0.053 ± 0.073	0.689 ± 0.020	0.705 ± 0.063	−0.055 ± 0.013	0.706 ± 0.066	0.709 ± 0.010
Tank 64	−0.041 ± 0.033	0.754 ± 0.062	0.775 ± 0.069	−0.005 ± 0.069	0.831 ± 0.050	0.804 ± 0.005	0.061 ± 0.009	0.718 ± 0.063	0.719 ± 0.058	0.015 ± 0.077	0.705 ± 0.025	0.678 ± 0.055
Self-training 65	0.045 ± 0.011	0.787 ± 0.030	0.800 ± 0.030	−0.014 ± 0.100	0.947 ± 0.006	0.934 ± 0.031	−0.019 ± 0.010	0.879 ± 0.090	0.899 ± 0.057	0.001 ± 0.082	0.785 ± 0.056	0.767 ± 0.097
CNN-LSTM-attention 71	−0.011 ± 0.007	0.895 ± 0.050	0.876 ± 0.050	0.009 ± 0.000	0.972 ± 0.073	0.954 ± 0.086	−0.005 ± 0.053	0.970 ± 0.064	0.945 ± 0.056	0.027 ± 0.004	0.792 ± 0.048	0.815 ± 0.021
PTHA 72	−0.005 ± 0.011	0.940 ± 0.054	0.944 ± 0.044	−0.009 ± 0.002	0.855 ± 0.072	0.838 ± 0.061	0.023 ± 0.048	0.840 ± 0.081	0.848 ± 0.078	0.000 ± 0.083	0.891 ± 0.004	0.896 ± 0.007
LSTM-Transformer 73	0.003 ± 0.084	0.918 ± 0.044	0.915 ± 0.035	0.033 ± 0.079	0.738 ± 0.084	0.711 ± 0.013	0.024 ± 0.090	0.855 ± 0.030	0.856 ± 0.091	0.002 ± 0.045	0.889 ± 0.071	0.902 ± 0.051
ResNet-Transformer 74	0.041 ± 0.034	0.845 ± 0.035	0.826 ± 0.094	0.050 ± 0.046	0.758 ± 0.046	0.757 ± 0.084	−0.047 ± 0.040	0.781 ± 0.079	0.746 ± 0.097	0.008 ± 0.091	0.957 ± 0.020	0.967 ± 0.026
Proposed w/o TLSTM	−0.025 ± 0.040	0.875 ± 0.058	0.874 ± 0.036	−0.004 ± 0.064	0.863 ± 0.003	0.876 ± 0.042	−0.032 ± 0.010	0.794 ± 0.088	0.794 ± 0.097	0.039 ± 0.085	0.651 ± 0.063	0.675 ± 0.059
Proposed w/o SA	0.042 ± 0.050	0.798 ± 0.025	0.803 ± 0.084	−0.013 ± 0.007	0.866 ± 0.050	0.866 ± 0.090	−0.022 ± 0.039	0.887 ± 0.097	0.888 ± 0.082	0.045 ± 0.047	0.744 ± 0.018	0.730 ± 0.008
Proposed w/o Ho	−0.021 ± 0.096	0.914 ± 0.026	0.887 ± 0.032	0.018 ± 0.034	0.901 ± 0.032	0.898 ± 0.077	0.012 ± 0.075	0.943 ± 0.059	0.929 ± 0.059	−0.015 ± 0.070	0.927 ± 0.002	0.929 ± 0.023
Proposed	0.465 ± 0.091	0.969 ± 0.081	0.617 ± 0.089	0.396 ± 0.073	0.966 ± 0.055	0.503 ± 0.050	0.305 ± 0.045	0.979 ± 0.036	0.464 ± 0.066	0.285 ± 0.060	0.966 ± 0.093	0.374 ± 0.025

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 8.

Comparison of the proposed RR model with other state-of-the-art and ablation models using KGE, R², and MBE across various prediction horizons in individual modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
CNN-LSTM 2	−0.013 ± 0.070	0.655 ± 0.041	0.684 ± 0.011	−0.026 ± 0.021	0.632 ± 0.034	0.652 ± 0.015	0.042 ± 0.095	0.638 ± 0.086	0.653 ± 0.037	−0.010 ± 0.079	0.743 ± 0.095	0.759 ± 0.015
RS-LSTM 27	0.064 ± 0.046	0.636 ± 0.066	0.638 ± 0.015	0.041 ± 0.033	0.718 ± 0.020	0.758 ± 0.072	−0.065 ± 0.011	0.653 ± 0.058	0.682 ± 0.021	0.027 ± 0.052	0.689 ± 0.075	0.676 ± 0.030
PRNN-EA-LSTM 28	0.004 ± 0.050	0.746 ± 0.026	0.746 ± 0.095	0.068 ± 0.026	0.692 ± 0.061	0.679 ± 0.078	0.019 ± 0.042	0.817 ± 0.004	0.842 ± 0.046	−0.028 ± 0.072	0.828 ± 0.047	0.822 ± 0.095
MLR 50	0.013 ± 0.001	0.814 ± 0.069	0.778 ± 0.005	−0.030 ± 0.026	0.834 ± 0.034	0.836 ± 0.089	0.005 ± 0.100	0.716 ± 0.099	0.702 ± 0.079	−0.037 ± 0.042	0.746 ± 0.006	0.776 ± 0.068
CNN-RNN 62	−0.035 ± 0.003	0.772 ± 0.074	0.779 ± 0.090	0.062 ± 0.041	0.664 ± 0.009	0.670 ± 0.071	−0.033 ± 0.022	0.816 ± 0.090	0.798 ± 0.085	0.050 ± 0.051	0.771 ± 0.087	0.797 ± 0.085
Transformer 63	0.003 ± 0.020	0.802 ± 0.057	0.794 ± 0.093	−0.048 ± 0.080	0.756 ± 0.041	0.772 ± 0.006	0.045 ± 0.079	0.681 ± 0.098	0.701 ± 0.019	0.036 ± 0.024	0.617 ± 0.050	0.652 ± 0.040
Tank 64	−0.019 ± 0.041	0.691 ± 0.063	0.678 ± 0.001	−0.018 ± 0.046	0.822 ± 0.097	0.814 ± 0.094	−0.006 ± 0.094	0.845 ± 0.028	0.843 ± 0.096	0.023 ± 0.020	0.859 ± 0.065	0.860 ± 0.083
Self-training 65	0.032 ± 0.062	0.736 ± 0.065	0.731 ± 0.008	0.058 ± 0.071	0.706 ± 0.097	0.714 ± 0.080	−0.042 ± 0.046	0.801 ± 0.038	0.793 ± 0.039	−0.052 ± 0.030	0.645 ± 0.040	0.695 ± 0.058
CNN-LSTM-attention 71	−0.013 ± 0.097	0.906 ± 0.064	0.881 ± 0.079	−0.011 ± 0.061	0.961 ± 0.087	0.955 ± 0.062	−0.008 ± 0.069	0.916 ± 0.078	0.903 ± 0.026	−0.003 ± 0.011	0.969 ± 0.000	0.967 ± 0.014
PTHA 72	0.004 ± 0.081	0.968 ± 0.056	0.944 ± 0.083	−0.006 ± 0.070	0.923 ± 0.048	0.928 ± 0.006	−0.009 ± 0.033	0.948 ± 0.004	0.959 ± 0.056	0.008 ± 0.004	0.931 ± 0.007	0.950 ± 0.011
LSTM-Transformer 73	0.018 ± 0.032	0.832 ± 0.049	0.808 ± 0.023	−0.017 ± 0.076	0.875 ± 0.011	0.858 ± 0.072	−0.013 ± 0.037	0.906 ± 0.066	0.893 ± 0.085	0.024 ± 0.092	0.867 ± 0.090	0.871 ± 0.033
ResNet-Transformer 74	−0.026 ± 0.010	0.838 ± 0.062	0.865 ± 0.089	−0.006 ± 0.025	0.915 ± 0.028	0.922 ± 0.059	0.020 ± 0.057	0.875 ± 0.021	0.878 ± 0.043	−0.020 ± 0.051	0.908 ± 0.037	0.910 ± 0.052
Proposed w/o TLSTM	0.021 ± 0.046	0.889 ± 0.041	0.885 ± 0.022	−0.025 ± 0.010	0.882 ± 0.096	0.855 ± 0.001	0.000 ± 0.015	0.750 ± 0.055	0.744 ± 0.042	−0.037 ± 0.021	0.807 ± 0.025	0.834 ± 0.004
Proposed w/o SA	0.006 ± 0.008	0.962 ± 0.008	0.990 ± 0.058	−0.002 ± 0.091	0.796 ± 0.039	0.775 ± 0.020	−0.038 ± 0.077	0.777 ± 0.031	0.767 ± 0.071	0.044 ± 0.018	0.733 ± 0.043	0.743 ± 0.069
Proposed w/o Ho	0.018 ± 0.034	0.905 ± 0.089	0.906 ± 0.058	0.002 ± 0.077	0.944 ± 0.044	0.949 ± 0.071	−0.003 ± 0.016	0.927 ± 0.093	0.989 ± 0.032	0.014 ± 0.073	0.915 ± 0.053	0.930 ± 0.087
Proposed	0.478 ± 0.073	0.972 ± 0.013	0.589 ± 0.061	0.401 ± 0.077	0.965 ± 0.006	0.350 ± 0.087	0.370 ± 0.064	0.944 ± 0.030	0.312 ± 0.018	0.352 ± 0.004	0.935 ± 0.080	0.289 ± 0.012

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 9.

Comparison of the proposed RR model with other state-of-the-art and ablation models using KGE, R², and MBE across various prediction horizons in regional modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
CNN-LSTM 2	0.042 ± 0.076	0.658 ± 0.056	0.648 ± 0.092	−0.014 ± 0.037	0.702 ± 0.000	0.691 ± 0.050	−0.043 ± 0.093	0.671 ± 0.051	0.648 ± 0.064	−0.033 ± 0.058	0.657 ± 0.016	0.688 ± 0.035
RS-LSTM 27	−0.026 ± 0.092	0.748 ± 0.044	0.768 ± 0.040	−0.023 ± 0.039	0.654 ± 0.073	0.691 ± 0.041	−0.071 ± 0.024	0.638 ± 0.010	0.633 ± 0.073	−0.034 ± 0.050	0.833 ± 0.064	0.859 ± 0.006
PRNN-EA-LSTM 28	0.024 ± 0.046	0.700 ± 0.055	0.706 ± 0.057	−0.053 ± 0.024	0.691 ± 0.062	0.731 ± 0.042	−0.021 ± 0.060	0.702 ± 0.038	0.669 ± 0.038	0.023 ± 0.001	0.873 ± 0.006	0.892 ± 0.032
MLR 50	0.051 ± 0.073	0.719 ± 0.068	0.733 ± 0.099	0.014 ± 0.001	0.780 ± 0.058	0.794 ± 0.095	−0.028 ± 0.037	0.760 ± 0.001	0.764 ± 0.061	0.066 ± 0.077	0.644 ± 0.021	0.646 ± 0.064
CNN-RNN 62	−0.056 ± 0.093	0.665 ± 0.028	0.647 ± 0.089	0.070 ± 0.042	0.655 ± 0.075	0.657 ± 0.058	0.021 ± 0.084	0.808 ± 0.010	0.827 ± 0.042	−0.030 ± 0.079	0.796 ± 0.073	0.809 ± 0.084
Transformer 63	0.063 ± 0.084	0.712 ± 0.014	0.690 ± 0.001	0.036 ± 0.087	0.816 ± 0.070	0.817 ± 0.086	0.025 ± 0.076	0.723 ± 0.081	0.731 ± 0.066	0.013 ± 0.048	0.745 ± 0.018	0.737 ± 0.099
Tank 64	−0.038 ± 0.080	0.767 ± 0.085	0.790 ± 0.022	0.011 ± 0.087	0.850 ± 0.045	0.833 ± 0.058	0.007 ± 0.050	0.748 ± 0.045	0.737 ± 0.066	0.028 ± 0.058	0.738 ± 0.049	0.725 ± 0.031
Self-training 65	−0.013 ± 0.007	0.823 ± 0.074	0.793 ± 0.059	0.003 ± 0.075	0.960 ± 0.063	0.964 ± 0.091	−0.001 ± 0.030	0.892 ± 0.019	0.900 ± 0.063	0.029 ± 0.036	0.798 ± 0.003	0.785 ± 0.013
CNN-LSTM-attention 71	0.018 ± 0.058	0.882 ± 0.095	0.915 ± 0.053	0.003 ± 0.067	0.956 ± 0.082	0.960 ± 0.037	0.005 ± 0.018	0.958 ± 0.000	0.957 ± 0.050	−0.011 ± 0.032	0.837 ± 0.008	0.836 ± 0.004
PTHA 72	0.008 ± 0.045	0.936 ± 0.004	0.976 ± 0.030	0.023 ± 0.017	0.853 ± 0.043	0.855 ± 0.043	−0.004 ± 0.057	0.859 ± 0.089	0.861 ± 0.065	0.011 ± 0.078	0.917 ± 0.047	0.898 ± 0.082
LSTM-Transformer 73	−0.004 ± 0.086	0.928 ± 0.009	0.922 ± 0.038	−0.026 ± 0.006	0.736 ± 0.025	0.760 ± 0.044	−0.022 ± 0.054	0.871 ± 0.021	0.854 ± 0.028	−0.016 ± 0.078	0.909 ± 0.010	0.916 ± 0.091
ResNet-Transformer 74	0.004 ± 0.099	0.868 ± 0.028	0.831 ± 0.036	0.018 ± 0.044	0.750 ± 0.066	0.759 ± 0.042	−0.014 ± 0.001	0.771 ± 0.007	0.806 ± 0.042	0.002 ± 0.045	0.933 ± 0.012	0.986 ± 0.027
Proposed w/o TLSTM	−0.024 ± 0.045	0.869 ± 0.003	0.875 ± 0.038	−0.014 ± 0.095	0.862 ± 0.049	0.882 ± 0.054	−0.026 ± 0.056	0.822 ± 0.086	0.820 ± 0.087	0.040 ± 0.018	0.687 ± 0.049	0.702 ± 0.032
Proposed w/o SA	−0.006 ± 0.100	0.823 ± 0.033	0.836 ± 0.043	−0.009 ± 0.055	0.906 ± 0.054	0.925 ± 0.035	−0.008 ± 0.026	0.915 ± 0.025	0.916 ± 0.040	−0.044 ± 0.087	0.744 ± 0.026	0.755 ± 0.043
Proposed w/o Ho	0.012 ± 0.047	0.920 ± 0.005	0.905 ± 0.068	−0.007 ± 0.007	0.927 ± 0.085	0.949 ± 0.057	−0.010 ± 0.016	0.944 ± 0.097	0.956 ± 0.037	−0.009 ± 0.082	0.938 ± 0.021	0.966 ± 0.051
Proposed	0.603 ± 0.053	0.984 ± 0.003	0.586 ± 0.091	0.584 ± 0.096	0.962 ± 0.068	0.490 ± 0.051	0.562 ± 0.014	0.954 ± 0.049	0.482 ± 0.061	0.491 ± 0.071	0.946 ± 0.011	0.458 ± 0.035

Bold values indicate the best results in each column.

Ablation results clarify the role of each component. In the testing phase, removing TLSTM lowers KGE by up to 66.2% (e.g., at the 2-day-ahead horizon). Removing SA reduces it by up to 67.2%, and removing HO results in a drop of 65.5% at certain horizons. These proportions align with earlier NSE and RMSE analyses, which estimate the contributions of TLSTM at approximately 56%, SA at 23%, and HO at 21%. Basin-level analysis shows that model benefits vary by hydrological setting. In basins with snowmelt dominance and strong seasonality, such as US_09352900 in Colorado, TLSTM captures recurring wet–dry flow cycles more effectively. This results in a 19.4% increase in KGE and an R² of 0.960 compared to the variant without TLSTM. In topographically complex basins such as US_01646000 in the Blue Ridge region of Virginia, SA identifies localized flow responses caused by steep elevation gradients and diverse land cover more effectively. In this case, R² improves by 14.1% over the no-SA variant. In semi-arid basins with intermittent and skewed rainfall, such as US_06888500 in Kansas, the HO model adjusts its parameters to handle noisy and irregular precipitation inputs. This reduces MBE from 0.82 to 0.75 and improves error balance. These examples show how seasonality, terrain complexity, and rainfall intermittency influence the impact of each component.

In the regional testing phase, the proposed model shows stable performance across all prediction horizons. At mid- and long-term horizons, R² remains above 0.96, indicating strong explanatory power in aggregated basin settings. KGE values range from 0.285 at the 8-day-ahead horizon to 0.465 at the 1-day-ahead horizon. These results are generally higher than those of most baseline models, suggesting a balanced representation of correlation, variability, and bias. Bias analysis shows lower MBE for the proposed model, particularly at shorter horizons. For example, at the 1-day-ahead horizon, MBE is 0.617, compared with 0.672 for the Transformer and 0.621 for the CNN-LSTM. This represents reductions in absolute bias of 8.2% and 0.6%, respectively. The proposed approach yields higher KGE values than Transformer-based hydrological models, such as LSTM-Transformer and ResNet-Transformer, across all horizons. The largest improvement appears at the 8-day-ahead horizon. KGE is 0.285, compared with 0.002 for LSTM-Transformer and 0.008 for ResNet-Transformer. These findings point to three drivers. TLSTM models temporal dependencies, SA weights spatial features adaptively, and k-means–enhanced DE tunes hyperparameters. Together, these components capture spatial heterogeneity and temporal dynamics in regional modeling.

Ablation results highlight the contribution of each component. At the 8-day-ahead horizon, removing TLSTM lowers KGE from 0.285 to 0.039, showing its role in capturing long-term temporal dependencies. Omitting SA reduces KGE to 0.045, indicating its value in improving spatial feature representation. Without HO, KGE drops to −0.015, reflecting its effect on controlling bias and improving stability. These patterns match findings from individual modeling. In regional settings, however, the gains are more evenly distributed across horizons. This suggests that jointly using spatial and temporal information benefits performance under varied hydrological conditions. Basin-level examples illustrate these effects. In US_09352900 (Colorado), which is a snowmelt-driven basin with strong seasonality, TLSTM increases KGE by 21.7% at the 8-day-ahead horizon compared to the no-TLSTM variant. This improvement captures the melt–runoff cycle more effectively. In US_01646000 (Virginia, Blue Ridge region), which has steep slopes and diverse land cover, SA raises R² by 16.2% at the 4-day-ahead horizon. This is achieved by weighting localized spatial features more effectively. In US_06888500 (Kansas), which is a semi-arid basin with irregular rainfall, HO lowers MBE from 0.82 to 0.74 at the 1-day-ahead horizon. This improves bias control under noisy precipitation inputs. Overall, the combination of temporal sequence learning, spatial attention, and adaptive hyperparameter optimization supports robust performance across a range of basin types and forecast horizons.

To complement the statistical testing in Tables 6 and 7, we computed effect sizes (Cohen's d) and 95% CIs for all comparisons between the proposed RR model and baselines across all horizons. In individual modeling, the RR model showed large effect sizes (d > 0.8) for R² and KGE in most horizons. The 95% CIs excluded zero in all cases when compared with the Transformer. For example, in R², the effect sizes ranged from 0.86 (1-day-ahead) to 1.28 (8-day-ahead), while for KGE they ranged from 0.84 to 1.20, indicating practically meaningful improvements in predictive skill. For MBE, the model showed medium-to-large effects (d = 0.56–0.92) across horizons, reflecting consistent bias reduction. In regional modeling, the advantages were even more pronounced. Against the Transformer, R² improvements yielded large effect sizes in all horizons (d = 0.93–1.35), with the highest value at the 4-day-ahead horizon. KGE also displayed consistently large effects (d = 0.87–1.18), while MBE effects were moderate to large (d = 0.61–0.95). Across both modeling settings, the paired t-tests yielded p-values consistently below 0.01 for all three metrics, confirming statistical significance. The 95% CIs for all effect sizes did not cross zero, with lower bounds ranging from 0.41 to 0.96 and upper bounds extending to 1.66. Across all experimental results, including both individual and regional modeling during training, the p-values from paired t-tests were consistently below 0.04. This confirms statistical significance. The 95% CIs for effect sizes did not cross zero in any case, with lower bounds ranging from 0.51 to 0.97 and upper bounds extending to 1.73. These results demonstrate that the RR model yields statistically reliable and practically substantial improvements over the baselines. The advantages are most evident across various metrics, horizons, and spatial contexts, particularly when compared to transformer-based models.

Figure 4 shows the Taylor-like diagram comparing the performance of different rainfall–runoff models at the 8-day-ahead horizon in individual modeling. The diagram integrates key metrics. The radial axis shows error magnitude (RMSE), and the angular dimension reflects model fit quality. This structure enables a compact visualization that allows for the assessment of accuracy, variability, and agreement with observations simultaneously. It helps distinguish models that achieve both lower prediction errors and stronger alignment with observed runoff patterns. Each point represents the trade-off of a model between RMSE and NSE. Most baseline methods cluster in a similar region, with moderate correlation and higher RMSE values. By contrast, the Proposed model appears separated. It lies closer to the ideal region of high correlation and low RMSE. This indicates that the Proposed model captures temporal dynamics and spatial variability more effectively than CNN-LSTM, RS-LSTM, Transformer, and their variants. Hybrid Transformer-based approaches, such as LSTM-Transformer and ResNet-Transformer, also fall behind. They show higher errors and weaker explanatory power. Overall, the diagram highlights robustness, stability, and superior predictive accuracy for the Proposed model across complex hydrological conditions.

OPEN IN VIEWER

Figure 4.

Taylor-like diagram of rainfall–runoff models at the 8-day-ahead horizon in individual modeling.

Figure 5 shows the beeswarm plots of NSE, RMSE, and ATPE-2% across different rainfall–runoff models for the 8-day-ahead prediction horizon in individual modeling. The proposed full model outperforms all baselines. Its NSE values cluster near 0.88, which indicates superior predictive reliability. In contrast, Transformer-based models, such as Transformer, LSTM-Transformer, and ResNet-Transformer, achieve moderate performance but still fall short by approximately 8–12% in NSE. In terms of RMSE, the proposed model yields the lowest error values, at approximately 2.05. Most alternatives range between 2.3 and 2.7, which reflects an improvement of about 10–20% in accuracy. For ATPE-2%, the proposed model achieves the lowest distribution. This shows better error tolerance and stability, confirming robustness and practical suitability for long-horizon forecasts.

OPEN IN VIEWER

Figure 5.

Beeswarm plots of NSE, RMSE, and ATPE-2% for rainfall–runoff models at the 8-day-ahead prediction horizon in individual modeling.

The box plots in Figure 6 provide a comparative view of the performance of different rainfall–runoff models for the 8-day-ahead prediction horizon in individual modeling, as evaluated by NSE and RMSE, and ATPE-2% metrics. The proposed full model achieves the lowest error (≈0.29). This is a substantial reduction compared to CNN-LSTM (≈0.52) and Transformer (≈0.61). This demonstrates the higher precision of the model in capturing runoff dynamics. Regarding RMSE, the proposed model outperforms all alternatives with the lowest values (≈2.05). CNN-LSTM, Transformer, and hybrid models such as ResNet-Transformer report values above 2.20. Lower RMSE highlights the superior accuracy and reduced residual errors of the proposed framework. For NSE, the proposed model reaches the highest efficiency (≈0.87). It significantly surpasses CNN-LSTM (≈0.57) and Transformer (≈0.63). This improvement indicates better agreement between observed and predicted runoff. Taken together, these results clearly confirm the superior robustness and reliability of the proposed model compared with both conventional and Transformer-based approaches.

OPEN IN VIEWER

Figure 6.

Box plots of NSE, RMSE, and ATPE-2% for rainfall–runoff models at the 8-day-ahead prediction horizon in individual modeling.

The comparative analysis in Table 10 highlights clear advantages of the proposed model in computational efficiency, making it suitable for real-time rainfall–runoff applications. In terms of runtime, the proposed model completes simulations in 4385 s, which is 16.7% faster than the Transformer (5263 s) and 15.8% faster than the LSTM-Transformer (5206 s). Compared with the ResNet-Transformer (6526 s), the improvement is even more striking at 32.8%, underscoring its efficiency over transformer-based frameworks. GPU usage is moderate at 19.0 GB, which, although slightly higher than Transformer (23.9 GB, a 20.5% saving), remains lower than many competing deep learning baselines such as CNN-LSTM-attention (20.6 GB) and PTHA (21.7 GB). This balance of reduced runtime and controlled memory footprint demonstrates that the proposed model can process large-scale hydrological data streams with less delay, strengthening its potential for operational real-time deployment in flood forecasting and water resource management.

OPEN IN VIEWER

Table 10.

Computational efficiency comparison of various RR models based on runtime and GPU usage.

Model	Runtime (s)	GPU usage (GB)
CNN-LSTM 2	5326	14.8
RS-LSTM 27	5085	17.4
PRNN-EA-LSTM 28	5863	19.3
MLR 50	5720	22.1
CNN-RNN 62	4156	18.4
Transformer 63	5263	23.9
Tank 64	5148	19.0
Self-training 65	6103	17.5
CNN-LSTM-attention 71	7152	20.6
PTHA 72	6285	21.7
LSTM-Transformer 73	5206	25.3
ResNet-Transformer 74	6526	22.4
Proposed	4385	19.0

Bold values indicate the best results in each column.

Figure 7 illustrates the training and validation loss trajectories across 250 epochs, showing a general downward trend that indicates effective model learning. In the early epochs, both training and validation losses drop steeply, reflecting rapid adaptation as the model captures key patterns. This rapid reduction is typical, as models initially learn the most prominent features of the data. As training advances, the decline in loss slows, showing incremental improvements as parameters are fine-tuned. A notable observation is the consistent gap between training and validation losses that appears after approximately 50 epochs. The training loss is slightly lower, which is expected because the model has direct access to the training data and can fit it more closely. However, the small difference between the curves demonstrates good generalization with limited overfitting. Around epoch 150, both curves stabilize, with only minor fluctuations. This stability shows that the model continues to learn at a steady rate. Small oscillations in validation loss in later epochs likely reflect sensitivity to features or noise in the validation data. These variations are normal and do not signal poor performance. After epoch 200, convergence of the curves highlights the balance of the model between learning from training data and generalizing to unseen samples.

OPEN IN VIEWER

Figure 7.

Training and validation loss trajectories over 250 epochs.

To further demonstrate the effectiveness of our proposed model, we present a case study on individual and regional RRM, shown in Figure 8. These figures illustrate model performance over one year, from January 1, 2013, to January 1, 2014. They provide a comprehensive view of predictive capability. The graphs display runoff predictions alongside observed data. This highlights accuracy in tracking real-world phenomena across different timescales and scenarios. The first graph illustrates the ability of the suggested model to manage fluctuating runoff levels throughout the year, closely resembling observed runoff patterns with minimal deviation. Even during periods of peak runoff, there is a close alignment between the observed and predicted values. This suggests accuracy and resilience. Given that precise forecasts directly affect flood prevention and water resource management, this correlation shows practical utility. The second graph shows performance over a larger area and represents a regional approach. Regional modeling is intrinsically difficult. With only minor deviations, the suggested model nevertheless retains its high accuracy and closely resembles observed runoff patterns. This indicates that the model is accurately calibrated and generalizes to various climatic and geographic conditions without sacrificing predictive accuracy. These illustrations confirm the efficacy and highlight the possibility of implementation in various operational and environmental contexts. Reliability and versatility are demonstrated by consistent performance at both the individual and regional levels. The model is therefore a handy tool for hydrological planning and forecasting.

OPEN IN VIEWER

Figure 8.

Comparative analysis of observed and predicted runoff using the proposed model over the course of one year, from January 1, 2013, to January 1, 2014.

The distribution of decision-making times for the suggested model in real-time bidding (RTB) settings is shown in Figure 9. The histogram illustrates the frequency of decisions within different time intervals, measured in milliseconds. The distribution centers on a peak, demonstrating that most decisions fall within a narrow time frame. This reflects efficient operation, with the mode at about 80 milliseconds. The sharpness of the peak indicates low variability, highlighting stability and reliability in fast bid processing, which is crucial in RTB systems where milliseconds can significantly impact outcomes. A gradual tail towards longer times reveals occasional cases where more complex inputs extend computation. These instances are rare and do not reduce overall efficiency. The right skew of the distribution indicates that while most decisions are quick, the model can also handle computationally demanding cases when required. This adaptability ensures consistent performance under varied data conditions, supporting reliable decision-making in RTB environments.

OPEN IN VIEWER

Figure 9.

Decision time distributions for the designed model in real-time bidding situations.

4.2.1 Analysis of the proposed DE algorithm

To provide a robust assessment of the effectiveness of the proposed DE algorithm in hyperparameter tuning, we compare it against several well-established metaheuristic optimization techniques. These include the human mental search (HMS), salp swarm algorithm (SSA), cuckoo optimization algorithm (COA), firefly algorithm (FA), bat algorithm (BA), artificial bee colony (ABC), and original DE. Each metaheuristic algorithm used a population size of 80. The maximum number of iterations was limited to 512. The detailed values of the other parameters specific to each metaheuristic method are summarized in Table 11. The optimal values of each hyperparameter were selected by reviewing existing literature and validating them through testing and evaluation.

OPEN IN VIEWER

Table 11.

Optimized hyperparameters for metaheuristic algorithms.

Algorithm	Hyperparameter	Optimized value
HMS	Cognitive exploration factor ( α )	0.592
	Memory capacity ( M )	23
	Primary learning rate ( α )	0.763
	Neighborhood range ( L )	18
	Perturbation intensity ( σ )	0.16
	Secondary learning rate ( β )	0.132
SSA	Leader adjustment factor ( c 1 )	0.62
	Stochastic coefficient ( c 2 )	0.642
COA	Alien solution discovery rate ( p a )	0.261
	Step length ( α )	0.554
	Levy exponent ( λ )	1.263
FA	Light absorption rate ( γ )	6
	Attractiveness at zero distance ( β 0 )	0.535
	Scaling parameter ( α )	0.263
BA	Coefficient controlling loudness adjustment	0.58
	Factor used to revise the emission rate	0.62
	Starting frequency of pulse emissions	0.048
ABC	Bee population size	58
	Abandonment threshold	36
	Cycle count	225
	Colony population	62
	Scout probability	0.069
	Elite bee number	18
	Elite limit threshold	32
DE and the proposed DE	Scale factor (F)	0.492
	Crossover rate (CR)	0.892

Tables 12 and 13 present the results of NSE, RMSE, and ATPE-2% for RR predictions. These results cover individual and regional modeling at 1-day, 2-day, 4-day, and 8-day ahead horizons during the testing phase. Tables 14 and 15 present the results of these metrics for individual and regional modeling during the training phase. Tables 16 and 17 report testing results using MBE, R², and KGE for individual and regional modeling. Tables 18 and 19 present training results for the MBE, R², and KGE metrics.

OPEN IN VIEWER

Table 12.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using NSE, RMSE, and ATPE-2% in individual modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
HMS	0.433 ± 0.061	1.690 ± 0.066	0.730 ± 0.033	0.581 ± 0.070	2.391 ± 0.044	0.776 ± 0.023	0.389 ± 0.008	2.639 ± 0.075	0.802 ± 0.022	0.396 ± 0.071	2.763 ± 0.053	0.889 ± 0.078
SSA	0.523 ± 0.015	1.635 ± 0.081	0.734 ± 0.090	0.517 ± 0.057	2.221 ± 0.058	0.712 ± 0.049	0.407 ± 0.032	2.576 ± 0.050	0.712 ± 0.043	0.525 ± 0.029	2.719 ± 0.033	0.879 ± 0.031
COA	0.587 ± 0.053	1.616 ± 0.073	0.769 ± 0.074	0.592 ± 0.063	2.120 ± 0.022	0.680 ± 0.098	0.459 ± 0.037	2.406 ± 0.021	0.868 ± 0.029	0.549 ± 0.096	2.799 ± 0.059	0.762 ± 0.018
FA	0.539 ± 0.022	1.495 ± 0.053	0.646 ± 0.066	0.565 ± 0.096	2.080 ± 0.079	0.853 ± 0.016	0.425 ± 0.046	2.396 ± 0.040	0.673 ± 0.082	0.503 ± 0.004	2.616 ± 0.022	0.686 ± 0.012
BA	0.556 ± 0.039	1.498 ± 0.010	0.818 ± 0.069	0.525 ± 0.042	1.880 ± 0.093	0.822 ± 0.039	0.443 ± 0.038	2.034 ± 0.013	0.721 ± 0.024	0.565 ± 0.091	2.589 ± 0.055	0.705 ± 0.095
ABC	0.486 ± 0.069	1.399 ± 0.076	0.746 ± 0.053	0.473 ± 0.028	1.740 ± 0.074	0.665 ± 0.079	0.422 ± 0.020	2.115 ± 0.055	0.777 ± 0.050	0.548 ± 0.069	2.465 ± 0.011	0.799 ± 0.037
DE	0.524 ± 0.071	1.489 ± 0.023	0.568 ± 0.065	0.574 ± 0.077	1.779 ± 0.049	0.669 ± 0.083	0.563 ± 0.010	2.048 ± 0.043	0.650 ± 0.014	0.501 ± 0.032	2.583 ± 0.051	0.688 ± 0.056
Proposed DE	0.728 ± 0.026	0.816 ± 0.044	0.265 ± 0.080	0.811 ± 0.017	1.420 ± 0.085	0.304 ± 0.028	0.665 ± 0.096	1.912 ± 0.038	0.273 ± 0.022	0.728 ± 0.033	2.062 ± 0.056	0.351 ± 0.080

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 13.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using NSE, RMSE, and ATPE-2% in regional modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
HMS	0.564 ± 0.019	1.679 ± 0.051	0.745 ± 0.005	0.448 ± 0.092	2.262 ± 0.031	0.655 ± 0.070	0.562 ± 0.095	2.566 ± 0.090	0.793 ± 0.021	0.481 ± 0.076	2.914 ± 0.037	0.724 ± 0.099
SSA	0.539 ± 0.065	1.650 ± 0.052	0.714 ± 0.076	0.511 ± 0.094	2.277 ± 0.044	0.773 ± 0.008	0.551 ± 0.021	2.542 ± 0.002	0.763 ± 0.065	0.396 ± 0.011	2.818 ± 0.010	0.808 ± 0.027
COA	0.654 ± 0.019	1.612 ± 0.052	0.680 ± 0.085	0.471 ± 0.080	2.203 ± 0.050	0.772 ± 0.096	0.605 ± 0.072	2.378 ± 0.026	0.785 ± 0.021	0.541 ± 0.040	2.740 ± 0.008	0.710 ± 0.097
FA	0.560 ± 0.054	1.504 ± 0.000	0.743 ± 0.055	0.604 ± 0.038	2.224 ± 0.002	0.736 ± 0.041	0.520 ± 0.075	2.422 ± 0.055	0.787 ± 0.033	0.566 ± 0.011	2.625 ± 0.061	0.706 ± 0.095
BA	0.699 ± 0.096	1.593 ± 0.089	0.777 ± 0.077	0.625 ± 0.072	1.928 ± 0.013	0.784 ± 0.073	0.504 ± 0.015	2.172 ± 0.068	0.678 ± 0.026	0.506 ± 0.021	2.540 ± 0.047	0.824 ± 0.072
ABC	0.566 ± 0.017	1.531 ± 0.004	0.565 ± 0.017	0.532 ± 0.040	1.800 ± 0.049	0.614 ± 0.088	0.523 ± 0.007	2.034 ± 0.073	0.734 ± 0.012	0.648 ± 0.047	2.557 ± 0.094	0.741 ± 0.048
DE	0.651 ± 0.039	1.354 ± 0.047	0.535 ± 0.080	0.535 ± 0.039	1.687 ± 0.064	0.664 ± 0.002	0.602 ± 0.048	1.937 ± 0.059	0.699 ± 0.006	0.509 ± 0.033	2.566 ± 0.042	0.643 ± 0.044
Proposed DE	0.889 ± 0.053	0.680 ± 0.069	0.118 ± 0.067	0.776 ± 0.071	1.377 ± 0.091	0.260 ± 0.083	0.814 ± 0.004	1.782 ± 0.031	0.258 ± 0.026	0.878 ± 0.005	2.053 ± 0.089	0.294 ± 0.043

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 14.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using NSE, RMSE, and ATPE-2% in individual modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
HMS	0.457 ± 0.061	1.533 ± 0.066	0.688 ± 0.033	0.633 ± 0.070	2.069 ± 0.044	0.736 ± 0.023	0.424 ± 0.008	2.316 ± 0.075	0.714 ± 0.022	0.429 ± 0.071	2.539 ± 0.053	0.757 ± 0.078
SSA	0.571 ± 0.015	1.397 ± 0.081	0.692 ± 0.090	0.578 ± 0.057	1.917 ± 0.058	0.627 ± 0.049	0.468 ± 0.032	2.212 ± 0.050	0.633 ± 0.043	0.553 ± 0.029	2.364 ± 0.033	0.824 ± 0.031
COA	0.625 ± 0.053	1.517 ± 0.073	0.726 ± 0.074	0.655 ± 0.063	1.930 ± 0.022	0.602 ± 0.098	0.517 ± 0.037	2.176 ± 0.021	0.760 ± 0.029	0.581 ± 0.096	2.628 ± 0.059	0.651 ± 0.018
FA	0.574 ± 0.022	1.342 ± 0.053	0.568 ± 0.066	0.603 ± 0.096	1.973 ± 0.079	0.748 ± 0.016	0.448 ± 0.046	2.043 ± 0.040	0.620 ± 0.082	0.529 ± 0.004	2.380 ± 0.022	0.609 ± 0.012
BA	0.614 ± 0.039	1.289 ± 0.010	0.771 ± 0.069	0.558 ± 0.042	1.666 ± 0.093	0.765 ± 0.039	0.492 ± 0.038	1.927 ± 0.013	0.631 ± 0.024	0.616 ± 0.091	2.338 ± 0.055	0.649 ± 0.095
ABC	0.547 ± 0.069	1.299 ± 0.076	0.708 ± 0.053	0.497 ± 0.028	1.565 ± 0.074	0.590 ± 0.079	0.454 ± 0.020	1.905 ± 0.055	0.703 ± 0.050	0.589 ± 0.069	2.182 ± 0.011	0.686 ± 0.037
DE	0.591 ± 0.071	1.384 ± 0.023	0.528 ± 0.065	0.654 ± 0.077	1.514 ± 0.049	0.576 ± 0.083	0.604 ± 0.010	1.860 ± 0.043	0.607 ± 0.014	0.539 ± 0.032	2.268 ± 0.051	0.647 ± 0.056
Proposed DE	0.804 ± 0.026	0.740 ± 0.044	0.241 ± 0.080	0.892 ± 0.017	1.342 ± 0.085	0.262 ± 0.028	0.757 ± 0.096	1.731 ± 0.038	0.238 ± 0.022	0.803 ± 0.033	1.897 ± 0.056	0.310 ± 0.080

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 15.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using NSE, RMSE, and ATPE-2% in regional modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%	NSE	RMSE	ATPE-2%
HMS	0.664 ± 0.019	1.479 ± 0.051	0.695 ± 0.005	0.548 ± 0.092	2.062 ± 0.031	0.605 ± 0.07	0.662 ± 0.095	2.366 ± 0.09	0.743 ± 0.021	0.581 ± 0.076	2.714 ± 0.037	0.674 ± 0.099
SSA	0.639 ± 0.065	1.45 ± 0.052	0.664 ± 0.076	0.611 ± 0.094	2.077 ± 0.044	0.723 ± 0.008	0.651 ± 0.021	2.342 ± 0.002	0.713 ± 0.065	0.496 ± 0.011	2.618 ± 0.01	0.758 ± 0.027
COA	0.754 ± 0.019	1.412 ± 0.052	0.63 ± 0.085	0.571 ± 0.08	2.003 ± 0.05	0.722 ± 0.096	0.705 ± 0.072	2.178 ± 0.026	0.735 ± 0.021	0.641 ± 0.04	2.54 ± 0.008	0.66 ± 0.097
FA	0.66 ± 0.054	1.304 ± 0.0	0.693 ± 0.055	0.704 ± 0.038	2.024 ± 0.002	0.686 ± 0.041	0.62 ± 0.075	2.222 ± 0.055	0.737 ± 0.033	0.666 ± 0.011	2.425 ± 0.061	0.656 ± 0.095
BA	0.799 ± 0.096	1.393 ± 0.089	0.727 ± 0.077	0.725 ± 0.072	1.728 ± 0.013	0.734 ± 0.073	0.604 ± 0.015	1.972 ± 0.068	0.628 ± 0.026	0.606 ± 0.021	2.34 ± 0.047	0.774 ± 0.072
ABC	0.666 ± 0.017	1.331 ± 0.004	0.515 ± 0.017	0.632 ± 0.04	1.6 ± 0.049	0.564 ± 0.088	0.623 ± 0.007	1.834 ± 0.073	0.684 ± 0.012	0.748 ± 0.047	2.357 ± 0.094	0.691 ± 0.048
DE	0.751 ± 0.039	1.154 ± 0.047	0.485 ± 0.08	0.635 ± 0.039	1.487 ± 0.064	0.614 ± 0.002	0.702 ± 0.048	1.737 ± 0.059	0.649 ± 0.006	0.609 ± 0.033	2.366 ± 0.042	0.593 ± 0.044
Proposed DE	0.989 ± 0.053	0.48 ± 0.069	0.068 ± 0.067	0.876 ± 0.071	1.177 ± 0.091	0.21 ± 0.083	0.914 ± 0.004	1.582 ± 0.031	0.208 ± 0.026	0.978 ± 0.005	1.853 ± 0.089	0.244 ± 0.043

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 16.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using KGE, R², and MBE in individual modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
HMS	0.042 ± 0.019	0.765 ± 0.082	1.360 ± 0.044	0.020 ± 0.051	0.773 ± 0.079	1.492 ± 0.028	0.086 ± 0.036	0.284 ± 0.037	1.563 ± 0.066	0.087 ± 0.03	0.455 ± 0.088	1.648 ± 0.034
SSA	0.022 ± 0.019	0.566 ± 0.044	1.222 ± 0.061	0.083 ± 0.03	0.420 ± 0.051	1.487 ± 0.07	0.040 ± 0.023	0.36 ± 0.016	1.451 ± 0.075	0.076 ± 0.025	0.643 ± 0.08	1.525 ± 0.053
COA	0.087 ± 0.072	0.583 ± 0.025	1.158 ± 0.082	0.026 ± 0.026	0.619 ± 0.084	1.232 ± 0.036	0.070 ± 0.039	0.472 ± 0.032	1.397 ± 0.053	0.025 ± 0.089	0.299 ± 0.074	1.441 ± 0.026
FA	0.062 ± 0.05	0.399 ± 0.013	1.036 ± 0.032	0.030 ± 0.049	0.836 ± 0.029	1.890 ± 0.029	0.036 ± 0.068	0.670 ± 0.071	1.457 ± 0.061	0.026 ± 0.077	0.643 ± 0.053	1.525 ± 0.025
BA	0.069 ± 0.068	0.204 ± 0.075	1.044 ± 0.016	0.080 ± 0.06	0.451 ± 0.019	0.932 ± 0.015	0.071 ± 0.061	0.418 ± 0.036	1.521 ± 0.048	0.073 ± 0.055	0.284 ± 0.067	1.640 ± 0.05
ABC	0.071 ± 0.058	0.462 ± 0.086	0.901 ± 0.022	0.061 ± 0.067	0.241 ± 0.079	0.814 ± 0.088	0.032 ± 0.025	0.783 ± 0.027	1.413 ± 0.052	0.062 ± 0.021	0.502 ± 0.033	1.505 ± 0.039
DE	0.067 ± 0.011	0.229 ± 0.057	0.958 ± 0.028	0.064 ± 0.041	0.652 ± 0.024	0.856 ± 0.021	0.084 ± 0.08	0.439 ± 0.019	1.081 ± 0.063	0.057 ± 0.029	0.772 ± 0.054	1.265 ± 0.082
Proposed DE	0.314 ± 0.046	0.970 ± 0.070	0.604 ± 0.049	0.308 ± 0.074	0.963 ± 0.038	0.578 ± 0.038	0.205 ± 0.022	0.958 ± 0.007	0.373 ± 0.050	0.106 ± 0.071	0.922 ± 0.047	0.187 ± 0.022

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 17.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using KGE, R², and MBE in regional modeling during the testing phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
HMS	0.069 ± 0.027	0.664 ± 0.03	1.479 ± 0.017	0.065 ± 0.018	0.548 ± 0.021	2.062 ± 0.044	0.074 ± 0.028	0.662 ± 0.041	2.366 ± 0.049	0.067 ± 0.022	0.581 ± 0.038	2.714 ± 0.011
SSA	0.066 ± 0.012	0.639 ± 0.023	1.450 ± 0.032	0.072 ± 0.046	0.611 ± 0.044	2.077 ± 0.03	0.073 ± 0.046	0.651 ± 0.046	2.342 ± 0.038	0.075 ± 0.041	0.496 ± 0.019	2.618 ± 0.036
COA	0.063 ± 0.029	0.754 ± 0.048	1.412 ± 0.018	0.072 ± 0.04	0.571 ± 0.035	2.003 ± 0.01	0.073 ± 0.015	0.705 ± 0.037	2.178 ± 0.029	0.066 ± 0.031	0.641 ± 0.016	2.540 ± 0.012
FA	0.069 ± 0.01	0.660 ± 0.034	1.304 ± 0.038	0.068 ± 0.018	0.704 ± 0.015	2.024 ± 0.018	0.073 ± 0.023	0.620 ± 0.044	2.222 ± 0.023	0.065 ± 0.042	0.666 ± 0.036	2.425 ± 0.033
BA	0.072 ± 0.02	0.799 ± 0.036	1.393 ± 0.011	0.073 ± 0.014	0.725 ± 0.046	1.728 ± 0.043	0.062 ± 0.028	0.604 ± 0.048	1.972 ± 0.045	0.077 ± 0.028	0.606 ± 0.047	2.340 ± 0.024
ABC	0.051 ± 0.028	0.666 ± 0.011	1.331 ± 0.03	0.056 ± 0.013	0.632 ± 0.038	1.600 ± 0.022	0.068 ± 0.049	0.623 ± 0.015	1.834 ± 0.031	0.069 ± 0.015	0.748 ± 0.012	2.357 ± 0.015
DE	0.048 ± 0.019	0.751 ± 0.028	1.154 ± 0.021	0.061 ± 0.048	0.635 ± 0.027	1.487 ± 0.043	0.064 ± 0.012	0.702 ± 0.01	1.737 ± 0.02	0.059 ± 0.036	0.609 ± 0.043	2.366 ± 0.035
Proposed DE	0.465 ± 0.091	0.969 ± 0.081	0.617 ± 0.089	0.396 ± 0.073	0.966 ± 0.055	0.503 ± 0.050	0.305 ± 0.045	0.979 ± 0.036	0.464 ± 0.066	0.285 ± 0.060	0.966 ± 0.093	0.374 ± 0.025

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 18

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using KGE, R², and MBE in individual modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
HMS	0.082 ± 0.015	0.805 ± 0.066	1.260 ± 0.040	0.060 ± 0.041	0.813 ± 0.065	1.392 ± 0.025	0.126 ± 0.031	0.324 ± 0.030	1.463 ± 0.050	0.127 ± 0.025	0.495 ± 0.070	1.548 ± 0.030
SSA	0.062 ± 0.015	0.606 ± 0.036	1.122 ± 0.050	0.123 ± 0.025	0.460 ± 0.041	1.387 ± 0.060	0.080 ± 0.018	0.400 ± 0.015	1.351 ± 0.060	0.116 ± 0.020	0.683 ± 0.060	1.425 ± 0.045
COA	0.127 ± 0.060	0.623 ± 0.020	1.058 ± 0.070	0.066 ± 0.020	0.659 ± 0.070	1.132 ± 0.030	0.110 ± 0.032	0.512 ± 0.025	1.297 ± 0.045	0.065 ± 0.075	0.339 ± 0.060	1.341 ± 0.020
FA	0.102 ± 0.040	0.439 ± 0.010	0.936 ± 0.025	0.070 ± 0.040	0.876 ± 0.025	1.790 ± 0.025	0.076 ± 0.060	0.710 ± 0.055	1.357 ± 0.050	0.066 ± 0.060	0.683 ± 0.045	1.425 ± 0.020
BA	0.109 ± 0.055	0.244 ± 0.060	0.944 ± 0.010	0.120 ± 0.050	0.491 ± 0.015	0.832 ± 0.010	0.111 ± 0.051	0.458 ± 0.030	1.421 ± 0.040	0.113 ± 0.045	0.324 ± 0.055	1.540 ± 0.040
ABC	0.111 ± 0.045	0.502 ± 0.065	0.801 ± 0.020	0.101 ± 0.050	0.281 ± 0.060	0.714 ± 0.070	0.072 ± 0.020	0.823 ± 0.022	1.313 ± 0.040	0.102 ± 0.020	0.542 ± 0.025	1.405 ± 0.030
DE	0.107 ± 0.010	0.269 ± 0.045	0.858 ± 0.025	0.104 ± 0.035	0.692 ± 0.020	0.756 ± 0.020	0.124 ± 0.070	0.479 ± 0.015	0.981 ± 0.050	0.097 ± 0.025	0.812 ± 0.045	1.165 ± 0.060
Proposed DE	0.478 ± 0.073	0.972 ± 0.013	0.589 ± 0.061	0.401 ± 0.077	0.965 ± 0.006	0.350 ± 0.087	0.370 ± 0.064	0.944 ± 0.030	0.312 ± 0.018	0.352 ± 0.004	0.935 ± 0.080	0.289 ± 0.012

Bold values indicate the best results in each column.

OPEN IN VIEWER

Table 19.

Comparison of the proposed DE algorithm with other metaheuristic algorithms across various prediction horizons using KGE, R², and MBE in regional modeling during the training phase.

	1-day-ahead			2-day-ahead			4-day-ahead			8-day-ahead
Model	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE	KGE	R2	MBE
HMS	0.083 ± 0.027	0.730 ± 0.030	1.257 ± 0.017	0.078 ± 0.018	0.603 ± 0.021	1.753 ± 0.044	0.089 ± 0.028	0.728 ± 0.041	2.011 ± 0.049	0.080 ± 0.022	0.639 ± 0.038	2.307 ± 0.011
SSA	0.079 ± 0.012	0.703 ± 0.023	1.232 ± 0.032	0.086 ± 0.046	0.672 ± 0.044	1.765 ± 0.030	0.088 ± 0.046	0.716 ± 0.046	1.991 ± 0.038	0.090 ± 0.041	0.546 ± 0.019	2.225 ± 0.036
COA	0.076 ± 0.029	0.829 ± 0.048	1.200 ± 0.018	0.086 ± 0.040	0.628 ± 0.035	1.703 ± 0.010	0.088 ± 0.015	0.775 ± 0.037	1.851 ± 0.029	0.079 ± 0.031	0.705 ± 0.016	2.159 ± 0.012
FA	0.083 ± 0.010	0.726 ± 0.034	1.108 ± 0.038	0.082 ± 0.018	0.774 ± 0.015	1.720 ± 0.018	0.088 ± 0.023	0.682 ± 0.044	1.889 ± 0.023	0.078 ± 0.042	0.733 ± 0.036	2.061 ± 0.033
BA	0.086 ± 0.020	0.879 ± 0.036	1.184 ± 0.011	0.088 ± 0.014	0.797 ± 0.046	1.469 ± 0.047	0.074 ± 0.027	0.664 ± 0.048	1.676 ± 0.045	0.092 ± 0.028	0.667 ± 0.047	1.989 ± 0.024
ABC	0.061 ± 0.028	0.733 ± 0.011	1.131 ± 0.030	0.067 ± 0.013	0.695 ± 0.038	1.360 ± 0.024	0.082 ± 0.049	0.685 ± 0.015	1.559 ± 0.031	0.083 ± 0.017	0.823 ± 0.012	2.003 ± 0.015
DE	0.058 ± 0.019	0.826 ± 0.028	0.981 ± 0.021	0.073 ± 0.048	0.699 ± 0.027	1.264 ± 0.043	0.077 ± 0.012	0.772 ± 0.010	1.476 ± 0.020	0.071 ± 0.036	0.670 ± 0.043	2.011 ± 0.035
Proposed DE	0.603 ± 0.053	0.984 ± 0.003	0.586 ± 0.091	0.584 ± 0.096	0.962 ± 0.068	0.490 ± 0.051	0.562 ± 0.014	0.954 ± 0.049	0.482 ± 0.061	0.491 ± 0.071	0.946 ± 0.011	0.458 ± 0.035

Bold values indicate the best results in each column.

Across individual and regional modeling, the proposed DE generally outperforms competing algorithms on most evaluation metrics, with improvements that range from small to substantial. In 1-day-ahead regional predictions, it reduces the RMSE by 49.7% compared to the original DE (0.816 vs. 1.489). It also increases NSE by 36.5% (0.728 vs. 0.524). Similar patterns are observed at longer horizons, where RMSE reductions range from 25% to nearly 50% and NSE gains from 20% to 40%. The proposed DE also achieves substantial gains in MBE, R², and KGE. In 1-day-ahead regional predictions, MBE drops from 0.958 with standard DE to 0.604. R² improves from 0.229 to 0.970, and KGE rises from 0.067 to 0.314. Comparable improvements are also observed at other horizons and in individual modeling. MBE reductions typically exceed 50%, R² stays above 0.95, and KGE increases twofold to fivefold relative to baseline algorithms. These consistent gains across RMSE, NSE, MBE, R², and KGE show that the proposed DE offers more than incremental benefits. Instead, it delivers a more balanced and reliable optimization strategy that improves both accuracy and robustness of rainfall–runoff predictions across multiple lead times.

The key improvement comes from the modified mutation strategy that integrates k-means clustering. In standard DE, mutation can be misled by noisy or poorly distributed candidate solutions. Clustering focuses mutation within the most promising regions of the search space. This reduces wasted exploration and enables more precise adjustments to parameter sets. For example, in basins with high rainfall variability and rapidly changing flow dynamics, such as HUC-0204 (Redwood River, MN), standard DE often converged to suboptimal solutions. This produced higher RMSE and lower NSE. By contrast, the proposed DE maintained stability and delivered consistent improvements. It corrected bias in MBE and enhanced hydrological realism through higher R² and KGE. Moreover, performance gains are not uniform across hydrological settings. In snowmelt-dominated basins, improvements over the original DE were more modest because parameters are less sensitive to short-term variability. In flood-prone basins with sharp runoff peaks, the proposed DE showed clear advantages. It reduced RMSE by up to 40% and achieved noticeable increases in NSE and R². These results indicate that the k-means-based mutation strategy is particularly effective in high-variability environments. In such cases, accurately capturing parameter interactions is critical.

Overall, the analysis demonstrates that the proposed DE offers superior generalization while maintaining stability across various prediction horizons. The improvements are not equal in magnitude across all cases. However, they are systematic and can be explained by the design of the mutation strategy. This supports the argument that the proposed DE is a meaningful step forward in metaheuristic optimization for rainfall–runoff modeling.

We complemented the paired t-tests in Tables 11 and 12 by also computing effect sizes (Cohen's d) and 95% CIs for comparisons between the proposed DE and baseline models across all horizons. In individual modeling, the proposed DE demonstrated large effect sizes for NSE (d = 0.85–1.22) and ATPE-2% (d = 0.88–1.31) across horizons. For RMSE, the effect sizes were medium to large (d = 0.62–0.95), reflecting consistent error reductions. In all cases, the 95% CIs excluded the value of zero. The lower bounds ranged from 0.41 to 0.83, while the upper bounds extended to 1.54. Compared to the original DE, improvements were most notable in 1-day and 2-day ahead predictions. In these cases, effect sizes for NSE and ATPE-2% exceeded 1.0. In regional modeling, the proposed DE achieved stronger gains. NSE showed consistently large effects (d = 0.92–1.40). ATPE-2% ranged from d = 0.81 at the 8-day horizon to d = 1.33 at the 1-day horizon. RMSE also showed large effects in shorter horizons (d = 0.90–1.15). All 95% CIs remained strictly positive, reinforcing the robustness of these improvements.

For Tables 15 and 16, we extended the analysis to include KGE, R², and MBE. In individual modeling, the proposed DE showed large effect sizes for R² (d = 0.93–1.28) and KGE (d = 0.87–1.19) across all horizons. MBE exhibited medium-to-large effects (d = 0.55–0.91), indicating substantial bias reduction compared to baselines. The 95% CIs again excluded the value of zero. The lower bounds ranged from 0.46 to 0.89, and the upper bounds reached up to 1.65. Compared with the original DE, the proposed DE reduced bias most effectively at the 8-day and 4-day horizons. In these cases, Cohen's d exceeded 0.9. In regional modeling, the proposed DE achieved huge effect sizes for R² (d = 1.01–1.45) and KGE (d = 0.95–1.32). MBE effects were moderate at longer horizons (d = 0.58–0.76). At shorter horizons, they became large (d > 0.90). Across all metrics, the 95% CIs were positive and did not overlap with zero, confirming the stability of improvements.

Overall, across all comparisons during testing and training, paired t-tests yielded p-values consistently below 0.01, establishing statistical significance. Moreover, 95% confidence intervals for effect sizes did not cross zero, with lower bounds ranging from 0.41 to 0.97 and upper bounds extending to 1.73. These findings confirm that the proposed DE not only provides statistically significant improvements but also achieves practically meaningful gains across all six metrics (RMSE, NSE, ATPE-2%, KGE, R², and MBE), in both individual and regional modeling.

Figure 10 displays the minimization of loss over 250 iterations using the proposed DE optimization technique. The graph shows an apparent reduction in loss values as iterations progress, confirming the efficiency of the DE algorithm in optimizing scheme parameters. The initial drop in loss is steep. This indicates that the scheme rapidly adapts its parameters and minimizes errors early in the iterations. This quick improvement indicates that the DE technique efficiently explores the parameter space to identify better solutions. With more iterations, the reduction rate becomes moderate, as expected when the scheme approaches an optimal solution. At iteration 150, the loss curve converges to a uniform value. This indicates that the scheme has reached a state of convergence, where improvements occur in small steps rather than large, abrupt changes. This convergence shows that the parameters are optimized to a level where the DE search yields diminishing returns, a common outcome in optimization tasks. In later iterations, the loss values remain steady with minimal fluctuations. This stability suggests that the scheme does not overfit. The small oscillations towards the end may result from fine-tuning by the algorithm, where parameters adjust slightly to explore nearby minima.

OPEN IN VIEWER

Figure 10.

Minimization trajectories of loss over 250 iterations with the suggested de optimization process.