A physically based distributed karst hydrological model (QMG 1 model-V1.0) for flood simulations

Karst trough and valleys landforms are prone to flooding, primarily because of 13 the unique hydrogeological features of karst landform, which are conducive to the spread of 14 rapid runoff. Hydrological models that represent the complicated hydrological processes in 15 karst regions are effective for predicting karst flooding, but their application has been 16 hampered by their complex model structures and associated parameter set, especially so for 17 distributed hydrological models, which require large amounts of hydrogeological data. 18 Distributed hydrological models for predicting the Karst flooding is highly dependent on 19 distributed structrues modeling, complicated boundary parameters setting, and tremendous 20 hydrogeological data processing that is both time and computational power consuming. 21 Proposed here is a distributed physically-based karst hydrological model, known as the 22


Introduction
Karst trough and valley landformsKarst trough valleys are very common in China, especially in the southwest.In general, these karst areas are water scarce during most of the year because their surfaces store very little rainfall, but they are also potential birthplaces for floods because trough and valley landforms and topographic features facilitate the formation and propagation of floods (White, 2002;Li et al., 2021).The coexistence of drought and flood is a typical phenomenon in these karst trough and valley areas.Taking the example of the present study area, i.e. the Qingmuguan karst trough valley, floods used to happen here constantly during the rainy season.In recent years, with more extreme rainfall events and the increased area of construction land in the region, rainfall infiltration has decreased and rapid runoff over impervious surfaces has increased, resulting in frequent catastrophic flooding in the basin (Liu et al., 2009).Excess water overflows from karst sinkholes and underground river outlets often occur during floods (Jourde et al., 2007(Jourde et al., , 2014;;Martinotti et al., 2017), flooding large areas of farmland and residential areas and causing serious economic losses (Gutierrez, 2010;Parise, 2010;Yu et al., 2020).Therefore, it is both important and urgent to simulate and predict karst flooding events in karst trough and valleys such as the study area.
Hydrological models can be effective for forecasting floods and evaluating water resources in karst areas (Bonacci et al., 2006;Ford and Williams, 2007;Williams, 2008Williams, , 2009)).However, modelling floods in karst regions is extremely difficult because of the complex hydrogeological structure.Karst water-bearing systems consist of multiple media under the influence of complex karst development dynamics (Worthington et al., 2000;Ková cs and Perrochet, 2008;Gutierrez, 2010), such as karst caves, conduits, fissures and pores, and are usually highly spatially heterogeneous (Chang and Liu, 2015;TeixeiraparenteMario et al., 2019).In addition, the intricate surface hydrogeological conditions and the hydrodynamic conditions inside the karst water-bearing medium result in significant temporal and spatial differences in the hydrological processes in karst areas (Geyer et al., 2008;Bittner et al., 2020).
In early studies of flood forecasting in karst regions, simplified lumped hydrological models were commonly used to describe the rainfall-discharge relationship (e.g.Ková cs and Sauter, 2007;Fleury et al., 2007b;Jukić and Denić, 2009;Hartmann et al., 2014a).With the development of physical exploration technology and the progress made in mathematics, computing and other interdisciplinary disciplines, the level of modelling has gradually improved (Hartmann and Baker, 2017;Hartmann, 2018;Petrie et al., 2021), and distributed hydrological models have subsequently become widely used in karst areas.The main difference between lumped and distributed hydrological models is that the latter divide the entire basin into many sub-basins to calculate the runoff generation and confluence, thereby better describing the physical properties of the hydrological processes inside the karst water-bearing system (Jourde et al., 2007;Hartmann, 2018;Epting et al., 2018).
Because of their simple structure and little demand for modelling data, lumped hydrological models have been used widely in karst areas (Kurtulus and Razack, 2007;Ladouche et al., 2014).In a lumped model, the river basin is considered as a whole to calculate the runoff generation and confluence, and there is no division running into sub-basins (Dewandel et al., 2003;Bittner et al., 2020).Lumped models usually consider the inputs and outputs of the model (Liedl and Sauter, 2003;Hartmann andBake, 2013, 2017).
In addition, most of the model parameters are not optimized in a lumped model, and the physical meaning of each parameter is unclear (Chen, 2009;Bittner et al., 2020).
Distributed hydrological models are of active interest in flood simulation and forecasting research (Ambroise et al., 1996;Beven and Binley, 2006;Zhu and Li, 2014).
Compared with a lumped model, a distributed model has a more definite physical significance for the model structure in terms of its mechanism (Meng and Wang, 2010;Epting et al., 2018).In a distributed hydrological model, an entire karst basin can be divided into many sub-basins (Birk et al., 2005) using high-resolution digital elevation map (DEM) data.In the rainfall-runoff algorithm of the model, the hydrogeological conditions and karst aquifer characteristics can be considered fully to simulate precisely the runoff generation and confluence (Martinotti et al.,2017;Gang et al., 2019).The commonly used basin distributed hydrological models (i.e.not a special groundwater numerical model such as MODFLOW) have also been applied widely in karst areas and include the SHE/MIKE SHE model (Abbott et al., 1986a,b;Doummar et al., 2012), SWMM model (Peterson and Wicks, 2006;Blansett and Hamlett, 2010;Blansett, 2011), TOPMODEL (Ambroise et al., 1996;Suo et al., 2007;Lu et al., 2013;Pan, 2014) and SWAT model (Peterson and Hamlett, 1998;Ren, 2006).
The commonly used distributed hydrological models have multiple structures and numerous parameters (Lu et al., 2013;Pan, 2014), which means that a distributed model may need vast amounts of data to build its framework in karst regions.For example, the distributed groundwater model MODFLOW-CFPM1 requires detailed data regarding the distribution of karst conduits in a study area (Reimann et al., 2009).Another example is the Karst-Liuxihe model (Li et al., 2019), which has there are fifteen parameters and five underground vertical structureslayers in the model.structure and has 15 parameters, thereby Such a complex structure makes the modeling data demand is large, and the modeling in karst area is extremely difficult.makingit difficult to model in karst areas.In addition, a special borehole pumping test may be required to obtain the rock permeability coefficient.
To overcome the difficulty of the large modelling-data demands for distributed hydrological models in karst areas, a new physically based distributed hydrological model-known as the QMG (Qingmuguan) model-V1.0-wasdeveloped in the present study.Other commonly used karst groundwater models with complex structure and parameters-such as the aforementioned MODFLOW-CFPM1 model-require a lot of hydrogeological data for modelling in karst areas (Qin and Jiang, 2014).The new QMG model has a high potential for application in karst hydrological simulation and forecasting.It has certain advantages in its framework and structural design, having a double-layer structure and fewer parameters.The horizontal structure is divided into river channel units and slope units, and the vertical structure below the surface is divided into a shallow karst aquifer and a deep karst aquifer system.This relatively simple model structure reduces the demand for modelling data in karst areas, and only a small amount of hydrogeological data is needed for modelling.
To ensure that the QMG model work well in karst flood simulation and prediction in the case of relatively simple structure and parameters.We carefully designed the algorithms of runoff generation and confluence in the model.To ensure that the QMG model works well in karst flood simulation and prediction despite its relatively simple structure and parameters, we carefully designed the algorithms for runoff generation and confluence in the model.Also, to verify the applicability of the QMG model to flood simulation in karst basins, we selected the Qingmuguan karst trough valley in Chongqing, China as the study area for a flood simulation and uncertainty analysis.In particular, we analysed the sensitivity of the model parameters.

