Analysis of parameter uncertainty in SWAT model using a Bayesian Box–Cox transformation three-level factorial analysis method: a case of Naryn River Basin

Acronyms

SWAT

Soil and water assessment tool

DREAM

Differential evolutionary adaptive metropolis

MCMC

Markov chain Monte Carlo

PDF

Probability density function

LH-OAT

The Latin hypercube one factor at a time

TFA

Three-level factorial analysis

ArcGIS

Geographic information system

HRU

Hydrologically comparable unit

R²

Coefficient of determination

BCT

Box–Cox (BC) transformation

NSE

Nash–Sutcliffe model efficiency

PBIAS

Percentage bias

Notations

θ

The set of aggregate parameters in SWAT

Q_obs

Set of observe discharge values

Q_sim

Set of simulation discharge values

P(Q_obs|θ)

Prior parameter probability distribution

P(θ)

Prior information for the parameters set θ

L(Q_obs|θ)

Probability density of the observed discharge values

P(Q_obs)

Normalizing constant

f(Q_obs|θ)

Probability density of the likelihood

ε_i(θ)

The residual on day i

σ_ε,i(θ)

Residual standard deviation

λ

Transformation parameters

ε_ijk

A portion of uncertainty due to probability

y_ijk

Measured reaction for A in ith level and B in jth level for kth replicate

(τβ)_ji

Effect of the interaction between τ_i and β_j

τ_i

Effect of factor A in i level

β_j

Effect of factor B in j level

μ

Total average effect

SS

The sums of squares

SS_AB

Sum of squares for the interaction between factors A and B

SS_T

Total sum of squares

SS_E

Sum of squares error

Central Asia is a typical region that is characterized by the influence from the dominant Westerlies with the contributions from the East Asian and Indian monsoons. The combination of the unique atmospheric circulation and its distinctive geographical location make it one of the most ecologically and environmentally fragile regions in the world (e.g., the Qinghai-Tibet Plateau). Except for high-altitude regions, most of Central Asia is characterized by arid and semi-arid climates. Hydrological modeling in these regions is challenging due to the unique hydroclimatic variables and the sparse distribution of meteorological stations (Kan et al. 2017; Didovets et al. 2021). As the two largest rivers in Central Asia, the Amu Darya and the Syr Darya both terminate in the Aral Sea. Since the dissolution of the Soviet Union in 1991, unreasonable management and allocation of water resources between upstream countries (e.g., Kyrgyzstan and Tajikistan) and downstream countries (e.g., Kazakhstan, Uzbekistan, and Turkmenistan) have caused conflicts and highlighted water scarcity and water security issues in the downstream areas (Schöne et al. 2018; Zheng et al. 2019). The Naryn River originates in the southern foothills of the Tien Shan Mountains in eastern Kyrgyzstan and flows from east to west through Kyrgyzstan into Uzbekistan before joining the Kara River to form the Syr Darya. However, it is still unclear how the water flow responds to natural and anthropogenic factors. Hydrological modeling of the Naryn River Basin is crucial for analyzing water balance, forecasting long-range streamflow, and predicting flood risk. By analyzing the hydrological processes of the basin, decision makers can formulate appropriate policies to alleviate the pressure of water scarcity in such an arid and information-deficient region.

Hydrological models, classified into physically based and machine learning (ML) models, were extensively used for analyzing water balance, forecasting long-range streamflow, predicting real-time flood, and investigating climate-change impact in watershed management (Yuan et al. 2018; Cheng et al. 2020; Ansell & Valle 2021; Adnan et al. 2022). Adnan et al. (2021) proposed a novel hybrid method, ELM-PSOGWO, which integrates particle swarm optimization (PSO) and gray wolf optimization (GWO) with extreme learning machine (ELM) for predicting monthly runoff in the Mangla watershed in northern Pakistan. Ikram et al. (2022a) developed a hybrid ML method based on the extended marine predators algorithm (EMPA) and artificial neural network (ANN-EMPA) for flow estimation in the upstream of the Indus River Basin. Hussainzada & Lee (2022) used the Soil and Water Assessment Tool (SWAT) model to determine the hydrological processes and assess the water resource capacity of the Barhab River Basin (BRB) in northern Afghanistan. Generally, the ML models can directly explore the relationship between hydrological cycle variables and runoff in theoretical system models without considering the involved physical processes. The majority of the parameters in the SWAT model are not readily accessible. They can only be derived by calibration against input and output data of the watershed response (Vrugt et al. 2008; Mengistu et al. 2019; Liu et al. 2022; Schuster-Wallace et al. 2022). Besides, hydrological models often encounter substantial uncertainties with respect to the input data, initial and boundary conditions, model structure, and parameters due to insufficient observation data and difference in spatiotemporal scale (Sidle 2021; Wang et al. 2021). These uncertainties can limit the reliability and practical value of the model simulation results while they have not been evaluated in practice (Baig et al. 2022; Han & Morrison 2022).

