One-month-ahead river discharge forecasting using machine learning and limited hydrometeorological data: the Karatoya-Atrai River basin, Bangladesh

1 Introduction

There has been an increase in the intensity of the hydrological cycle due to climate change and global warming (IPCC, 2013), leading to changes in the timing and distribution of water resources (Kundzewicz et al., 2008). The increasing shortage of fresh and uncontaminated water during the dry season is a paramount challenge confronting humanity in the 21st century (Raihan et al., 2021; Simonovic, 2003; Zhu et al., 2021). On the other hand, flooding is one of the most challenging disasters during the wet season. Climate change is projected to intensify the frequency, duration, and variability of floods and droughts, necessitating urgent advances in water resource management (Hanak and Lund, 2012; Hong et al., 2021). The climatic system of our planet relies heavily on rivers, which play a crucial role in maintaining the connection between water and the atmosphere (Meybeck, 2003; Vörösmarty et al., 2000). In many regions, rivers are the only accessible water source, hindering socioeconomic development where monitoring infrastructure is lacking (Sullivan, 2002). Therefore, monitoring and quantifying river flow is vital for sustainability and future projections (Zakharova et al., 2020).

Climate change poses an added challenge for water resource management in South Asian nations, affecting various sectors of the economy in both the short and long-term (Shaw et al., 2013). While these challenges are global, Bangladesh exemplifies the acute impacts of these challenges due to its unique geographical and socioeconomic vulnerabilities. Bangladesh is highly prone to disasters, ranking fifth globally in disaster occurrence (Dastagir, 2015; Islam et al., 2021). Floods and droughts are the most challenging and common disasters in the country (Ahmed et al., 2021; Barua et al., 2016; Rahman, 2019; Zhang et al., 2019), due to being in a tropical monsoon climate with a high-density population and heavy reliance on natural resources (Zaman and Mondal, 2020). There has been a rise in salinity levels and a decrease in water flow during the dry season (Agrawala et al., 2003; Raihan et al., 2020). These disasters lead to significant crop losses, resulting in food insecurity that can vary by region and season, loss of livelihoods, crises related to migration, psychological issues, and environmental harm (Azad et al., 2013; Barua et al., 2016; Haque et al., 2019; Islam et al., 2014; Mamun et al., 2021; Mardy et al., 2018). Every year, millions of individuals lose their jobs, property, and social standing because of disasters that happen regularly (Alam et al., 2018; Habiba et al., 2012). For developing countries such as Bangladesh, managing transboundary water distribution is crucial for flood and drought mitigation. Northern and northwestern Bangladesh face annual challenges from floods and droughts of varying severity. (Rahaman et al., 2016). The Karatoya-Atrai River plays a pivotal role in draining a significant portion of the northwestern region of Bangladesh. Although the river consistently provides fishing opportunities, it frequently contributes to flooding in numerous regions during the monsoon season and droughts in the dry season. Therefore, it is important to properly manage the water of this river basin in all seasons.

Hydrology and water resource management rely on accurate streamflow/discharge data. The flow of streams is indicative of the overall condition of a watershed, which encompasses factors such as weather patterns, climate conditions, vegetation, soil properties, and land utilization (Rudra and Alam, 2023). This information is crucial for predicting floods and droughts, managing water supply, and displaying watershed characteristics (Brakenridge et al., 2012). However, obtaining reliable and readily available streamflow data is challenging. Currently, streamflow is measured at river gauging stations; however, access to this data is often intermittent or unavailable, particularly in developing and undeveloped nations like Bangladesh. Previously, many methods and techniques have been developed to predict discharge. Among them, traditional physical models and data-driven models are common. Models that predict discharge through physical methods rely on specific hydrological assumptions and necessitate extensive hydrological data for their calibration (Nearing et al., 2021). The hydrological data input (Rainfall, volume, intensity, dispersion) in these physically based models need to be accurate. However, acquiring such information can be challenging due to time constraints and higher costs. Especially in regions like Bangladesh, where there are limited monitoring systems, expensive data collection, or complicated logistics. In that case, the availability and reliability of hydrological data might limit the application of these models (Liang et al., 2019). This is why these models’ calibration and validation procedures are hindered by a lack of accurate and thorough hydrological data, leading to decreased forecasting accuracy (Dinic et al., 2021).

Data-driven models, such as ML models, require less input compared to physically based models and are efficient in data-scarce settings but do not represent physical processes explicitly. By avoiding physical characteristics and connections, ML models have the capability of finding accurate patterns and relations with few parameters (Wang et al., 2022). Due to this ability to precisely capture non-linear relationships among different variables, ML models are more flexible and reliable. ML models outperform traditional and other statistical models in flow prediction (Khosravi et al., 2022; Kumar et al., 2023; Rahimzad et al., 2021; Yaseen et al., 2019). The importance of ML techniques has proven to be powerful tools in this context, showcasing their ability to navigate complex relationships and find various applications in water resources (Van et al., 2020). However, there are debates on the optimum ML model. According to Wang et al. (2009), a comparison of different models, including Autoregressive Moving Average (ARMA), Artificial Neural Network (ANN), Adaptive Neural-based Fuzzy Inference System (ANFIS), Genetic Programming (GP), and Support Vector Machine (SVM) found that SVM, ANFIS, and GP performed the best based on training and validation results. In addition, previous studies concluded that some feature parameters unsuited for other algorithms can be used with SVM because of their distinctive properties (Hao et al., 2006). This enables the support vector machine approach to be extremely promising. However, (Nivesh et al., 2022) identified that, when contrasted with other models, such as the Random Forest (RF), Regression Tree (RT), Reduced Error Pruning (REP), Gaussian Process (GP), and SVM–based models, the performance of the RF models was discovered to be the most effective.

