We analyzed the characteristics of main karstic/non-karst reaches of the Lijiang River to uncover the causes behind different flood behaviors by providing a better understanding of the flood formation. Having 63 years of rainfall-runoff data and applying the HEC-HMS model, geo/hydrological features were investigated. The available reservoir capacity of karts (ARCK) was included through soil moisture accounting loss data to assess its impact. In particular, the expected instantaneous peak discharge rates/times were found largely imbalanced with generated unit hydrographs. Moreover, significant gaps among the floods’ features for different subbasins in terms of required peak modifications (2–4 times larger for mid-upstream, respectively) were mainly associated with the unique karst structure and initial condition due to various ARCK in rainy/dry seasons. Besides, notable dissimilarities between the wedge/prism storage volumes and the hydrograph’s wave traveling/receding time were observed owing to the geomorphological conditions. Although the contribution rates of drivers in karst flood formation cannot be quantitively modeled, based on our results the ARCK emerged to play a substantial role on the forecasted results, comparatively. Our results suggest that since ARCK varies, taking it into account (as initial abstraction) results in a more reliable estimation. This was underpinned by the results in which the unmodified simulations had a qualified rate of 52% accuracy on average and increased to 67.5% after the ARCK inclusion. This work adds to the body of evidence illustrating that in karst hydrology, ignoring the situational circumstances in modeling might lead to inaccuracies in flood forecasting for such dynamic watersheds.

  • Hydrological models inaccurately forecast flood features in karst basins.

  • The seasonality of available karst reservoir capacity drives flood peaks.

  • Initial conditions must be considered in model calibration for karstic areas.

Graphical Abstract

Graphical Abstract
Graphical Abstract

Natural resource availability is one of the key factors in the southern and eastern parts of China hosting a dense population and providing the basis for rapid development. Huge rivers meander through the semitropical areas of south and southeast China due to considerable amounts of annual rain. However, flood has been among the main environmental disasters in these regions, which is responsible for the loss of life, properties, and other spillover effects. According to statistics around 7000 people died in China due to floods in a decade between 2006 and 2015 (Liu & Shi 2017). Hence, studying and analyzing the flood drivers is vital for prevention measures. Although the influencing degree of different variables (i.e., land use/land cover change) on flood formation in different watersheds is a controversial topic, nevertheless, flood is basically an environmental depended phenomenon and can be driven by the surrounding's situation in watersheds (Dai et al. 2021).

Along with the intricacy of the methods, came the necessity for a more thorough understanding of the hydrological processes. Due to the nature of the hydrological phenomenon and its dependency on the variety of environmental indices, complex models have been developed and a vast range of parameters are involved (Zhou et al. 2019). However, with remarkable differences between the prediction results magnitudes when applying different hydrological formulas, it is imperative to select the right approaches for any sensible attempt to predict future requirements such as flood mitigation planning (Sharafati et al. 2020). As the impact of environmental variables on flood formation in watersheds can drastically change the subject to the ambient conditions, therefore, the model's output calibration is vital in accordance with the basin features and actual recorded data (Elmoustafa & Moussa 2010).

Although the Soil Conservation Service Unit Hydrograph (SCS-UH) and Snyder's synthetic Unit Hydrograph (SUH) are widely applicable methods, nonetheless, they both might sometimes fail to provide accurate values in terms of the hydrograph peak discharge and time estimation versus expected values in certain watersheds such as in our case. This is due to the generalization of values and coefficients used in hydrological models without taking the specific characteristics of each basin into consideration (Pan 2018). Hence, the modification of the produced hydrographs based on the observed data is obligatory (Tassew et al. 2019). Employing computer methods as lumped or distributed hydrological models in different environments (Brirhet & Benaabidate 2016) and circumstances in hydrology have been exemplified in numerous studies. Among which, the Hydrologic Engineering Canter's Hydrologic Modeling System (HEC-HMS) is a well-known flood prediction model, which has been widely used by hydrologists all around the world (Choudhury et al. 2002; Kabeja et al. 2020; Sharafati et al. 2020), especially for prediction purposes. Besides, due to its comprehensiveness and accuracy, the application of the Muskingum technique in flood routing is reported in a vast number of literatures (Tewolde & Smithers 2006; Delphi et al. 2010; Mo et al. 2021). To minimize the uncertainties and predict the temporal/spatial variations of flood hydrograph, the method takes into account the prism and wedge storages beside the hydrograph waive traveling time, which both factors are highly dependent on the watershed features and its physical descriptions.