Over the past decades, a number of research works were conducted for investigating the effect of parameter uncertainty on modeling performance. Farsi & Mahjouri (2019) analyzed the uncertainty of the parameters of the models using the differential evolutionary adaptive metropolis (DREAM) algorithm to decrease the uncertainty of model estimations. Diao et al. (2021) elaborated the solution mechanism of the standard DREAM algorithm, and the algorithm was applied for the optimal operation model of the reservoir group in Jialing River. Xu et al. (2022) used DREAM algorithm to analyze the parameters’ posterior distributions and the nutrient loadings’ simulated uncertainty ranges for agriculture in SWAT model. Wang et al. (2022) used DREAM to obtain the posterior distribution of SWAT model parameters to consider parameter uncertainty. In general, the DREAM algorithm can obtain the accurate posterior distribution of parameters through large-scale simulations (Yang et al. 2008; Zhang et al. 2020; Liu et al. 2021). However, it still uses the traditional statistical assumption of residual homoscedasticity. In the hydrologic modeling, errors are not homoscedastic but heteroskedastic since the size of errors tends to rise with the observed flow magnitude.

Box–Cox transformation (BCT) is effective for eliminating the heteroscedasticity of residuals in the simulation, which can improve the normality and linearity of the data series. BCT has achieved good normality in hydrological frequency analysis and calculations (Qiao et al. 2022). Guo et al. (2022) utilized residual-based modeling to generate and evaluate probabilistic hydrologic and water quality predictions at 18 sites in the United States; where the BCT scheme provided the best predictions of uncertainty for all case studies. Zhang et al. (2022) proposed a new hybrid prediction model called variational mode decomposition (VMD)-BC-Elman by introducing the VMD and BCT into the Elman neural network; the results demonstrated that BCT could remove the skewness of the original streamflow sequence, transform skewed candidate input variables into normal distribution, and enhance the correlation between input variables and output variables. The BCT can be employed to eliminate the heteroskedasticity of model residuals, and the residuals can be assumed to adhere to a generalized error distribution (Cheng et al. 2014, 2018). However, even after translation, model data have trouble showing how several uncertain factors influence model results.

In the uncertainty analysis, it is very important to identify how much each parameter affects the final result. Sensitivity analysis is often used for quantifying the influence of various parameters on the model output, specifically assessing the contribution of parameters to overall model output uncertainty (Uusitalo et al. 2015; Koo et al. 2020). Van Griensven & Meixner (2003) proposed an improved sensitivity analysis approach, named as Latin Hypercube One factor At a Time (LH-OAT), for the sensitivity analyses for a large number of parameters. However, the LH-OAT approach cannot quantify the effect of parameter interactions on model output, it can only offer a sensitivity ranking of parameters (Liu et al. 2017). In comparison, factorial analysis (FA) can quantify the main effect of a single parameter and different-level interactions of multiple input parameters on the system performance (Gao et al. 2021). As an extension of the traditional FA, three-level factorial analysis (TFA) is useful for detecting the curvilinear relationship between the parameters and the response (Zhang et al. 2016; Liu et al. 2019; Jia et al. 2020). Unfortunately, no research work has been reported on the integration of DREAM, BCT, and TFA techniques into SWAT hydrological model to disclose the influence of parameter uncertainty and their interactions on the runoff for the river basin in Central Asia.

