Lijiang flood characteristics and implication of karst storage through Muskingum flood routing via HEC-HMS, S. China

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).

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 km² 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 km² (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 km² (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 km² 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 m³/s in the dry season and up to 12,000 m³/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).

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.

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.

where CN is the empirical coefficient related to land use and soil type, P is the total rainfall, I_a 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 m³/s, t_p 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.

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.

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 m³/s initial calculated peak flood (via SUH) was modified to 2843 m³/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 m³/s, which was supposed to be 6400 m³/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 m³/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 m³/s was modified to 6947 m³/s (based on the 8004 m³/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 m³/s discharge versus 11,455 m³/s calculated via Hyfran-Plus (605 m³/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 m³/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.

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.

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).

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 m³, 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 m³/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 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.