Predicting floods in a large karst river basin by coupling PERSIANN-CCS QPEs with a physically based distributed hydrological model

In general, there are no long-term meteorological or hydrological data available for karst river basins. The lack of rainfall data is a great challenge that hinders the development of hydrological models. Quantitative precipitation estimates (QPEs) based on weather satellites offer a potential method by which rainfall data in karst areas could be obtained. Furthermore, coupling QPEs with a distributed hydrological model has the potential to improve the precision of flood predictions in large karst watersheds. Estimating precipitation from remotely sensed information using an artificial neural network-cloud classification system (PERSIANN-CCS) is a type of QPE technology based on satellites that has achieved broad research results worldwide. However, only a few studies on PERSIANN-CCS QPEs have occurred in large karst basins, and the accuracy is generally poor in terms of practical applications. This paper studied the feasibility of coupling a fully physically based distributed hydrological model, i.e., the Liuxihe model, with PERSIANN-CCS QPEs for predicting floods in a large river basin, i.e., the Liujiang karst river basin, which has a watershed area of 58 270 km², in southern China. The model structure and function require further refinement to suit the karst basins. For instance, the sub-basins in this paper are divided into many karst hydrology response units (KHRUs) to ensure that the model structure is adequately refined for karst areas. In addition, the convergence of the underground runoff calculation method within the original Liuxihe model is changed to suit the karst water-bearing media, and the Muskingum routing method is used in the model to calculate the underground runoff in this study. Additionally, the epikarst zone, as a distinctive structure of the KHRU, is carefully considered in the model. The result of the QPEs shows that compared with the observed precipitation measured by a rain gauge, the distribution of precipitation predicted by the PERSIANN-CCS QPEs was very similar. However, the quantity of precipitation predicted by the PERSIANN-CCS QPEs was smaller. A post-processing method is proposed to revise the products of the PERSIANN-CCS QPEs. The karst flood simulation results show that coupling the post-processed PERSIANN-CCS QPEs with the Liuxihe model has a better performance relative to the result based on the initial PERSIANN-CCS QPEs. Moreover, the performance of the coupled model largely improves with parameter re-optimization via the post-processed PERSIANN-CCS QPEs. The average values of the six evaluation indices change as follows: the Nash–Sutcliffe coefficient increases by 14 %, the correlation coefficient increases by 15 %, the process relative error decreases by 8 %, the peak flow relative error decreases by 18 %, the water balance coefficient increases by 8 %, and the peak flow time error displays a 5 h decrease. Among these parameters, the peak flow relative error shows the greatest improvement; thus, these parameters are of the greatest concern for flood prediction. The rational flood simulation results from the coupled model provide a great practical application prospect for flood prediction in large karst river basins.

1 Introduction

The highly anisotropic karst water-bearing media and intricate hydraulic conditions cause karst flood processes to exhibit significant differences in time and space, which leads to laminar flow and turbulent flow transmutation in karst areas; thus, flood events in karst river basins are more complicated than those in non-karst areas (Ford and Williams, 1989; Goldscheider and Drew, 2007). This difference makes it difficult to precisely simulate and forecast the karst flood process using a hydrological model. It is common practice to simplify karst water-bearing media before building a model. For example, the karst river basin could be made into a multiple and nested spatial structure, the underground river could be made into an intelligible river system in the model, and the cave could be an anisotropic medium with a large vertical infiltration coefficient and porosity but a small specific yield. Even so, it is still hard to quantify the spatial structure of karst water-bearing media with a physics–mathematics model. Karst flood simulation results usually have some errors that cannot be ignored, and these errors represent the main problem in flood prediction in karst river basins (Kovacs and Perrochet, 2011).

Because the dynamic changes in karst hydrological processes and the hydraulic conditions of the underlying surface are complicated and nonlinear in karst areas, obtaining hydrogeological parameters, such as specific yield, hydraulic conductivity and aquifer transmissivity, is difficult. With the rapid development of remote sensing, GIS technology and hydrogeology, the technology used in field work, including tracer tests (Birk et al., 2005; Doummar et al., 2012) and infiltration tests, has made significant progress. However, accurately simulating the laws of motion of karst hydrological processes in karst water-bearing media based on these experimental tests remains difficult. Therefore, traditional methods, such as lumped hydrological models, are not suitable for flood prediction in karst areas (Hartmann et al., 2013). Compared with the performance of lumped hydrological models, physically based distributed hydrological models (PBDHMs) have some advantages in terms of generating karst flood predictions. PBDHMs divide the entire karst river basin into a series of small grid units named karst sub-streams, which precisely reflect the real rules of hydrological processes and karst development characteristics. Therefore, the PBDHM approach has great application potential in terms of improving karst flood simulation and prediction capabilities (Ambroise et al., 1996). Many PBDHMs have been proposed since the blueprint of the PBDHM was published by Freeze and Harlan (1969). The first full PBDHM, called the SHE model, was published in 1987 (Abbott et al., 1986a, b). Shustert and White (1971) attempted to use the PBDHM in karst areas. In their research, the dissolved carbonate species were analyzed in the waters of 14 carbonate springs in the central Appalachians. These springs were classified into diffuse-flow feeder-system types and conduit feeder-system types. PBDHMs have obtained several good research results in karst areas (Atkinson, 1977; Quinlan and Ewers, 1985; Quinlan et al., 1996; Duan and Miller, 1997; Ren, 2006; Liu et al., 2013; Zhang et al., 2007).

The PBDHM used in this paper is the Liuxihe model (Chen, 2009), which is a fully distributed model with 14 physically based parameters. After the karst mechanisms were added, the number of parameters was 20. Unlike other distributed hydrological models, there are some special structural designs in this model. For instance, the whole model structure is divided into eight independent parts, which are called sub-models. These sub-models include the (1) watershed division and data mining sub-model, (2) unit classification and river section estimation sub-model, (3) rainfall fusion computational sub-model, (4) evapotranspiration calculation sub-model, (5) runoff calculation sub-model, (6) confluence calculation sub-model, (7) parametric sensitivity analysis sub-model, and (8) parameter optimization sub-model. Unlike other distributed models, separate parameter uncertainty analysis calculations must be performed outside the model. However, the parametric sensitivity analysis is a fixed module in the Liuxihe model, which means that when the model is built for flood prediction, parametric uncertainty analysis has already been carried out. The parametric uncertainty analysis in the Liuxihe model is based on a multi-parameter sensitivity analysis that was presented by Choi et al. (1999).

In actual flood predictions, people may pay more attention to the flood process at specific points of the river section. For example, focus may be directed at the mouth of the river or the outlet of the basin. These points have special significance in relation to procedures such as flood warnings and evacuations. Therefore, extracting the flood processes at these points is important and should be given special consideration. In the Liuxihe model, these points are named early warning points, and flood prediction, which is urgently needed in karst areas, can be performed separately at these points. For example, the confluence of underground rivers could be established through a field survey and a geological borehole test and set as an early warning point because this is a point at which the influence of karst may dominate the runoff processes.

In addition, the catchment property data for the Liuxihe model, which primarily include the digital elevation model (DEM), land use and soil types, can be easily downloaded from open-access databases for free. Therefore, the Liuxihe model can be built in any area. Though it is not easy to obtain the basic data needed to build a distributed hydrological model in karst areas, only a very small amount of data must be downloaded from the web to build the Liuxihe model, making it a feasible option for flood simulation and prediction in karst basins.

The regulation and storage capacity of karst water-bearing media are weak. When the accumulated rainfall exceeds the maximum drainage capacity of the channel during a heavy rain storm, a karst immersion-waterlogging hazard is much more likely to occur. The hazard will become increasingly serious with the intensification of extreme global weather events. Therefore, some effective measures need to be taken to reduce losses caused by floods. For example, effectively and reliably simulating and predicting the karst flood process using a PBDHM are important non-project measures for flood control. However, there are insufficient rain gauges and long-term meteorological or hydrogeological data available to build a PBDHM in karst river basins classified as ungauged basins. Predictions in ungauged basins (PUBs) are the theme of the international hydrological decade, at the core of which is runoff calculation (Li and Ren, 2009). Therefore, it is more difficult to forecast flood events in karst river basins than in non-karst areas. How to solve the problem of rainfall sources is a key factor in the current karst flood prediction challenge. Quantitative precipitation estimates (QPEs) and, particularly, satellite QPE technology, make it possible to obtain reasonable rainfall data in karst areas. However, the current application of QPEs is immature, which results in poor QPE accuracy, and the effect of the karst flood simulation and prediction is also poor.