Although ML models have been applied to streamflow forecasting worldwide, their use in Bangladesh’s data-scarce, transboundary river systems remains limited, and existing applications leave several gaps. First, many studies estimate discharge from same-month meteorological inputs yet describe the result as “forecast” even though no forecast lead time is defined and such contemporaneous estimation is of limited value for early warning. Second, in data-scarce basins, it is seldom established which predictors carry forecast skill, whether empirically simulated runoff contributes information beyond rainfall and temperature, or whether antecedent hydrological conditions dominate the response. Third, evaluation is frequently based on random data splits that leak temporal information, on aggregate error metrics that obscure performance during floods and droughts, and without comparison to the simple baselines that any useful model must exceed. This study addresses these gaps by developing an operational, one-month-ahead discharge forecasting framework at the gauging-station scale, using only the variables available in the basin (rainfall, temperature, and discharge). For the Karatoya–Atrai River, a transboundary river in northwestern Bangladesh, discharge in each month is forecast from information available in preceding months, including lagged rainfall and temperature, accumulated rainfall, empirically simulated (Khosla) runoff, and antecedent discharge, all derived from the existing record without additional data collection. A range of ML models spanning linear, regularized, kernel-based and ensemble methods is compared against persistence and seasonal-climatology baselines under a strictly chronological validation scheme, and an ablation experiment quantifies the contribution of each predictor group, including the empirically simulated runoff whose role in the runoff–discharge relationship has not previously been examined in this setting. Model performance is further assessed separately for high- and low-flow months and for the monsoon and dry seasons, in keeping with the flood- and drought-management context that motivates the work. The contribution of this study is therefore threefold. Methodologically, it establishes a leakage-free, baseline-benchmarked and regime-aware forecasting protocol suited to the short, incomplete records typical of data-scarce basins. Scientifically, it clarifies which of the few available predictors govern discharge forecast skill in a monsoonal, transboundary catchment, providing new insight into the runoff–discharge relationship in a region where such ML applications have not previously been explored. Practically, it yields a parsimonious, low-parameter model that relies only on routinely available meteorological and discharge data. Because the analysis rests on a single gauging station, its transferability to other stations and basins remains a hypothesis for future work rather than a demonstrated result; if confirmed, such a parsimonious framework could offer a basis for discharge forecasting and for flood- and drought-related climate-resilience efforts in data-scarce communities.

2 Study area

The Karatoya-Atrai River (Figure 1), a major transboundary river flowing through India and Bangladesh (Anik et al., 2022; Chaki et al., 2014), holds historical significance as one of the main tributaries of the Teesta River. The river originates in the depressed area of Baikunthapur, India, before crossing into Bangladesh at Bardeshwari in the Panchagarh district (26°28′ latitude, 88°36′ longitude). After flowing through Bangladesh, it re-enters India through West Bengal and eventually returns to Bangladesh (Nowreen et al., 2020).

The Karatoya-Atrai River spans approximately 50 km within India and has a total length of around 340 km (Quader, 1995). It can be divided into two sections: the upper and lower Karatoya-Atrai.

The upper section stretches from Panchagarh to Dinajpur, while the lower section extends from Naogaon to the Karatoya-Atrai Railway Bridge (Quader, 1995). The deepest point of the river reaches 30 meters, the river basin spans almost 200 meters in width, and the yearly average precipitation totals 1,500 millimeters, with 80% of this rainfall occurring from June to mid-October (Ghosh et al., 2024; Saha and Pal, 2019). Flooding happens nearly every year in the watershed due to brief but intense rainfall events and the absence of any flood prevention strategies. High temperatures (exceeding 35 °C) and minimal rainfall from March to May lead to decreased water levels, causing many areas of the river to become entirely dry (Pal and Sarda, 2020).

Nine stations monitor water levels, and four stations measure the discharge along this river. Among these four discharge stations, Station SW144 (Chak Hariharpur) is particularly important, as it is located near the India-Bangladesh border and monitors water flow throughout the rest of the region.

3 Materials and methods

3.1 Materials

This study uses ML models to forecast the discharge using various secondary data sources from different organizations (Figure 2). The secondary data sources include hydrological and meteorological data, reports, studies, and other publications from government agencies, academic institutions, and non-governmental organizations (NGOs) to evaluate the model. Monthly temperature data were collected from 1964 to 1994 from the Bangladesh Meteorological Department (BMD). Monthly rainfall and discharge data were collected from the Bangladesh Water Development Board (BWDB) for the same period. These data sources provide valuable information about the region’s hydrology and ecology, as well as the potential impact on the floodplain of the Karatoya-Atrai River. Runoff was then calculated based on rainfall and temperature data.

3.2 Methods

The methodology proceeds through a sequence of linked phases, summarized in Figure 2. Hydro-meteorological data were collected and pre-processed, with an empirical runoff series derived from precipitation and temperature as an additional predictor. Lagged, antecedent, and calendar-based predictors were then constructed and used to train and tune multiple models, alongside two reference forecasts, persistence and seasonal climatology, which set the minimum skill any model must exceed. Models were evaluated using multiple metrics, with an ablation analysis quantifying each predictor group’s contribution and a regime-specific assessment covering high/low flow, seasons, and peaks. The best-performing configuration was selected as the final one-month-ahead forecast. All the analysis was done using python 3.

3.2.1 Forecasting framework

The task was formulated as one-month-ahead forecasting of monthly discharge. For the target month t, discharge Q t was predicted from information available up to the month − 1, together with the calendar month of t, which is known in advance:

Q t = f ( x ≤ t − 1 , m t ) (1)

Where x ≤ t − 1 denotes predictors derived from observations up to the month t − 1 and m tis the calendar month. Restricting the predictors in this way ensures that no contemporaneous meteorology enters the model, so that the reported skill reflects operational forecasting at a one-month horizon.

3.2.2 Data pre-processing and predictor construction

Runoff was calculated using precipitation and temperature data. Here, runoff refers to the volume of water generated at the surface due to precipitation, while discharge represents the flow rate measured at a gauging station. Khosla’s runoff calculation formula was employed for this purpose. Khosla (1949) examined the precipitation, overflow, and temperature information for different catchments in India and the USA to build an exact connection between spillover and precipitation (Archana et al., 2023). The time span must be at least a month. His relationship with month-to-month overflow for the T > 4.5 °C is as follows:

