https://orcid.org/0000-0002-8136-0299
Abstract
Hydrological simulation in karst areas is of great importance and challenge. It is a practical way to enhance the performance of existing hydrological models in karst areas by coupling karst modules that represent hydrological processes in these areas. The near-surface critical zone structure affects runoff generation in karst areas significantly and its complex hydrological processes could be simplified with threshold behaviors. This study proposed a three-thresholds-based karst runoff generation module (3T-KRGM), which used three reservoirs to represent water storage in the soil zone, soil–epikarst interface, and epikarst zone. The 3T-KRGM is coupled with the Xinanjiang (XAJ) model to extend the applicability of the model to karst areas. Both the improved XAJ model and the original XAJ model were used in the Shibantang watershed, which is a typical karst watershed located in southwest China. The results indicate that the performance of daily discharge simulations was obviously improved by introducing the 3T-KRGM. In addition, both the parameter sensitivity analysis and baseflow simulation demonstrate that the 3T-KRGM is rational in structure. The 3T-KRGM could also be easily coupled into other hydrological models, thus benefiting the hydrological simulation in karst areas.
HIGHLIGHTS
We proposed a three-thresholds-based karst runoff generation module to improve model performance in karst areas.
The theoretical basis of the module is a modification of a three-thresholds mechanism revealed by hillslope experiments.
We coupled the module into the Xinanjiang model to extend its applicability in karst areas.
The advantages of the module for hydrological simulation were verified from various aspects.
INTRODUCTION
Karst is the geological process in which water causes the dissolution of soluble rocks such as carbonate and sulfate (Hartmann et al. 2014; Malagò et al. 2016). Karst forms distinctive surface and subsurface features, including karren, sinkholes, karst conduits, and underground caves, which makes the hydrological process in karst areas highly complex (Bailly-Comte et al. 2009; Zhou et al. 2019; Chen et al. 2022). Globally, about 15% of the non-ice continental surface is karst areas, with 16.5% of the population living in these areas (Goldscheider et al. 2020). China has one of the largest continuous karst areas (located in the southwest of China) and its karst percentage is as high as 26.5% (Zhang et al. 2011; Goldscheider et al. 2020). Karst waters, which include carbonate aquifers, provide fresh water to around 25% of the population globally (Ford & Williams 2007). Hydrological simulations in karst catchments are critical not only for water resource management but also for flash flood prevention (Bonacci et al. 2006; Yang et al. 2022) and ecosystem services evaluation (Tian et al. 2016).
Hydrological models are essential tools for hydrological simulations. The hydrological models employed in karst areas can be split into two groups, data-driven and process-driven models, in line with other regions (Hartmann et al. 2014; Chen et al. 2022). The data-driven models, such as long short-term memory (LSTM) (Fang & Shao 2022), convert rainfall and other input into runoff directly based on transfer functions. The advantages of data-driven models are few parameters and considerable adaptability to data situations and regions (Zhou et al. 2019). The lack of physical process representation is where data-driven models have long been criticized (Kuczera & Mroczkowski 1998). The process-driven models are further classified into distributed, semi-distributed, and lumped models (Paul et al. 2021). The distributed and semi-distributed models discretize the karst watershed into several basic units (e.g., grids or sub-basins) and assign certain governing equations, parameters, and state variables to each basic unit according to the internal karst landform (Hartmann et al. 2014). The advantage of distributed and semi-distributed models is that they offer a thorough insight into hydrological processes within a karst watershed while also considering temporal and spatial variability (Zhang et al. 2011; Xu et al. 2020; Li et al. 2021). However, the high dataset requirements and considerable computing expense restrict the application of distributed and semi-distributed models in practice (Fatichi et al. 2016; Simmons et al. 2020). The lumped models treat the whole karst watershed as a basic unit and further conceptualize its general physical processes (Hartmann et al. 2014; Zhou et al. 2019). Because of their simplicity and effectiveness, lumped models are still widely employed in practical applications (Dubois et al. 2020; Duran et al. 2020).
Lumped hydrological models use reservoirs and linkage functions between these reservoirs to provide a simplified interpretation of hydrological processes in areas of consideration (Tritz et al. 2011). The most widely used lumped hydrological model in China is the Xinanjiang (XAJ) model (Zhao 1992). Instead of building a lumped model for karst areas from scratch, a practical alternative is to adapt an existing lumped model by coupling modules that depict karst hydrological processes (Yang et al. 2022). One reason is the existence of non-karst areas in a karst-dominated watershed. The runoff generation process in karst and non-karst areas should be considered separately in a karst watershed (Zhou et al. 2019), which could be described by the coupled karst modules and original modules of the existing lumped model correspondingly. Another reason is that existing code and expertise in the model setup, parameter calibration, and uncertainty quantification can be reused.
The near-surface critical zone structure in karst areas denotes the soil–epikarst system and is usually characterized by high permeability and considerable water storage capacity (Perrin et al. 2003b). It plays an important role in runoff generation as it controls the conversion of precipitation into fast-flow and slow-flow components (Jukić & Denić-Jukić 2009). The soil and epikarst zone were treated separately in the Moisture Balance (MB) model (Jukić & Denić-Jukić 2009), and effective rainfalls were calculated based on the moisture balance of these two zones. However, the linkage between soil and epikarst zone is no more than the direct recharge to epikarst after the saturation of the soil zone in this model. Tritz et al. (2011) conceptualized the soil and epikarst zone as a whole and used two thresholds to represent the activation and deactivation of fast-flow paths in the soil–epikarst system. The interflow generated in the soil–epikarst system is a critical hydrological mechanism in subtropical karst areas (Fu et al. 2015), which could partly be attributed to the infiltration rate of the soil–epikarst interface (SEI) (Wang et al. 2020). The infiltration rate of the SEI functions as a threshold because it could reflect the change in the hydrological flux pattern. The complex hydrological processes in the soil–epikarst system could be simplified by identifying threshold behaviors (Wang et al. 2022). Thus, the near-surface critical zone structure in karst areas and related threshold behaviors should be considered in karst modules.
The objective of this study was to develop a new karst runoff generation module capable of representing hydrological processes in the soil–epikarst system explicitly. The proposed karst runoff generation module was then coupled into the XAJ model to improve model performance in karst areas. The Shibantang watershed in southwest China was selected as a study area. The parameter sensitivity of the improved XAJ model was performed to diagnose the proposed module and provide guidance for model calibration. The improved XAJ model and original XAJ model were used to simulate daily runoff in the study area, and a comparison was conducted based on the observed discharge series.
KARST RUNOFF GENERATION MODULE
Perceptual model
The typical subsurface critical zone structure in karst areas (GS refers to the ground surface; SEI denotes soil–epikarst interface; ETI refers to the interface of epikarst and transmission zone; GWL denotes groundwater level; SFI denotes the interface of saturation zone and fresh bedrock).
The typical subsurface critical zone structure in karst areas (GS refers to the ground surface; SEI denotes soil–epikarst interface; ETI refers to the interface of epikarst and transmission zone; GWL denotes groundwater level; SFI denotes the interface of saturation zone and fresh bedrock).
The infiltrated water quickly reaches the SEI due to the shallow soil depth and the activation of the preferential pathway after the onset of precipitation (Wang et al. 2020). The infiltration-excess surface flow is rarely seen in karst areas as the infiltration rate exceeds the precipitation intensity in most cases. The infiltration rate of the SEI controls the amount of water that enters the epikarst. It is usually smaller than that of the soil zone, which causes water accumulation at the SEI (Fu et al. 2015). According to the fill-spill hypothesis (Tromp-van Meerveld & McDonnell 2006), the interflow would begin once the depressions at the SEI were filled and the saturated areas at the SEI were connected (Wang et al. 2020). With the rise of temporary water tables originating from the SEI, the fast and concentrated recharge to the saturation zone through karst fracture would be activated gradually (Tritz et al. 2011). Meanwhile, the saturated surface runoff would begin as the temporary water tables reached the ground surface (Wang et al. 2020, 2022). Rapid groundwater runoff would result from the fast and concentrated recharge to the saturation zone, which is characterized by rapid speed and huge quantity (Chen et al. 2022). The saturated surface runoff and rapid groundwater runoff could be treated together as rapid runoff. The water infiltrated into the epikarst would produce epikarst seepage runoff in the epikarst zone (Hartmann et al. 2014).
Three thresholds exist in the above-mentioned processes (Wang et al. 2020). The first threshold (T1) is the infiltration rate of the SEI, and the water would accumulate at the interface once the threshold is satisfied. The second threshold (T2) is the quantity of precipitation required to fill the depression and form connectivity at the SEI, after which the interflow runoff would be generated. The last threshold (T3) is the quantity of precipitation required to both saturate the soil zone and activate the fast and concentrated recharge path, after which it would generate rapid runoff.
Mathematical model
The structure of the proposed three-thresholds-based karst runoff generation module (3T-KRGM). The module consists of three reservoirs: soil–epikarst interface (SEI), epikarst, and soil reservoir.
The structure of the proposed three-thresholds-based karst runoff generation module (3T-KRGM). The module consists of three reservoirs: soil–epikarst interface (SEI), epikarst, and soil reservoir.
SEI reservoir
denotes the storage of the SEI reservoir (mm), 
is net precipitation (mm), 
denotes the initial storage of the SEI reservoir (mm), 
is the water percolated into the epikarst reservoir (mm), and 
is the water run into the soil reservoir (mm). The portion of precipitation that exceeds evapotranspiration capacity forms net precipitation
. It could be expressed as:
denotes evapotranspiration capacity (mm), and 
denotes net evapotranspiration (mm). To minimize the number of parameters, the way to calculate 
is set to the same as the XAJ model, which is expressed as:
is evapotranspiration coefficient (–), and 
is observed pan evaporation (mm).
is calculated with the GR4J type function (Perrin et al. 2003a) and expressed as:
denotes the storage of epikarst reservoir (mm), 
denotes the capacity of epikarst reservoir (mm), and tanh is the hyperbolic tangent function. According to the mechanism of T1, only the part of 
that exceeds 
could be retained in the SEI reservoir.
exceeds a threshold 
. The equation of 
could be expressed as:
Epikarst reservoir
denotes the initial value of epikarst reservoir (mm), 
is epikarst seepage runoff (mm). The equation of 
could be expressed as:
Soil reservoir
denotes the storage of soil reservoir (mm), 
denotes the initial storage of soil reservoir (mm), 
and 
are the interflow runoff and rapid runoff (mm), and 
is the evapotranspiration from soil reservoir (mm). The equation of 
could be expressed as:
denotes the interflow coefficient on karst pervious areas (–). The equation of 
could be expressed as:
is the capacity of the soil reservoir (mm), 
denotes the watershed-averaged amount of water stored in the soil zone (mm), 
denotes the saturation water storage capacity on the karst area curve exponent (–). The equation of 
could be expressed as:
could be expressed as:
METHODS
The XAJ model
The XAJ model (Zhao 1992) is a conceptual hydrological model developed in 1973 and perfected in the 1980s. The critical assumption of the XAJ model is that no runoff is generated until soil moisture reaches field capacity, which is appropriate in the humid and semi-humid watersheds of China (Zang et al. 2021). After the total amount of runoff is evaluated, the runoff separation procedure separates generated runoff into surface, interflow, and groundwater runoff to employ different hillslope concentration methods. The XAJ model's major characteristic is interpreting the spatial distribution of tension and free water storage using two curves (Zhao 1992), which are also the basis of runoff generation and separation calculation.
The XAJ model's inputs are watershed-averaged precipitation and observed pan evaporation, while its outputs are actual evapotranspiration from the watershed and discharge at the watershed outlet. The basic structure of the XAJ model contains four modules: evapotranspiration calculation, runoff generation, runoff separation, and flow concentration. The whole watershed is subdivided into pervious and impervious areas according to the permeability of the underlying surface. The first three modules are combined to determine the generated surface, interflow, and groundwater runoff in the pervious areas. The portion of watershed-averaged precipitation that exceeds actual evapotranspiration directly forms the surface runoff in the impervious areas. The surface runoff on the whole watershed and the interflow and groundwater runoff in the pervious areas are routed from the hillslope to the watershed outlet using the flow concentration module.
The coupling of the 3T-KRGM and the XAJ model
denoting how much of the watershed is made up of karst pervious areas. The proposed 3T-KRGM is solely applied to the karst pervious areas, whereas the corresponding modules of the XAJ model are applied to the non-karst pervious areas and impervious areas. The 3T-KRGM shares the same inputs as the XAJ model, producing 
, 
, and 
as output. The corresponding modules of the XAJ model produce surface runoff on the non-karst pervious areas (
), interflow on the non-karst pervious areas (
), groundwater on the non-karst pervious areas (
), and surface runoff on the impervious areas (
). 
, 
, and 
formed the total surface runoff (
). 
and 
formed the total interflow runoff (
). 
and 
formed the total groundwater runoff (
). 
, 
, and 
were routed to the watershed outlet and formed discharge Q by using the flow concentration module of the XAJ model. The XAJ models with and without 3T-KRGM are referred to as the improved XAJ model and the original XAJ model, respectively. The parameters of the improved XAJ model and original XAJ model are listed in Table 1. Parameters of the original and improved XAJ model (with signification and range)
| Parameter | Signification | Range |
|---|---|---|
 | The ratio of potential evapotranspiration to pan evaporation | [0.6,1.5] |
 | Evapotranspiration coefficient of deep soil layer | [0.1,0.2] |
 | Averaged tension water storage capacity of upper soil layer | [5,30] |
 | Averaged tension water storage capacity of lower soil layer | [60,90] |
 | Averaged tension water storage capacity of deep soil layer | [15,40] |
 | The ratio of impervious watershed area | [0.01,0.2] |
 | The exponent of the tension water storage capacity curve | [0.1,0.4] |
 | Averaged soil layer free water storage capacity | [10,30] |
 | The exponent of the free water storage capacity curve | [1,1.5] |
 | Outflow coefficient of interflow from free water storage | [0.1,0.55] |
 | Outflow coefficient of groundwater from free water storage | 0.7 −  |
 | Recession coefficient of interflow storage | [0.5,0.9] |
 | Recession coefficient of groundwater storage | [0.95,0.998] |
 | Recession coefficient of water storage in river network | [0.1,0.9] |
 | Lag time | [1,10] |
 | The ratio of karst pervious area | [0.1,0.9] |
 | The capacity of the epikarst reservoir | [10,100] |
 | The capacity of the soil reservoir | [10,80] |
 | Saturation water storage capacity on the karst area curve exponent | [0.1,0.4] |
 | Interflow coefficient on karst pervious areas | [0.01,0.8] |
 | Threshold of the SEI | [0,10] |