The development of numerical weather prediction models in recent decades has provided a reasonable and accurate QPE product that can be used in karst areas. The current mainstream QPEs include weather radar QPEs (Delrieu et al., 2014; Rafieei et al., 2014; Faure et al., 2015), satellite QPEs and radar-merging satellite QPEs (Stenz, 2014; Bartsotas et al., 2017; Goudenhoofdt and Delobbe, 2009; Wardhana et al., 2017). Additionally, precipitation can be estimated from remotely sensed information using artificial neural networks/PERSIANN QPEs (Soroosh et al., 2000; Hirpa et al., 2010; Romilly, 2011; Yang et al., 2007), the dPERSIANN-climate data record/PERSIANN-CDR (Ashouri et al., 2014; Liu et al., 2017; Tan and Santo, 2018; Hussain et al., 2018), and the PERSIANN-cloud classification system/PERSIANN-CCS (Yang et al., 2004, 2007; Mekonnen and Hossain, 2010). Studying the QPE products from meteorological satellites has become a popular topic in rainfall prediction research (Hu et al., 2013).

Many scholars at home and abroad have performed considerable research using QPE technology, and they have achieved many acceptable results. However, considerable uncertainty exists in the application of these results, which causes the precision of the QPEs to be low; thus, the precipitation result generated from the QPEs may be unsatisfactory. Two effective measures could reduce the uncertainty of QPE results in the karst area. One measure is to match the appropriate resolution of the model. The resolution can directly affect the results of the QPEs; thus, if the resolution is too low, then the division of the grid units is too coarse, which causes a considerable error in the rainfall estimates. However, if the resolution is too high, then the meteorological model structure is complicated and unstable. Furthermore, the required computational resources will increase exponentially as the model spatial resolution increases (Chen et al., 2017), which leads to a large number of calculations and low efficiency. Therefore, using the appropriate model spatial resolution is extremely important in terms of the QPE results. The other measure that affects uncertainty is that the current technology of QPEs still has some systematic errors due to uncertainties in the structure and mathematical algorithms. For this reason, when compared with the precipitation observed using rain gauges, the results of QPEs have some relative errors, and these errors cause the karst flood simulation results from the coupled model (i.e., those from coupling the QPEs with a PBDHM) to have uncertainties that largely affect the model's performance. Therefore, the results of the initial QPEs could not be directly used to build the coupled model. In this study, a post-processing method was employed to revise the productions of the PERSIANN-CCS QPEs products, which caused the QPE results to be more credible and receivable.

There have been many studies of PERSIANN-CCS QPEs (Yang et al., 2007). However, most of these studies have been conducted in small non-karst watersheds. In this study, the PERSIANN-CCS QPEs were employed in an attempt to estimate the rainfall data in a large karst river basin, i.e., the Liujiang karst river basin (LKRB), which has an area of 5.8×10⁴ km² and is located in Guangxi Province, China. Watershed flood prediction relies on a PBDHM as a computation tool, while precipitation is the driving force behind the model (Li et al., 2017). This method has the potential to improve the accuracy of karst flood predictions by coupling PERSIANN-CCS QPEs with a PBDHM. The PBDHM in this study is the Liuxihe model (Chen, 2009). This study is the first time the Liuxihe model for flood simulation and prediction in karst basins has been used.

Therefore, the model structure and function have been improved to suit the requirements of the karst basin. For instance, in this study, the entire river basin will be divided into many small sub-basins using the DEM data, and this process is adequate when considering non-karst basins. However, to ensure the effect and accuracy of the model in karst areas, the model structure must be more refined. Thus, in this paper, the sub-basins will be further divided into many karst hydrology response units (KHRUs). The entire karst hydrological process, including the storage and regulation processes of the epikarst zone and the spatial interpolation of the precipitation, evapotranspiration and rainfall–runoff, are all calculated based on the KHRUs. Furthermore, in the original Liuxihe model, the underground layer is treated as an integral unit, and a linear reservoir method is adopted to calculate the amount of underground runoff. However, because the structure of the karst underground layer is nonlinear, the original linear reservoir method of the Liuxihe model is not appropriate. Therefore, in this study, the Muskingum routing method is used to improve the convergence of the underground runoff calculations. Additionally, the epikarst zone, as a distinctive structure of the KHRU, is carefully considered in the model. An exponential decay equation is used to calculate the regulation and storage processes in the epikarst zone.

The spatial resolution of the Liuxihe model for the LKRB is 200 m ×200 m. The PERSIANN-CCS QPE products, which have a spatial resolution of 0.04^∘ ×0.04^∘ and a time interval of 30 min, are employed to estimate the precipitation results for the LKRB. The resolution of the PERSIANN-CCS QPEs must be downscaled to the same size as the Liuxihe model before the coupled model can be built. After post-processing, the PERSIANN-CCS QPE products could offer high-precision precipitation results for the LKRB in locations where there is an inadequate number of rain gauges. Additionally, the model performance can be greatly improved by coupling the post-processed PERSIANN-CCS QPEs with the Liuxihe model. A modified PSO algorithm (Chen et al., 2016) is used to optimize the coupled model parameters in this paper, and this method could control the uncertainty of parameterization.

2 Study area and data

2.1 Study area

The LKRB in southern China was selected as the study area for this research. The LKRB is the second largest tributary of the Pearl River and covers three provinces, namely Guizhou, Guangxi and Hunan. The LKRB is the most developed karst area of China, with a drainage area of 58 270 km² and a channel length of 1121 km. Moreover, the LKRB is a typical karst-mountainous catchment that has experienced frequent flash flooding in past centuries. The peak forest-plain area is the main karst landform on the ground, while the karst conduit and fissure are well developed underground. There are also many complicated underground rivers and springs with large flows (Li, 1996). The karst water-bearing media are highly nonlinear and heterogeneous, which makes it very difficult to simulate and forecast the karst hydrological process.

The LKRB is in the sub-tropical monsoon climate zone, with an average annual precipitation between 1400 and 1700 mm, and the precipitation distribution is highly uneven on spatial and temporal scales. The precipitation from April to September accounts for 75 %–80 % of the annual precipitation. A sketch map of the LKRB is shown in Fig. 1a.

Figure 1 Sketch map of Liujiang and the Beijiang catchment.

The most developed karst area in LKRB is the Beijiang catchment, where the influence of karst features highly dominates the rainfall–runoff processes. The Beijiang catchment is a tributary of the middle and upper reaches of the Liujiang River, lying at 25^∘06^′–25^∘27^′ N and 108^∘38^′–109^∘18^′ E. The drainage area of the Beijiang catchment is 1790 km², and the length is 130 km. The catchment has a dense river system (Fig. 1b) and is surrounded by high mountains with peak elevation at 1000–1800 m (Fig. 1c), in which the peak-cluster depression covers most of the area. The average valley slope gradient is 0.143.

2.2 Landform, tectonics and hydrogeology information

The LKRB is located in the central part of Guangxi Province, China. The terrain is high on all sides and low in the middle. The cross-strait terraces of the Liujiang River are well developed, especially near the Liuzhou River gauge (as shown in Fig. 1), which is located at the outlet of the LKRB. The northern part of the basin has transmeridional arc-like folded belts, where the soluble rock forms syncline and the sand shale forms anticline. Sand shale formations and carbonate and carbonate clastic rocks are widely distributed here. The karst valley is the main landform in the southern part of the basin, and the overlying lithology is clay and gravel with poor water permeability. The underlying bedrock is mainly carbonate and dolomite, and the karst fissures are well developed, in which a large amount of water is stored (He, 2017).

The western part of the basin has a large area of limestone in a continuous distribution, and a peak-cluster depression covers most of the area. The landform of the eastern basin is mainly hilly, where the rocks are soft-hard due to their different anti-erosion abilities. The hard rocks form low mountains that move towards the gentle slope and then back to the steep slope. The landforms of the central part of the basin are mainly the isolated peak plain and the peak forest plain. Overall, the main landforms of the LKRB are the peak forest plain and the peak-cluster depression.

The Liujiang River is located in the karst valley basin, which is covered by Quaternary loose deposits. The underlying surface is dominated by alluvium, diluvium and katatectic layers due to the fluviraption of the Liujiang River and the karst geological background, and the thickness is approximately 10–20 m. Carbonate, sandstone, shale and carbonate clastic rocks are widely distributed in the basin. Among them, the area of the carbonate rocks is about 19 230 km², which accounts for 33 % of the entire watershed. The outcrops in the basin mainly include Upper Devonian limestone (D₃), Lower Carboniferous Datangpo formation limestone (C₁d, C₁d³), Middle (C₂d) and Upper Carboniferous (C₃) limestone, Upper Permian carbonate and clastic rocks (P₂d, P_2 h), Lower Triassic clastic and carbonate rocks (T₁), Lower Cretaceous clastic and carbonate rocks, and loose rock groups of the Quaternary Pleistocene (Q, Q_p) and Holocene (Q_h).

After studying the karst geomorphology of the LKRB, Williams (1987) believed that the peak-cluster depression had developed into turreted peak-forest landforms after a long evolutionary process, which is equivalent to the late prime of life, i.e., entering old age in terms of geomorphologic evolution. Allogeneic water, especially from the Liujiang River, is the main driving force behind the development of peak-forest landforms. Therefore, the peak-forest plains and valleys are often distributed in contiguous areas near the main trunk stream of the Liujiang River. The main karst landform of the LKRB is peak-forest plain, and there are also some peak-cluster depressions and peak-forest valleys. Figure 2 shows the DEM and three-dimensional topographical map of the LKRB.