R = P − L (2)

Here, R (greater than 0) represents monthly runoff in mm; P is the monthly rainfall; and L is the monthly loss, equal to 0.48*Temperature. This formula, a well-established method in hydrology, enabled the estimation of runoff from these meteorological variables. The equation has been tested on various catchments in India (Dahe and Deshmukh, 2018; Dalavi et al., 2024; Ghose et al., 2018; Madhusoodhanan et al., 2012; Samantaray and Ghose, 2020) and is found to offer great outcomes for the yearly period to be utilized in preliminary examinations. However, this equation is a generic empirical relationship that was not calibrated or validated for the Karatoya-Atrai Basin, and does not represent infiltration, evapotranspiration partitioning, groundwater or routing.

Before modeling, the monthly records were checked for physically implausible values. At this station, temperature and discharge are always above zero (minimum observed values were 23.0 °C and 5.7 m³ s⁻¹). In contrast, zero rainfall and zero Khosla-derived runoff can naturally occur during the dry season, so these values were kept. In total, only three missing values for rainfall (April 1965 and 1972, May 1970, 0.20% of the hydro-meteorological data) were corrected using linear interpolation, which is suitable for evenly spaced monthly data. A final consistency check confirmed that there were no months with zero rainfall but non-zero runoff.

To forecast the discharge of a given month one month in advance, three kinds of information were used. The first is the time of year. Because the basin has a strong annual cycle, the calendar month was encoded using a sine and a cosine term. This treats the months as a repeating cycle, so that December and January are recognized as neighbors rather than as the two ends of a numerical scale. The second is the recent weather. This comprised the rainfall, temperature, and Khosla-derived runoff of the preceding month, the rainfall and temperature of the month before that, and the total rainfall accumulated over the three preceding months. Together, these describe how wet and how warm the basin has been in the period leading up to the month being forecast. The third is the river’s own recent behavior. The discharge observed in each of the three preceding months was included, since a catchment retains water in soil, groundwater, and channel storage, and flow in 1 month therefore carries over into the next.

The data were split chronologically after the lead-aware feature engineering (369 from 372 months; the first 3 months, 1964–01 to 1964–03, are consumed by the three-month discharge lags), the first 80% (1964–04 to 1988–10, 295 months) serving as the training set and the last 20% (1988–11 to 1994–12, 74 months), retained for testing. Because several predictors are lagged, a chronological partition, combined with blocked five-fold time-series cross-validation on the training period, was adopted to preserve the temporal ordering of the record and to prevent information from future months from entering either model fitting or model selection. Hyperparameters were tuned by cross-validated search over the time-series folds, using grid search for low-dimensional parameter spaces and randomized search otherwise, with the root-mean-square error as the selection criterion. Standardization and target transforms were fitted separately within each fold.

Within the validation framework, the contribution of each group of predictors (Table 1) was quantified through an ablation experiment in which the groups were added cumulatively, and a consistently tuned RF was re-evaluated at each step. The first set contained only the calendar terms together with the previous month’s rainfall and temperature. The second added Khosla-derived runoff from the preceding month. The third instead extended the first set with the lagged meteorology, namely the rainfall and temperature of 2 months earlier and the three-month accumulated rainfall. The fourth added the discharge-memory terms, that is, the discharge of the three preceding months. A fifth configuration appended Khosla-derived runoff to this full set.

Table 1. Predictors group.

The comparisons (B – A) and (E – D) measure the contribution of Khosla-derived runoff, while (C – A) and (D – C) measure the contributions of lagged meteorology and antecedent discharge, respectively.

Because the intended application concerns floods and droughts, performance was assessed not only in aggregate but separately across flow regimes. High-flow months were defined as those whose observed discharge exceeded the 75th percentile of the record, low-flow months as those falling below the 25th percentile, and extreme flows as those above the 90th percentile; performance was additionally evaluated for the monsoon season, from June to October, and for the dry season, from November to May.

3.2.3 Machine learning models

Multiple ML methods were compared, chosen to span the main families used in data-driven streamflow modeling. Ordinary linear regression served as a simple and transparent reference. Ridge regression is a linear model that penalizes large coefficients and so guards against overfitting when predictors are correlated (Hoerl and Kennard, 1970). Polynomial regression is a type of regression analysis that models the relationship between independent and dependent variables using an nth-degree polynomial (Sinha, 2013). Second-degree polynomial regression allows the relationship between predictors and discharge to curve rather than follow a straight line and was implemented by explicitly forming the squared and interaction terms of the predictors. Support vector regression is a supervised ML technique. It operates within a theoretical framework of linear functions in a multi-dimensional feature space (Bouzgou, 2006; Melki, 2018; Ojo et al., 2022). It aims to minimize the error between predicted and actual values while allowing for a tolerance margin. It has the capability to capture both linear and non-linear associations between independent and dependent variables, enhancing its effectiveness and efficiency (Chen et al., 2017; Oguntunde et al., 2018).

Finally, two tree-based ensemble methods built from decision trees were included: random forest (RF) and extreme gradient boosting (XGBoost). RF is a supervised ML algorithm. It is one of the most used algorithms due to its accuracy, simplicity, and flexibility (Abd Algani et al., 2022; González-Sopeña et al., 2021; Singh et al., 2016; Teles et al., 2021; Tyralis et al., 2019). This technique combines multiple decision trees to make predictions. Each tree is built using a subset of data and features, and the final prediction is obtained by aggregating the predictions of individual trees. XGBoost is an improved version of the gradient boosting algorithm that uses regularization and efficient tree-building techniques to improve prediction accuracy and reduce overfitting (Chen and Guestrin, 2016). It builds trees in sequence so that each corrects the errors of its predecessors. Two simple reference forecasts established the level of skill that any model must exceed. Persistence forecasting, also known as naive forecasting, where the latest observation is used as the prediction for the next step, is commonly used as a baseline model for assessing forecast skill (Ghimire and Krajewski, 2020; Ghimire et al., 2020; Mittermaier, 2008) It assumes that “tomorrow will look like today” (Ghimire and Krajewski, 2020). While Seasonal climatology (seasonal naive forecast) uses the historical average for the corresponding month (or season) as the prediction. It is a standard benchmark in hydrology and climate forecasting because many environmental variables exhibit a strong seasonal cycle (Seibert, 2001). Because the basin has a strong annual cycle, neither reference is trivial, and comparison shows how much a model adds beyond repeating the previous month’s flow or the long-term seasonal average.