Landform and topography
The Qingmuguan karst trough valley is located in the southeastern part of the Sichuan Basin, China at the junction of the Beibei and Shapingba districts in Chongqing, with the coordinates of 29°40′N-29°47′N, 106°17′E-106°20′E.The basin covers an area of 13.4 km 2 and is part of the southern extension of the anticline at Wintang Gorge in the Jinyun Mountains, with the anticlinal axis of Qingmuguan located in a parallel valley in eastern Sichuan (Yang et al., 2008).The surface of the anticline is heavily fragmented, and faults are extremely well developed with large areas of Triassic carbonate rocks exposed.Under the long-term erosion of karst water, a typical karst trough landforms pattern of 'three mountains and two troughs' has formed (Liu et al., 2009).This karst trough landform provides convenient conditions for flood propagation, and the development of karst landforms is extremely common in the karst region of southwest China, especially in the karst region of Chongqing.Similar regions include the karst trough valley of the Zhongliang Mountains and the Laolongdong karst basin in Nanshan, Chongqing.
The basin is oriented north-north-east and south-south-west in a narrow band of slightly curved arcs and is ~12 km long from north to south.The direction of the mountains in the region is basically the same as that of the tectonic line.The direction of the mountains in the region is generally consistent with the direction of the tectonic line.The difference in relative elevation is 200-300 m.The map in Figure . 1 gives an overview of the Qingmuguan karst basin.

Hydrogeological conditions
The Qingmuguan basin is located within the subtropical humid monsoon climate zone, with an average temperature of 16.5°C and an average precipitation of 1250 mm that is concentrated mainly in May-September.An underground river system has developed in the karst trough valley, with a length of 7.4 km, and the water supply of the underground river is mainly rainfall recharge (Zhang, 2012).Most of the precipitation is collected along the hill slope into the karst depressions at the bottom of the trough valley, where it is recharged to the underground river through the dispersed infiltration of surface karst fissures and concentrated injection from sinkholes (Fig. 1a).An upstream surface river collects in a gentle valley and enters the underground river through the Yankou sinkhole (elevation 524 m).Surface water in the middle and lower reaches of the river system enters the underground river system mainly through catenuliform cover collapse sinkholes (Gutierrez et al., 2014) or fissures.
The stratigraphic and lithologic characteristics of the basin are dominated largely by carbonate rocks of the Lower Triassic Jialingjiang Group (T 1j ) and Middle Triassic Leikou Slope Group (T 2l ) on both sides of the slope, with some quartz sandstone and mudstone outcrops of the Upper Triassic Xujiahe Group (T 3xj ) (Zhang, 2012).The topography of the basin presents a general anticline (Fig. 1b), where carbonate rocks on the surface are corroded and fragmented, high permeabilitywith a large permeability coefficient.Compared with the core of the anticline, the shale rocks of the two wings of the anticline are less eroded and form a good waterproof layer.
To investigate the distribution of karst conduits in the underground river system, we conducted a tracer test in the study area.The tracer was placed into the Yankou sinkhole and recovered in the Jiangjia spring (Fig. 1a,c).According to the tracer test results (Gou et al., 2010), the karst water-bearing medium in the aquifer was anisotropic, whereas the soluble carbonate rocks were extremely permeable.and The the karst conduits in the underground river were extremely well developed, and there was a large single-channel underground river about five meters wide.The response of the underground river to rainfall was very fast, with the peak flow observed at the outlet of Jiangjia spring 6-8 h after rainfall based on the tracer test results.The flood peak rose quickly and the duration of the peak flow was short.The underground river system in the study area is dominated by large karst conduits, which is not conducive to water storage in water-bearing media, but is very conducive to the propagation of floods.