Figure 2 The DEM and three-dimensional topographical map of the LKRB.

2.3 Rain gauges and karst flood process

There are 68 rain gauges and 131 grid points for the PERSIANN-CCS QPEs within the LKRB, and data from 30 karst flood events that occurred between 1982 and 2013 were collected. There was one flood event each year. Among them, five karst flood events between 2008 and 2013 were used to test the effect of coupling PERSIANN-CCS QPEs with the Liuxihe model. The karst flood process in the LKRB has typical characteristics: the flood peak flows usually exceed 10 000 m³ s⁻¹, and there is an expression of a multi-peak flood process. A flood process usually lasts approximately 10 days, and the shortest flood event duration was only approximately 3 days, while the longest was 25 days. Hourly precipitation data were collected from the rain gauges in this study, and these results were compared with the results from the PERSIANN-CCS QPEs. The rain gauges, the grid points of the PERSIANN-CCS QPEs and the Liuzhou River gauge that is located close to the outlet of the LKRB are shown in Fig. 1a.

There are 11 early warning points set in the Beijiang catchment (Fig. 1b), and 10 karst flood events at the Goutan warning point were collected to validate the flood simulation effect based on the Liuxihe model, in which the Goutan point is the outlet of the Beijiang catchment. In fact, the Beijiang catchment is in the center of the storm area of Guangxi Province, China. According to field observation data, the observed maximum 24 h accumulated precipitation is 779.11 mm in the Beijiang catchment, and the maximum 3-day accumulated precipitation is 1335.15 mm. Karst floods are typical flash floods with rapid discharge and water level fluctuation, mainly caused by storms, and the developed karst landform plays an important role in flood propagation. For instance, the karst depressions can store some water content during heavy rain. Additionally, the regulation functions of the karst fissure system can slow the flood propagation process.

2.4 Property data

The catchment property data for the distributed hydrological models mainly include the DEM, land use and soil types. These data were downloaded from open-access databases. The DEM was downloaded from the Shuttle Radar Topography Mission database at http://srtm.csi.cgiar.org (last access: 12 June 2018) (Falorni et al., 2005; Sharma et al., 2014). The downloaded DEM had an initial spatial resolution of 90 m ×90 m, and after many model resolution tests, the most appropriate resolution of the Liuxihe model in the LKRB was confirmed to be 200 m ×200 m. Therefore, the spatial resolution of the initial DEM was rescaled to 200 m ×200 m in this study, and this value represents the high resolution for the Liuxihe model in the LKRB. The DEM is shown in Fig. 2a. The land use-type data were downloaded from http://landcover.usgs.gov (last access: 12 June 2018) (Loveland et al., 1991, 2000), and the soil-type data were downloaded from http://www.isric.org (last access: 12 June 2018). The initial spatial resolutions of the land use-type and soil-type data were both 1000 m ×1000 m. However, both resolutions had to be rescaled to 200 m ×200 m in this study. Figure 3a shows the land use types, and b shows the soil types.

Figure 3 The property data for the Liuxihe model in the LKRB.

3 PERSIANN-CCS QPEs and post-processing results

3.1 PERSIANN-CCS QPEs

The original PERSIANN system (Hsu et al., 1999) was based on geostationary infrared imagery and was later extended to include the use of both infrared and daytime visible imagery. This method represents an automated system for estimating precipitation from remotely sensed information through the use of artificial neural networks. The method for rainfall estimation that is under development at the University of Arizona is continuously improving as technology advances (Soroosh et al., 2000). The fundamental algorithm of the PERSIANN system is based on a neural network. The network parameters could be optimized by an adaptive training characteristic, which can estimate the precipitation from a geosynchronous satellite at any time and place.

The PERSIANN-CCS (Yang et al., 2004; Hsu et al., 2007) is a patch-based cloud classification and rainfall estimation system from low Earth orbit and geostationary satellites that uses pattern recognition technology and computer imaging technology (Yang et al., 2007). Satellite-based precipitation retrieval algorithms use information ranging from visible (VIS) to infrared (IR) spectral bands of geostationary earth orbiting (GEO) satellites and microwave (MW) spectral bands (Hsu et al., 2007).

The QPE products of PERSIANN-CCS have generated precipitation estimates at a resolution of 0.04^∘ ×0.04^∘ scale and at a time interval of 30 min since 2000. The output of PERSIANN-CCS QPEs was downscaled at 200 m × 200 m to achieve the same spatial resolution as that of the Liuxihe model in the LKRB. The down-scaling method used in this paper was based on statistical relationships between the meteorological variables and DEM data using the LOO (leave-one-out) cross evaluation method and spatial autocorrelation analysis methods (Fan et al., 2017).

The hourly precipitation data from the PERSIANN-CCS QPEs were collected and compared with the precipitation observed by the rain gauges.

The estimation of rainfall from the PERSIANN-CCS consists of the following steps (Hsu, 2007): (1) IR cloud image segmentation, (2) characteristic extraction from IR cloud patches, (3) patch characteristic classification, (4) obtaining the rainfall estimation results of the QPE products, and (5) evaluating and revising the results of the QPE products.

In this paper, the PERSIANN-CCS QPEs real-time data used in the LKRB from the current version of PERSIANN-CCS are available and downloadable online (http://cics.umd.edu/ipwg/us_web.html, last access: 28 June 2018).

3.2 Precipitation estimation results

The QPE product of the PERSIANN-CCS generated precipitation results for the LKRB. There were 131 grid points of PERSIANN-CCS QPEs within the LKRB, and these points were representative and completely covered the entire watershed (as shown in Fig. 1). The spatial resolution was 200 m ×200 m, and the time interval was 1 h. The respective QPE products of the PERSIANN-CCS in 2008, 2009, 2011, 2012 and 2013 were produced, and the results indicated that five rainfall events corresponded to the five karst flood processes. Figures 4–8 show the average precipitation pattern comparisons of the two precipitation products of the 5 years, where (a) is the average precipitation based on data from the rain gauges, (b) is the average precipitation based on the data from the PERSIANN-CCS QPEs, and (c) is the quantile–quantile plot, in which the 45^∘ line is used to compare two precipitation products.

Figure 4 Precipitation pattern comparison of two precipitation products (2008): (a) is the average precipitation of rain gauges, (b) is the average precipitation of PERSIANN-CCS QPEs, and (c) is the quantile–quantile plot, in which the 45^∘ line is used to compare the two precipitation products.

Figure 5 Precipitation pattern comparison of two precipitation products (2009): (a) is the average precipitation of rain gauges, (b) is the average precipitation of PERSIANN-CCS QPEs, and (c) is the quantile–quantile plot, in which the 45^∘ line is used to compare the two precipitation products.

Figure 6 Precipitation pattern comparison of two precipitation products (2011): (a) is the average precipitation of rain gauges, (b) is the average precipitation of PERSIANN-CCS QPEs, and (c) is the quantile–quantile plot, in which the 45^∘ line is used to compare the two precipitation products.

Figure 7 Precipitation pattern comparison of two precipitation products (2012): (a) is the average precipitation of rain gauges, (b) is the average precipitation of PERSIANN-CCS QPEs, and (c) is the quantile–quantile plot, in which the 45^∘ line is used to compare the two precipitation products.

Figure 8 Precipitation pattern comparison of two precipitation products (2013): (a) is the average precipitation of rain gauges, (b) is the average precipitation of PERSIANN-CCS QPEs, and (c) is the quantile–quantile plot, in which the 45^∘ line is used to compare the two precipitation products.

According to the results of Figs. 4–8, it appears that the temporal average precipitation patterns of both products are quite similar, especially in terms of the rainfall distribution, while there are some differences in the quantitative values. The results from the PERSIANN-CCS QPEs are smaller than those from the rain gauges, which means that a relative error exists between the two products. From the quantile–quantile plot, the two rainfall scatter plots are closely distributed on both sides of the 45^∘ line, which means that the rainfall distributions of both products are close to each other.

3.3 Evaluation of PERSIANN-CCS QPEs

To quantitatively evaluate the results of the PERSIANN-CCS QPEs, the precipitation from the PERSIANN-CCS QPEs and the precipitation from the rain gauges were compared in this study. The rainfall distribution of both products is shown in Figs. 4–8. For further comparison, the average precipitation of the five karst flood events was calculated, and the results are shown in Table 1.

Table 1 Precipitation pattern comparison of the two precipitation products.

According to the results of Table 1, there are obvious relative errors between the two precipitation products. The average precipitation values of the PERSIANN-CCS QPEs were lower than those from the rain gauges. For the five karst flood events from 2008 to 2013, the relative errors between the two products were −11 %, −16 %, −7 %, −19 % and −20 %, respectively. The average relative error was −14 %, and the maximum error was −20 %, which means that these relative errors cannot be ignored. Therefore, the precipitation results generated by the PERSIANN QPEs must be revised effectively, and the precipitation data observed by the rain gauges can be used to revise the results of the PERSIANN QPEs in this study.

3.4 The post-processed PERSIANN-CCS QPEs

To make the results of the PERSIANN QPEs more credible and receivable, the precipitation results were revised using the observed precipitation measured by the rain gauges. First, it was necessary to locate the grid points of the PERSIANN-CCS QPEs that were closest to the rain gauges (as shown in Fig. 1). There were 23 grid points in the LKRB. Second, the average precipitation values of the PERSIANN-CCS QPEs and the rain gauges were calculated, and the average precipitation from the rain gauges was used as the true precipitation value. Third, the process of revising the results of the PERSIANN QPEs based on the average precipitation observed by the rain gauges is summarized as follows.

1. The average precipitation of these 23 grid points based on the PERSIANN-CCS QPEs was calculated with the following equation:

   $$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{P}}_{\mathrm{PERSIANN}-\mathrm{CCS}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N{P}_i{F}_i}{N}},\end{array}$$

   where ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{PERSIANN}-\mathrm{CCS}}$ is the average precipitation of the 23 grid points based on the PERSIANN-CCS QPEs, P_i is the precipitation based on the PERSIANN-CCS QPEs at the i grid point, F_i is the catchment area of the i grid point, and N is the number of grid points.

