Comparison of Bayesian and frequentist quantile regressions in studying the trend of discharge changes in several hydrometric stations of the Gorganroud basin in Iran

River discharge is a crucial factor in hydrology and water resources, and it is affected by various climatic factors like temperature and precipitation. Any change in these factors can lead to alterations in the discharge trend (Koutalakis et al. 2019; Wei et al. 2020). Reducing runoff poses a significant challenge in managing water resources, especially in arid and semi-arid regions such as Iran (Kanani et al. 2019). Lower rainfall can decrease runoff and subsequently reduce water resources (Sabbaghi et al. 2020). Conversely, heavy rainfall can result in a sudden increase in discharge, leading to economic and social damage in sectors like agriculture, livelihood, and infrastructure. This can seriously disrupt regional development plans (Baba'eeyan et al. 2013). Therefore, studying the changes in river discharge and runoff volume due to rainfall can help predict flood time and reduce damages caused by it. Recently, many basins in Iran have undergone significant changes, with a decrease in the amount of runoff and river flow, or in some cases, an increase in the form of floods (Sharafati et al. 2020; H. Ahmadi et al. 2020a; M. Ahmadi et al. 2020; Kousali et al. 2022). Trend analysis is a suitable model to evaluate hydro-climatological conditions in basins, as it helps study changes in a variable over time (Binesh et al. 2017). Most studies on trend detection use the Mann–Kendall non-parametric test, which was developed by Kendall (1975) and Mann (1945). This method has been used to investigate trends in temperature by Ahmadi et al. (2018) and Alhaji et al. (2018), precipitation by Basarir et al. (2018), Nyikadzino et al. (2020) and Salehi et al. (2020), flow by Ghorbani et al. (2016), Tehrani et al. (2019), Galkate et al. (2020), Norouzi (2020), Kuriqi et al. (2020), and Owolabi et al. (2021), and droughts by Kalisa et al. (2020) and Ghorbani et al. (2020). However, these studies focus on the average and median of the data, which might not capture trends that occur within a specific range of the data. Additionally, the mean alone is not sufficient to fully understand the distribution pattern of the variable being studied at different levels of the independent variable, as noted by Bararkhanpour et al. (2020). The pattern of discharge is not consistent and can be affected by changes in weather conditions, resulting in abnormalities and asymmetry. Traditional parametric models may not accurately capture these changes over time. Meanwhile, non-parametric methods can be problematic when dealing with small data samples with short lengths. To strike a balance between efficiency and accuracy, Ker & Coble (2003) proposed a semi-parametric model based on quantile regression. This model can plot regression curves that match different percentile points without the constraints of normal linear regression. Many studies have since applied this model to analyze changes in climatic and hydrological parameters, such as rainfall intensity and trends (Fan & Chen 2016; Dunn et al. 2019; Gao & Franzke 2017; Salarijazi 2017; Salman et al. 2018). For example, Dhakal & Tharu (2018) observed changes in heavy rainfall distribution and their relationship to tropical storms in the southeastern United States. Abbas et al. (2019) found that the quantile regression model can provide patterns for very dry and very wet conditions. Bararkhanpour et al. (2020) studied seasonal changes in precipitation in Gorgan, Iran, and found that the lower quantiles in autumn and the middle quantiles in spring and winter have a decreasing trend slope. During summer and winter, there is an upward trend in the upper quantiles. Mohsenipour et al. (2020) conducted a study on the distribution of precipitation during the monsoon months in Bangladesh (June to September). The findings revealed that rainfall decreased in many quantiles from June to August and increased in September at most stations. Zhang et al. (2020) found that the Arctic Sea ice's quantiles from 0.01 to 0.99 had a statistically significant decrease over time. The study also confirmed the important role that climate patterns play in reducing sea ice, especially in the upper and lower quantiles of sea ice concentrations. Treppiedi et al. (2021) researched the trend of changes in precipitation in the Mediterranean region and found an increasing trend in less than an hour of rainfall, particularly in the upper quantiles of precipitation. Despite the quantile regression model's advantages, the quantile curves in this model are often estimated independently. As a result, the quantile curves are evaluated independently, and estimates at extreme quantiles do not consider information from estimates in near-distributed masses (Ramsey 2020). Due to their independent estimation, quantile curves can intersect, which invalidates the distribution of independent variables (Ramsey 2020).

Quantile regression Bayesian models are a solution to the limitations of frequentist models and provide accurate estimates by utilizing information from quantile levels and locations. This method surpasses the shortcomings of the frequentist model and has been proven to be more precise and accurate (Hu et al. 2021; Alahmadi et al. 2022). Reich et al. (2011) proposed Bayesian quantile regression, which is an effective model for measuring variable distribution over time. This model offers a smooth function for the relationship between all possible quantiles, unlike the frequentist model that provides point estimates. Additionally, the Bayesian model borrows information from all quantile levels, leading to more precise estimates. Despite limited research on the subject, some studies have utilized the Bayesian quantile regression model to analyze factors affecting various parameters and their trend over time. Tan et al. (2019) modeled the distribution of changes in winter rainfall in Canada and discovered an increase in the upper and lower quantiles in eastern Canada. Tharu & Dhakal (2020) used the Bayesian regression model to analyze trends in monthly and daily maximum precipitation indices across 1,108 locations in the United States over 65 years. The study revealed that the upper quantiles of heavy rainfall distribution changed at a much higher rate than the lower quantiles. Amini et al. (2020) investigated the relationship between oceanic–atmospheric indices and standardized precipitation index (SPI) in Iran and found that La Niña intensified drought conditions in several regions of the Caspian Sea, while the North Atlantic Oscillation (NAO) had the weakest effect on drought in different quantiles and regions of Iran. Yong et al. (2021) also conducted research on the dependence of heavy precipitation on dew point temperature and surface air temperature at 78 stations on the Tibetan Plateau. The study concluded that there is a stronger correlation between heavy precipitation and dew point temperature. This study aims to examine the trend in flow values, particularly low and high flows, in the Gorganroud River basin. The investigation will be conducted in three stations, namely Arazkouseh, Tamar, and Galikesh, located in Golestan province, Iran, during different seasons. Two statistical methods, frequentist quantile regression and Bayesian quantile regression, will be used to analyze the data. The study will also compare the results obtained from both methods to determine the significance of using the Bayesian approach in quantile regression.

The Gorganroud basin, situated in Golestan province in Iran, is a flood-prone area that causes significant damage to its economic and social sectors every year. Despite being located in northern Iran's green and forested belt, the region is susceptible to drought. Therefore, examining the discharge in its rivers at different time scales is crucial. Arazkouseh, Tamar, and Galikesh stations, the oldest hydrometric stations with daily discharge statistics spanning half a century, were selected as a case study in this basin.

The Gorganroud basin covers an area of approximately 11,380 km², which is equivalent to 48% of Golestan province. Its main branch, running through the Gorgan and Gonbad plains, spans 359.4 km and contains 67% of the surface water in the province, amounting to 828 million cubic metres. The region's average annual rainfall ranges from 200 to 700 mm. Forests dominate the basin's southern and eastern regions, while alluvial plains are used for agriculture and pasture in the north and west.

The study utilized a daily time series of discharge data provided by the Regional Water Authority of Golestan Province. The selection of stations was based on their sufficient statistics for over 30 years and not having a dam to collect and store water from upstream. Therefore, Arazkouseh, Tamar, and Galikesh stations were chosen as they had 50-year daily discharge data from 1966 to 2016 (1346–1395). The study calculated the average monthly discharges, formed their seasonal series, and conducted trend tests for analysis.

This study utilized two models, namely frequency-oriented quantile regression and Bayesian multiple regression, to analyze the seasonal trend of changes in various monthly discharge values.

Koenker (2006) developed Quantreg software in the R programming language for quantile regression problems. The rank inverse model is used to estimate the intercept and slope coefficient of the quantile regression model. The standard errors, confidence intervals, t-statistic, and significance for these coefficients are also estimated using the Quantreg software package.

Ramsey (2020) explains that the covariate effect varies in both quantile level and location. To ensure a valid quantile process, the quantiles must be non-decreasing for all predictor values. Estimating the quantile using this model has the advantage of sharing information or data across quantile levels. However, formulating Bayesian quantile regression is more complex than classical quantile regression because the effects of parameters must be theoretically determined over an infinite number of quantile levels.

This is Bernstein's basis polynomial equation, and Equation (4) is specifically a Bernstein polynomial of degree M. Bernstein's polynomials have proven to be commonplace for economic and statistical problems in which shape constraints must be applied (Wang & Ghosh 2012); As Bernstein's polynomial degree tends to infinity, it converges uniformly to any continuous function at intervals [0,1].

Determining as a spatial Gaussian process leads to smoothing across spatial locations. Two approaches may be considered for model estimation. The complete model can be estimated with a Monte Carlo Markov Chain Sampling Program (MCMC) in small and moderate-sized samples. Reich et al. (2011) developed an approximate model that first fits quantile regressions separately over some quantiles and then uses the Bayesian spatial model to construct the quantile process during these initial estimates. In terms of mean square error, coverage, and power, complete and approximate models perform better than classical quantile regression. MCMC algorithms are described in the Appendix to Reich et al. (2011). The Bayesian quantile regression model eventually provides a tool for examining changes in a conditional distribution while accepting a high degree of flexibility in periods of covariate effects. Covariates affect higher sequences and the distribution shape; these effects change across locations and multiple levels. The model is calculated using the BSquare package of R software (Smith et al. 2013).

Table 1 displays the descriptive statistics for three stations. The Arazkouseh station has a higher average monthly discharge compared with the other two stations. Additionally, the discharge data from this station have a larger deviation from the average data, as indicated by a larger standard deviation. All three stations have positively skewed data with positive kurtosis. The positive skewness is due to the presence of high and very high discharge data in the data series, which is more pronounced in the Arazkouseh station. On the other hand, positive kurtosis implies an excessive concentration of data around the mean, with the highest kurtosis observed in the Tamar station.

This research is based on seasonal data, and its results for different seasons are as follows.

For Tamar and Arazkouseh stations, the slope of changes in the extreme upper and upper quantiles was much higher than in the lower quantiles, indicating a noticeable decrease in very high discharge values. However, the trend slopes of different quantiles in Galikesh station showed values close to each other in almost all studied quantiles.

The covariate effect in different quantiles was measured with a lower degree of uncertainty for the Tamar station and the lower and median quantiles at the Galikesh station. Comparing the results of the frequentist quantile regression and the Bayesian quantile regression models based on the slope and intercept of the quantiles, it can be concluded that in the Bayesian model, the trend slope values in different quantiles are estimated more uniformly compared with the frequentist approach. Additionally, the Bayesian approach shows smoother skewness and kurtosis curves than the frequentist quantile model.

Table 2 shows the regression coefficients (slope and intercept) and statistical significance of the trends. Arazkouseh and Galikesh stations have a significant decreasing trend in most selected quantiles, while the Tamar station has only had a significant trend in a limited number of quantiles. The highest significant trends are found in the upper quantiles of the Arazkouseh and Galikesh stations. For example, the extreme upper quantile of 0.95 decreased with a slope of −0.31 (equivalent to a decrease of 0.31 m³/s per year) at Arazkouseh station and with a slope of −0.17 (equivalent to a decrease of 0.17 m³/s per year) at Galikesh station.

Table 3 provides statistical significance and regression coefficients for summer trends. Most selected quantiles across various stations did not have statistically significant trends, except for extremely lower quantiles (quantile < 0.5). At Arazkouseh, the 0.05, 0.1, and 0.25 quantiles showed a significant increasing trend, indicating that discharge values equivalent to those quantiles in summer have increased by 0.0017, 0.002, and 0.0026 m³/s per year. Conversely, lower quantiles (0.1 and 0.25) at Galikesh and Tamar stations showed a significant decreasing trend with negative slopes. Therefore, it can be stated that the 0.1 and 0.25 quantiles of discharge at Tamar station have decreased by 0.0008 and 0.0032 m³/s per year, respectively, and the 0.25 quantile of discharge at Galikesh station has decreased by 0.0093 m³/s per year. The most significant and noticeable estimated slope was found at the 0.95 quantile in the Tamar station, showing an increasing positive slope of 0.12 m³/s per year.

Table 4 displays the statistical significance of the trend during autumn. The findings indicate that most of the quantiles for all three stations showed no significant trend. Specifically, the low discharge values in the Araskouseh (quantile 0.05), Tamar (quantile 0.1), and Galikash (quantile 0.1) stations have significantly decreased by 0.0045, 0.0055, and 0.0063 m³/s per year, respectively. On the other hand, the high discharge values in Arazkouseh (quantile 0.9) and Galikesh (quantile 0.75) stations have significantly decreased by 0.069 and 0.008 m³/s per year, respectively, while it has increased by 0.038 m³/s per year in Tamar station (quantile 0.95).

Based on the data presented in Table 5, it appears that there is no significant trend in any of the selected quantiles at the Galikesh station. However, at the Tamar station, there was a notable decrease in the lower values of discharge, specifically in quantiles 0.05 and 0.1, with trend slopes of 0.0048 and 0.0064 m³/s per year, respectively. Similarly, at the Arazkouseh station, the upper 0.75 quantile exhibited a significant decreasing trend with a slope of −0.09 m³/s per year. It is worth noting that most of the slopes observed were relatively minor and insignificant, particularly during the autumn season, suggesting limited changes in discharge during this period.

The research findings suggest that quantile regression models offer a more comprehensive understanding of hydrological and meteorological variables than non-parametric tests. In particular, the Bayesian quantile regression model has been found to detect trends that were previously unnoticed by other trend detection methods (Tharu & Dhakal 2020).

A comparison of the results obtained from quantile regression models and Bayesian quantile regression reveals that all studied stations exhibited a decreasing trend across all quantiles during the spring season. This trend was more pronounced in the upper and extreme upper quantiles. Conversely, during the summer season, lower quantiles showed almost negligible slopes for all stations, while upper quantiles, particularly the extreme upper ones, exhibited a significant negative slope at the Arazkouseh station and a positive slope at Tamar and Galikesh stations. In the autumn season, quantile regression showed a negative trend in the upper quantiles for the Arazkouseh station and a significant positive trend for the Tamar station. However, the Galikesh station only showed a noteworthy positive trend in extreme upper quantiles. Bayesian quantile regression showed similar trends in all quantiles of Tamar and Galikesh stations. Finally, during the winter season, the upper quantiles in Arazkouseh and Galikesh stations showed a negative trend, while the upper quantiles of the Tamar station showed a significant positive trend. However, Bayesian quantile regression indicated that the trends in all three stations were approximately equal and negative.

Recent studies have highlighted that the construction of the Golestan Dam has resulted in a decrease in average monthly discharge and an increase in maximum flow indicators in Gonbad and Qazaqali stations in the Gorganroud basin (Daiechini et al. 2020). Further, the maximum flow indicators in the Aq Qhala station have increased after the construction of the Voshmgir Dam. H. Ahmadi et al. (2020) and (M. Ahmadi et al. 2020) suggest that rainfall-runoff in the Gorganroud region, under the SSP8.5 scenario, displays an increasing trend for the future, with maximum discharge occurring in April and May. They further found that an increase in extreme rainfall could lead to more frequent flood events in the region. Additionally, Norouzi Nazar et al. (2023) have projected that base flow in the Gorganroud watershed would increase by 166.7% and 77.2% during the periods of 2021–2040 and 2040–2060, respectively, due to climate and land use changes. They also note that runoff is expected to increase by 54.2% and 41% during the same periods, respectively.

The analysis of flow data has revealed that changes in discharge distribution have affected several quantiles above river flow, particularly during monsoon months, in stations such as Tamar and Galikesh, leading to more flooding. Therefore, quantile regression is preferred as it can reveal the temporal dependencies of the variable in question for the mean value and its quantiles (Abbas et al. 2019). Saadi et al. (2017) used the quantile regression model to examine the distributional changes of annual monsoon rainfall and river flow for 31 stations in Sarawak and found that upward trends were mainly in the lower quantiles and in January and December in the upper quantiles. Similarly, Yves et al. (2019) found that severe floods had an increasing trend only in a limited number of basins, while most areas exhibited a decreasing trend in their study of flood trends and their drivers in 171 basins in the south of France. According to Pasch & Engelke (2022), the flood risk assessment requires an examination of changes in high quantiles. To predict floods and very high discharges more effectively, authorities can use the quantile regression method. In some cases, there has been a decreasing trend in high discharge values for certain stations, such as Arazkouseh station in all seasons, which may lead to problems with available water supply. Additionally, changes in flow patterns can pose significant threats to animal and plant species and have adverse effects on the environment. Nourali (2023) suggested that accurate flood forecasting and flood control planning must be considered in response to sudden floods in the Tamar basin area located in the Gorganrood basin in Iran. The results of Donyaii & Sarraf (2021) showed that discharge and temperature values in the Gorganroud basin will experience decreasing and increasing trends, respectively, at the end of the 21st century. Meanwhile, precipitation will not change much. Yves et al. (2019) indicated that the decrease in flood events could be due to the decrease in soil moisture caused by an increase in temperature, evaporation, and transpiration, and a decrease in rainfall in the region. Daiechini et al. (2020) stated that the construction of the dam and changes in climatic variables mostly influenced the change in hydrological indicators of the flow. Khoie et al. (2023) demonstrated that more than 60% of the changes in streamflow in the Gorganroud basin were due to land use changes.

When comparing how frequentist quantile and Bayesian quantile curves are fitted, we find that the estimates of the slope of the quantile curve in the frequentist and Bayesian models are different in some cases. These differences mostly appear when the slope of the quantile curve in the desired quantile is estimated with frequentist quantile regression, causing the crossing of adjacent quantile lines or fitting different slopes with adjacent polynomials. Such estimates have lower accuracy and validity. Tan & Shao (2017) showed that the Bayesian quantile regression model produces smoother confidence intervals of trend coefficients at extremely high and low quantile levels. It also estimates a parameter of uncertainty more accurately than the frequentist quantile regression approach (Kozumi & Kobayashi 2011). Therefore, the Bayesian regression model modifies the invalid measurements of frequentist quantile regression and provides valid estimates for all possible quantiles.

This study investigated changes in discharge values over 50 years for the Arazkouseh, Tamar, and Galikesh stations of the Gorganroud basin on a seasonal scale. Two new regression models, frequentist quantile regression and Bayesian quantile regression, were used to determine trends. The Bayesian quantile regression model provided more accurate results as it considered information in all discharge quantiles. These models are important for flood and drought management as they show the status of changes in different discharges.

The study found that high discharge values had different decreasing and increasing trends, with most cases showing negative and decreasing trends, especially in the higher extreme values of discharge. Tamar and Galikesh stations had an increasing trend in summer, which could result in summer floods. Although the discharge is decreasing in most quantiles, the flood months (upper quantiles) are increasing. This finding is consistent with recent reports of floods in the region. The decreasing trend in the Gorganroud basin, which will become more severe in the future due to global climate changes, along with the increase in the water demand of the region for various purposes due to population growth and expansion of industries, is a serious warning for politicians, planners, and local managers to prevent the possible occurrence of water crisis in the region in the future with appropriate planning. Future research should quantify the climatic changes in land use in the decreasing and increasing trends of discharge in the Gorganroud basin.

This research project has been supported by the Gorgan University of Agricultural Sciences and Natural Resources (grant no. 98-412-8). The authors are grateful to Gorgan University of Agricultural Sciences and Natural Resources for providing financial support for conducting this research.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.