Data
To build the QMG model to simulate the karst flood events, the necessary modelling baseline data had to be collected, including: 1) high-resolution DEM data and hydrogeological data (e.g., the thickness of the epikarst zone, rainfall infiltration coefficient on different karst landforms, and permeability coefficient of rock); 2) land-use and soil type data; and 3) rainfall data in the basin and water flow data of the underground river.The DEM data was downloaded from a free database on the public Internet, with an initial spatial resolution of 30  30 m.The spatial resolution of landuse and soil types were 1000  1000 m, and they were also downloaded from the Internet.After considering the applicability of modelling and computational strength, as well as the size of the basin in the study area (13.4 km 2 ), the spatial resolution of the three types of data was resampled uniformly in the QMG model and downscaled to 15  15 m based on a spatial discrete method by Berry et. (2010).
The hydrogeological data necessary for modelling was obtained in three simple ways. 1) A basin survey was conducted to obtain the thickness of the epikarst zone, which was achieved by observing the rock formations on hillsides following cutting for road construction.Information was collected regarding the location, general shape, and size of karst depressions and sinkholes, which had a significant impact on compiling the DEM data and determining the convergence process of surface runoff.And the sinkholes in the basin are cover collapse sinkholes (Gutierrez et al., 2014) according to the basin survey.There are 3 large sinkholes (more than 3 meters in diameter) and 12 small sinkholes less than 1 meter in diameter.The rest of the sinkholes between 1 and 3 meters in diameter are 5 in total.
The confluence calculation of these sinkholes in the model was based on the results of previous study (Meng et al., 2009).2) Empirical equations developed for similar basins were used to obtain the rainfall infiltration coefficient for different karst landforms and the permeability coefficient of rock.For example, the rock permeability coefficient was calculated based on an empirical equation from a pumping test in a coal mine in the study area (Li et al., 2019).3) A tracer experiment was conducted in the study area (Gou et al., 2010) to obtain information on the underground river direction and flow velocity., for instance, underground karst conduits are well developed in the area, which form an underground river about five meters wide.There is no hydraulic connection between the underground river system in the area and the adjacent basin, means there is no overflow recharge.
Rainfall and flood data are important model inputs, and represent the driving factors that allow hydrological models to operate.In the study area, rainfall data was acquired by two rain gauges located in the basin (Fig. 1a).Point rainfall was then spatially interpolated into basin-level rainfall (for such a small basin area rainfall results obtained from two rain gauges was considered representative).There were 18 karst flood events in the period of 14 April 2017 to 10 June 2019.We built a rectangular open channel at the underground river outlet and set up a river gauge on it (Fig. 1a) to record the water level and flow data every 15 minutes.

Hydrological model framework and algorithms
The hydrological model developed in this study was named the QMG model after the basin for which it was developed and to which it was first applied, i.e. the Qingmuguan basin.The QMG model proposed in this study has a two-layer structure, including a surface part and an underground part, .with theThe former surface structure mainly performing the calculation of runoff generation and the confluence of the surface river, while the underground structure latter performs the confluence calculation of the underground river system.
The structure of the QMG model is divided into a two-layer structure, both horizontally and vertically.The horizontal structure of the model is divided into river channel units and slope units.The vertical structure below the surface is divided into a shallow karst aquifer (including soil layers, karst fissures and conduit systems in the epikarst zone) and a deep karst aquifer system (bedrockrock stratum and underground river system).This relatively simple model structure means that only a small amount of hydrogeological data is needed when modelling in karst regions.Figure 2 shows a flowchart of the modelling and calculation procedures required for the QMG model.To describe accurately the runoff generation and confluence on a grid scale, these karst sub-basins are further divided into many karst hydrological response units (KHRUs) based on the high-resolution (15  15 m) DEM data in the model.The specific steps involved in the division were adopted by referring to studies of hydrological response units (HRUs) in TOPMODEL by Pan (2014).As the smallest basin computing units, the KHRUs can effectively ignore the spatial differences of karst development within the units and reduce the uncertainty in the classification of model units.Figure 3 shows the spatial structure of the KHRUs.The modelling and operation of the QMG model consists of three main stages: 1) spatial interpolation, and the retention of rainfall and evaporation calculations; 2) runoff generation and confluence calculation for the surface river; and 3) confluence calculation for the underground runoff, including the confluence in the shallow karst aquifer and the underground river system.

Rainfall and evaporation calculation
In the QMG model, the spatial interpolation of rainfall is accomplished by a kriging method using the ArcGIS 10.2 software.The Tyson polygon method may be a simpler method for rainfall interpolation if the number of rainfall gauges in the basin is sufficient.The point rainfall observed by the two rainfall gauges in the basin (Fig. 1a) was interpolated spatially into areal rainfall for the entire basin.
Basin evapotranspiration in the KHRUs was mainly vegetal, soil evaporation and water surface evaporation.They were calculated using the following equations (modified from Li et al., 2020): is the wind speed 150 m above the water surface.

