The effect of karst system occurrence on flood peaks in small watersheds, southwest China

The increasing magnitude and frequency of flood disasters result in an increasing number of’ human deaths and economic losses. Accurate flood simulation has been a hot topic in hydrological sciences as it is very helpful for flood warning, flood control, and flood disaster reduction. Commonly used hydrological models are often divided into two categories: the lumped hydrological model and the distributed hydrological model. The lumped hydrological model uses an approximate method to describe the runoff process without considering the physical mechanism of the hydrological phenomenon (Pan 2018). Therefore, distributed hydrological models with clear physical meaning parameters and accurate descriptions of the hydrological process have been developed and are widely used in hydrology (Gumindoga et al. 2017; Wang 2018a; Busico et al. 2020; Ferguson & Fenner 2020; Khaleghi et al. 2020). It is well known that karst fissures, sinkholes, and skylights are widely distributed in the karst basin, constituting a unique three-dimensional spatial system and making the hydrological process more complex. However, certain approaches for flood simulation in the karst basin remain limited. To improve the accuracy of flood simulation in the karst basin, researchers usually modify the original lumped or distributed hydrological model, as summarized in Table 1. As can be seen from Table 1, the improved hydrological models are only suitable for specific basins and have no universality. In addition, previous studies have shown that the HEC-HMS model has good accuracy for flood simulation in a small watershed (Kamali et al. 2013; Silva et al. 2014). In China, some studies have indicated that the HEC-HMS model is the optimum model for flood simulation in a small basin (Zhao et al. 2018; Zhang et al. 2019).

For southwest China, floods are the main disaster during the flood season. However, most of the basins are karst basins, with complex topography and geological conditions (Pu 2009). Karst peak clusters, sinkholes, and skylights are widely distributed, forming a unique three-dimensional spatial system. Thus, the mechanism of flood formation is more complex, and the simulations of the flooding are more difficult than those in non-karst basins. Therefore, the objective of this study is to add a karst reservoir unit to the HEC-HMS model to form a flood simulation tool for a small karst basin, the Xiajia basin in southwest China. Furthermore, this study tries to determine whether the added karst reservoir unit has a certain effect on the HEC-HMS model performance. The findings of this study can provide a useful insight for flood warning and management in karst areas.

The Xiajia basin is located in Lingyun County, Baise city, Guangxi, China, with a total catchment area of 799.20 km² and an average annual rainfall of 1,708.3 mm. Information for each station is shown in Table 2, and the location of the Xiajia basin is shown in Figure 1. There is a natural karst cave – water source cave – located on the southeast side of Lingyun County, and most of the runoff in the upper reaches of the cave is the underground undercurrent. The middle part of the Xiajia basin is hilly landform, while the northeast side has a karst peak-cluster depression topography. The basin is a typical karst basin in southwest China (Huang 1995).

First, the DEM data with a resolution of 30 m was downloaded from www.gscloud.cn/. Then, HEC-GeoHMS was applied to obtain the sub-basins and river channel of the Xiajia basin, as shown in Figure 2 and Table 3. W01–W12 represent 12 sub-basin units, and R01–R05 represent four river channel units. In addition, soil data and land cover data with resolutions of 1 × 1 km were downloaded from the HWSD and IGBP_LUCC database, respectively, and the land cover and soil type of the Xiajia basin were obtained (Figure 3). The common methods for runoff calculation are the SCS curve method (Fang et al. 2007) and the initial losses-lateral losses method (Zhang et al. 2017). In this study, it was found that the SCS curve method had better results in the sub-basins W02–W04 and W06–W12, while the initial losses-lateral losses method had better results only in the sub-basins W01, W03 and W05. Therefore, the SCS curve method was adopted. As suggested by Ge (2016) and Wang (2018b), the Snyder unit line method was used for direct runoff confluence calculation, and the exponential attenuation method was applied to the base-flow separation. Finally, the Muskingum method was employed to calculate the river confluence due to its better applicability, fewer parameters, and easy calibration (Li 2020). The principle flow of the model is shown in Figure 4 (Li & He 2015).

A sub-basin is divided into karst areas and non-karst areas, and for karst areas, the karst features (natural karst cave, etc.) are considered as karst reservoirs, and the structure is shown in Figure 5. Suppose there is a sub-basin W_a, the sub-basin adjacent to W_a is W_b, and the direction of water flow is from W_a to W_b. Then, the flow yield in the karst area of W_a can be estimated by the karst reservoir (VB₁), and the water flowing from W_a to W_b is calculated by using a linear reservoir (VB₂) and the unit line method. The principle of the improved HEC-HMS model (Cheng 1988) is shown in Figure 6, and the constructions for karst reservoirs and linear reservoirs are as follows.