Hyperparameters were selected by cross-validated search over the training period, using blocked time-series cross-validation with five expanding-window folds; the held-out test set was not involved in tuning. For models with small parameter spaces (Ridge, polynomial regression, SVR) an exhaustive grid search was used, and for the ensemble models (RF, XGBoost) a randomized search of 25 iterations was used. In all cases the selection criterion was the mean cross-validated root mean square error, and a fixed random seed (42) ensured reproducibility. Table 2 reports, for each model, the range searched and the value selected. The analysis was implemented in Python 3 using scikit-learn (linear, ridge, polynomial and support vector regression, RF, and the cross-validation and search procedures) and the XGBoost library; the seasonal statistical benchmark used statsmodels.

Table 2. Hyperparameters of the models used in the study.

Model selection used blocked time-series cross-validation with five expanding-window folds on the training period (1964–04 to 1988–10, 295 months). The training block grew successively the earliest 50, 99, 148, 197 and 246 months and each fold was validated on the immediately following 49-month block (1968–06 to 1972–06; 1972–07 to 1976–07; 1976–08 to 1980–08; 1980–09 to 1984–09; and 1984–10 to 1988–10, respectively). Every validation block is therefore strictly later than the data used to train it, and the held-out test set (1988–11 to 1994–12, 74 months) was excluded from all tuning.

3.2.4 Evaluation and analysis tools-metrics

Eight different evaluation tools (Equations 3–13) are used to evaluate the performance of ML models. These include Mean Absolute Error (MAE), Mean Absolute Percentage Error (MAPE), Root Mean Square Error (RMSE), Normalized RMSE (NRMSE), Nash-Sutcliffe Efficiency (NSE), Kling–Gupta Efficiency (KGE), and Percent Bias (PBIAS). These metrics are used to evaluate both the model training and testing phases.

Mean Absolute Error (MAE): This metric quantifies the average magnitude of errors between the predicted and actual discharge values, providing a clear understanding of the model’s overall precision. In simple terms, it provides the arithmetic mean of the absolute error between the predicted and observed values.

MAE = 1 N ∑ i = 1 N ∣ X − Y ∣ (3)

Here, X represents the actual or observed value, Y represents the predicted value, and N stands for the number of observations.

Mean Absolute Percentage Error (MAPE): The MAPE measures the average percentage discrepancy between predicted and observed discharge, providing insight into the magnitude of relative error. It lies within the same category of measures of MAE and MSE (Mean Squared Error) (De Myttenaere et al., 2016). However, it focuses on the percentage of error (De Myttenaere et al., 2016). Many scholars have considered it unsuitable for evaluation due to large errors.

MAPE = 1 N ∑ i = 1 N ∣ X − Y X ∣ × 100 (4)

Mean Absolute Scaled Error (MASE): The MASE measures the average absolute forecasting error scaled by the average absolute error of a naive forecasting model. Unlike MAPE, MASE is scale-independent and remains well-defined when observed values are zero or close to zero, making it suitable for comparing forecasting performance across different datasets and time series (Hyndman and Koehler, 2006). A MASE value less than 1 indicates that the forecasting model performs better than the naïve benchmark, whereas a value greater than 1 indicates poorer performance.

MASE = 1 N ∑ i = 1 N ∣ X − Y ∣ 1 N − 1 ∑ i = 2 N ∣ X − Y t − 1 ∣ (5)

Here, 1 N − 1 ∑ i = 2 N ∣ X − Y t − 1 ∣ represents the mean absolute error of the naive one-step-ahead forecast.

Root Mean Square Error (RMSE): RMSE calculates the square root of the average squared differences between predicted and actual discharge, providing a measure of the model’s ability to capture variations. It is used to quantify the error in prediction in terms of the units of the variables. The values extend from 0 to infinity, with an RMSE of 0 signifying an ideal fit. The unit of RMSE is the same as the variable being measured.

RMSE = 1 N ∑ i = 1 N ( X − Y ) 2 (6)

Normalized Root Mean Squared Error (NRMSE): The NRMSE normalizes the RMSE to the scale of the observed discharge, facilitating the comparison of models across different scales and units. In simple words, RMSE describes the mean error, while NRMSE determines the significance of that error relative to the observed data. Two normalizations are reported: NRMSE_mean, which scales the error by the mean observed discharge, and NRMSE_range, which scales it by the observed range.

NRMSE mean = 1 X ˉ 1 N ∑ i = 1 N ( X − Y ) 2 (7)

NRMSE range = 1 X max − X min 1 N ∑ i = 1 N ( X − Y ) 2 (8)

Where X̄ is the average observed value. Both are dimensionless, with lower values indicating closer agreement; for NRMSE_range, values below approximately 0.1 denote excellent and 0.1–0.2 strong agreement. Occasionally, writers must standardize the MAE and RMSE to quantitatively assess the prediction accuracy of certain models. Their standardization is represented by the normalized root mean squared error (Piotrowski et al., 2022).

Nash-Sutcliffe Efficiency (NSE): The NSE is a standardized measure that assesses the ratio of the residual variance to the variance of the observed data (Nash and Sutcliffe, 1970). NSE assesses the model’s ability to capture discharge variability by comparing predicted values to observed values, reflecting the model’s goodness of fitness.

NSE = 1 − [ ∑ i = 1 N ( X − Y ) 2 ∑ i = 1 N ( X i − X ) 2 ] (9)