| Parameter | Signification | Range |
|---|---|---|
| | The ratio of potential evapotranspiration to pan evaporation | [0.6,1.5] |
| | Evapotranspiration coefficient of deep soil layer | [0.1,0.2] |
| | Averaged tension water storage capacity of upper soil layer | [5,30] |
| | Averaged tension water storage capacity of lower soil layer | [60,90] |
| | Averaged tension water storage capacity of deep soil layer | [15,40] |
| | The ratio of impervious watershed area | [0.01,0.2] |
| | The exponent of the tension water storage capacity curve | [0.1,0.4] |
| | Averaged soil layer free water storage capacity | [10,30] |
| | The exponent of the free water storage capacity curve | [1,1.5] |
| | Outflow coefficient of interflow from free water storage | [0.1,0.55] |
| | Outflow coefficient of groundwater from free water storage | 0.7 − |
| | Recession coefficient of interflow storage | [0.5,0.9] |
| | Recession coefficient of groundwater storage | [0.95,0.998] |
| | Recession coefficient of water storage in river network | [0.1,0.9] |
| | Lag time | [1,10] |
| | The ratio of karst pervious area | [0.1,0.9] |
| | The capacity of the epikarst reservoir | [10,100] |
| | The capacity of the soil reservoir | [10,80] |
| | Saturation water storage capacity on the karst area curve exponent | [0.1,0.4] |
| | Interflow coefficient on karst pervious areas | [0.01,0.8] |
| | Threshold of the SEI | [0,10] |
The schematic diagram of the coupling of 3T-KRGM and the XAJ model. The dashed arrow denotes runoff components generated from the 3T-KRGM. The gray square denotes the summation of two runoff components.
The schematic diagram of the coupling of 3T-KRGM and the XAJ model. The dashed arrow denotes runoff components generated from the 3T-KRGM. The gray square denotes the summation of two runoff components.
Study area and data
The study area's location, topography, and river system with hydrological and precipitation stations.
The study area's location, topography, and river system with hydrological and precipitation stations.
This study includes daily precipitation series from seven precipitation stations (Baina, Gantang, Hongjiadu, Liulong, Shachang, Shawo, and Shibantang), daily pan evaporation, and daily discharge series from Shibantang hydrological station. All data series used are from January 2006 to December 2018. The Annual Hydrological Report of China (Volume VI, Book X) is the data source for the daily precipitation and discharge, while the National Meteorological Science Data Center of China (http://data.cma.cn/) is the source for the daily pan evaporation. The Thiessen polygon method was used to calculate the watershed-averaged daily precipitation from all seven daily precipitation series. The elevation data of the watershed come from the SRTM 90 m DEM dataset (Geospatial Data Cloud of China, http://www.gscloud.cn/).
Model setup and parameter calibration
is simulated flood peak (m3/s), 
is observed flood peak (m3/s), 
is the index of datum, N is the total amount of the data, 
and 
represent the i-th simulated and observed discharge, 
is the average value of observed discharge. The Shuffled Complex Evolution (SCE-UA) algorithm (Duan et al. 1992) was used to optimize model parameters in the calibration period, and its objective function is 
.Sensitivity analysis of the improved XAJ model
The sensitivity analysis of the improved XAJ model is to identify sensitive model parameters that can be used for model calibration and diagnostic evaluation (Pianosi et al. 2016). This study adopted the widely used PAWN method (Pianosi & Wagener 2015; Puy et al. 2020) for global sensitivity analysis. The PAWN method treats the sensitivity of an input factor as the distance between the unconditional cumulative distribution function (CDF) of output obtained when all inputs vary simultaneously, and the conditional CDF obtained when varying all inputs but that input. The Kolmogorov–Smirnov (KS) statistic, i.e., the maximum absolute difference, is used to calculate the distance between the unconditional and conditional CDF. Given the complexity of the model input–output relationship, the KS statistic is calculated approximately rather than analytically. The approximation strategy detailed by Pianosi & Wagener (2018) was utilized, with the range of each input variable divided into n equally spaced intervals. The median of KS statistics across intervals was used to indicate the sensitivity of model output to model parameters. In this study, PAWN was conducted with 50,000 randomly generated parameter sets using the Latin hypercube sampling methods. The KGE performance metric was used to summarize the discharge series produced by the improved XAJ model into a scalar for sensitivity analysis. The value of n was set to 10, 20, and 40 to validate the robustness of the sensitivity analysis result according to the suggestion provided by Pianosi & Wagener (2018). The sensitivity analysis was conducted using the SALib Python library (Iwanaga et al. 2022).
RESULTS AND DISCUSSION
Sensitivity analysis of the improved model
The sensitivity of the improved XAJ model's parameters by the PAWN method. The parameters are in descending order using the result of n = 20.
The sensitivity of the improved XAJ model's parameters by the PAWN method. The parameters are in descending order using the result of n = 20.
The most sensitive parameter is 
, which regulates the evapotranspiration capacity and, consequently, the quantity of total runoff. The parameters L and 
, which are used to route the runoff (originating from both karst and non-karst area) stored in the river channel to the watershed outlet, have ranks of 2 and 6, respectively. The three parameters mentioned above come from the original XAJ model and affect the improved XAJ model. 
, L, and 
are all sensitive model parameters according to previous studies (Zhao 1992). The fourth sensitive parameter is 
, which represents the ratio of impervious area within a watershed. In impervious areas, surface runoff is identical to net precipitation and is an important source of flood peaks. Thus, 
should be a sensitive parameter according to existing model practice. 
is the third sensitive parameter, which represents the ratio of the karst pervious areas of the watershed. It supports the importance of the 3T-KRGM because this module is only used on the karst pervious areas. 
and 
mainly controls the interflow and rapid runoff generated on the karst area (
and 
), with ranks of 5 and 7, respectively. It indicates 
and 
play an important role in the simulation of the Shibantang watershed.
Comparison of performances and parameters of the improved and original model
The optimized parameters of the original and improved XAJ model
| Parameter |  |  |  |  |  |  |  |  |  |  |  |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Original | 0.86 | 0.14 | 9.9 | 76.8 | 27.5 | 0.10 | 0.23 | 14.9 | 1.14 | 0.4 | 0.3 |
| Improved | 1.10 | 0.14 | 23.1 | 75.2 | 27.0 | 0.11 | 0.21 | 22.1 | 1.28 | 0.3 | 0.4 |
| Parameter |  |  |  |  |  |  |  |  |  |  | |
| Original | 0.62 | 0.96 | 0.57 | 1 | – | – | – | – | – | – | |
| Improved | 0.73 | 0.97 | 0.58 | 1 | 0.81 | 36.1 | 41.6 | 0.26 | 0.03 | 0.85 |
| Parameter | | | | | | | | | | | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Original | 0.86 | 0.14 | 9.9 | 76.8 | 27.5 | 0.10 | 0.23 | 14.9 | 1.14 | 0.4 | 0.3 |
| Improved | 1.10 | 0.14 | 23.1 | 75.2 | 27.0 | 0.11 | 0.21 | 22.1 | 1.28 | 0.3 | 0.4 |
| Parameter | | | | | | | | | | | |
| Original | 0.62 | 0.96 | 0.57 | 1 | – | – | – | – | – | – | |
| Improved | 0.73 | 0.97 | 0.58 | 1 | 0.81 | 36.1 | 41.6 | 0.26 | 0.03 | 0.85 |
The performance metrics of the original and improved XAJ model in the calibration and validation period
| Model | Period | NSE | KGE | FPRE (%) | FVRE (%) | RMSE |
|---|---|---|---|---|---|---|
| Original | Calibration | 0.86 | 0.79 | −21.9 | 15.3 | 14.9 |
| Validation | 0.76 | 0.83 | 2.4 | 14.0 | 18.8 | |
| Improved | Calibration | 0.90 | 0.94 | −8.7 | −2.0 | 12.2 |
| Validation | 0.86 | 0.88 | 14.0 | −6.0 | 14.1 |
| Model | Period | NSE | KGE | FPRE (%) | FVRE (%) | RMSE |
|---|---|---|---|---|---|---|
| Original | Calibration | 0.86 | 0.79 | −21.9 | 15.3 | 14.9 |
| Validation | 0.76 | 0.83 | 2.4 | 14.0 | 18.8 | |
| Improved | Calibration | 0.90 | 0.94 | −8.7 | −2.0 | 12.2 |
| Validation | 0.86 | 0.88 | 14.0 | −6.0 | 14.1 |
The comparison of simulation results of the original and improved XAJ model. (a) Simulated daily stream discharge hydrographs accompanied by observations; (b,c) the annual performance metrics of simulation results of these two models.
The comparison of simulation results of the original and improved XAJ model. (a) Simulated daily stream discharge hydrographs accompanied by observations; (b,c) the annual performance metrics of simulation results of these two models.
The parameter 
decides the evapotranspiration capacity according to Equation (4). The parameter 
controls the water storage for direct evapotranspiration. A smaller 
and smaller 
combination means more runoff would be generated, which is revealed in Figure 6(b), as the FVRE of the original XAJ model is bigger than the improved XAJ model for most cases. The optimized values of 
, 
, and 
for the original XAJ model were different from those of the improved XAJ model, which indicates more surface runoff, more interflow runoff, and faster interflow concentration velocity. The comparison of model parameter values demonstrates that the rapid flow runoff component in karst areas was compensated with more surface runoff, more and faster interflow in the original XAJ model. The compensation is at the expense of a higher FVRE, slower hydrograph metrics (NSE and KGE), and, most importantly, the right result for the wrong reasons (Goswami et al. 2010). The SCE-UA optimized parameter value of the 
was 0.81, which is quite similar to the karst landscape ratio of the Shibantang watershed (Xu et al. 2020). It explains the rationality of model parameters of the improved XAJ model to a certain extent.
As can be seen from Table 3, the improved XAJ model was better than the original XAJ model for both calibration and validation periods in terms of the NSE, KGE, FVRE, and RMSE performance metrics. The boxplots of Figure 6(b) represent the distribution of annual performance metric values with 25%, 50% (median), and 75% quantile, as well as outliers. Compared to the original XAJ model, the NSE and KGE values of the improved XAJ model are close to 1, and the FVRE value of the improved XAJ model is close to 0, as shown in Figure 6(b). According to Figure 6(c), the RMSE values of the improved XAJ model are lower than the original XAJ model over the majority of years. Meanwhile, these two models each have advantages and disadvantages in the calibration and validation period for the FPRE metric, as shown in Table 3. The absolute mean values of the annually evaluated FPRE of the improved XAJ model for the calibration period and the validation period are 16.2 and 15.4%, respectively. In comparison, the corresponding values for the original XAJ model are 22.0 and 17.9%. It revealed that the improved XAJ model was better than the original XAJ model in terms of the FPRE metric. To conclude, the improved XAJ model performs better than the original XAJ model in the Shibantang watershed during the simulation period, indicating the model performance improvement of the proposed 3T-KRGM.
Comparison of the simulation result of the improved and original model in flood events
Simulated daily discharge series of the original and improved XAJ model for three selected typical flood events.
Simulated daily discharge series of the original and improved XAJ model for three selected typical flood events.
Comparison of simulation-based baseflow and UKIH-estimated baseflow
and 
after slope concentration with linear reservoirs denoted as 
and 
was treated as simulation-based baseflow produced by the improved XAJ model. The simulation-based baseflow did not include the fast groundwater flow component in karst areas because it was a part of 
. Based on the observed daily discharge series, the estimated baseflow was generated using the United Kingdom Institute of Hydrology (UKIH) method (Piggott et al. 2005). This method separated the base flow by identifying the turning points of the hydrograph and piecewise linear interpolation. Figure 8 depicts comparisons between simulation-based baseflow and UKIH-estimated baseflow. 
The comparisons of simulation-based baseflow of the improved XAJ model and the UKIH-estimated baseflow for three selected flood events. R in the chart on the right denotes the Pearson correlation coefficient of simulation-based baseflow and UKIH-estimated baseflow.
The comparisons of simulation-based baseflow of the improved XAJ model and the UKIH-estimated baseflow for three selected flood events. R in the chart on the right denotes the Pearson correlation coefficient of simulation-based baseflow and UKIH-estimated baseflow.
As shown in Figure 8, the simulation-based baseflow agreed well with the UKIH-estimated baseflow in trends. The Pearson correlation coefficient between two baseflow series in flood events 20090619, 20120711, and 20180618 were 0.88, 0.92, and 0.72, respectively. Although there were some differences between simulation-based baseflow and UKIH-estimated baseflow, especially in the 20180618 flood event, we believed that the improved XAJ model could reproduce the baseflow well, considering the unavoidable uncertainty that existed in baseflow separation methods (Partington et al. 2012). The capacity of the improved XAJ model to reproduce the baseflow in karst areas further validates the rationality of the proposed 3T-KRGM.
CONCLUSIONS
A conceptual runoff generation module, called the three-thresholds-based karst runoff generation module (3T-KRGM), has been proposed to improve the performance of an existing hydrological model in karst areas. The 3T-KRGM was coupled into the XAJ model, a widely used lumped hydrological model in China. The XAJ model coupled with 3T-KRGM was referred to as the improved XAJ model, while the XAJ model not coupled with 3T-KRGM was referred to as the original XAJ model. The Shibantang watershed, with a karst ratio of 86%, was selected as the study area. The parameter sensitivity analysis of the improved XAJ model was conducted for diagnostic evaluation and model calibration. The improved XAJ model and the original XAJ model were used for the daily discharge simulation in the study area to evaluate the proposed 3T-KRGM. The main conclusions are given as follows:
- (1)
The proposed 3T-KRGM considers the near-surface karst zone structure explicitly and uses three reservoirs to represent the water storage in the soil zone, soil–epikarst, and epikarst zone. The theoretical basis of 3T-KRGM is a modification of three-threshold mechanism revealed by hillslope experiments (Wang et al. 2020) with consideration of rapid groundwater runoff. The 3T-KRGM has six parameters and produces rapid, interflow, and epikarst seepage runoff on karst pervious areas.
- (2)
The parameter sensitivity analysis result obtained using the PAWN method reveals that three parameters from the 3T-KRGM are sensitive. Among these parameters, the ratio of karst pervious area
is found to be the most sensitive. The calibrated value of 
is consistent with a reported value in the study area, which indicates the rationality of the 3T-KRGM's structure. - (3)
The improved XAJ model outperforms the original model with respect to performance metrics for the simulation period and the flood hygrograph for selected events in the study area. The improved XAJ model simulated slow-flow component matches UKIH-estimated baseflow during these events. The results show that the proposed 3T-KRGM improves the performance of existing hydrological model in karst areas and is reasonable.
Overall, the proposed three-thresholds-based karst runoff generation module (3T-KRGM) helps improve the XAJ model's performance in a karst-dominated watershed. The 3T-KRGM could also be coupled with other hydrological models, hence benefiting hydrological simulation in karst areas. Further research should be conducted to validate the proposed module with more watersheds.
ACKNOWLEDGEMENTS
This study was supported by the National Natural Science Foundation of China (NSFC) (grants, 41730750; 41877147). The authors appreciate the editors and anonymous reviewers for their critical comments and constructive suggestions, which helped improve the manuscript.
AUTHOR CONTRIBUTION STATEMENT
J.Z., G.L., and Z.L. conceptualized the study; J.Z., G.L., and Y.D. studied methodology and did formal analysis; J.Z. did the investigation and wrote the original draft; Y.D. did data curation; J.Z. and Y.H. visualized the study; G.L., Y.D., Y.H., B.L., and Z.L. wrote, reviewed, and edited the article; Z.L. and B.L. conducted funding acquisition; Z.L. supervised the work. All authors have read and agreed to the submitted version of the manuscript.
DATA AVAILABILITY STATEMENT
Data cannot be made publicly available; readers should contact the corresponding author for details.
CONFLICT OF INTEREST
The authors declare there is no conflict.