The discharge and precipitation data (2002–2018) used in this study were recorded by Chaoli, Donghe, Lingyun and Xiajia stations and provided by the Chengbi River Reservoir Management Bureau. Wang et al. (2019) showed that the determination of hydrological model parameters is closely related to the flood peak and the flood volume of the flood used during the calibration. A calibration using variable size flood samples is helpful to obtain more stable and representative model parameters. Therefore, 31 different size flood samples were applied to calibrate the model in this study. Among the flood samples, the maximum flood peak was 358 m³/s (number 20100628), the minimum flood peak was 101 m³/s (number 20050614), the average value was 179 m³/s, and the standard deviation was 69 m³/s. Take the maximum measured flood peak (358 m³/s) as the reference and take one-third (119 m³/s) and two-thirds (239 m³/s) as the bases for flood classification. Finally, the samples can be divided into three categories: category (i) (Q_m < 119 m³/s), category (ii) (119 m³/s < Q_m < 239 m³/s), and category (iii) (Q_m > 239 m³/s). The flood and their corresponding rainfall samples are shown in the appendix (Supplementary material, Table A1). In addition, to study the influence of the rainfall data with different centers on the model simulation results, the watershed was divided into four different rainfall centers and the division results are shown in Figure 7 (Sun et al. 2014).

The relative error (Re) of the flood peak, the of the flood volume, the deviation of the peak time ⁠, and the Nash efficiency coefficient (N_SE) (Moriasi et al. 2007) were selected as the evaluation indices of the model simulation accuracy. A simulation with Re ≤ 20%, ⁠, and N_SE ≥ 0.5 can be considered as a qualified simulation. Additionally, the square of R (R²) and the mean absolute error (MAE) were applied to give a better understanding of the overall goodness-of-fitness. The optimization parameters obtained by the Nelder–Mead method (Zettel & Hitzmann 2017) were substituted into the initial HEC-HMS model to simulate 31 floods in the calibration period, and the parameters are shown in the appendix (Supplementary material, Table A3). Finally, it is shown that the simulation result was poor, and the qualified rate (⁠⁠, n is the number of qualified simulations, N is the number of total simulations) was only 38.71%. Due to space limitations, only the simulations with three different size flood samples numbered 20080926 (category i), 20050620 (category ii), and 20100628 (category iii) were analysed, and the results are shown in Table 5 and Figure 8. For the simulation with flood sample numbered 20080926, the simulated curve and measured curve were close, and the peak time was the same, but the simulated flood peak and flood volume were slightly larger than the measured values. After the flood peak, a small amount of precipitation also caused obvious fluctuation in the simulated flow. The values of the flood peak and flood volume were 23.53 and 40.65%; the was 0; and the N_SE, R², and MAE were 0.38, 0.96 and 18.9, respectively. For the simulation with the flood sample numbered 20050620, the simulated flood peak was obviously larger than the measured value, and the simulated flood volume was slightly larger than the measured value. After the flood peak, there was also a fluctuation in the large flow caused by small precipitation. The values of the flood peak and flood volume were 76.09 and 10.14%, the was –11, and the N_SE, R², and MAE were –1.5, 0.82 and 26.3. For the simulation with the flood sample numbered 20100628, the simulated flood peak and flood volume were far beyond the measured values. The values of the flood peak and flood volumes were 149.44 and 182.17%, the was 4, and the N_SE, R², and MAE were –15.6, 0.69 and 284.2, respectively.