NSE is mainly used for the evaluation of hydrological models, specifically the predictive ability of hydrological models (Guinot et al., 2011; Ritter and Muñoz-Carpena, 2013). The NSE value ranges between negative infinity and 1, where 1 indicates the highest level of match of the predicted and observed value.

Kling–Gupta Efficiency (KGE): KGE (Gupta et al., 2009) is widely used to evaluate the performance of hydrological models. KGE is a performance metric that decomposes the error into three components: correlation, variability, and bias.

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

Here, 𝑟 is the linear correlation between observations and simulations, 𝛼 is a measure of the flow variability error, and 𝛽 is a bias term. The terms αand βare defined as follows:

β = μ s μ o (11)

α = σ s σ o (12)

Where, σ sis the standard deviation of the simulation, σ o is the standard deviation of the observation, μ sis the mean of the simulated time series (e.g., flows predicted by the model) and μ ois the mean of the observed time series.

Percent Bias (PBIAS): The PBIAS measures the average tendency of the predicted discharge to be larger or smaller than the observed discharge, expressed as a percentage of the total observed flow. Unlike MAE and RMSE, which quantify the magnitude of the error irrespective of its direction, PBIAS retains the sign of the error and therefore reveals systematic over- or under-prediction that error magnitudes alone would conceal.

PBIAS = [ ∑ i = 1 N ( Y − X ) ∑ i = 1 N X ] × 100 (13)

The optimal value is 0%, with values of larger magnitude indicating greater systematic bias. A negative PBIAS indicates that the model underestimates the observed discharge, whereas a positive PBIAS indicates over-estimation.

4 Results

4.1 Descriptive statistics and trend analysis of hydrometeorological parameters

Descriptive statistics of collected parameters, including rainfall, temperature, runoff, and discharge, are analyzed, as shown in Figure 3. Rainfall values range from 0 mm to 1379.33 mm, with a median of 132.83 mm. For temperature, the median value is 31.93 degrees Celsius (°C), with a range of 23 to 41 °C. Discharge exhibits a median value of 70.98 cubic meters per second (m³s⁻¹), with values ranging from 5.7 to 734 m³s⁻¹; most values fall below 150 m³s⁻¹. Calculated runoff ranges from 0 to 1,218 mm, with a median value of 0, as negative runoff values are adjusted to 0(Figure 3).

The descriptive statistics (Table 3) for the dataset reveal notable characteristics for each variable. The ranges for each variable indicate significant variation, especially in Rainfall (1379.33) and Runoff (1218.1), highlighting the large disparities in the data. The skewness values indicate that Rainfall and Runoff exhibit positive skewness, suggesting a rightward distribution, whereas Temperature and Discharge display more balanced or slightly skewed distributions, suggesting different underlying patterns in these variables.

Table 3. Descriptive statistics of the parameters.

A heatmap analysis assesses the correlations among the parameters used for predicting discharge (Figure 4). The values in the figure represent the strength and direction of the relationships among the variables, ranging from +1 (strong positive correlation) to −1 (strong negative correlation). Rainfall and runoff exhibit a high correlation value of 0.97, indicating that rainfall strongly influences runoff. These two parameters have also played a significant role in predicting discharge, with positive correlation values of 0.49 and 0.45, respectively.

In contrast, temperature shows the weakest correlation with discharge (coefficient of +0.20) and is most strongly correlated with rainfall (coefficient of +0.42). The absence of negative correlations in the dataset suggests that all variables contribute positively to discharge prediction.

To assess whether the hydro-meteorological regime is stationary over the study period, the annual series of rainfall, temperature and discharge were tested for monotonic trends and abrupt shifts (Table 4). The analysis was performed on annual values of total rainfall, mean temperature, and discharge so that the strong seasonal cycle does not confound the trend signal. The Mann–Kendall test evaluates the presence and direction of a monotonic trend through the statistic Z, whose sign indicates the direction and whose magnitude, compared against the standard normal distribution, determines significance; the Sen’s slope estimator gives the median rate of change per year; and the Pettitt test locates the most probable point of an abrupt shift in the mean, testing whether it is statistically significant.

Table 4. Trend and change-point analysis of the annual hydro-meteorological series (1964–1994).

The results show a divergence between the meteorological drivers and the discharge response. Neither annual rainfall (Z = −0.75, p = 0.455) nor annual mean temperature (Z = +1.50, p = 0.135) exhibits a statistically significant trend over 1964–1994. The slight negative slope in rainfall and positive slope in temperature is within the range expected from interannual variability. Annual discharge, by contrast, shows a statistically significant increasing trend (Z = +2.46, p = 0.014), rising at approximately 2 m³ s⁻¹ per year a cumulative increase of the order of 60 m³ s⁻¹ across the record. As the evaluation uses a chronological split, the increasing trend places the highest flows in the later test period, partly beyond the range represented in training which, given the inability of tree ensembles to extrapolate, contributes to the underestimation of high and extreme flows. The presence of a discharge trend without corresponding trends in local rainfall or temperature is notable and is consistent with the transboundary setting of the Karatoya–Atrai, in which streamflow is influenced by upstream conditions in India that are not captured by the local meteorological record. The Pettitt test detects no statistically significant change point in any of the three series, indicating that the discharge increase is a gradual trend rather than an abrupt regime shift.

4.2 Model comparisons

This study evaluates ML models for forecasting discharge and identifies the best-performing model, which is Random Forest (RF). The RF model achieves the highest accuracy and closely follows the ideal line observed versus forecasted discharge values, as shown in Figure 5.

However, Figure 5 provides only a visual representation of the models’ performance and does not quantify their accuracy. To validate the forecasted performance of the models, forecasts were made for the test data. Evaluation metrics, including NSE, KGE, MAE, MAPE, RMSE, and NRMSE, are calculated based on the forecast of each model and are compared to one another to identify the optimal model. These values for RF, SVR, LR, PR, RR and XGBoost are shown in Figure 6.