2. The average precipitation of the 23 rain gauges was calculated using the following equation:

   $$\begin{array}{}\text{(2)}& {\stackrel{\mathrm{‾}}{P}}_{\mathrm{2}}={\displaystyle \frac{\sum _{j=\mathrm{1}}^M{P}_j}{M}},\end{array}$$

   where ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{2}}$ is the average precipitation observed by the 23 rain gauges, P_j is the precipitation observed at the j rain gauge, and M is the number of rain gauges.

3. The precipitation values observed by the adjacent rain gauges were used to revise the results of the PERSIANN-CCS QPEs with the following equation:

   $$\begin{array}{}\text{(3)}& P_i^{\prime }=P_i{\displaystyle \frac{{\stackrel{\mathrm{‾}}{P}}_{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{P}}_{\mathrm{PERSIANN}-\mathrm{CCS}}}},\end{array}$$

   where $P_i^{\prime }$ is the value of precipitation based on the PERSIANN-CCS QPEs after revision on the i grid point, and ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{2}}/{\stackrel{\mathrm{‾}}{P}}_{\mathrm{PERSIANN}-\mathrm{CCS}}$ is the revised factor.

4. After revision, the precipitation results based on the PERSIANN-CCS QPEs were used as input data for the Liuxihe model to test its feasibility for use in the flood simulation.

After running the post-processing procedure for the PERSIANN-CCS QPEs described above, it was determined that the revised factor ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{2}}/{\stackrel{\mathrm{‾}}{P}}_{\mathrm{PERSIANN}-\mathrm{CCS}}$ was a key factor that made the results of the PERSIANN-CCS QPEs much closer to the value of observed precipitation recorded by the rain gauges, indicating that the systematic errors of the PERSIANN-CCS QPEs could be corrected effectively. Therefore, the post-processing method described in this paper is both feasible and necessary. Additionally, it could greatly improve the accuracy of the coupled model in the simulation and prediction of karst floods. Furthermore, the revised factor could be preserved as an empirical value for future flood prediction in the LKRB.

4 Hydrological model

4.1 Liuxihe model

The Liuxihe model proposed by Yangbo Chen (Chen, 2009) of Sun Yat-Sen University, China, is employed as the fully distributed hydrological model in this study, which is a physically based distributed hydrological model (PBDHM) mainly for catchment flood simulation and prediction (Chen et al., 2016, 2017; Li et al., 2017). The Liuxihe model earned its name by being the first successful application in the Liuxihe catchment, Guangdong Province, China. There are three layers vertically, including the canopy layer, the soil layer and the underground layer in the model, and the whole catchment is divided into a great number of grid cells horizontally using the high-resolution DEM data, with the divisions called sub-basins. Each grid is considered a uniform basin, and the elevation, land cover type, soil type, and other model elements including rainfall–runoff, evapotranspiration, etc. are calculated in the uniform basin. All cells are categorized into three types, namely hillslope cell, river cell and reservoir cell.

An improved PSO algorithm (Chen et al., 2016) is employed to optimize the model parameters in this study, which can make the model's performance much better in flood prediction in karst river basins. The observed meteorological and hydrological data and the development conditions of the karst underground river are used to optimize the model parameters. The terrain property data, such as the DEM, land use type and soil type, can be downloaded freely from an open-access database online. The model is validated against observed karst flood events. These factors of the model are physically based and rational to truly reflect the underlying surface of the karst basin. Therefore, this implies that the Liuxihe model could be used for real-time flood prediction in karst river basins. Figure 9 shows the structure of the Liuxihe model.

Figure 9 The structure of the Liuxihe model.

4.2 Improvement of the Liuxihe model

The Liuxihe model has been successfully applied for flood predictions in many river basins. However, none of these basins were karst areas. This study is the first time the model has been used in a karst river basin. The structure of the model should be improved to suit the needs of the karst basin in question. Therefore, some effective measures should be taken before building the model. First, the karst water-bearing media should be simplified, and this process could include making the karst basin a multiple and nested spatial structure. The underground river could be included as the intelligible channel system in the model, and the cave could be used as the anisotropic medium with a large vertical infiltration coefficient and porosity but a small specific yield. Finally, the fault could be used as the anisotropic medium with a large vertical infiltration coefficient and a specific yield. Second, the entire karst river basin can be divided into many small karst sub-basins using high-resolution DEM data. Furthermore, to suit the karst area, the karst sub-basins can be divided into many KHRUs, which are generally independent of each other. The entire karst hydrological process, including the storage and regulation processes of the epikarst zone, the spatial interpolation of precipitation, the evapotranspiration and the rainfall–runoff, are all calculated based on this KHRU. Then, these hydrological processes can be summarized for each of the karst sub-basins. Additionally, the outlet flow is formed through the river confluence among each karst sub-basin from the upstream region to the downstream region. This type of multi-structure distributed hydrological model could utilize variously scaled information effectively and optimize the use of observed meteorological, hydrological and geological data.

In this study, the KHRUs were divided by GIS technology combined with karst topography, land use type and soil type (Ren, 2006). Each KHRU in this study had its own model characteristics, such as meteorological and hydrological characteristics, as well as the karst developmental characteristics. The KHRU was proposed to describe the spatial variation of the karst sub-basins. The differences within the KHRUs were smaller than those among the KHRUs. Then, each KHRU was vertically divided into five layers: the canopy, the soil, the epikarst zone, the bedrock and the underground river. A sketch map of the KHRU is shown in Fig. 10.

In Fig. 10b, the three-dimensional model of the KHRU in the LKRB was built in the laboratory to better understand how groundwater moves in the karst media and converges with the surface river. Then, the hydrological model could be built and visualized in this way.

To satisfy the applicability of the model in karst areas, the epikarst zone, which is a distinctive structure of the KHRU, was carefully considered in the model. The epikarst zone is composed of karst rocks with macro cracks and tiny fissures. When rain falls on the ground, it is intercepted by plants, held in depressions and experiences some evapotranspiration. Then, the rainfall infiltrates into the soil and rock layer and satisfies the water shortage of the unsaturated zone. Part of the water in the epikarst zone may form karst springs that emerge from the surface. Another part will enter the superficial karst water system of the epikarst zone. When the rainfall intensity is heavy enough to form surface runoff on the exposed bedrock, part of the water will enter the karst conduit through sinkholes.

Figure 10 Sketch map of the KHRU.

The karst hydrological process of the epikarst zone could be divided into rapid fissure flow and slow fissure flow. After heavy rain, a large amount of water in the epikarst zone is stagnant and can form a surface karst aquifer with a temporary water table. If there are large cracks or fractures under the water table, a precipitation funnel will form and be associated with a drop in the water table. Rapid fissure flow refers to rainfall that infiltrates into the karst conduit through the precipitation funnel, and this flow occurs in the macro cracks and has high speeds. When rainfall enters the superficial karst water system of the epikarst zone, the macro cracks will fill first. This part of the saturated water content, named rapid fissure flow, will move directly into the karst conduit through the macro crack. Because this rapid fissure flow will pass quickly through the karst conduit system without stopping, and because the water regulation and storage functions are weak, the regulation and storage of the rapid fissure flow were ignored in this study. The rest of the water content in the epikarst zone infiltrates through tiny fissures. This part of the water, named slow fissure flow, plays an important role in the process of rainfall regulation. The water content of the slow fissure flow can be described by the following equation:

$$\begin{array}{}\text{(4)}& {\mathrm{SW}}_{\mathrm{epi}}=Q_{\inf }-V_{\mathrm{crk}},\end{array}$$

where SW_epi is the water content of the slow fissure flow in the epikarst zone.

Q_inf is the infiltration water content of the rainfall, and V_crk is the water content of the rapid fissure flow in the macro crack.