Karst that is, in fact, residual carbonate hills (Jiang et al. 2020) is widely spread in some regions and accounted for about 15% of the lands in the world (Zhou et al. 2019), such as China, Malaysia (Sarawak), and Vietnam. In China, around one-third of the country is covered by karst, especially concentrated in the south and southwest; however, the karst flow process is yet the complicated research areas in hydrology (Cao et al. 2018; Jiang et al. 2019a, 2020; Mo et al. 2021). Among the difficulties of flood routing in karst hydrology are unknown inflow–outflow currents (to and from the study areas) through karst pipe networks that cause studying an open basin (not closed), karst cave water storage with retention/stagnation effects on stormwater, as well as buried flows, just to mention few (Mo et al. 2021). Hence, the conceptualization of karst flow is challenging. The flood's different behavior within a relatively small distance in the Lijiang catchment is a clear manifestation of such impacts in the karst environment. The bare and covered karst (via the soil or loose sediment) underground water movements in the Lijiang watershed was studied by Jiang et al. (2019a). Applying remote sensing imagery (Jia et al. 2019) and hydrological models in hydro-meteorological studies of the Lijiang River watershed has been exemplified by researchers during the past two decades (Cao et al. 2018; Xu & Dai 2018); however, the quantitative karst condition impacts, weightage, and its influencing degree on the flood characteristics have been rarely reported. In fact, the karstic watershed is a dynamic environment that has different initial conditions prior to each flood event, and therefore, no single method could effectively predict floods in such an environment for all the different scenarios. In other words, since the initial condition in such watershed is constantly changing, hence the model would require to be recalibrated for each flood. This is due to the change in the volume of the retained water in karst cracks and underground lakes and its discharge rates, which is a function of rainfall intervals, evaporation, and other ambient conditions (Dai et al. 2022).

Generally, there are three methods (empirical, conceptual, and distributed methods) for more accurately including the karstic feature's effect in the hydrological models. These methods use linear and nonlinear rainfall-runoff functions. These three categories are either using lesser parameters for ease of applicability where no much data are available (i.e., empirical), or are complex with so many parameters, data, and calibrations required (i.e., distributed), which make it difficult to apply. The conceptual models, as the third group, however, are to some extent known to be hybrid methods as an integration of the previous two (Zhou et al. 2019). One of the well-known karst models is XAJ, which is a conceptual model for humid regions. It was modified by Zhou et al. where they offered a conceptual coupled hydrologic model (K-XAJ) to simulate the rainfall-runoff in karst-dominated areas (Zhou et al. 2019). Mo et al. studied the quantitative karst effect on HEC-HMS simulation in the Lijiang area through different size floods to calibrate the model and enhanced the model outcomes via taking the karst reservoir involved as a variable, which resulted in a 13–24% improvement in model forecasting accuracy (Mo et al. 2021). However, what is missing in these models is that karst is an unsteady watershed in which if the same rainfall happens in two different times, different amounts of it will be converted to runoff as the ARCK varies over time.

This study has focused on the flood characteristics in the upstream and midstream of the Lijiang River (Guilin and Yangshuo, respectively) and aims to uncover the main drivers of the stormwater different responses in these sub-watersheds through analyzing the runoff formation process and exploring its influencing elements such as fluctuations in the available reservoir capacity of karst. The work is vital for future river management, water consumption, and decision-making in this area. Moreover, the results are important since the main income source in this karst landform is through the Lijiang River-based tourism industry, which its peak time overlaps the monsoon and the flood season (Yao & Mallik 2020).

Study area

The Lijiang River lies between 109°45′–110°56′E and 24°38′–25°55′N latitudes and longitudes, and flows from the northeast of the Guangxi autonomous region (in south China) toward the south as a part of the well-known Zhujiang River system (or the Pearl River as the third biggest river of China) (Figure 1). Located at the Yunnan-Guizhou Plateau, more than one-third of Guangxi is covered by karst landforms (Jiang et al. 2019b) with a developed karst system which is mainly around the middle and lower reaches of the Lijiang River watershed. Here, the carbonate rocks are mainly made of limestone, which are formed during the middle-upper Cambrian to the middle Triassic geological periods (Wei et al. 2018).
Figure 1

The study area with the Lijiang River situation (a) and its three subbasins (b).

Figure 1

The study area with the Lijiang River situation (a) and its three subbasins (b).

Close modal

The river originates from the hilly areas at the eastern side of Maoer mountain at its upper reaches and crosses through Guilin City, then flat areas of Yangshuo county, and finally, the flow path approaches the Pingle County and continues further south under another name (Guijiang River, as it merges other rivers) before being released into the Zhujiang River. With a total of 214 km in length and 13 major tributaries (Jia et al. 2019), Lijiang has more than 11,3000 km2 catchment area comprised of three distinguishable sub-watersheds.

First, the upstream section, which is the smallest subbasin, is Guilin and has an area of 2388 km2 (21% of the catchment area). The fan-shaped Guilin subwatershed is marked with a prominent magnitude yearly precipitation of 75 inches (1890 mm), averagely, due to its subtropical monsoon climate. Distinguishably, most parts of this region are steep and hilly, which are composed of jungles, woodlands, and timber land covers. Impacted by humid systems entering this subbasin through the South China Sea and Guangdong province, the elevated Maoer mountain causes more annual precipitations than mid- and downstream (Liu et al. 2008). Second, the midstream section of the Lijiang catchment is Yangshuo county with an area of 3425 km2 (30% of the catchment). It is more plateaued and karstic and has a yearly rainfall of 61 inches (1550 mm). Third, the downstream section of the catchment is the mountain Pingle subbasin that is the greatest among the three. The area of this section is 5484 km2 as 49% of the Lijiang catchment, which has a poor urbanization rate (comparatively) and 54 inches (1370 mm) of average yearly precipitation. Seasonality of rainfall in the Lijiang watershed has caused dramatic fluctuations in the river flow rate (as low as 12 m3/s in the dry season and up to 12,000 m3/s in the raining season) (Chen et al. 2011). The average annual temperature in the Lijiang watershed is between 17 and 20 °C. As the majority of Lijiang's scenic spots and tourist sites are located near the river path or at the confluences, therefore, when flood arises, they will be the most impacted zones. Furthermore, the residential areas (counties and cities) are majorly placed at the subbasins' pour points. With 20 AAA-level scenic spots, the classic karst landscapes of this watershed are attracting up to 50 million tourists every year (Yao & Mallik 2020).