According to Figure 6 and Table 5, the performance of the evaluated models varies considerably, indicating that some approaches are more effective than others in forecasting discharge. Overall, the RF model demonstrates the strongest predictive capability, achieving the lowest error values and the highest overall efficiency among the tested models. RF obtains the lowest MAE (43.97 m³s⁻¹), RMSE (69.72 m³s⁻¹), NRMSE_mean (0.379), and NRMSE_range (0.113), along with the highest NSE value (0.857), indicating its superior accuracy and ability to capture discharge variability. The RF model is followed closely by XGBoost, which achieves comparable performance with MAE of 46.39 m³s⁻¹, RMSE of 71.26 m³s⁻¹, NRMSE_mean of 0.388, NRMSE_range of 0.116, and NSE of 0.851. Although XGBoost provides the highest KGE value (0.902) and a lower PBIAS (−6.19%) compared with RF, RF offers the best overall balance across multiple evaluation criteria.

Table 5. Evaluation of metric values for different models.

The regression-based models, including Linear Regression, Polynomial Regression, SVR, and Ridge Regression, show moderate performance but generally exhibit higher errors and lower efficiency compared with the ensemble learning approaches. Polynomial Regression achieves relatively competitive results with MAE of 49.10 m³s⁻¹ and MAPE of 40.88%, whereas SVR provides the lowest MAPE (39.74%) among all models but suffers from higher RMSE (79.36 m³s⁻¹), lower NSE (0.815), and lower KGE (0.825). The Persistence model shows the weakest forecasting capability, with the highest MAE (88.16 m³s⁻¹), RMSE (128.55 m³s⁻¹), MAPE (70.79%), and lowest NSE (0.514), demonstrating limited ability to reproduce discharge variations. Overall, the results highlight the superiority of machine-learning ensemble models, particularly RF and XGBoost, for accurate discharge prediction compared with conventional statistical and single regression-based approaches.

To visually represent the correlation between actual and predicted values by the top three (based on RMSE) models in this study, a scatter plot of prediction errors is presented in Figure 7. Four different error categories are used to determine which model has the least prediction errors. Most of the error margins in RF lie within the 0 to 20% range, while LR and XGBoost have error margins in the range of 20 to 35% and above 35%. Additionally, most of the values in the RF model align with the linear line of perfect prediction in RF mode, indicating the model’s accuracy.

4.3 Performance across flow regimes

Disaggregating the RF forecasts by flow regime (Table 6) shows that the model’s overall skill is concentrated in the mid-range and seasonal conditions. Monsoon and dry-season forecasts exhibit good performance (NSE = 0.87 and 0.85) with MASE values below 1 (0.56 and 0.69), indicating that the model outperforms a naïve benchmark in both seasons. Performance deteriorates for the flows most relevant to hazard management. High-flow months (above the 75th percentile) and extreme-flow months (above the 90th percentile) exhibit large absolute errors (RMSE = 123 and 149 m³ s⁻¹), systematic underestimation (PBIAS = −11 and −23%), and MASE values greater than 1 (1.49 and 1.85), indicating that the model performs worse than a naïve forecasting benchmark under these conditions. Although the negative NSE values within these narrow flow ranges are partly attributable to the limited variance within each subset, the consistently high RMSE, negative PBIAS, and MASE > 1 collectively confirm reduced reliability in predicting flood peaks. In contrast, low-flow months have small absolute errors (RMSE = 16 m³ s⁻¹; MAE = 11 m³ s⁻¹) and a low MASE (0.16), demonstrating strong performance relative to the naïve benchmark despite the very high MAPE (116%). The inflated MAPE reflects the well-known instability of percentage-based error metrics when observed discharge approaches zero rather than poor predictive accuracy. Overall, these results indicate that the primary limitation of the model lies in forecasting high and extreme flows, while low-flow predictions remain accurate in absolute terms.

Table 6. Regime-specific test performance of the RF.

4.4 Contribution of predictor groups

To quantify the independent contribution of each predictor group and to determine whether the Khosla-derived runoff carries information beyond rainfall and temperature an ablation experiment was conducted in which predictor groups were added cumulatively (sets A–E; Table 1). Two points regarding the use of data should be emphasized. First, predictor-set selection was performed entirely on the training data, using the blocked five-fold time-series cross-validation described in Section 3.2.2; the held-out test set was not used to select predictors, tune hyperparameters, or decide the final model. Second, for completeness and transparency, the performance of each predictor set is also reported on the test set, so that the training-based selection and the out-of-sample behavior can be compared directly. Table 7 presents both: the cross-validated scores on the training folds, which form the basis of selection, and the corresponding single-evaluation scores on the held-out test set.

Table 7. Ablation of predictor groups (random forest).

Cross-validated scores are the mean ± standard deviation over the five training-period folds and form the basis for predictor selection; test-set scores are reported for confirmation and were not used in selection.

The cross-validation results identify antecedent discharge as the decisive predictor group: adding discharge memory raises the cross-validated NSE from 0.745 (set C) to 0.779 (set D), an improvement of +0.034, by far the largest of any group. In contrast, the Khosla-derived runoff term provides no benefit on the selection metric: adding it changes the cross-validated NSE by −0.005 (B versus A) and by −0.004 (E versus D) that is, it slightly reduces skill in both comparisons. This is consistent with the near collinearity of runoff and rainfall (r = 0.97; Figure 4): because the runoff term is an algebraic transform of rainfall, it introduces no independent information and marginally degrades the model through added redundancy. The lagged-meteorology group likewise contributes little on the cross-validation metric (C versus A, −0.001), indicating that most of the predictive skill beyond seasonality is carried by antecedent discharge rather than by antecedent rainfall and temperature.

On this basis, set D seasonality, antecedent meteorology and antecedent-discharge memory, without the runoff term was selected as the final predictor set on the grounds of parsimony and non-redundancy. The test-set scores confirm that this choice generalizes: set D attains a test NSE of 0.857, and although set E (D + runoff) shows a marginally higher test NSE (0.866), this apparent advantage is small (+0.009), lies within the seed-to-seed variability of the RF, and is contradicted by the cross-validation result and by the B-versus-A comparison. Retaining a redundant predictor for a marginal, non-robust gain on a single test partition would not be justified and would be precisely the kind of test-set-influenced decision that a rigorous evaluation should avoid.