The slow fissure flow in the epikarst zone is calculated by an exponential decay equation (Ren, 2006) as follows:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{W}_{sep}=W_{\mathrm{epi}}\left(\mathrm{1}-\exp \left({\displaystyle \frac{-\mathrm{\Delta }T}{{\mathrm{TT}}_{perc}}}\right)\right),\\ W_{\mathrm{epi},t+\mathrm{1}}=W_{\mathrm{epi},t}+{\mathrm{SW}}_{\mathrm{epi},t+\mathrm{1}}-W_{sep,t+\mathrm{1}},\\ {\mathrm{TT}}_{perc}={\displaystyle \frac{{\mathrm{SAT}}_{\mathrm{epi}}-{\mathrm{FC}}_{\mathrm{epi}}}{K_{\mathrm{epi}}}},\end{array}\right.\end{array}$$

where W_sep is the water content that flows from the epikarst zone to the underground river. Because the regulation and storage functions of the rapid fissure flow are ignored in this study, W_sep refers to the slow fissure flow, W_epi is the current water content of the slow fissure flow in the epikarst zone (where subscripted epi stands for epikarst, and the same applies below), ΔT is the simulation time step, TT_perc is the attenuation coefficient, SAT_epi is the saturation water content of the slow fissure flow, FC_epi is the field capacity, and K_epi is the saturated hydraulic conductivity of the slow fissure flow.

The linear reservoir model is employed to calculate the regulation process of the superficial karst fissure system in the epikarst zone, and the base discharge is calculated by the hydraulic gradient of the KHRU (Neitsch et al., 2002) as follows:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}{Q}_{gw}=\mathrm{8000}{\displaystyle \frac{K_{\mathrm{epi}}{h}_{wtbl}}{{\left(L_{gw}\right)}^{\mathrm{2}}}},\\ Q_{gw,i}=Q_{gw,i-\mathrm{1}}\exp \left(-a_{gw}\mathrm{\Delta }t\right),\\ +W_{rchrg}\left[\mathrm{1}-\exp \left(-a_{gw}\mathrm{\Delta }t\right)\right],\\ W_{rchrg,i}=W_{\mathrm{seep}}\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{gw}}}\right)\right],\\ +W_{rchrg,i-\mathrm{1}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{gw}}}\right),\end{array}\right.\end{array}$$

where Q_gw is the base discharge, Q_gw, i and $Q_{gw,i-\mathrm{1}}$ are the supply quantities of the base discharge that converges into the karst conduit or underground river on the i day and the (i−1) day, respectively, K_epi is the saturated hydraulic conductivity of the epikarst zone, h_wtbl is the hydraulic gradient, L_gw is the length of the KHRU, a_gw is the depletion coefficient of the base discharge, ΔT is the simulation time step (day), W_rchrg, i is the supply quantity of the aquifer on the i day (mm d⁻¹), W_seep is the water flux through the bottom of the soil profile into the underground aquifer on the i day (mm d⁻¹), and δ_gw is the delay time of the supply (day).

In the original Liuxihe model, the underground layer is treated as an integral unit, and a linear reservoir method is used to calculate the underground runoff. However, the structure of the karst underground layer is nonlinear; thus, the linear reservoir method is obviously not appropriate here. Therefore, in this study, the Muskingum routing method was used to calculate the convergence process of the karst underground river, and the equation is as follows:

$$\begin{array}{}\text{(7)}& W=K[xI+(\mathrm{1}-x\left)O\right]=KO{}^{\prime },\end{array}$$

where O^′ is the water storage content, O is the outlet flow of the river reach, x is the dimensionless proportion factor, I is the inflow discharge of the river reach, and K is the slope of the correlation curve of the water storage content and the discharge.

The finite difference method is used to calculate the water balance equation and the Muskingum routing method:

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{l}{O}_{\mathrm{2}}=C_{\mathrm{0}}{I}_{\mathrm{2}}+C_{\mathrm{1}}{I}_{\mathrm{1}}+C_{\mathrm{2}}{O}_{\mathrm{1}},\\ C_{\mathrm{0}}+C_{\mathrm{1}}+C_{\mathrm{2}}=\mathrm{1},\end{array}\right.\end{array}$$

where

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}{C}_{\mathrm{0}}={\displaystyle \frac{\mathrm{0.5}\mathrm{\Delta }t-Kx}{\mathrm{0.5}\mathrm{\Delta }t+K-Kx}},\\ C_{\mathrm{1}}={\displaystyle \frac{\mathrm{0.5}\mathrm{\Delta }t+Kx}{\mathrm{0.5}\mathrm{\Delta }t+K-Kx}},\\ C_{\mathrm{2}}={\displaystyle \frac{-\mathrm{0.5}\mathrm{\Delta }t+K-Kx}{\mathrm{0.5}\mathrm{\Delta }t+K-Kx}}.\end{array}\right.\end{array}$$

If the Muskingum routing method parameters of K and x can be determined for a karst underground river reach, then the values of C₀, C₁ and C₂ can be calculated by Eq. (6). When Δt=2Kx, C₀=0, which means that the karst flood prediction lead time will be 2Kx. Under this condition, the Muskingum routing method can be simplified as follows:

$$\begin{array}{}\text{(10)}& O_{\mathrm{2}}=C_{\mathrm{1}}{I}_{\mathrm{1}}+C_{\mathrm{2}}{O}_{\mathrm{1}}.\end{array}$$

One of the key problems of the Muskingum routing method involves determining how to optimize the parameters K and x in practical applications. It is hard to generalize the parameters K and x in flood simulation and prediction due to their variability with flow conditions. Ahilan et al. (2012) used the generalized extreme value (GEV) to analyze the flood frequency distributions in Irish rivers, and the result showed that a Type II distribution appears in a single cluster in the karst area, which reflects the finite nature of karst storage and the effects of saturation when storage is no longer available. In this study, 30 karst flood events are collected to validate the performance of the Muskingum model in study area. The least squares method is used to optimize the parameters K and x in this study as follows:

$$\begin{array}{}\text{(11)}& min\left\{\left.E=\sum _{j=\mathrm{1}}^n{\left\{W_{\mathrm{0}}\left(j\right)-W_{\mathrm{1}}\left(j\right)-C\right\}}^{\mathrm{2}}\right\},\right.\end{array}$$

where E is the objective function between the observed water storage content and the simulated water storage content, which requires only the least squares approximation with regard to the functional value; W₀(j) and W₁(j) are the observed and simulated water storage contents within the j period, respectively; $W_{\mathrm{1}}\left(j\right)=K[xI+(\mathrm{1}-x\left)O\right]$; n is the total number of observation periods; and C is the absolute value of the water storage content.

To simplify the calculation, $A=K\cdot x$ and $B=K\cdot (\mathrm{1}-x)$; then, the partials can be taken with respect to A,B, and C, respectively:

$$\begin{array}{}\text{(12)}& \left\{\begin{array}{l}\sum W_{\mathrm{0}}I=A\sum I^{\mathrm{2}}+B\sum \left(OI\right)+C\sum I^{\mathrm{2}},\\ \sum W_{\mathrm{0}}O=A\sum \left(OI\right)+B\sum O^{\mathrm{2}}+C\sum O,\\ \sum W_{\mathrm{0}}I=A\sum I+B\sum O+Cn.\end{array}\right.\end{array}$$

Then, the values of $A,B,$ and C can be calculated as follows:

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{l}A={\displaystyle \frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}}-{\displaystyle \frac{y_{\mathrm{3}}}{y_{\mathrm{2}}}},\\ B={\displaystyle \frac{y_{\mathrm{1}}{z}_{\mathrm{2}}}{y_{\mathrm{3}}{z}_{\mathrm{2}}}}-{\displaystyle \frac{y_{\mathrm{2}}{z}_{\mathrm{1}}}{y_{\mathrm{2}}{z}_{\mathrm{3}}}},\\ C=\sum {\displaystyle \frac{W_{\mathrm{0}}-A\sum I-B\sum O}{n}},\end{array}\right.\end{array}$$

where

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{\mathrm{y}}_{\mathrm{1}}=\sum \left(W_{\mathrm{0}}I\right)-{\displaystyle \frac{\sum W_{\mathrm{0}}\sum I}{n}},\\ {\mathrm{y}}_{\mathrm{2}}=\sum I^{\mathrm{2}}-{\displaystyle \frac{{\left(\sum I\right)}^{\mathrm{2}}}{n}},\\ {\mathrm{y}}_{\mathrm{3}}=\sum \left(IO\right)-{\displaystyle \frac{\sum O\sum I}{n}},\\ z_{\mathrm{1}}=\sum \left(W_{\mathrm{0}}O\right)-{\displaystyle \frac{\sum W_{\mathrm{0}}\sum O}{n}},\\ z_{\mathrm{2}}=\sum O^{\mathrm{2}}-{\displaystyle \frac{{\left(\sum O\right)}^{\mathrm{2}}}{n}},\\ z_{\mathrm{3}}=\sum IO-\left({\displaystyle \frac{\sum O\sum I}{n}}\right),\\ K=A+B,\\ x={\displaystyle \frac{K}{A}}.\end{array}\right.\end{array}$$