Runoff generation
In the QMG model, the surface runoff generation in river channel units means the rainfall in the river system after deducting evaporation losses.This portion of the runoff will participate in the confluence process directly through the river system, rather than undergoing infiltration.In contrast, the process of runoff generation in slope units is more complex, and its classification is related to the developmental characteristics of surface karst in the basin, rainfall intensity and soil moisture.For example, when the soil moisture content is already saturated, there is the potential for excess infiltration surface runoff in exposed karst slope units.The surface runoff generation of the KHRUs in the river channel units and slope units can be described by the following equations (modified from Chen, 2009Chen, , 2018;;Li et al., 2020): Here, P r (t) [mm]  In the KHRUs (Fig. 3), underground runoff is generated primarily from the infiltration of rainwater and direct confluence recharge from sinkholes or skylights.In the QMG model, the underground runoff is calculated by the following equations (modified from Chen, 2018): Here, R g [mm] is the underground runoff depth (this part of the underground runoff is mainly from the direct confluence supply of the karst sinkholes or karst windowsskylights in the study area), R 0 [mm]  then the vadose zone is not yet full, there will be no underground runoff generation, and rainfall infiltration at this time will continue to compensate for the lack of water in the vadose zone until it is full and before runoff is generated.

Channel routing and confluence
In the QMG model, the calculation of runoff confluence on the KHRUs includes the confluence of surface river channel and underground runoff.There are already many mature and classical algorithms available for calculating the runoff confluence in river channel units and slope units, such as the Saint-Venant equations and Muskingum convergence model.In this study, the Saint-Venant equations were adopted to describe the confluence in the surface river and hill slope units, for which a wave movement equation was adopted to calculate confluence in slope units (Chen, 2009): Here, we customized two variables a and b: Equation ( 7) was substituted into Eq.( 5) and discretized by a finite-difference method, giving ( 1) The Newton-Raphson method was used for the iterative calculation using Eq. ( 8): where Q [L/s] is the confluence of water flow in slope units, L [dm] is its runoff width, h [dm] is the runoff depth and q [dm 2 /s] is the lateral inflow on the KHRUs.Here, the friction slope f S equals the hill slope 0 S , and the inertia term and the pressure term in the motion equation of the Saint-Venant equations were ignored.The term v [dm/s] is the flow velocity of surface runoff in the slope units as calculated by the Manning equation, n is the roughness coefficient of the slope units, is the slope inflow in the KHRU at time t+1 and is the slope discharge in the upper adjacent KHRU at time t+1.
Similarly, the surface river channel confluence was described based on the Saint-Venant equation, where a diffusion wave movement equation was adopted, meaning that the inertia term in the motion equation was ignored: A finite-difference method and the Newton-Raphson method were used for the iterative calculation of the above equation: where Q [L/s] is the water flow in surface river channel units, A [dm 2 ] is the discharge section area, c is a custom intermediate variable and  [dm] is the wetted perimeter of the discharge section area.
The underground runoff in the model includes the confluence of the epikarst zone and underground river.In the epikarst zone, the karst water-bearing media are highly heterogeneous (Williams, 2008).For example, the anisotropiccrisscrossed karst fissure systems and conduit systems consist of large the corrosion fractures.When rainfall infiltrates into the epikarst zone, water moves slowly through the small (less than 10cm in this study) karst fissure systems, while it flows rapidly in larger (more than 10cm) conduits.
The key to determining the confluence velocity lies in the width of karst fractures.In the KHRUs (Fig. 3), the 10-cm width of the fracture was used as a threshold value (Atkinson, 1977) based on the borehole pumping test in the basin, meaning that if the fracture width exceeded 10 cm, then the water movement into it was defined as rapid flow; otherwise, it was defined as slow flow.The confluence in the epikarst zone was calculated by the following equation (modified from Beven and Binley, 2006): The distinction between rapid and slow flows in the epikarst zone is not absolute.The 10-cm width of a karst fracture as the dividing threshold is underrepresented due to the only five limited boreholes have been tested for pumping in the region.alsohas some subjectivity.
In fact, there is usually water exchange between the rapid and slow flows at the junction of large and small fissures in karst aquifers.In the QMG model, this water exchange can be described with this equation (modified form Li et al., 2021): Here, ,, [dm] is the water head difference between the rapid and slow flows at the junction of large and small fissures in KHRUs, np is the number of fissure systems connected to the adjacent conduit systems,   ,, The confluence of the underground river system plays an important role for the confluence at the basin outlet.To facilitate the calculation of confluence in the QMG model, the underground river systems can be generalized into large multiple conduit systems.
During floods, these conduit systems are mostly under pressure.Whether the water flow is laminar or turbulent depends on the flow regime at that time.The water flow into these conduits is calculated by the Hagen-Poiseuille equation and the Darcy-Weisbach equation (Shoemaker et al., 2008) Here,