In this work, the HEC-HMS rainfall-runoff model was applied in order to analyze the flood and its influencing factors in different subbasins of the Lijiang River. The required information and primary input data to run the model are obtained via using different software, functions, sources, and methods as presented in the following framework (Figure 2). After the model establishment and calibration, considering the ARCK, initial abstractions were introduced for validation using real data.
Figure 2

The flowchart of the applied method in the current study.

Figure 2

The flowchart of the applied method in the current study.

Close modal

Digital Elevation Model, Universal Transverse Mercator, and river network layers

The two basic essential layers of each watershed in GIS are the DEM (Digital Elevation Model) layer and the river network. These two data are provided in Arc-Hydro which is a GIS-based extension to be used later in the HEC-GeoHMS extension of GIS and finally in the HEC-HMS model. First, the DEM layer of the watershed (which is a raster data) was obtained from the US Geological Survey (USGS – available at https://earthexplorer.usgs.gov/) with an acceptable resolution of 30 × 30 m cell dimensions. After being spliced and cropped, this layer was masked in Arc-Hydro according to the watershed boundaries and based on the region's UTM zone (Universal Transverse Mercator), which specifies the study area for coordination purposes in GIS. This layer helps to classify the slopes, flow direction, and identify the conical karst spots through the Landsat imagery in Hec-geoHMS to quantitatively calculate how many percent of each subwatershed is covered by bare karst. The UTM map projection system, in fact, identifies different locations on the ground (in meters) in which the earth is divided into 60 longitudinal zones (every 6° width is considered as one zone), and Lijiang falls in the 49th zone.

Then, the river network shape file (which is vector data) for the basin was obtained from the open street map website (available at https://www.openstreetmap.org) and embedded into the Arc-Hydro extension after clipping to match the watershed boundaries. Next, for calibration, this vector data were post-processed via actual river network data as well as compared it with different maps including local maps and Google earth maps through direct observation to match the topographical situation of the river channel, geometrically.

Although GIS can automatically specify the river network from the DEM layer through considering the topography of the basin, however, in low elevation areas of the basin, any road or valley might be taken as a river; therefore, the outcome of this method must be calibrated observationally as well as comparatively. Finally, the basin and its subbasins boundaries were delineated via HEC-GeoHMS software (geospatial hydrologic modeling extensions of the GIS model) considering the river network data. The areas and the river length were also computed through this extension.

Landuse data, soil type, flow types, and CN number

Raster landuse data (the latest available images in 2013) were derived from the China Geographical information Monitoring Cloud Platform (GIM Cloud). The soil type data layer of the watershed was downloaded (https://www.fao.org/soils-portal/data) from the FAO portal. This layer would help identify the soil groups, hydrologically, to apply the infiltration rates for different groups in CN number calculations via GIS.

The remote sensing information together with the soil type and the DEM layer was used to calculate the Soil Conservation Service Curve Number (SCS-CN) grid of each subwatershed in GIS (in HEC-GeoHMS extension). The CN calculations were done taking into consideration the different landuse types of the watershed such as agricultural, industrial, water, and timberland specified through color code. Using Arc-Hydro, the geographical and situational characteristics (length of the longest river, area, the distance to the basin gravity center, CN, and slope) of the three subbasins' hydrometric locations were pulled out from GIS.

According to the identified landuse, area size, and CN numbers, the flow types in each subbasin were defined for the model proportionally. In bare karst, three types of flow can happen such as rapid direct runoff on the impermeable surface, rapid detention in underground lakes through the large size cracks and shafts, and slow leaching through smaller cracks. In covered karst, the flow types before and after saturation would be different. Here, based on the soil infiltration rate, the flow will be divided into two parts, the surface runoff will form in accordance with the CN number, and the influx part will join the underground current through the underlying karst network. Therefore, the total runoff is the summation of the runoff in these three parts (bare karst, cover karst, and non-karst), considering the area size in the sub-watershed, proportionally, in which coefficients are applied to convert the flow into the runoff for each part (Zhou et al. 2019). Lastly, the initial abstraction would be calculated based on these values, as well as ARCK that is a function of the intervals and the current amount of soil moisture. The soil moisture accounting loss data was acquired for a period of 1 year, and the fluctuation of soil water content in different depths was analyzed versus the rainfall and river flow to obtain the ARCK for different seasons with various interval possibilities according to our previous work (Dai et al. 2022). This will yield the ARCK that can result in the initial abstraction in the model for more accurate flood simulation in the karstic watershed.

Rainfall data

Daily and average monthly rainfall data for the three subbasins of the Lijiang River (Guilin, Yangshuo, and Pingle) during 63 years (from 1954 to 2017) were obtained from their respective meteorological stations. In order to extend the precipitation values of these three stations to their respective subbasin in the entire study area, zonal distribution (in GIS) using Radial Basis Function (RBF) was applied (Figure 3(b)). It is clear that the northern areas of the watershed have the highest amounts of precipitation and it decreases gradually to the southeast and further lower when it goes to southwestern parts of the Lijiang basin due to elevation changes (Liu et al. 2008). Moreover, the available runoff data (daily, average monthly, and maximum monthly in m3/s), as well as daily and average monthly water level (in m) were gathered from the three stations that are located at the subbasin outlets.
Figure 3

Zonal distribution map of rainfall using RBF (a) and Hyfran-Plus computed peak flood values in m3/s for different frequency design storms (% confident level) in the Guilin subbasin (b).

Figure 3

Zonal distribution map of rainfall using RBF (a) and Hyfran-Plus computed peak flood values in m3/s for different frequency design storms (% confident level) in the Guilin subbasin (b).

Close modal

Unit hydrographs

SCS and Snyder's SUH as commonly used methods were employed to estimate the runoff before being applied in the HEC-HMS model. In our study, the SUH was used to estimate the flood hydrograph, while the SCS technique was applied for comparison and control to check if the obtained data via SUH and SCS are different. Considering our watershed area scale, the Snyder Unit Hydrograph is applied which is one of the most frequently used methods. To obtain the unit hydrograph, initially, the peak instantaneous discharge value for each year was calculated from the existing daily and average monthly discharge data. Then, via having long-term yearly peak instantaneous discharge data, Hydrological Frequency Analysis Plus DSS (Hyfran-Plus) software was applied to calculate the expected maximum discharge for different return periods (2–10,000 years return periods) in all the three sub-watersheds. Hyfran-Plus is a developed mathematical process to fit statistical distributions, especially for extreme events. The DSS stands for Decision Support System that helps to choose the most appropriate type among various distributions of data. Theoretical and empirical distributions for fitting the given maximum discharge data by Hyfran-Plus were analyzed based on the null hypothesis (H0) of accuracy tests. Three distribution tests of Weibul, Gamma, and Gumbel were run to check if the adequacy and accuracy of the data are statistically correct at different significant levels via an acceptable P-value (less than 0.05 α-value for a 95% confidence level) to reject the null hypothesis (H0) or accept the alternate hypothesis (H1). Lastly, the unit hydrographs of the three subbasins were generated using both methods and modified according to the expected maximum discharge (with different return periods). The SCS equations (Equations (1)–(3)) and Snyder formula (Equation (4)) are as follows:
formula
(1)
formula
(2)
formula
(3)
formula
(4)

where CN is the empirical coefficient related to land use and soil type, P is the total rainfall, Ia and F are the values of the primitive and total loss for precipitation both in mm, S is the holding capacity (retention) of the watershed, Q is the calculated flow discharge rate in m3/s, tp is the lag time, L is the length of the main river, and Y is the slope in percentage.

The computed SUH offers certain discharge values at specific times (including the peak time). Therefore, the obtained SUH was converted to an hourly hydrograph (via regressions for each limb of the graph). And then, the dimensions of the SUH can be multiplied in any rainfall value (net amount of precipitation) to generate the runoff hydrograph of that single rain in its respective subbasin. The net amount of rain resulted from the individual rains in particular return periods, hence, the 6 h. Rainfall duration with different frequencies (50 years in our case) for the three stations was calculated.

SUH sometimes is not accurate enough in calculating the hydrological components such as the base and/or peak time of the hydrograph. Comparing the expected peak discharge for a given frequency (obtained via Hyfran-Plus) with the calculated values from SUH (and SCS) will result in a modified hydrograph to match the real observed data from the watershed. The modifications are mainly to the peak discharge and the peak time for a total constant flood volume. And lastly, the modified hydrographs can be used for calibration purposes in the HEC-HMS model. Since models are applying simplified equations and uncertainties compared to the actual conditions, hence the calculations will always require calibration and modifications through comparing the outcomes with real rainfall-runoff scenarios, as well as the calculated values (via Snyder's method in our case) for any watershed.

Muskingum flood routing

The HEC-HMS model uses the Muskingum method for flood routing. Unlike in SCS and Snyder unit hydrographs, the Muskingum in the flood routing process is taking into account the form and characteristics impacts of the watershed (slope, river path shape, topography, and geomorphological characteristics of the watershed) in prism and wedge storage. Prism is the steady part of the flow, while the wedge is the additional part. The X and K coefficients are defined for greater degrees of sophistication and reliability in the peak flood estimation. In fact, the final shape of a predicted hydrograph is a function of a complex set of criteria. The relationship between X (dimensionless) and K in the Muskingum formula is as follows Equation (5):
formula
(5)
where X represents the weighting factor to distribute the storage volume between the prism and wedge (flow-specific gravity coefficient), while K is the flow travel time (or propagation time) in hour in this method. For instance, if X= 0.3, then the wedge is 30% and the prism would be 70% of the storage total volume. As for the K, the longer it takes for the peak to move toward downstream (the greater K value), the more attenuation takes place which means the wave will have a lower peak discharge at the downstream. Comparing X and K in the Muskingum method, the latest one (K) has more impact on the hydrograph shape and peak, while higher values of X are translated to exponential magnitudes of wedges. The equations used in Muskingum flood routing are as follows (Equations (6)–(9)) (Mo et al. 2021):
formula
(6)
formula
(7)
formula
(8)
formula
(9)
where Qdo, Qdi, Quo, and Qui are denoted as the downstream and upstream outflow and inflow, respectively, t is the computation interval, and .

After the calibration of the model, six different floods that occurred in this watershed were used for validation and analysis. Real data for floods at the beginning or end of the rainfall season at three subbasins were used with respect to the seasonality of the ARCK.

The calculated results through both techniques (SCS versus Snyder) for peak discharge and the peak time of the unit hydrograph in each one of the three subbasins were adjacent (with approximately 7% difference in peak discharge rates). However, the expected peak discharge via Hyfran-Plus according to the actual observed data found to have a remarkable gap with the two methods' obtained values. For instance, the 50 years frequency rain's hydrograph was simulated with SUH which for Guilin a 134.7 mm rain did produce 67,224 m3 of total runoff volume, a peak discharge of 1266 m3/s during the base and peak times of 182 and 26 h, respectively. The expected peak discharge (via Hyfran-Plus from observed data) for the same rain was 4930 m3/s, which is approximately four times higher. For the same total runoff volume, this will result in a much shorter base time (47 h) and a peak time (7 h). In other words, we must calibrate the SUH hydrograph with four times greater peak during almost a quarter of the SUH computed base/peak times for Guilin (Figure 4(a)).
Figure 4

Direct and calibrated hydrographs via SUH for Guilin subbasin (a) and SUH computed-amended values for 50 years frequency rain in the Lijiang basin (b).

Figure 4

Direct and calibrated hydrographs via SUH for Guilin subbasin (a) and SUH computed-amended values for 50 years frequency rain in the Lijiang basin (b).

Close modal

Then, a comparison between two of the subbasins (Guilin and Yangshuo at its downstream) was performed. Although with 30% larger size compared to the Guilin subbasin, and approximately the same amount of rain, the flood hydrograph in Yangshuo has a much lower peak. Here also, a gap between the actual and simulated hydrographs was found for this subbasin (Yangshuo); however, the required modification to match the observed data was much lower than Guilin, comparatively. Indeed, a 1424 m3/s initial calculated peak flood (via SUH) was modified to 2843 m3/s expected peak (doubled), in the same way as explained above, and the 230 h calculated base time was reduced to 106 h modified base time. After adding the Guilin flood hydrograph with a lag time of 9 h to the Yangshuo hydrograph, the total calculated peak discharge rate was 7856 m3/s, which was supposed to be 6400 m3/s based on Hyfran-Plus computed data. For this total value, the difference between the expected and the modified peak is due to the fact that the expected peak resulted from Hyfran-Plus with 95% confidence level is somewhere between the two blue lines margin (Figure 3(b)) and the best possibility of estimation or the average value falls on the red line in the below graph; therefore, the primary calculated value (7856 m3/s) by SUH was accepted and no modifications were required.

Lastly, for Pingle which is the biggest subbasin at downstream, an initial obtained peak discharge rate of 2323 m3/s was modified to 6947 m3/s (based on the 8004 m3/s expected peak discharge computed by Hyfran-Plus), and the 185 h calculated base time via SUH reduced to 54 h. Yangshuo (containing Guilin flood) hydrograph was added to it with a 2 h lag time, which accounted for a total of 12,060 m3/s discharge versus 11,455 m3/s calculated via Hyfran-Plus (605 m3/s difference) with a double-peak discharge shape of the hydrograph. Figure 4(a) and 4(b) shows these differences in the Guilin subbasin graphically as well as the final Lijiang watershed hydrograph in Pingle in m3/s during the time. After the SUH graph modifications, we could use the calculated lag parts and the obtained hydrograph to run the HEC-HMS model to fit the final results.

In the HEC-HMS model, the 0.00001 tolerance, 3000 iteration, and a high Nash–Sutcliffe Efficiency or NSE value of 0.972 resulted, indicating that the model can be applied for simulation purpose in the Lijiang River. The NSE equal to 1 is considered a perfect fit, >0.75 is very good, 0.64–0.74 is good, 0.5–0.64 is satisfactory, and <0.5 is considered as unsatisfactory fit (Moriasi et al. 2007). The model-generated outlet discharge graph for a 50 years design storm with 12,223 m3/s was obtained as below in Figure 5, which represents an obvious matching with the SUH amended hydrograph as shown in Figure 4(b). The calibrated output of the model was only reliable after modifying the unit hydrograph peak discharge rate and peak/base time.
Figure 5

HEC-HMS-generated hydrograph at Pingle for the entire Lijiang watershed.

Figure 5

HEC-HMS-generated hydrograph at Pingle for the entire Lijiang watershed.

Close modal

Moreover, the average values of X and K in the Muskingum flood routing formula for different reaches of the three subbasins are tabulated and comparatively presented below in Table 1.

Table 1

Fifty years frequency flood hydrological parameters in three subbasins and geological comparison via GIS/HEC-HMS

Area (km2)Slope (%)Qmax.cal (m3/s)Qmax.exp (m3/s)Qmax.mod (m3/s)Tb (h)Tp (h)P (mm)Qtotal (m3)K (h)X
Guilin 2387.7 13.6 1266 4930 5013 47 134.7 67,200 1.56 0.24 
Yangshuo 3425.2 6.7 1424 3103 2843 106 19 130 95,500 3.78 0.38 
Pingle 5484.3 12.7 2323 8004 6947 54 12 117 123,000 1.65 0.21 
Area (km2)Slope (%)Qmax.cal (m3/s)Qmax.exp (m3/s)Qmax.mod (m3/s)Tb (h)Tp (h)P (mm)Qtotal (m3)K (h)X
Guilin 2387.7 13.6 1266 4930 5013 47 134.7 67,200 1.56 0.24 
Yangshuo 3425.2 6.7 1424 3103 2843 106 19 130 95,500 3.78 0.38 
Pingle 5484.3 12.7 2323 8004 6947 54 12 117 123,000 1.65 0.21 

Qmax.cal: calculated maximum discharge; Qmax.exp: expected maximum discharge; Qmax.mod: modified maximum discharge.

Lastly, the established model was validated using actual data. In fact, in terms of qualified rate, the simulations for six different actual flood events during 2017, 2018, 2019, and 2020 had only 52% accuracy, on average, in different subbasins prior to modification, with some unsatisfactory NSE value (<0.5) while increased to 67.5% after the inclusion of ARCK (Table 2).

Table 2

Average NSE values for the six simulated flood events in three subbasins before and after ARCK modifications

NSE values2017 May2018 Apr2018 Aug2019 Apr2019 Aug2020 AprAverage
Initial 0.54 0.74 0.41 0.67 0.33 0.41 0.52 
After ARCK 0.66 0.85 0.62 0.92 0.51 0.49 0.675 
NSE values2017 May2018 Apr2018 Aug2019 Apr2019 Aug2020 AprAverage
Initial 0.54 0.74 0.41 0.67 0.33 0.41 0.52 
After ARCK 0.66 0.85 0.62 0.92 0.51 0.49 0.675 

The huge gap between the observed and calculated peak discharges in the individual subbasins and also different required modification rates for different subbasins are probably because of the three justifications discussed here. These three are including slope, karst reservoir capacity, and river cross-section. For a 50-year frequency design storm in Guilin versus Yangshuo (134.7 and 130 mm precipitation, respectively), the total produced runoff volumes are 67,200 versus 95,500 m3, respectively, which the difference is obviously due to the watershed sizes as the Guilin subbasin area is 30% smaller than Yangshuo. However, the flood hydrograph peak in Guilin was remarkably higher than Yangshuo (5013 compared to 2843 m3/s) with a sudden onset (9 h in Guilin versus 19 h in Yangshuo), which means Guilin has 76% higher peak during almost half of the peak time in Yangshuo. Although the slope variation might cause this gap to some extent, however, the main trigger lies between the karst condition and prism to wedge storage volume. As evidence, comparing these two subbasins, the calculated average watershed slopes via GIS are 13.6 versus 7.6% for Guilin and Yangshuo, respectively. If we decrease the Guilin slope (in the SCS equation) from 13.6 to 7.6%, the calculated peak in the unit hydrograph will be degraded by 29%, only which justifies a part of the peak difference between Guilin and Yangshuo. Moreover, the impact of CN number variation was also negligible. Therefore, the other two important drivers are more likely to cause the high peak in Guilin compared to Yangshuo.

The second reason is the ARCK in Yangshuo. Geomorphologically, Devonian limestone bedrocks are the dominant parent material in the karst regions of the Lijiang catchment. Mountains here are mainly vertical-sided limestone tower karsts (peak forest or Fenglin karst type), raised individually from the alluvial plains of the Lijiang watershed. Besides, the connected conical shape hills (peak cluster or Fengcong karst type) of this basin, which stand on a common base and are connected, form a continuous terrain like an egg-box topography in the northeast of the province (Figure 6).
Figure 6

The dramatic landscape of Fenglin karst and conical towers surrounding the Lijiang River (a), the karst map of Yangshuo (Zhao et al. 2020) (b), a covered karst area in the Yangshuo subbasin being uncovered (c), the schematic process of water flow in the karstic network structure from surface flow, to stone leaching, and retention in underground lakes (d).

Figure 6

The dramatic landscape of Fenglin karst and conical towers surrounding the Lijiang River (a), the karst map of Yangshuo (Zhao et al. 2020) (b), a covered karst area in the Yangshuo subbasin being uncovered (c), the schematic process of water flow in the karstic network structure from surface flow, to stone leaching, and retention in underground lakes (d).

Close modal

The towers and conical hills are generally believed that being left from coalescence dolines due to the dissolution of limestone rocks during long evolution progress as they are fragile and vulnerable to environment changes (Jiang et al. 2019a). Fenglin karst type is mostly formed by stormwater erosion, while the Fengcong karst type has been mainly formed through the solution process due to rainfall. Unlike in Guilin, most parts of its downstream (Yangshuo) are covered by karst as shown in Figure 6(b).

From the lithology point of view, these hills and towers drain the rainwater completely to the surface/subsurface water through their steep surface, cave networks, shafts, and skylights, which work the same as a 3D drainage system (Dai et al. 2021; Mo et al. 2021). For hydrologists' concern, these conical rocky towers should transform the rainfall into runoff quickly and form a sharp hydrograph with a short peak time due to the absence of soil in these karst towers together with a high slope (Cao et al. 2018). However, in karst watersheds, besides the surface flow, the stormwater partially converts to undercurrent through the epikarst and infiltration zone owing to the developed karrens, sinkholes, foot caves, and a natural drainage network structurally (Figure 6(d)). This part (which is magnitude is greatly dependent on the ARCK and its drainage network) will flow into the saturated zone and subterranean rivers and will avoid sudden raises in the produced peak floods (Jiang et al. 2019a). This will have an immediate attenuation impact on the expected hydrograph, especially in bare karst (peak cluster) due to its developed sinkholes as compared to covered karst.

Considering the rainfall intervals, the possible scenarios of flood formation are between the margin of minimum and maximum available reservoir capacity of karst. During the wet season, a worst-case scenario might be two heavy rainfalls in a row with a short interval in which the first will fill the ARCK and consequently the second one will entirely be converted to runoff. On the other hand, if the same scenario happens during the dry season, no flood may occur as both rainfalls might just fill the ARCK and saturate the soil (Dai et al. 2022). Any other scenario also could be defined in between these two margins; therefore, the magnitude of ARCK which is a function of other parameters (such as intensity and duration of the previous rains, evaporation, temperature, humidity, and season, karstic water discharge rate which depend on the morphology and shape of the underground network) can largely control the flood features. The inclusion of this variable into the model did enhance the prediction accuracy by an average of 15.5% for six simulated floods. The influencing degree of these drainage networks and ARCK on the flood formation is highly dependent on the development of these cracks and hidden pipe systems. Technically, the development of the drainage network, sinkholes, and reservoirs is an evolutionary progress, in which the geomorphological stage of the watershed shows how developed the network might be. The karst in Yangshuo is in its old age (matured or monadnock), which indicates that the underground system and networks are most likely well developed and in a stable stage (Jia et al. 2019). Moreover, the definite magnitude of the ARCK in detaining the stormwater or the amount of inflow–outflow to and from these open basins is still undeterminable accurately (Mo et al. 2021). This was underpinned by the obtained simulated hydrographs compared with the actual data in our work, especially at Yangshuo where a modified initial abstraction resulting from ARCK could calibrate the model outputs with higher accuracy.

On the other hand, as abovementioned, the model uses the Muskingum method for flood routing. As can be seen in Table 1, subbasins with higher slopes (Guilin and Pingle) are having smaller wedges and shorter travel times (or smaller X and K, respectively), while the plateau landforms (Yangshuo watershed) have longer travel time and higher X or ‘wedge to prism storage ratio’. Hence, besides the impact of the slope and the karst reservoir, the third probable cause behind the different flood behaviors (in peak flood magnitude, peak time, and receding time) between these two sub-watersheds might be the geomorphology of the river path in these two parts of the watershed, which is considered in Muskingum with more details unlike for SCS or Snyder.

In the Guilin subbasin, lower values of X represent smaller wedges and bigger prisms (compared to Yangshuo). This means that the river network in Guilin must have a bigger cross-sectional area (in total) to accommodate higher volumes of the prism than the wedge comparatively. However, bigger values of K, which represents the wave traveling time, cause a sharp hydrograph. In other words, although the magnitude of the wedge proportionally to the prism is smaller than in Yangshuo; however, this smaller volume travels downstream within a relatively shorter time due to the higher average slope of the subwatershed. Therefore, although the total river cross-section is larger, the flood wave peak is not getting any lower or attenuated, but typically a quick flood with a progressive peak happens in Guilin City at the subbasin's outlet. The bigger or smaller cross-section size is referred to as the bathymetry-measured river cross-section size as a function of the watershed area, in which vast areas would require bigger cross-section sizes accordingly. Guilin has more tributaries and hence a bigger ‘total cross-section size to watershed area ratio’ as compared to Yangshuo at its downstream with fewer tributaries. In contrast, having a larger X value which means less prism and more wedge, the total river cross-section in Yangshuo is smaller. It also has a higher flood wave traveling time (2.4 times higher than Guilin) due to its lower slope. Hence, unlike in Guilin, the cause behind the flood in Yangshuo is the longer traveling time (due to the lower slope) of an attenuated peak. And since flat watersheds produce shallow and widespread floods, therefore, this lower volume of the wedge (compared to Guilin) will overflow to the surrounding areas of the river path, forming a shallow flood, and recedes after a relatively longer period of time.

Altogether, although the flood formation driver's contribution rates in this watershed cannot be quantitatively modeled, however, the ARCK emerged to play the most important role as a decisive factor among the three justifications behind the lower peak hydrograph in Yangshuo, since the basic ingredients found uncontroversial and it is in agreement with other findings in this catchment (Bailly-Comte et al. 2008; Zhao et al. 2019).

This work has conducted an analysis of the flood formation influencing factors in karstic areas of the Lijiang watershed. Applying the HEC-HMS model and comparing the simulated versus actual data, our results suggest that besides slope impact, a higher prism to wedge ratios in upstream causes a sharp hydrograph with a relatively short peak/base time. In contrast, the same amount of rainfall in midstream causes a more attenuated peak with a longer receding time. Irrespective of the details, this is more likely due to the available karst reservoir in this subbasin as well as the lower slope, while the latest seems to have less influence. The different capacities of the karstic reservoir during dry or wet seasons and its fluctuations within the rainy season resulted in the model calibration accuracy to be highly dependent on the initial abstract values as output had a low qualified rate of 52% accuracy prior to calibration which increased to 67.5%. We may, therefore, consider the available reservoir capacity of karst (ARCK) as the main driver in the model calibration for reliable flood forecasting in this karstic watershed.

This work was funded by the Natural Science Foundation of China (Grant Nos. 52150410408 and 51979046), the National Key Research and Development Program of China (Grant No. 2019YFC0507500), and the Science and Technology Major Project of Guangxi, China (Grant No. AA20161004).

S.R. developed methodology, wrote the review, and edited the article. J.D. conceptualized the whole article, edited and supervised the article. J.X. conducted data curation and wrote the original draft. Z.W. validated, investigated, and visualized the article. L.L. also helped Z.L. and L.P. in formal analysis, administering the project, and brought resources.

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

The authors declare there is no conflict.

Bailly-Comte
V.
,
Jourde
H.
,
Roesch
A.
,
Pistre
S.
&
Batiot- Guilhe
C.
2008
Time series analyses for karst/river interactions assessment: case of the Coulazou river (southern France)
.
Journal of Hydrology
349
(
1–2
),
98
11
.
Brirhet
H.
&
Benaabidate
L.
2016
Comparison of two hydrological models (lumped and distributed) over a pilot area of the Issen Watershed in the Souss Basin, Morocco
.
European Scientific Journal, ESJ
12
(
18
),
347
.
https://doi.org/10.19044/esj.2016.v12n18p347
.
Choudhury
P.
,
Shrivastava
R. K.
&
Narulkar
S. M.
2002
Flood routing in river networks using equivalent Muskingum inflow
.
Journal of Hydrologic Engineering
7
,
413
419
.
Dai
J.
,
Rad
S.
,
Xu
J.
,
Pen
S.
,
Gan
L.
,
Chen
X.
,
Yu
C. W.
&
Zhang
S.
2021
Impacts of climate change versus landuse change on recent Lijiang River flood regime, South China
.
Tecnologia y Ciencias del Agua
12
(
3
),
257
303
.
Delphi
M.
,
Shooshtari
M. M.
&
Zadeh
H. H.
2010
Application of diffusion wave method for flood routing in Karun River
.
International Journal of Environmental Science and Development
1
,
432
.
Elmoustafa
A. M.
&
Moussa
A. M.
2010
Evaluating of Muskingum Hydrologic Model Parameters in Reach Four of Nile River Using 1D Model
.
Scientific Bulletin Cairo University
,
Cairo, Egypt
.
Jia
Z.
,
Liu
X.
,
Cai
X.
&
Xing
L.
2019
Quantitative remote sensing analysis of the geomorphological development of the Lijiang River Basin, Southern China
.
Journal of the Indian Society of Remote Sensing
47
(
5
),
737
747
.
Jiang
G.
,
Guo
F.
,
Lo
K. F. A.
,
Liu
F.
,
Guo
Y.
,
Wu
W.
&
Polk
J. S.
2019b
Utilization status of rainwater harvesting and its improvement techniques in bare karst areas for domestic use and ecological restoration
.
Carbonates and Evaporites
34
,
1381
1390
.
Kabeja
C.
,
Li
R.
,
Guo
J.
,
Rwatangabo
D. E. R.
,
Manyifika
M.
,
Gao
Z.
,
Wang
Y.
&
Zhang
Y.
2020
The impact of reforestation induced land cover change (1990–2017) on flood peak discharge using HEC-HMS hydrological model and satellite observations: a study in two mountain basins, China
.
Water
12
(
5
),
1
22
.
 https://doi.org/10.3390/w12051347
.
Liu
Y. Q.
,
Chen
C. M.
&
Zheng
Y.
2008
The characteristics of spatial distribution and types of Apr. ∼Sept
.
Rainfall in the Pearl River Basin
24
,
67
73
.
Mo
C.
,
Wang
Y.
,
Ruan
Y.
,
Qin
J.
,
Zhang
M.
,
Sun
G.
&
Jin
J.
2021
The effect of karst system occurrence on flood peaks in small watersheds, southwest China
.
Hydrology Research
52
(
1
),
305
322
.
Moriasi
D.
,
Arnold
J.
,
Van Liew
M.
,
Bingner
R.
,
Harmel
R.
&
Veith
T.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Transactions of the ASABE
50
,
885
.
Pan
Z.
2018
Distributed NAM Hydrological Model Research and System Implementation
.
Huazhong University of Science and Technology
,
Wuhan
,
China
(in Chinese with English abstract)
.
Sharafati
A.
,
Yaseen
Z. M.
&
Pezeshki
E.
2020
Strategic assessment of dam overtopping reliability using a stochastic process approach
.
Journal of Hydrologic Engineering
25
(
7
),
04020029
.
Tewolde
M. H.
&
Smithers
J.
2006
Flood routing in ungauged catchments using Muskingum methods
.
Water SA
32
,
379
388
.
Wei
Y. L.
,
Li
C. Z.
&
Chen
W. H.
2018
Characteristics and formation and evolution analysis of the karst landscape of Guangxi
.
Guangxi Sciences
25
,
465
504
.
Xu
J.
&
Dai
J.
2018
Variation characteristics of runoff in upper reaches of Lijiang River basin
.
Yangtze River
49
(
10
),
41
46
(in Chinese with English abstract)
.
Zhao
H.
,
Xiao
Q.
,
Zhang
C.
,
Zhang
Q.
,
Wu
X.
,
Yu
S.
,
Miao
Y.
&
Wang
Q.
2020
Transformation of DIC into POC in a karst river system: evidence from δ13CDIC and δ13CPOC in Lijiang, Southwest China
.
Environmental Earth Sciences
79
(
12
),
1
12
.
Zhou
Q.
,
Chen
L.
,
Singh
V. P.
,
Zhou
J.
,
Chen
X.
&
Xiong
L.
2019
Rainfall-runoff simulation in karst dominated areas based on a coupled conceptual hydrological model
.
Journal of Hydrology
573
,
524
533
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).