4.5 Model performance in training and testing phase

Figure 8 compares the training and testing RMSE of each model on the chronological split, with Δ denoting the difference between them. The comparison shows that training fit is a poor guide to forecast skill, and that the models differ markedly in how well they generalize. XGBoost reproduces the training record almost exactly (RMSE ≈ 0.2 m³ s⁻¹) yet attains a test RMSE of 71.26 m³ s⁻¹, a gap of 71.1 m³ s⁻¹ and a clear signature of overfitting on a record of this length. The RF also fits the training data closely (≈ 30.6 m³ s⁻¹) but generalizes better, achieving the lowest test RMSE of any model (69.72 m³ s⁻¹) despite a residual gap of 39.1 m³ s⁻¹; a gap of this kind is expected for bagged tree ensembles, which by construction fit the training sample closely. Support vector regression shows a smaller gap (11.3 m³ s⁻¹), and the polynomial model essentially none (0.2 m³ s⁻¹).

Notably, the linear and ridge models return negative differences (−7.3 and −10.6 m³ s⁻¹), performing slightly better on the unseen test period than on the training period. This indicates an absence of overfitting in these models and confirms that the chronological hold-out does not constitute an unrepresentatively difficult test period. Taken together, the results show that the ranking of models by training error (in which XGBoost would appear best) is unrelated to their ranking by test error, and that the RF was selected based on out-of-sample forecast skill rather than training fit.

4.6 Optimal ML model

The Random Forest (RF) model demonstrates the best overall predictive performance (Figure 9) compared with climatology, persistence, XGBoost, LR, PR, SVR, and RR models. RF achieves the lowest MAE (43.971 m³s⁻¹), RMSE (69.715 m³s⁻¹), NRMSE_mean (0.379), and NRMSE_range (0.113), indicating superior accuracy and reduced prediction errors. Furthermore, RF obtains the highest Nash–Sutcliffe efficiency (NSE = 0.857), followed by XGBoost (0.851), highlighting its strong capability to reproduce observed variations. RF also achieves a high Kling–Gupta efficiency (KGE = 0.887), although XGBoost provides the highest KGE value (0.902) among all models. Overall, the results indicate that RF provides the most balanced and reliable performance across multiple evaluation criteria, outperforming traditional regression approaches such as PR, SVR, LR, and RR.

The hydrograph comparison (Figure 10) indicates that the RF model generally reproduces the temporal and seasonal variability of observed monthly discharge during the test period at a one-month lead time. The model successfully captures the timing and overall shape of the major discharge peaks and recession periods. However, the comparison of the largest high-flow events shows a tendency to underestimate peak discharge, particularly for the most extreme events, while performance is comparatively closer during moderate-flow conditions. Thus, although the model demonstrates good skill in reproducing overall discharge dynamics, its reduced accuracy in representing extreme high flows should be considered when interpreting its suitability for flood-related forecasting.

Finally, to interpret the RF forecasts, a SHAP (SHapley Additive exPlanations) analysis was applied to the test set (Figure 11). The analysis shows that the forecasts are governed chiefly by seasonality and by the river’s recent flow. The cyclic month terms month_sin and month_cos exert the largest and most widely spread influence, reflecting the strong and systematic annual cycle of the basin, while the previous month’s discharge, Q_L, is the most influential hydrological predictor and displays a clear monotonic effect. High Q_L (red) shifts the forecast upward, and low Q_L (blue) shifts it downward, consistent with the catchment-memory mechanism by which storage in soils, groundwater and channels carries flow from 1 month into the next. The rainfall predictors Rain_L, Rain_acc3 and Rain_prev contribute moderately and in the expected direction, whereas the temperature terms Temp_L and Temp_prev and the more distant discharge lags Q_L1 and Q_L2 have only minor effects.

These attributions corroborate the ablation experiment (Section 4.4): the predictive skill of the model derives from the combination of seasonality (month_sin, month_cos), antecedent discharge (Q_L) and recent rainfall (Rain_L, Rain_acc3, Rain_prev), rather than from Khosla-derived runoff, which was excluded from the final predictor set for its negligible independent contribution and accordingly does not appear among the influential variables. It should be noted that SHAP values quantify each predictor’s contribution within the fitted model rather than a causal effect, so that influence may be shared among correlated predictors, such as month_sin and month_cos or the successive discharge lags Q_L, Q_L1 and Q_L2; the analysis therefore explains how the model uses the available information and reinforces, by an independent route, the conclusion that catchment memory and seasonality are the dominant controls on one-month-ahead discharge at this station.

5 Discussions

Among the models compared, RF produced the most accurate one-month-ahead forecasts, with the lowest RMSE, MAE, and NRMSE, and the highest NSE. Its NSE of 0.86 indicates that the model reproduces the great majority of the observed variance in monthly discharge (Rezaie-Balf et al., 2019; Schaefli and Gupta, 2007), while a KGE of 0.89 shows that this skill is achieved with a well-matched mean and variability rather than through error compensation (Gupta et al., 2009). A range-normalized RMSE of 0.11 places the forecasts in the range conventionally regarded as strong, and a PBIAS of − 8.4% indicates only a modest overall tendency to under-predict (Moriasi et al., 2007)

These values are, however, best judged against a reference rather than against fixed thresholds, since a performance score is interpretable only relative to what a trivial model already achieves (Knoben et al., 2019; Seibert et al., 2018). Seasonal climatology, which simply predicts the long-term mean discharge of the target calendar month, itself attains an NSE of 0.80, reflecting the strength and regularity of the monsoonal annual cycle. Persistence, which carries the previous month’s flow forward, achieves NSE of 0.51. The RF therefore improves on the climatological benchmark by a real but moderate margin (ΔNSE ≈ 0.06; ΔRMSE ≈ 13.3 m³ s⁻¹). This tempers, rather than negates, the practical value of the model: a substantial part of the predictability in this basin is seasonal and is available without any machine learning, and the contribution of the ML models lies in the additional skill they extract from antecedent conditions.