Due to the importance of water resources management, hydrological simulation has always been a topic of concern for researchers. The majority of the parameters in the SWAT model are not readily accessible. They can only be derived by calibration against input and output data of the watershed response (Vrugt et al. 2008; Mengistu et al. 2019). Parameter errors are always presented throughout the calibration and derivation processes. In addition, the previous hydrological studies in the Naryn River Basin mainly focus on using physical or ML models to simulate river runoff, and there is little research on revealing parameter uncertainty effecting hydrological processes through physical model (i.e., SWAT). Due to the scarcity of data in the Naryn River Basin, the availability of hydrological and meteorological observation data limits streamflow simulation performance. It is critical to undertake a thorough investigation of the uncertainties in the SWAT model simulation process for enhancing the accuracy of the hydrological simulation for the Naryn River Basin. In order to fill these research gaps, further studies are needed to explore the uncertainties associated with hydrological processes in the Naryn River Basin.

This study aims to develop a Bayesian Box–Cox transformation three-level factorial analysis (BBC-TFA) method for assessing parameter uncertainties and their interactions on hydrological model responses. BBC-TFA incorporates BCT, DREAM, and TFA techniques within a general framework. BBC-TFA will then be applied to the Naryn River Basin to quantify the uncertainty of hydrological parameters in the SWAT model. The main novelty and contribution are as follows: (i) this is the first attempt of integrating BCT, DREAM, and TFA techniques into SWAT model; BBC-TFA can remove heteroskedasticity and get a more trustworthy posterior distribution of the parameters by introducing BCT into the DREAM method; (ii) BBC-TFA can obtain the parameters’ contribution; (iii) BBC-TFA is capable of identifying the influence of sensitive parameters on runoff; (iv) BBC-TFA is first applied to the Naryn River Basin (an information-poor area) to reveal the effects of parameter uncertainty. The results can assist relevant researchers in establishing effective hydrological models and help decision makers understand the physical factors associated with the watershed and provide theoretical support for studies in similar watersheds.

The organization of the remaining sections of this paper is as follows: Section 2 describes the framework of BBC-TFA and the evaluation criteria used in this study. Section 3 describes the study area and the data used in this paper. Section 4 discusses the results of this study. Section 5 concludes the study and discusses the limitations and potential extensions of the BBC-TFA method.

The residuals of the original data are of particular importance, not the converted data, thus we need to add a fraction in the likelihood function to account for this.

By dividing the sum of squares for each input factor (SS_A, SS_B, and SS_AB) by the total sum of squares for all parameters (SS_T), the sensitivity index of each factor was obtained. The model's parameters and their values are selected using the DREAM method, which includes the posterior distribution's shape, statistical properties of the posterior parameters, and correlation coefficients. The TFA has a wider range of parameter values, which allows for a more accurate assessment of parameters’ fuzziness and the links between them.