The precipitation and flood samples used in the calibration for the improved HEC-HMS model are the same as for the initial HEC-HMS model. The simulation results are summarized in Table 6. Overall, 16 of the 31 simulations are qualified, and the qualified rate is 51.61%. Among the unqualified simulations, six simulations are unqualified for the of the flood peak, seven simulations are unqualified for the ⁠, eight simulations are unqualified for the of the flood volume, and seven simulations are unqualified for the Nash efficiency coefficient. In addition, the spatial variation of the evaluation indices is shown in Figure 9. Six precipitation centers fall in zone I, and the corresponding qualified rate of the flood simulation is 66.67%, the NSE value is 0.93, and R² is 0.94. Five precipitation centers fall in zone II, and the corresponding qualified rate of the flood simulation is 40.00%, the NSE value is 0.44, and R² is 0.91. Ten precipitation centers fall in zone III, and the corresponding qualified rate of the flood simulation is 70.00%, the NSE value is 0.88, and the R² value is 0.90. Ten precipitation centers fall in zone IV, and the corresponding qualified rate of the flood simulation is 60.00%, the NSE value is 0.83, and R² is 0.88. Thus, the correlation between the flood simulation effectiveness and precipitation center is not significant, and the difference in the precipitation center has little effect on the flood simulation result. Additionally, according to precipitation forcing of the People's Republic of China (GB/T 28592–2012), 24-hour rainfall can be divided into six grades: light rain, moderate rain, heavy rain, rainstorms, heavy rainstorms, and extraordinary rainstorms. It was found that among the precipitation samples, one sample was heavy rain, and the corresponding flood simulation result was unqualified. Eight samples were rainstorms, and the maximum precipitation was 96.74 mm, the minimum precipitation was 54.18 mm, the average was 76.03 mm, the standard deviation was 12.62, and the corresponding flood simulation qualified rate was 87.5%. Twenty samples were heavy rainstorms, and the maximum and minimum precipitation were 239.46 and 102.71 mm, respectively. The average was 160.17 mm, the standard deviation was 41.38, and the corresponding flood simulation qualified rate was 55.0%. Two samples were extraordinary rainstorms, and the maximum precipitation was 263.63 mm, the minimum precipitation was 253.18 mm, the average was 258.41 mm, the standard deviation was 7.39, and the corresponding flood simulation qualified rate was 50.0%. For flood peak type, 19 floods were single-peak, with concentrated corresponding precipitation processes and large flood peak values, and the qualified rate was 73.68%. Twelve floods were multi-peak, with relatively scattered corresponding precipitation processes, and the qualified rate was 41.67%. The flood simulation effectiveness was related to the flood peak type, and the flood simulation effectiveness of the single-peak floods was better than that of the multi-peak floods.

Based on the comprehensive analysis of the precipitation level, precipitation center and flood peak type, it can be concluded that the HEC-HMS model constructed in this study has better simulation results under the conditions of heavy rainfall and concentrated precipitation processes.

The precipitation and flood samples used in the validation for the improved HEC-HMS model included 18 different size flood samples and their corresponding precipitations, and the simulation results are summarized in Table 7. Overall, 11 simulations were qualified, and the qualified rate was 61.11%. Three floods were selected as the representative, and the variation of flood hydrograph and karst storage is shown in Figure 10. It can be seen from Figure 10 that the karst storage was variable, which shows that the karst reservoir unit can simulate the stagnation of the runoff process in the karst area well. In addition, the spatial variation of the evaluation indices is shown in Figure 11. One precipitation center falls in zone I, and the corresponding qualified rate of the flood simulation was 100% with the NSE and R² values of 1.00. Two precipitation centers fall in zone II, and the corresponding qualified rate of the flood simulation was 100.00% with the NSE and R² values of 1.00. Seven precipitation centers fall in zone III, and the corresponding qualified rate of the flood simulation was 71.43% with the NSE value of 0.94 and R² of 0.96. Eight precipitation centers fall in zone IV, and the corresponding qualified rate of the flood simulation was 37.50% with the NSE value of 0.55, and R² of 0.85. For precipitation level, three precipitation samples were rainstorms, the maximum precipitation was 95.49 mm, the minimum precipitation was 76.63 mm, the average was 84.90 mm, the standard deviation was 10.00, and the simulation qualified rate was 66.67%. Thirteen precipitation samples were heavy rainstorms with the maximum and minimum values of 203.06 and 109.58 mm. The average was 141.45 mm, the standard deviation was 23.83, and the simulation qualified rate was 61.54%. Two precipitation samples were extraordinary rainstorms with maximum and minimum values of 301.38 and 273.51 mm, respectively. The standard deviation was 19.71, and the simulation qualified rate was 50.00%. For flood peak type, nine flood samples were single-peak with concentrated corresponding precipitation processes and large flood peak values, and the qualified rate was 77.78%. Nine flood samples were multi-peak with relatively scattered corresponding precipitation processes, and the qualified rate was 44.44%. The flood simulation effectiveness is related to the flood peak type, and the flood simulation effectiveness of the single-peak floods is better than that of the multi-peak floods.

The calibration results of the original HEC-HMS model are generally poor, and the simulated values overestimate the measured values. To determine the reasons why the simulation results are poor, three hypotheses are proposed: (1) there are incorrect parameter values; (2) the studied basin is not a closed basin due to the influence of karst landforms, leading to a large number of flows flowing down to other basins through karst pipelines (Xu et al. 2017); and (3) the karst cave, called a water source cave, has a detention effect on floods (Chen 2018).