Parameter optimization
In total, the QMG model has 12 parameters, of which flow direction and slope are topographic parameters that can be determined from the DEM without parametric optimization, while the remaining 10 parameters require calibration.Other distributed hydrological models with multiple structures usually have many parameters.For example, the Karst-Liuxihe model (Li et al., 2021) has 15 parameters that must be calibrated.In the QMG model, each parameter is normalized as where i x is the dimensionless parameter value i after it is normalized, * i x is the parameter value i in actual physical units, and 0 i x is the initial or final value of i x .Through the processing of Eq. ( 16), the value range of the model parameters is limited to a hypercube ignores the influence of the spatiotemporal variation of the underlying surface attributes on the parameters, while also simplifying the classification and number of the model parameters to a certain extent.Accordingly, the model parameters can be divided further into rainfall-evaporation ones, epikarst-zone ones and underground-river ones.Table 1 lists the parameters of the QMG model.
Because the QMG model has relatively few parameters, it is possible to calibrate them manually, which has the advantage that the operation is easy to implement and does not require a special program for parameter optimization.However, the disadvantage is that it is subjective, which can lead to great uncertainty in the manual parameter calibration process.
To compare the effects of parameter optimization on model performance, this study used both manual parameter calibration and the improved chaotic particle swarm optimization algorithm (IPSO) for the automatic calibration of model parameters, and compared the effects of both on flood simulation.
In general, the structure and parameters of a standard particle swarm optimization algorithm (PSO) are simple, with the initial parameter values obtained at random.For parameter optimization in high-dimensional multi-peak hydrological models, the standard PSO is easily limited to a local convergence and cannot achieve the optimal effect, while the late evolution of the algorithm may also cause problems, such as precocity and stagnant evolution, due to the 'inert' aggregation of particles, which seriously affects the efficiency of parameter selection.It is necessary to overcome the above problems and make the algorithm converge to the global optimal solution with a high probability.In parameter optimization for the QMG model, we improved the standard PSO algorithm by adding chaos theory, and developed the IPSO, where 10 cycles of chaotic disturbances were added to improve the activity of the particles.The inverse mapping equation of the chaotic variable is where X ij is the optimization variable for the model parameters, max min () XX  is the difference between its maximum and its minimum, Z ij is the variable before the disturbance is added and Z ' ij represents the chaotic variables after a disturbance is added,  is a variable determined by the adaptive algorithm, 0 ≤ α ≤ 1, and Z * is the chaotic variable formed when

Uncertainty analysis
Uncertainties in hydrological model simulation results usually originate from three aspects: input data, model structure and model parameters (Krzysztofowicz, 2014).In the present study, the input data (e.g.rainfall, flood events and some hydrogeological data) were first validated and pre-processed through observations to reduce their uncertainties.
Second, we simplified the structure of the QMG model to reduce the structural uncertainty.As a mathematical and physical model, a hydrological model has some uncertainty in flood simulation and forecasting because of the errors in system structure and the algorithm (Krzysztofowicz and Kelly, 2000).The model was designed with full consideration of the relationship between the amount of data required to build the model and its performance for flood simulation and forecasting in karst regions, and the model's entire framework was integrated through simple structures and easy-to-implement algorithms, using the concept of distributed hydrological modelling.Conventionally, the extent of uncertainty is increased with the growing complexity of the model structure.We therefore ensured that the structure of the QMG model was simple when it was designed, and the model was divided into surface and underground double-layer structures to reduce its structural uncertainty.
Third, we focus on analysing the uncertainty and sensitivity of the model parameters and their optimization method, for which a multi-parametric sensitivity analysis method (Choi et al., 1999;Li et al., 2020) was used to analyse the sensitivity of the parameters in the QMG model.The steps in the parameter sensitivity analysis are as follows.

1) Selection of appropriate objective function
The Nash-Sutcliffe coefficient is widely used as the objective function to evaluate the performance of hydrological models (Li et al., 2020(Li et al., , 2021)).It was therefore used to assess the QMG model.Because the most important factor in flood forecasting is the peak discharge, it is used in the Nash coefficient equation: where NSC is the Nash-Sutcliffe coefficient, 2) Parameter sequence sampling The Monte Carlo sampling method was used to sample 8000 groups of parameter sequences.The parametric sensitivity of the QMG model was analysed and evaluated by comparing the differences between the a priori and a posteriori distributions of the parameters.

3) Parameter sensitivity assessment
The a priori distribution of a model parameter means its probability distribution, while the a posteriori distribution refers to the conditional distribution calculated after sample sampling, and it can be calculated based on the simulation result of the parametric optimization.If there is a significant difference between the priori distribution and its posteriori distributionthemof the parameter, then the parameter being tested has a high sensitivity, whereas if there is no obvious difference, then the parameter is insensitive.The parametric priori distribution is calculated as ,, 1 ( 0.85) 100 1 where P i,j is the a priori distribution's probability when , 0.85 ij NSC  .We used a simulated Nash coefficient of 0.85 as the threshold value, and n was the number of occurrences of a Nash coefficient greater than 0.85 in flood simulations.In each simulation, only a certain parameter was changed, while the remaining parameters remained unchanged.
If the Nash coefficient of this simulation exceeded 0.85, then the flood simulation results were considered acceptable.The term i  is the difference between the acceptable value and its mean, which represents the parametric sensitivity (0 < i  < 1).The higher the i  value, the more sensitive the parameter.N is the 8000 parameter sequences, and is the average value of the a priori distribution.

Model Setting
Once the model was built, some of the initial conditions had to be set before running it to simulate and forecast floods, such as basin division, the setting of initial soil moisture, and the assumption of the initial parameter range.1) In the study area, the entire Qingmuguan karst basin was divided into 893 KHRUs, including 65 surface river units, 466 hill slope units, and 362 underground river units.The division of these units formed the basis for calculating the process of runoff generation and convergence.
2) The initial soil moisture was set to 0-100% of the saturation moisture content in the basin, and the specific soil moisture before each flood had to be determined by a trial calculation.
3) The waterhead boundary conditions of the groundwater were determined by a tracer test in the basin, where a perennial stable water level adjacent the groundwater-divide was used as the fixed waterhead boundary.The base flow of the underground river was determined to be 35 L/s from the perennial average dry season runoff.4) The range of initial parameters and convergence conditions were assumed before parameter optimization (Figure 4).5) Parameter optimization and flood simulation validated the performance of the QMG model in karst basins.