The parameters of the Muskingum routing method can be optimized using the equations shown above. Then, the convergence process of the karst underground river can be calculated by the Muskingum routing method in the Liuxihe model.

5 Model setup

5.1 Hydrological model setup

The method that combines a DEM with a stream network leads to a more accurate drainage network in terms of surface runoff modeling (Li and Tao, 2000), especially in karst areas. In this study, based on the high resolution of 200 m ×200 m used for the Liuxihe model in the LKRB, the entire studied area was divided into 1 469 900 grid cells, which were named the karst sub-basins, using the DEM. The grid cells included 1 463 204 hillslope cells and 6696 river cells. Then, the karst sub-basins were further divided into many KHRUs. The river system was divided into 3 orders as shown in Fig. 1.

Because of the sinkholes and karst depressions in the karst watershed, as well as the systematic error of the DEM itself, there are many pits, including true and false pits, in the LKRB. Among them, the true pits include karst depressions and sinkholes, and they usually have a certain scale and elevational differences. The false pits were represented only by a few points with low elevation, which was due to the systematic errors of the DEM. Therefore, the true and false pits should be reliably distinguished before using the DEM data to divide the area into the karst sub-basins. First, we identified all of the pits with low elevation and connected them on a plane. Then, we distinguished the true pits from the false pits based on the on-site topographic survey. Finally, the model retained the true pits such as the sinkholes and karst depressions, but the false pits were filled (i.e., removed).

The KHRU was introduced in this study to reasonably describe the spatial variability of the karst water-bearing media (as shown in Fig. 10). The spatial characteristics of every KHRU have a definite physical meaning. Therefore, the calculation of the evapotranspiration, rainfall runoff and parameter optimization of the KHRU was physically based, which could truly reflect the differences of the underlying surface. After the division of the karst sub-basins and the KHRUs, the post-processed PERSIANN-CCS QPE results can be used as the input data for the Liuxihe model to simulate and forecast the karst flood process. The performance of the coupled model was reliably improved in this way.

In the Liuxihe model, the flood process of specific points, named the early warning points of some critical river sections, could be simulated and predicted. Figure 1 shows that there are few rain gauges located upstream of the Liujiang River (which is why the PERSIANN-CCS QPEs were used here). However, the karst is very developed here, and the influence of the karst dominates the runoff processes. Therefore, an early warning point was established at the Goutan River gauge (Fig. 1b) to extract the most developed karst area in the LKRB, Beijiang catchment, where the influence of karst features highly dominates the rainfall–runoff processes. There are 11 early warning points set in the Beijiang catchment (Fig. 1b).

5.2 Parameter optimization of the coupled model

There were 14 parameters that needed to be optimized for the original Liuxihe model, and after adding the karst mechanism, the number of parameters increased to 20, as shown in Table 2. The parameters of the epikarst zone were the most complicated due to the anisotropy of the karst water-bearing media, which made it difficult to measure and calculate the hydraulic characteristics.

Table 2 The parameters of the model.

The hydrogeology parameters used in this study, including the permeability coefficient of the rock mass, the rainfall infiltration coefficient, the specific yield of the aquifer, and the storage coefficient, were calculated by the field test and the experience function. For instance, the permeability coefficient K was calculated by an experience function according to the water inrush prediction of a coal mine in the study area:

$$\begin{array}{}\text{(15)}& {\displaystyle }\left\{\begin{array}{l}Q=\mathrm{1.366}K{\displaystyle \frac{(\mathrm{2}H-M)\cdot M-h^{\mathrm{2}}}{\lg R_{\mathrm{0}}-\lg r_{\mathrm{0}}}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{24}}},\\ R_{\mathrm{0}}=r_{\mathrm{0}}+\mathrm{10}\cdot S\sqrt{K},\\ r_{\mathrm{0}}=\sqrt{{\displaystyle \frac{a\cdot b}{\mathit{\pi }}}},\end{array}\right.\end{array}$$

where Q is the mine inflow, m³ h⁻¹; K is the permeability coefficient, m d⁻¹; H is the distance from the water-resisting floor to the water level of the confined aquifer, m; M is the aquifer thickness, m; h is the height of the dynamic water level, m; R₀ is the substitute influence radius, m; r₀ is the substitute radius, m; S is the drawdown value, m; and a⋅b is the area of the mine, m².

In the water inrush test of the coal mine, the other parameters in Eq. (15) were given, and the permeability coefficient K was calculated by anti-Eq. (15).

The parameters of the epikarst zone, including the thickness, saturated water content, field capacity and macro crack volume ratio, were obtained based on the field survey, geological borehole test and pumping test as well as on the empirical value observed in the study area.

The epikarst zone was mainly developed on the hard surface of pure carbonate rock, especially on Paleozoic limestone. The thicknesses and characteristics of the epikarst zone differ due to different climates, topography and landforms. The parameters of the coupled model and the epikarst zone are listed in Table 2a and b, and the rainfall infiltration coefficients of the different karst landforms are calculated based on the empirical values shown in Table 2c.

The soil type parameters, such as the saturated water content and the field capacity, were calculated using a software tool (Ren, 2006). The statistical relationship between the soil texture and the soil water can be easily queried in the software tool. In addition, this method has been effectively proven by many experiments (Servat and Sakho, 1995), and the calculated value of this method has a good fitting relationship with the measured value.

The Liuxihe model has been deployed on a supercomputer system with parallel computation technology (Chen et al., 2016). An improved PSO algorithm (Chen et al., 2017) was employed to optimize the parameters of the coupled model in this study. There are 30 karst flood events from 1982 to 2013 in the LKRB, and among them, 3 flood events – Floods 2004070300, 2009060908, and 2011010100 – were used for parameter optimization simulations in this paper. The flood simulation results are shown in Fig. 11 and Table 3.

Figure 11 The flood simulation results obtained through parameter optimization by the improved PSO algorithm.

From the flood simulation results in Fig. 11, it can be seen that the Flood 2009060908 simulated result is the best. The simulated process for this flood is closest to the observed process, and the valuation indices of flood simulation results including the Nash–Sutcliffe coefficient, C; correlation coefficient, R; process relative error, P %; peak flow relative error, E %; coefficient of water balance, W; and peak time error, T(h), are also the best. Table 3 shows the valuation indices of flood simulation results from the improved PSO algorithm. Therefore, Flood 2009060908 is finally adopted for the Liuxihe model parameter optimization. Other floods will be used to verify the model performance.

Table 3 The evaluation indices of flood simulation results obtained through parameter optimization by the improved PSO algorithm.

The parameter optimization results from the improved PSO algorithm are shown in Fig. 12 as follows: (a) the objective function evolution result, (b) the parameter evolution result, and (c) the simulated flood process using the optimized model parameters.

Figure 12 Parameter optimization results with the improved PSO algorithm.

To test the parameter optimization effect with different precipitation sources, both the precipitation of the rain gauge and the precipitation of the PERSIANN-CCS QPEs were used to optimize the parameters of the coupled model. For comparison, the simulated flood process of the coupled model with the same parameter from the rain gauges and the re-optimized parameter from the post-processed PERSIANN-CCS QPEs are drawn in Fig. 12c.

5.3 Parametric uncertainty analysis

In this study, parametric uncertainty analysis refers to sensitivity analysis, and this process is conducted using a fixed module called the parametric sensitivity analysis sub-model in the Liuxihe model. It is a parameter sensitivity analysis method that was developed based on the GLUE method, and it was named multi-parameter sensitivity analysis (MPSA) by Choi et al. (1999). Monte Carlo sampling was used to obtain the value of the parameter spatial variation. The sensitivity of each parameter could be obtained by running the model multiple times.

In this study, the Nash–Sutcliffe coefficient was used as the objective function value for the parametric sensitivity analysis, and the formula is as follows:

$$\begin{array}{}\text{(16)}& \mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left(Q_i-Q_i^{\prime }\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left(Q_i-\stackrel{\mathrm{‾}}{Q}\right)}^{\mathrm{2}}}},\end{array}$$

where NSE is the objective function value of the Nash–Sutcliffe coefficient, Q_i and $Q_i^{\prime }$ are the observed streamflow and the simulated streamflow, respectively, in m³ s⁻¹, $\stackrel{\mathrm{‾}}{Q}$ is the average value of the observed flows in m³ s⁻¹, and n is the number of observation periods in h.

First, the initial value range of the parameter was determined to be [0.5, 2.5]. Second, 6000 groups of parameter sequences were obtained by the Monte Carlo sampling method. Third, the Liuxihe model was run to simulate the objective function values of the Nash–Sutcliffe coefficient, and the karst flood processes were the three flood events also used for parameter optimization. In this study, the critical value of the Nash–Sutcliffe coefficient was 0.85, and the objective function values below this threshold were considered to be unacceptable values; otherwise, they were considered to be acceptable values. The degree of separation between these values indicates the sensitivity of the parameters. This degree of separation was calculated according to the Nash–Sutcliffe coefficient (NSD). To analyze parameter sensitivity more easily, a factor SI is given here, and SI $=\mathrm{1}-\left|\mathrm{NSD}\right|$ – the closer the value of SI is to 0, the less sensitive the parameter. Table 4 shows the SI values, which represent the sensitivity of the parameters in the Liuxihe model.