For assumption (1), the Morse method (Song et al. 2015a) is applied to analyze the sensitivity of the parameters, and the results are shown in Table 8 and Figure 12. Taking the floods numbered 20080926, 20050620, and 20100628 as examples, for flood peak, CN is the most sensitive parameter (e ranges from 5.818 to 13.451), followed by Cp, t_p and I_a. For peak time, the sensitivity of each parameter is low, among which, C_p is the most sensitive parameter (e ranges from 0.050 to 0.100), followed by K and t_p. For flood volume, CN is the most sensitive parameter (e ranges from 1.596 to 2.438), followed by C_p, t_p, k and r. The above sensitivity analysis results imply that adjusting the sensitive parameters can change the simulation output and therefore improve the simulation accuracy of the model to a certain extent. However, an over-adjustment of the model parameters may make it lose its actual physical meaning. For example, when the CN value in this study drops to approximately one-third of its initial value, the simulated flood peak value of the large flood is equivalent to the measured value, but the CN value will be far from that of the land condition represented by the original land cover and soil hydrological group. Chen et al. (2016) and Liao (2017) also proved that it is unreasonable to significantly reduce the CN value. Therefore, it is believed that the values of the parameters in this study are reasonable and that they are not the main cause of the low accuracy simulation.

For assumption (2), an analysis of water balance is applied to find out whether the studied basin is a closed basin or not, and the results are shown in Table 9. In terms of multi-year averages, the water of the studied basin is unbalanced, which implies that the studied basin is not a closed basin. However, the value of implies that some water from the adjacent basin may flow into the studied basin through karst pipes. Thus, assumption (2) is not true.

For assumption (3), as stated above, a natural karst cave, called a water source cave, is located in sub-basin W06. The calibration results show that the simulated values with the flood samples that belong to category i are similar to the measured values, while the simulated values with the flood samples that belong to category iii are far beyond the measured values. This phenomenon can be explained from the perspective of detention storage, that is, when the incoming flow is small, it is the free outflow of the open channel, and as the discharge increases, it gradually changes to pressure orifice outflow, and therefore, some water cannot be discharged and stagnate in the karst area during the flooding time, and finally this water will gradually discharge after the flood recedes. This suggests that the drainage basin in the karst area above the water source cave has a stagnant effect on water flow. Thus, assumption (3) is reasonable.

Additionally, it should be noted that although the simulation result has been improved after adding a reservoir unit to the original HEC-HEM model, the simulation accuracy of this study is still low compared to that in non-karst basins (Feng 2016). The reason may be that the hydrological processes in karst basins are extremely complex, and not all the karst features and their effects can be considered in the model improvement, resulting in some errors. In addition, it is found that the model will have better simulation results when the rainfall is heavy and the precipitation process is concentrated. The reason may lie in that when the rainfall runoff is heavy, the rapid flow increases. Meanwhile, the evaporation and the karst underground pipe network consumption are relatively small. Therefore, the simulation accuracy is relatively high.

An accurate flood simulation is crucial for flood warning, flood control, and flood disaster reduction, especially in karst areas with more complex flood features. Therefore, the objective of this study was to add a karst reservoir unit to the HEC-HMS model to analyze the impact of karst system occurrence on flood peaks and form a more accurate model for flood simulation. Some major findings were as follows: The HEC-HMS model with a karst reservoir unit outperformed the original model. After adding the karst reservoir unit, the qualified rate was increased by 22.40 and 22.58% during calibration and validation, respectively. Furthermore, for flood peaks and flood volumes, the parameter CN was the most sensitive parameter, followed by C_p and t_p. Also, a comprehensive analysis of the flood sample data (Tables A1 and A2 in the Supplementary material) and simulation results (Tables 6 and 7) showed that the results were better for the simulation with a rainstorm and a single-peak type flood than that with a non-rainstorm and a multi-peak type flood. Additionally, the simulation results were better when the gravity center of precipitation fell in zone III, zone I, and zone IV.

The novelty of this study lay in adding a karst reservoir unit to the HEC-HMS model to quantify the effect of the karst system on the surface flow in a karst basin. Finally, the results implied that karst systems had an impact on hydrology by effectively restraining flood peaks and increasing underground runoff, which agreed with the findings of Bailly-Comte et al. (2008), Palanisamy & Workman (2014), and Zhao et al. (2019). However, some shortcomings remain to be addressed. Only one hydrological station is located in the studied basin, meaning the model has to be calibrated and validated by the data from one hydrological station. Further hydrological research is required to study the application of a reservoir unit in flood modelling at karst basins. Additionally, the differences in runoff parameters between karst areas and non-karst areas, the alternation of open flow and buried flow, and the mutual conversion between surface water and groundwater need to be further studied, and some useful advice can be found in Sun (2015). Besides, only one basin in Guangxi Province was considered in this study, and it is not realistic to extrapolate the findings of this study to larger spatial scales. Thus, more universal results need to be studied.

The authors are grateful for the support of the Natural Science Foundation of China (51969004 and 51579059), the National Key Research and Development Program of China (2017YFC1502405 and 2016YFC0401303), the Guangxi Natural Science Foundation of China (2017GXNSFAA198361), and the Innovation Project of Guangxi Graduate Education (YCBZ2019022).

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