Parameter Sensitivity Results
The number of parameters in a distributed hydrological model is generally large, and it is important to perform a sensitivity analysis of each parameter to quantitatively assess the impact of the different parameters on model performance.In the QMG model, each parameter was divided into four categories according to its sensitivity: (i) highly sensitive, (ii) sensitive, (iii) moderately sensitive, and (v) insensitive.In the calibration of model parameters, insensitive ones do not need to be calibrated, which can greatly reduce the amount of calculation and improve the efficiency of model operation.
The flow process in the calibration period (14 April to 10 May 2017) was adopted to calculate the sensitivity of the model parameters, for which the calculation principle was equation ( 19), and the parameter sensitivity results are calculated in Table 2.
Table 2 Parametric sensitivity results in QMG model.
In Table 2, the value of i  [equation ( 19)] represents a parameter's sensitivity, and the higher the value, the more sensitive the parameter is.From the results in Table 2, it was found that the rainfall infiltration coefficient, rock permeability coefficient, rock porosity, and the related parameters of soil water content, such as the saturated water content, and field capacity, were sensitive parameters.The order of parameter sensitivity was as follows: infiltration coefficient > permeability coefficient > rock porosity > specific yield > saturated water content > field capacity > flow direction > thickness > slope > Soil coefficient > channel roughness > evaporation coefficient.
In the QMG model, parameters are classified as highly sensitive, sensitive, moderately sensitive, and insensitive according to their influence on the flood simulation results.In Table 4, we divided the sensitivity of model parameters into four levels based on the i  value: 1) highly sensitive parameters, 0.8 < i  < 1; 2) sensitive parameters, 0.65 < i  < 0.8; 3) moderately sensitive parameters, 0.45 < i  < 0.65; and 4) insensitive parameters, 0 < i  < 0.45.The highly sensitive parameters were the infiltration coefficient, permeability coefficient, rock porosity, and specific yield.The sensitive parameters were the saturated water content, field capacity, and thickness of the epikarst zone.The moderately sensitive parameters were flow direction, slope, and soil coefficient.The insensitive parameters were channel roughness and the evaporation coefficient.

Parametric Optimization
In total, the QMG model has 12 parameters, of which only eight need to be optimized, which is relatively few for distributed models.The parameters of flow direction and slope as well as the insensitive parameters of channel roughness and the evaporation coefficient need not be calibrated, which can improve the convergence efficiency of the model parameter optimization.
In the study area, 18 karst floods during the period of 14 April 2017 to 10 June 2019 were recorded at the underground river outlet to validate the effects of the QMG model in karst hydrological simulations.The calibration period was 14 April to 10 May 2017 at the beginning of the flow process, with the remainder of the time being the validation period.In the QMG model, the IPSO algorithm was used to optimize the model parameters.To show the necessity of parameter optimization for the distributed hydrological model, the study specifically compared the flood simulations obtained using the initial parameters of the model (without parameter calibration) and the optimized parameters.automatic parameter optimization iterations were required to reach convergence (Li et al., 2021), demonstrating the effectiveness of the IPSO algorithm.
To evaluate the effect of parameter optimization, the convergence efficiency of the algorithm, and more importantly, the parameters after calibration were used to simulate floods.parameters is necessary and that there was an improvement in parameter optimization through the use of the IPSO algorithm in this study.In addition, it was found that the flow simulation effect was better in the calibration periods than in the validation periods (Fig. 6).
To compare the results of the flow processes simulation with the initial model parameters and the optimized parameters, six evaluation indices (Nash-Sutcliffe coefficient, correlation coefficient, relative flow process error, flood peak error, water balance coefficient, and peak time error) were applied in this study, and the results are presented in Table 3.
Table 3 Flood simulation evaluation index through parametric optimization.
Table 3 shows that the evaluation indices of the flood simulations after parametric optimization were better than those of the initial model parameters.The average values of the initial parameters for these six indices 0.81, 0.74, 27%, 31%, 0.80, and 5 h, respectively.
For the optimized parameters, the average values were 0.90, 0.91, 16%, 14%, 0.94, and 3 h, respectively.The flood simulation effects after parameter optimization clearly improved, implying that parameter optimization for the QMG model is necessary, and the IPSO algorithm for parameter optimization is an effective approach that can greatly improve the convergence efficiency of parameter optimization, and also ensure that the model performs well in flood simulations.

Model Validation in Flood Simulations
Following parameter optimization, we simulated the whole flow process (14 April 2017 to 10 June 2019 ) based on the optimized and initial parameters of the QMG model (Fig. 6), which enabled a visual reflection of the model used in the simulation of a long series of flow processes.To reflect the simulation effect of the model for different flood events, we divided the whole flow process into 18 flood events, then used the initial parameters of the model and the optimized parameters, respectively, to verify the model performance in flood simulations.Figure .7 and Table 4 show the flood simulation effects and their evaluation indices using both the initial and the optimized parameters.
Figure 7 Flood simulation effects based on initial and optimized parameters.
Table 4 Flood simulation indices for model validation.4, the average water balance coefficient based on the initial parameters was 0.69, i.e., much less than 1, indicating that the simulated water in the model was unbalanced.After parameter optimization, the average value was 0.92, indicating that parameter optimization had a significant impact on the model water balance calculation.
Table 4 shows that the average values of the six indices (Nash-Sutcliffe coefficient, correlation coefficient, relative flow process error, flood peak error, water balance coefficient, and peak time error) for the initial parameters were 0.79, 0.74, 26%, 25%, 0.69, and 5 h, respectively, while for the optimized parameters the average values were 0.92, 0.90, 10%, 11%, 0.92, and 2 h, respectively.All evaluation indices improved after parameter optimization, with the average values of the Nash coefficient, correlation coefficient, and water balance coefficient increasing by 0.13, 0.16, and 0.23, respectively.The average values of the relative flow process error, flood peak error, and peak time error decreased by 15%, 14%, and 3 h, respectively.These reasonable flood simulation results confirmed that parameter optimization by the IPSO algorithm was necessary and effective for the QMG model.
Compared with the overall flow process simulation shown in Figure 6, each flood process was better simulated by the QMG model (Fig. 7).This was because in the function of the QMG model and its algorithm design, the main consideration was the calculation of the flood process, but the correlation algorithm of the dry season runoff was not described well enough.For example, equations ( 12)-( 15) are the flood convergence algorithm.As a result, the model is not good at simulating other flow processes, such as dry season runoff, leading to a low accuracy in the overall flow process.The next phase of our research will focus on refining the algorithm related to dry season runoff and improving the comprehensive performance of the model.