Table 4 The calculation results of the parameter sensitivity in the Liuxihe model.

6 Results and discussion

6.1 Results of parameter optimization and sensitivity analysis

The results of the parameter optimization are shown in Fig. 12 as follows: (a) the objective function evolution result and (b) the parameter evolution result. From the results of Fig. 12a and b, it can be seen that the evolution number of the objective function for the parameter was 50, and the computation time of the parameter optimization based on the improved PSO algorithm was approximately 8 h, which means that convergence of the parameter optimization was achieved after only 50 cycles. In comparison, the computation time of the initial model parameters that were not optimized was approximately 55 h. This result implies that the improved PSO algorithm had high efficiency in terms of parameter optimization.

To test the parameter optimization effect using the improved PSO algorithm (Chen et al., 2017), the flood process simulated results achieved from the improved PSO algorithm, as well as the initial model parameter values, are shown in Fig. 12c. From the results shown in Fig. 12c, it can be seen that the coupled model does not simulate the observed karst flood process well when the initial model parameter values are used. Additionally, the simulated flood process obtained from using the improved PSO algorithm was very close to that from the observed process, which means that the improved PSO algorithm (Chen et al., 2017) in this study was effective and could largely improve the performance of the coupled model.

In this study, the sensitivity of the parameters in the Liuxihe model was calculated according to the Nash–Sutcliffe coefficient, as shown in Eq. (16). The values of SI $=\mathrm{1}-\left|\mathrm{NSD}\right|$, which represent the sensitivity of the parameters, and the results in Table 4 indicate that the SI values of the saturated water content parameter, θ_sat, were maximized, which means that the degree of separation between the unacceptable values and the acceptable values (NSD) was minimal. This parameter, θ_sat, was the most sensitive parameter in the Liuxihe model. When the SI value of a parameter is greater than 0.7, this parameter is identified as a highly sensitive parameter in the Liuxihe model, and SI values between 0.2 and 0.7 indicate that a parameter has medium sensitivity. When the SI value is less than 0.2, the parameter is insensitive. From Table 4, the SI values of the different parameters, from largest to smallest, are the saturated water content, θ_sat> saturation permeability coefficient, θ_s> field capacity, θ_fc> saturated hydraulic conductivity, K_s> macro crack volume ratio, V> Muskingum routing method (the slope of the water storage content and flow curve), K> Muskingum routing method (the proportion of the flow), χ> soil layer thickness, z> soil coefficient, b> bottom width, S_w> bottom slope, S_p> slope roughness, n> channel roughness, n₁> depletion coefficient, ω> evaporation coefficient, λ> potential evaporation, E_p> wilting percentage, and C_wl. Additionally, the θ_sat, θ_s, θ_fc, K_s, V, K, and χ parameters were highly sensitive; the z, b, S_w, S_p, n, n₁ and ω parameters had medium sensitivity; and the λ, E_p, and C_wl parameters were insensitive.

The flow direction, slope and thickness parameters of the epikarst zone could not be adjusted. Among them, the flow direction and the slope were directly calculated by the DEM data, and the thickness of the epikarst zone was a fixed value for a particular region. It was approximately 3–10 m of the study area according to the field survey.

6.2 Model validation results

To better test the effect of the Liuxihe model in flood simulation and prediction and to increase the results' acceptability, 30 karst flood events from 1982 to 2013 in LKRB are simulated by the Liuxihe model, and the evaluation indices of the simulated flood results are listed in Table 5. Table 5 shows that the six evaluation indices of the flood simulation results for the 30 flood events are credible and reasonable. The average value of the Nash–Sutcliffe coefficient (C) is 0.82, the correlation coefficient (R) is 0.83, the process relative error (P) is 0.22, the peak flow relative error (E) is 0.05, the water balance coefficient (W) is 0.87, and the peak flow time error (T) is −6 h. Among these results, the peak flow relative error (E) is minimal. The applicability of the Liuxihe model is proven through these accepted flood simulation effects in the LKRB.

Table 5 The evaluation indices of the simulated flood results based on the Liuxihe model in the LKRB.

Table 6 The evaluation indices of the simulated flood results based on the Liuxihe model in the Beijiang catchment.

To further validate the performance of the Liuxihe model in flood simulation and prediction, simulations are performed in a very developed karst area, where the influence of karst landforms plays an important role in hydrological processes. The most developed karst area in the whole basin examined in this study is the Beijiang catchment, and it is divided by the early warning point Goutan set in the Liuxihe model (Fig. 1b). In total, 10 karst flood events are simulated to test the performance of the Liuxihe model, and the evaluation indices of the simulated flood results are shown in Table 6. From these results, four karst flood simulation results are shown in Fig. 13.

Figure 13 Karst flood simulation results from the Liuxihe model in the Beijiang catchment.

From the results in Table 6, the evaluation indices of the simulated karst flood results produced by the Liuxihe model are quite good in the Beijiang catchment. The average value of the Nash–Sutcliffe coefficient (C) is 0.92, the correlation coefficient (R) is 0.91, the process relative error (P) is 0.11, the peak flow relative error (E) is 0.08, the water balance coefficient (W) is 0.94, and the peak flow time error (T) is 3 h. It is obvious that the evaluation indices of the simulated karst flood events based on the Liuxihe model are satisfying, and the accuracy is very high.

Additionally, from the flood simulation results in Fig. 13, the four reasonable karst flood simulation results, including those for Floods 2008071311, 2012080310, 2014061015, and 2016091501, prove the performance of the Liuxihe model in karst areas. The simulated flood discharge processes are very close to the observed values, especially for the peak flows. This finding implies that the Liuxihe model is feasible and effective in flood simulation and prediction in areas where karst is very well developed, as in the Beijiang catchment.

6.3 Results of flood simulation with the post-processed PERSIANN-CCS QPEs

After the correction was made, the post-processed PERSIANN-CCS QPE precipitation became much closer to the precipitation observed at the rain gauge. To analyze the effects of flood simulation with the initial PERSIANN-CCS QPEs and the post-processed QPEs, five karst flood events, including Floods 200806090200, 200906090800, 201106010900, 201206022000 and 201306011400, were simulated and compared; the results are shown in Fig. 14. In this simulation, the coupled model parameters remained unchanged; i.e., the original coupled model parameters based on the rain gauge precipitation were employed, while the re-optimized model parameters based on the precipitation of the post-processed PERSIANN-CCS QPEs were not.

Figure 14 The flood simulation results of the coupled model with the two precipitation products.

Table 7 Evaluation indices of simulated flood events using the initial PERSIANN-CCS QPEs and the post-processed values.

Figure 14 shows that the karst flood simulation results from the initial PERSIANN-CCS QPEs were not satisfactory, and the performance of the model was worse than that of the rain gauge precipitation. For instance, the simulated peak flows from the PERSIANN-CCS QPEs were lower than the observed peak flows. The performance of the coupled model with the post-processed PERSIANN-CCS QPEs was much better, and the evaluation indices of the flood simulation were largely improved (as shown in Table 7). The average value of the Nash–Sutcliffe coefficient (C) increased by 7 %, the correlation coefficient (R) increased by 8 %, the process relative error (P) decreased by 6 %, the peak flow relative error (E) decreased by 14 %, the water balance coefficient (W) increased by 5 %, and the peak flow time error (T) had a decrease of 2 h. Among these parameters, the peak flow relative error had the largest improvement, making it the most important factor in flood prediction. It was obvious that the evaluation indices improved substantially when the post-processed QPEs were used. Therefore, the post-processing method for PERSIANN-CCS QPEs in this paper was feasible and effective. In addition, coupling the post-processed PERSIANN-CCS QPEs with the Liuxihe model has the potential to improve the model performance in flood simulation and prediction in the LKRB.

6.4 Comparisons of different model parameters

The model parameters that were optimized using the precipitation from the rain gauge and those optimized using the PERSIANN-CCS QPEs were different, and the performance of the coupled model using the different parameters made a large difference in the flood simulation and prediction. To analyze this effect, the flood simulation results from two different sets of model parameters are shown in Fig. 15. One set used the parameters of the coupled model that was optimized by the precipitation from the rain gauge; i.e., the coupled flood simulation results used the same parameter as the rain gauge precipitation. The other used the parameters that were re-optimized by the post-processed PERSIANN-CCS QPEs. The flood process used for parameter re-optimization was also Flood 2009060908, and the other four flood events were used to validate the performance of the coupled model.

Figure 15 Coupled flood simulation results using the same parameter as the rain gauge precipitation and using the re-optimized parameter from the post-processed PERSIANN-CCS QPEs. (a) Flood event 200806090200. (b) Flood event 201106010900. (c) Flood event 201206022000. (d) Flood event 201306011400.

Table 8 The effect of recalibrating the coupling model parameters.