The CGIAR Consortium for Spatial Information (CGIAR-CSI) provided digital elevation model (DEM) grid sizes of 90 m × 90 m (https://srtm.csi.cgiar.org). Digital topographic data were utilized to divide watersheds into sub-basins, build river networks, and identify water bodies and outlets. From the Harmonized World Soil Database (HWSD) of the Food and Agriculture Organization of the United Nations, soil attributes are retrieved; the database has a spatial resolution of about 1,000 m (http://www.fao.org/soils-portal). Data on land use characteristics were taken straight from SWAT model database. Only 1951–1974 data of runoff from the Kurgan hydrographic station have been found because of the scarcity of data (https://portal.grdc.bafg.de). The historical (1951–1974) daily temperature data are downloaded from the website of Princeton Global Meteorological Forcing Dataset (PGMFD, http://hydrology.princeton.edu). PGMFD is derived by merging a number of observation-based datasets with the National Center for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) reanalysis. The temperature in the dataset is a perfect reflection of the actual value (Ahmed et al. 2019). The daily precipitation grids for the period 1951–1974 are taken from the Asian Precipitation-Highly Resolved Observational Data Integration Towards Evaluation (APHRODITE) dataset (http://aphrodite.st.hirosaki-u.ac.jp). APHRODITE is a long-term continental daily precipitation product that has been shown to be superior to other models and can be used directly (Ghimire et al. 2018).

This study uses the ArcSWAT2012, which is an ArcGIS extension and interface for SWAT, for model configuration and parameterization. A 100,000-ha threshold was used for the DEM data to split the study region into 49 sub-basins, which were then segmented into 541 hydrologically comparable units (HRUs). The thresholds of 10% for land use, 10% for soil, and 5% for slope were established, and the multiple slope option was used in the creation of each hydrologically comparable unit (HRU), with three classes of 0–10%, 10–20%, and >20%. There is a distinct set of land use types, soil types, and slope classes represented by each HRU in the SWAT model. Based on available data from 1951 to 1974, the runoff was simulated. The study ensures that the climate data used for calibration and validation do not differ significantly, with wet, mild, and dry years occurring in both periods (Arnold et al. 2012; Polong et al. 2023). In this study, the model is calibrated and validated using available data from 1951 to 1974. The calibration period is from 1954 to 1969, during which the average monthly flow is 403.1 mm/s. The validation period uses monthly data sets from 1970 to 1974, during which the average monthly flow is 375.4 mm/s.

In the SWAT model, there are several HRUs in an area, which might result in thousands of distinct parameter changes during model calibration, making it impossible for the researcher to alter them one by one. As a result, distinct parameters must be aggregated. To aggregate parameters in this article, two sorts of prefixes, ‘V_’ and ‘R_’, are employed. v__ means the default parameter is replaced by a given value, r__ means the existing parameter value is multiplied by (1 + a given value) (for the detailed parameter aggregation scheme, please refer to Abbaspour et al. 2007). Table 1 contains brief explanations of each parameter and its default values. According to the suggestion with Abbaspour et al. (2017), it is advisable to consider as few parameters as possible for model validation and calibration. In more than 200 parameter files, 19 parameters were chosen based on reviews of the literature, and the LH-OAT method was used to determine sensitivity rankings (Arnold et al. 2013; Abbaspour 2015).

Table 2 displays the correlation coefficients for the calculated parameters. Based on the data, the strongest positive correlations are seen between TIMP and SMFMX, with a value of 0.7. This is because both TIMP and SMFMX influence the rate at which snowmelt runoff is produced. SMFMX is the melt factor for snow on June 21. The Naryn River Basin is in the Northern Hemisphere, according to the file ‘SOIL AND WATER1 ASSESSMENT TOOL INPUT/OUTPUT FILE DOCUMENTATION’, the parameter SMFMX means the maximum melt factor in this study. It has affirmative effects on the generation of snowmelt runoff. TIMP also has a positive influence on snowmelt runoff production. It is therefore desirable to delve further into the interrelationship between the factors.

In this study, parameter uncertainty based on the SWAT model is also analyzed using the untransformed DREAM method. Figure 4(a) and 4(b) show the probability density histograms of the posterior distributions of the model parameters for the transformed DREAM method and the untransformed case, respectively, and the results show that there are big dissimilarities in the posterior distributions of the parameters considering the residual variance. Figure 4(b) shows that when applying the traditional DREAM method, the model parameters are distributed in a larger parameter interval. This is because the BCT used a similar logarithmic transformation to both observed and simulated data. It could greatly increase the data's normalcy, with the obvious result being that the bigger numbers make a bigger difference. This also made the chain of DREAM approaches need additional iterations to determine the optimal parameter distribution. The transformed DREAM method considering residual variance has a smaller distribution of parameter posteriors and more concentrated values of parameter posteriors than the untransformed DREAM method. Liu et al. (2017) used the DREAM algorithm to approximate the posterior distribution of model parameters through Bayesian inference. The results indicate that the posterior parameter distribution has a large range.

Figure 6 shows the ranges of model parameters and input uncertainty in discharge simulations with and without residual variance, respectively. Using the DREAM method, to calibrate and validate, PBIAS values are 11.8 and 11.7. Using the transformed DREAM method, the values of PBIAS for calibration and validation are −2 and −2.4. The results show that the range of uncertainty in the discharge simulations due to model input and parameter uncertainty is reduced considering the heteroskedasticity case, which may be caused by the reduction in the range of model parameters due to the introduction of the heteroskedasticity model. The BCT can improve the accuracy of the DREAM method for hydrological model uncertainty prediction and obtain more accurate parameter values and more representative uncertainty prediction intervals.

The parameter CH_K2 is effective hydraulic conductivity in main channel alluvium. It can affect the mutual recharge between groundwater and river (Zhang et al. 2016). Soil permeability coefficient is positively correlated with hydraulic conductivity, the higher the soil hydraulic conductivity, the more rainfall infiltration, and the smaller the runoff generated on the surface. The larger the value of CH_K2, the more water will permeate into groundwater and reduce the river discharge. In this study, it has a weak negative effect on SPF but has an insignificant effect on NPF and AAF. It means that the infiltration water in this area is relatively small when the water is abundant.

The parameter SMFMX can account for the impact of snowpack density on snowmelt like SMFMN. The SWAT model calculates the amount of snowmelt water in the watershed using the degree-day factor method, which works under the assumption that the potential snowmelt rate follows a date-dependent sinusoidal function, with the highest expected on June 21 and the lowest expected on December 21 (Liu et al. 2020a). The parameter has a positive slope for all responses and has a positive impact. This means that snowmelt contributes significantly to runoff in the basin.

The delay factor TIMP regulates how much the previous day's snow pile temperature affects the present snow temperature. The average daily temperature has a greater impact on snow pile temperature as TIMP gets closer to 1. Assuming the user sets a temperature threshold, the snow pile will only melt after it reaches or passes that point (Abbas et al. 2019). This parameter has a significant effect on all responses; it means that the amount of snow that melts daily fluctuates significantly from 1 day to the next. This extrapolation implies that the temperature of the basin varies significantly. These results demonstrate the significance of soil conditions and snowmelt, as represented by the model, for the study basin in response to the model. Studies on soil quality and snowmelt should be expanded if the model performance is to be enhanced.

In this study, a BBC-TFA method has been developed and applied to assess the uncertainty impacts of hydrological parameters on runoff in the Naryn River Basin (Central Asia). The BBC-TFA method has shown high accuracy in runoff simulation. The uncertainties of four sensitive parameters, such as snowpack temperature lag factor (TIMP), effective hydraulic conductivity in main channel alluvium (CH_K2), melt factor for snow on June 21 (SMFMX), and initial soil conservation service (SCS) runoff curve number for moisture condition II (CN2), have been investigated, and their effects and interactions on different hydrological responses have also been explored. SWAT provides the SCS curve number procedure (USDA SCS 1972) for estimating surface runoff, the SCS curve number procedure is an empirical model that came into common use in the 1950s.

Several findings are (i) the parameter CN2 has a positive effect on NPF, implying that the Naryn River Basin is mountainous (with steep topography and concentrated summer rainfall) and tends to trigger an infiltration-excess runoff; (ii) the parameter CH_K2 affects the mutual recharge between groundwater and rivers. The negative effect on all responses in this study is weak; this disclose that, in the Naryn River Basin, only a limited quantity of surface runoff infiltrates as groundwater; (iii) the parameter TIMP has a significant effect on all responses, which means that the temperature of the basin varies significantly; (iv) individual and interactive effects of parameter SMFMX and TIMP have a significant effect on the hydrological response, implying that snowmelt is a significant contributor to runoff in the Naryn River Basin; (v) the BCT can effectively eliminate the residual variance of the model output values from the observed values and obtain a lower range of uncertainty prediction intervals, enhancing the correlation between input variables and output variables; (vi) CN2 has the highest contribution to NPF with a rate of 88.22%, it means that the Naryn River Basin is at risk of flooding, but preparations can be made for water storage; TIMP has the highest contribution to both SPF and AAF with rates of 79.69 and 53.70%, respectively; it is indicated that the Naryn River has abundant glacier resources that can be appropriately developed and/or protected. Based on the above findings, the relationships between important parameters and the Naryn River Basin can be clarified, thus improving the simulation/prediction capability of the hydrological model for water resources in the Naryn River Basin.

This study is the first parametric uncertainty analysis of the hydrological modeling of the Naryn River Basin by the BBC-TFA method, and there are some limitations that can be further improved. First, BBC-TFA was only applied to the Naryn River Basin, and its application to more regions, especially for data-scarce areas, could be explored. Second, the Naryn River Basin is located in the Central Asian mountains, where the low density of observation points and data scarcity can increase the uncertainty of model performance and results. Therefore, to improve the simulation accuracy of the model, satellite remote sensing inversion can be used to increase the runoff sequence in future work. Third, this study used BCT to process data, resulting in a small range of posterior parameter distribution. In the future research work, new data preprocessing techniques can be attempted to eliminate homoscedasticity without changing the distribution of the data itself (Ikram et al. 2022b). Future research might take into account these extensions to enhance parameter and hydrologic modeling uncertainty estimations.

This research was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDA20060302), the Fund for Innovative Research Group of the National Natural Science Foundation of China (52221003) and the Start-up Project for Advanced Talent of Xiamen University of Technology in the First Half of 2022 (YKJ22011R). The authors are grateful to the editors and the anonymous reviewers for their insightful comments and suggestions.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.