Assessment and reduction of uncertainty
In general, the uncertainty in model simulation is due mainly to three aspects of the model: (i) the uncertainty of its input data, (ii) the uncertainty of its structure and algorithm and (iii) the uncertainty of its parameters.In the practical application of a hydrological model, these three uncertainties are usually interwoven, which leads to the overall uncertainty of the final simulation results (Krzysztofowicz, 2014).Therefore, the present study focused on the uncertainties in the input data, the model structure and the parameters to reduce the overall uncertainty of the simulation results.
First, the input data-mainly rainfall-runoff data and hydrogeological data-were pre-processed, which substantially reduced their uncertainty.Second, we simplified the structure of the QMG model, which is reflected in the fact that it has only two layers of spatial structure in the horizontal and vertical directions.This relatively simple structure reduced greatly the uncertainty due to the model structure.In contrast, the underground structure of our previous Karst-Liuxihe model (Li et al., 2021) has five layers, which leads to great uncertainty.Third, appropriate algorithms for runoff generation and confluence were selected.Different models were designed for different purposes, which leads to great differences in the algorithms used.In the QMG model, most of the rainfall-runoff algorithms used have been validated by the research results of others, and some of them were improved to suit karst flood simulation and forecasting by the QMG model.For example, the algorithm for the generation of excess infiltration runoff [Eq.( 2)] was an improvement of the version used in the Liuxihe model (Chen, 2009(Chen, , 2018;;Li et al., 2020).Finally, the algorithm for parameter optimization was improved.Considering the shortcomings of the standard PSO algorithm that tends to converge locally, this study developed the IPSO for parameter optimization by adding chaotic perturbation factors.The flood simulation results after parameter optimization were much better than those of the initial model parameters (Figs. 6 and 7 and Tables 2 and 3), which indicates that parameter optimization is necessary for a distributed hydrological model and can reduce the uncertainty of the model parameters.

Parameter sensitivity analysis
The parameter-sensitivity results in Table 2 show that the rainfall-infiltration coefficient in the QMG model was the most sensitive parameter.It was the key to determining the generation of excess infiltration surface runoff and separating surface runoff from subsurface runoff.If the rainfall infiltration coefficient was greater than the infiltration capacity, excess infiltration surface runoff was generated on the exposed karst landforms; otherwise, all rainfall would infiltrate to meet the water deficit in the vadose zone, and then continue to seep down into the underground river system, eventually flowing out of the basin through the underground river outlet.The confluence modes of surface runoff and underground runoff were completely different, resulting in a large difference in the simulated flow results.
Therefore, the rainfall infiltration coefficient had the greatest impact on the final flood simulation results.
Other highly sensitive parameters such as the rock permeability coefficient, rock large amount of data to build models in karst areas (Kraller et al., 2014).The QMG model has only a double-layer structure, with a clear physical meaning, and a small amount of basic data is needed to build the model in karst areas, such as some necessary hydrogeological data.For example, the distribution and flow direction of underground rivers is required, which can be inferred from a tracer test, leading to a low modelling cost.There were fewer parameters in the QMG model than in other distributed hydrological models, with only 10 parameters that needed to be calibrated.
The flood simulation after parameter optimization was much better than the simulation using the initial model parameters.After parameter optimization, the average values of the Nash coefficient, correlation coefficient and water balance coefficient increased by 0.13, 0.16 and 0.23, respectively, while the average relative flow process error, flood peak error and peak time error decreased by 15%, 14% and 3 h, respectively.Parameter optimization is necessary for a distributed hydrological model, and the improvement of the IPSO algorithm in this study was an effective way to achieve this.
In the QMG model, the rainfall infiltration coefficient I c , rock permeability coefficient K, rock porosity R p and the parameters related to the soil water content were sensitive parameters.The order of parameter sensitivity was infiltration coefficient > permeability coefficient > rock porosity > specific yield > saturated water content > field capacity > flow direction > thickness > slope > soil coefficient > channel roughness > evaporation coefficient.
This QMG model is suitable for karst trough and valley landform like this study areabasins, where the topography is conducive to the spread of flood water.Whether this model is applicable to other the karst areas of other landforms in non-trough valley regions still needs to be verified in the future studies.In addition, the basin area is very small, where the hydrological similarity between different small basin areas varies greatly (Kong and Rui, 2003).The size of the area to be modelled has a great influence on the choice of model spatial resolution (Chen et al., 2017).Therefore, whether the QMG model is suitable for flood forecasting in large karst basins needs to be determined.

Figure 3 .
Figure 3. Spatial structure of karst hydrological response units (KHRUs) (Li et al., 2021).The right-hand side of Figure.3 shows a three-dimensional spatial model of KHRUs established in the laboratory to reflect visually the storage and movement of water in the karst water-bearing medium with each spatial anisotropy, and to provide technical support for establishing the hydrological model.
is the rainfall variation by vegetation interception, P v [mm] is the vegetation interception of rainfall and E s [mm] is the actual soil evaporation.The term λ is the evaporation coefficient.The term E p [mm] is the evaporation capability, which can be measured experimentally or estimated by the water surface evaporation equation E w .The term F [mm] is the actual soil moisture, F sat [mm] is the saturation moisture content, F c [mm] is the field capacity, E w [mm/d] is the evaporation of the water surface and e  = e 0 −e 150 [hPa] is the draught head between the saturation vapour pressure of the water surface and the air vapour pressure 150 m above the water surface (150 m above the water surface was selected here because the altitude for temperature and humidity observations in the southwestern karst regions of China is usually set at 150-200 m).The term T the water surface and the temperature 150 m above the water surface,  is the relative humidity 150 m above the water surface and  [m/s] s] is the flow confluence in the epikarst zone at time t,  [g/L] is the density of the water flow, g [m/s 2 ] is gravitational acceleration, n is the valid computational units,[L] is the volume of the ijk-th KHRU, v is the kinematic viscosity coefficient, f ij is the attenuation coefficient in the epikarst zone, h ij[dm]  is the depth of shallow groundwater and z ij[dm]  is the thickness of the epikarst zone.
is the permeability coefficient at the junction of a fissure and conduit, ip d and ip r [dm] are the conduit diameter and radius, respectively, ip l [dm] is the length of the connection between conduits i and p, and ip  is the conduit curvature.Some of the parameters in this equation, such as   ,, conducting an infiltration test in the study area.
s] is the water flow of the laminar flow in the conduit systems, A [dm 2 ] is the conduit cross-sectional area, d [dm] is the conduit diameter,  [kg/dm 3 ] is the density of the underground river, =/    is the coefficient of kinematic viscosity, / hl   is the hydraulic slope of the conduits,  is the dimensionless conduit curvature, turbulent Q [L/s] is the turbulent flow in the conduit systems and H c [dm] is the average conduit wall height.
the optimal particle maps to the interval [0,1].In parameter optimization, the flowchart of the IPSO is shown in Figure.4.
Figure. 5 shows the iteration process of parameter optimization for the QMG model.