The ablation experiment clarifies where that additional skill originates. Antecedent discharge was the single most valuable group of predictors, and lagged and accumulated rainfall contributed further, whereas Khosla-derived runoff produced no consistent improvement once rainfall and temperature were available: adding it changed NSE by −0.005 when appended to the rainfall–temperature set and by +0.009 when appended to the full set. Because this runoff is computed from rainfall and temperature and correlates with rainfall at r = 0.97, it constitutes a transformation of information already supplied to the model rather than an independent hydrological observation, and set D was retained on grounds of parsimony. This finding is consistent with the SHAP analysis, which attributes the forecasts principally to the seasonal terms and to the previous month’s discharge (Lundberg and Lee, 2017), and it indicates that in this basin catchment memory and seasonality, rather than an empirical rainfall-loss relationship, govern monthly discharge at a one-month horizon.

Disaggregating performance by flow regime reveals a limitation of direct relevance to the intended application. Although monsoon and dry-season forecasts are both accurate (NSE 0.87 and 0.85), skill deteriorates markedly for the flows that matter most for hazard management: high flows above the 75th percentile and extreme flows above the 90th percentile carry RMSE values of 123 and 149 m³ s⁻¹ and a systematic negative bias (PBIAS −11 to −23%). Part of this behavior is intrinsic to tree ensembles, whose predictions are bounded by the range of the training data, and which therefore cannot extrapolate to flows exceeding those previously observed (Cai et al., 2023; Gnjato et al., 2026; Tyralis et al., 2019). Low-flow months show the converse pattern, with small absolute but large relative errors, a well-documented consequence of percentage-based metrics near zero (Hyndman and Koehler, 2006). Consequently, while the framework is informative for anticipating seasonal water availability, its reliability for flood-peak and drought-minimum forecasting is limited, and this qualification should temper its operational use.

The comparison among models is instructive in a way that aggregate rankings alone would obscure. Previous studies suggest that compared to traditional statistical models such as polynomial regression (PR) and autoregressive integrated moving average (ARIMA), RF demonstrates superior performance due to its ability to capture complex, nonlinear interactions among hydrological variables (Gacu et al., 2025; Papacharalampous et al., 2019; Tyralis et al., 2019). Additionally, when compared to ML models like support vector regression (SVR) and artificial neural networks (ANNs), RF provides better generalization and avoids overfitting, making it a preferred choice for discharge modeling (Shortridge et al., 2016; Yaseen et al., 2016). RF’s strong predictive ability makes it a valuable tool for water resource management, flood forecasting, and climate impact studies (Ahmed et al., 2024; Asif et al., 2025; Islam et al., 2025; Khan et al., 2026). By accurately modeling discharge levels, RF can help policymakers and hydrologists make informed decisions about reservoir operations, flood control measures, and drought mitigation strategies.

Though the study provides a robust and comprehensive understanding of ML models by finding the optimal model, certain limitations were encountered. The analysis rests on a single gauging station and a monthly time step; daily data would both enlarge the training sample and resolve the flood dynamics that monthly aggregation suppresses, and multi-station application would test the transferability claimed above. Despite its strengths, RF also has certain limitations. One drawback is its computational complexity, particularly when handling large datasets, as training multiple decision trees can be a time-consuming process. Additionally, while RF performs well under normal hydrological conditions, extreme events such as floods and droughts may require further enhancements, such as hybrid models or deep learning approaches, to improve forecasting accuracy (Gnjato et al., 2026; Mosavi et al., 2018; Pang et al., 2025; Vinokić et al., 2025). Future research could explore hybrid models that integrate RF with deep learning techniques, such as long short-term memory (LSTM) networks, to further improve RF’s forecasting performance. Finally, longer lead times than one month, and applying RF to real-time hydrological forecasting systems and testing its adaptability to different climate scenarios could provide deeper insights into its long-term reliability.

6 Conclusion

This study evaluated six machine-learning models for one-month-ahead forecasting of monthly discharge in the Karatoya–Atrai River, using only rainfall, temperature, and discharge records under chronological validation and benchmarked against persistence and seasonal-climatology references. The RF produced the most accurate forecasts (NSE 0.86, KGE 0.89, RMSE 69.71 m³ s⁻¹, MAE 43.97 m³ s⁻¹, MAPE 40.26%, NRMSE 0.11), followed closely by XGBoost. All models exceeded persistence (NSE 0.51) and seasonal climatology (NSE 0.80), though the moderate margin over seasonal climatology indicates that much of the monthly predictability arises from the basin’s pronounced annual cycle. Feature importance analysis showed that the previous month’s discharge and seasonality were the most important predictors, while Khosla-derived runoff added little improvement beyond rainfall and temperature, consistent with its near collinearity with rainfall (r = 0.97). A configuration using rainfall, temperature, and calendar information alone nonetheless retained most of the skill, which is relevant where real-time discharge is unavailable. Performance was not uniform across regimes; monsoon and dry-season forecasts were accurate (NSE 0.87 and 0.85), but high and extreme flows were systematically under-predicted, so the framework is better suited to anticipating seasonal water availability than to flood-peak or drought-minimum forecasting. Within these bounds, it offers a parsimonious, transferable basis for discharge forecasting in data-scarce basins. Further work should extend it to multiple stations, adopt a daily time step, examine architectures with explicit memory such as LSTM networks, and issue probabilistic rather than deterministic forecasts.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. The authors used Generative AI tools for language refinement and textual editing of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Mohammad Kamruzzaman, Bangladesh Rice Research Institute, Bangladesh

Reviewed by: Dilip Kumar Roy, Bangladesh Agricultural Research Institute, Bangladesh

K. Prashanth Kumar, NMAM Institute of Technology, India

Nayab Zahara, King Faisal University, Saudi Arabia