Figure 15 shows that the simulated flood results obtained using the re-optimized parameters from the post-processed PERSIANN-CCS QPEs were much better than those obtained using the same parameter as the rain gauge precipitation. The simulated flood discharge processes, especially the peak flows with the re-optimized parameter, were closer to the observed values. To further compare the flood simulation results, six evaluation indices were calculated and are shown in Table 8. The average value of the Nash–Sutcliffe coefficient increased by 7 %, the correlation coefficient increased by 7 %, the process relative error decreased by 2 %, the peak flow relative error decreased by 4 %, the water balance coefficient increased by 3 %, and the peak flow time error exhibited a 3 h decrease.

Moreover, compared with the simulated flood results from the initial PERSIANN-CCS QPEs in Table 8, the flood simulation results with the re-optimized parameters from the post-processed PERSIANN-CCS QPEs made great progress. The average value of the Nash–Sutcliffe coefficient increased by 14 %, the correlation coefficient increased by 15 %, the process relative error decreased by 8 %, the peak flow relative error decreased by 18 %, the water balance coefficient increased by 8 %, and the peak flow time error had a 5 h decrease (as shown in Tables 7 and 8). These results imply that the re-optimized parameters calculated using the post-processed PERSIANN-CCS QPEs are necessary and effective for the coupled model, and the model performance improved in terms of karst flood simulation and prediction.

6.5 Peak flow time error analysis

It is very important to accurately determine the flood peak flow time in karst areas, as this information could improve the response times of safe and rapid evacuations before a flood disaster appears. As shown in Figs. 14 and 15 and in Tables 7 and 8, all flood simulations had significant peak flow time errors, and all of the errors were negative, indicating that the simulated flood peaks appeared earlier than did the peaks in the observed values. Among them, the average peak flow time error from the precipitation from the rain gauge was −7 h, and this value was −10 h when the precipitation from the initial PERSIANN-CCS QPEs was used. This is an obvious error and cannot be ignored in flood prediction. The average peak flow time error of the coupled model that used the post-processed PERSIANN-CCS QPE precipitation and re-optimized parameters was −5 h. This result indicates that there is a great difference. It has been found that the average peak flow time errors of the Liuxihe model generated from the precipitation from the rain gauge and from the coupled model that used the precipitation from the post-processed PERSIANN-CCS QPEs and re-optimized parameters were −5 to −7 h (as shown in Tables 7 and 8). Therefore, the peak flow time error was −5 to −7 h for the coupled model in the LKRB, which means that the actual time of the flood peak may be 5–7 h later. This value is very important in flood prediction and is equivalent to a 5–7 h lead time in which safe evacuations can occur.

There are two reasons for the peak flow time errors. One reason is the systematic error of the coupled model itself. This error could be reduced by improving the model structure and function as well as by the reliable precipitation from the PERSIANN-CCS QPEs and parameter optimization. The other reason is due to the karst development laws and the characteristics of karst water-bearing media, which can regulate the rainfall process during floods. The karst depressions and other negative landforms in the upstream regions can hold back and store large amounts of floodwater. Furthermore, karst fissures can also slow the flood rate. These factors can play a crucial role in detaining natural floods and clipping the flood peaks. Therefore, the response times of the flood peak flow to the rainfall increased, and the observed flood peak times lagged behind. In comparison, the simulated flood peak flows appeared earlier.

As rainfall moves from the sky to the ground and, finally, to the point where the rainfall converges at the outlet of the basin, it has passed through the surface karst zone, the karst conduit and fissure as well as the underground river. The karst development laws and the characteristics of the karst water-bearing media have an obvious influence on the rainfall–runoff process during the entire hydrological process, which increases the response time of the flood peak flow to rainfall, and the simulated flood peak flow in the coupled model appears earlier. This result implies that there is a lead time that can be used for safe evacuation measures.

The flood peak flow time has a very close relationship with the flood rate, and the flood rate is very important in determining the key factors of the karst conduit, the underground river and the other hydrogeological parameters. The sensitive parameters in this paper, such as the underground river parameters (as shown in Table 2), could be estimated from the flood rate to build the coupled model in the karst area. According to the survey data and the tracing test in the study area, i.e., the LKRB, the flood flow rate is approximately 8.64–17.28 km d⁻¹ in the dry season, 17.28–43.2 km d⁻¹ in the normal season and 43.2–129.6 km d⁻¹ in the flood period. The extreme flow rate can reach 172.8 km d⁻¹, indicating that the karst conduit is highly developed in the LKRB.

7 Conclusions

Few reliable precipitation data from rain gauges are available in most karst river basins. How to obtain reasonable rainfall data for the development of a hydrological model that can be used for flood prediction is especially important. In this study, the PERSIANN-CCS QPEs offered effective precipitation results for the study area. After the correction, the post-processed PERSIANN-CCS QPEs coupled with a distributed hydrological model, i.e., the Liuxihe model, were proposed for karst flood simulation and prediction in the LKRB. The purpose of the study was not only to simulate the flood process well, but also to determine key information about how the karst hydrological process responds to the rainfall process in the coupled model. The coupled model employed in this paper had good performance in simulating flood events; thus, this method offers reasonable theoretical guidance for flood prediction, control and disaster reduction in karst river basins such as the LKRB. Based on the study results, the following conclusions can be drawn.

1. The quantitative precipitation estimates produced by the PERSIANN-CCS QPEs were very similar to the observed precipitation from the rain gauges, especially in terms of rainfall distribution. However, the PERSIANN-CCS QPEs underestimated the precipitation value. The average precipitation was 0.77 for the rain gauges and 0.66 for the PERSIANN-CCS QPEs. The average relative error was −14 % between the two precipitation products, and this relative error could be reasonably reduced by the post-processing method presented in this paper.

2. The applicability of the Liuxihe model is proven by 30 accepted flood simulation results in the LKRB and 10 in the Beijiang catchment. In particular, the simulated results are quite good for the 10 karst flood events in the Beijiang catchment, where the karst is very developed. The average value of the Nash–Sutcliffe coefficient (C) is 0.92, the correlation coefficient (R) is 0.91, the process relative error (P) is 0.11, the peak flow relative error (E) is 0.08, the water balance coefficient (W) is 0.94, and the peak flow time error (T) is 3 h.

   The parameter sensitivity analysis for the Liuxihe model shows that the parameters θ_sat, θ_s, θ_fc, K_s, V, K, and χ are highly sensitive; z, b, S_w, S_p, n, n₁ and ω have medium sensitivity; and λ, E_p, C_wl are insensitive parameters. The sequence of parameters sensitivity is as follows: saturated water content, θ_sat> saturation permeability coefficient, θ_s> field capacity, θ_fc> saturated hydraulic conductivity, K_s> macro crack volume ratio, V> Muskingum routing method (the slope of the water storage content and flow curve), K> Muskingum routing method (the proportion of the flow), χ> soil layer thickness, z> soil coefficient, b> bottom width, S_w> bottom slope, S_p> slope roughness, n> channel roughness, n₁> depletion coefficient, ω> evaporation coefficient, λ> potential evaporation, E_p> wilting percentage, C_wl.

3. The flood simulation results from the post-processed PERSIANN-CCS QPEs are better than that from the initial QPEs. The average values of the six evaluation indices, including the Nash–Sutcliffe coefficient (C), correlation coefficient (R), process relative error (P), peak flow relative error (E), water balance coefficient (W), and peak flow time error (T), with the initial PERSIANN-CCS QPEs were 0.66 %, 0.69 %, 0.28 %, 24 %, 0.81 and −10 h, respectively, while those from the post-processed QPEs were 0.73 %, 0.77 %, 0.22 %, 10 %, 0.86 and −8 h, respectively. This result indicates that the method used in this study for post-processing QPEs is effective and could improve the PERSIANN-CCS QPE capability.

4. The coupled model parameters should be re-optimized using the post-processed PERSIANN-CCS QPEs. This approach had better performance in the flood simulation than that when the model parameters were the same as those from the rain gauges. The average values of the Nash–Sutcliffe coefficient (C), correlation coefficient (R), process relative error (P), peak flow relative error (E), water balance coefficient (W), and peak flow time error (T) were 0.73 %, 0.77 %, 0.22 %, 10 %, 0.86 and −8 h, respectively, when the model parameters were the same as the rain gauge; however, those obtained from the re-optimized model parameters were 0.80 %, 0.84 %, 0.20 %, 6 %, 0.89 and −5 h, respectively. Thus, the proposed method significantly improves the model performance.

5. The simulated karst flood process based on the precipitation observed at the rain gauges was the best. In addition, the flood simulation results using the PERSIANN-CCS QPEs after post-processing and re-optimizing the model parameters improved the coupled model performance. The average value of the Nash–Sutcliffe coefficient increased by 14 %, the correlation coefficient increased by 15 %, the process relative error decreased by 8 %, the peak flow relative error decreased by 18 %, the water balance coefficient increased by 8 %, and the peak flow time error exhibited a 5 h decrease. Among these parameters, the peak flow relative error improved the most; thus, these parameters are the most important in terms of flood prediction in karst river basins.