Figure 5
Figure 5 Iteration process of parametric optimization.

Figure. 5
Figure. 5 shows that almost all parameters fluctuated widely at the beginning of the optimization, and then after about 15 iterations of the optimization calculation, most of the linear fluctuations become significantly less volatile, which indicated that the algorithm tended to converge (possibly only locally).When the number of iterations exceeded 25, all parameters remained essentially unchanged, meaning that the algorithm had converged (at this point there was global convergence).It took only 25 iterations to reach a definite Figure.6  shows the flood simulation effects.

Figure 6
Figure 6 Flow simulation results of QMG model based on parameter optimization.

Figure. 6
Figure.6  shows that the flows simulated by parameter optimization were better than those simulated by the initial model parameters.The simulated flow processes based on the initial parameters were relatively small, with the simulated peak flows in particular being smaller than the observed values, and there were large errors between the two values.In contrast, the simulated flows produced by the QMG model after parameter optimization were very similar to the observed values, which indicates that calibration of the model

Figure. 7
Figure.7  shows that the flood simulation results using the initial parameters were smaller than the observed values, and the model performance improved in flood simulations after parameter optimization.The simulated flood processes were in good agreement with observations, and were especially effective for simulating flood peak flows.From flood simulation indices in Table4, the average water balance coefficient based on the initial

Figure 4
Figure 4 Algorithm flow chart of the IPSO.

Figure 5
Figure 5 Iteration process of parametric optimization.

Figure 6
Figure 6 Flow simulation results of QMG model based on parameter optimization.

Figure 7
Figure 7 Flood simulation effects based on initial and optimized parameters.
is the net rainfall (deducting evaporation losses) in the river channel units at time t [h], P i (t) [mm] is the rainfall in the river channel units, L [m] is the length of the river channel, max W[m] is the maximum width of the river channel selected and A [m 2 ] is the cross-sectional area of the river channel.R si [mm] is termed the excess infiltration runoff in the QMG model, when the vadose zone is short of water and has not been filled.The infiltration capacity f max is different in different karst landform units, α, β are the parameters of the Holtan model and F s [mm] is the stable depth of soil water infiltration.
is the average depth of the underground runoff, p and m are attenuation coefficients calculated by conducting a tracer test in the study area, R e [L/s] is the underground runoff generated from rainfall infiltration in the epikarst zone, I w [mm] is the width of the underground runoff on the KHRUs, z [mm] is the thickness of the epikarst zone, R r [mm 2 /s] is the runoff recharge on the KHRUs during period t, R epi [mm 2 /s] is the water infiltration from rainfall, e v [mm/s] is the flow velocity of the underground runoff, K [mm/s] is the current permeability coefficient and  is the hydraulic gradient of the underground runoff.If the current soil moisture is less than the field capacity, i.e.
c FF  ,

Table 1
Parameters of the QMG model.

Table 2
Parametric sensitivity results in QMG model.

Table 3
Flood simulation evaluation index through parametric optimization.

Table 4
Flood simulation indices for model validation.