Simulation of Groundwater Flow in Fractured-Karst Aquifer with a Coupled Model in Maling Reservoir, China

: A coupled model has been developed to simulate groundwater ﬂow in fractured karst systems according to the complex geological and karst hydrogeological conditions of the dam site, where a 3D mathematical model based on Boussinesq equation was used to describe the movement of groundwater ﬂow in fractured medium, and a 1D conduit model for karst medium. The model was solved with the continuous hydraulic heads at the common boundaries. The hydraulic conductivities of karst medium were determined by geometrical parameters and ﬂux of pipes. Furthermore, the permeability parameters for fractured medium were calibrated by the measured and calculated groundwater levels. The calibrated model was employed to predict the variation of groundwater ﬂow ﬁeld and leakage from the karst pipes and underground powerhouse during the reservoir operation. The simulated results showed that the groundwater level of the powerhouse had decreased by about 2–5 m. The water level of conveyance pipeline had risen by 10–20 m, and the water level on both banks had risen by 15–25 m. The leakage of karst conduits for impervious failure was larger than that for normal seepage control. In addition, the leakage of the powerhouse was estimated to be about 1000–3000 m 3 /d, and the seepage control of karst pipes had little inﬂuence on the leakage of underground powerhouse.


Introduction
Many water conservancy and hydropower projects are located in fractured karst areas. The existence of karst pipes is a serious threat to the normal storage of reservoirs. The fractured karst media are usually divided into three basic types: pores and micro-fissures in rocks matrix, mesoscale fissures, and karst pipes and faults [1][2][3][4][5]. Pores, micro-fissures, and mesoscale fissures often account for the main part of the porosity of karst media, but their permeabilities are far less than those of karst pipes and faults. Therefore, these pores or fractures play a major role in water storage, where groundwater flow follows Darcy's law. However, karst pipes play a role in water conduction in the movement of groundwater; the groundwater flow follows non-Darcy's law or non-linear flow owing to faster flow velocity [6][7][8][9].
Due to the complexity of groundwater flow in fractured karst aquifer, the numerical simulation method is an effective tool for studying the process of groundwater flow [10][11][12].
In the past few decades, a variety of numerical simulation methods for studying groundwater flow have been developed. Many experts and scholars have explored models that simulate groundwater flow in karst media [13][14][15][16]. Wu and Malenica et al. [17,18] developed a generalized discrete-continuum model to simulate ground water flow in the karst aquifers, which took into account both quick conduit flow and diffusive fissure flow. Bauer et al. [19,20] built a hybrid continuum-discrete pipe flow model (CAVE) to study the effects of the coupling of two flow systems on the types and duration of early

Location
Maling Water Conservancy Project is located in the middle reaches of Mabie River in Xingyi City, Guizhou Province. The river is a large tributary of the Xijiang River system on the left bank of Nanpan River in the Pearl River Basin (Figure 1).
Its length is 142.5 km with a drop of 1462 m and drainage area of 2886 km 2 . This project is 3 km from Maling Town downstream, 16 km from Xingyi City, and 318 km from Guiyang City. The reservoir is a kind of karst canyon. The elevation of the dam is 90 m with a total reservoir capacity of 1.282 × 10 8 m 3 and regulating storage capacity of 1.072 × 10 8 m 3 . The project is located in the south of Yunnan-Guizhou Plateau and its topography is high in the north and low in the south. The highest elevation in the north is about 1600 m. It is gradually descending to the height of 1100-1200 m in the dam site. The main tasks of the project are a comprehensive utilization for water supply and irrigation in urban and rural areas, and power generation. The yearly average temperature is 16 Its length is 142.5 km with a drop of 1462 m and drainage area of 2886 km². This project is 3 km from Maling Town downstream, 16 km from Xingyi City, and 318 km from Guiyang City. The reservoir is a kind of karst canyon. The elevation of the dam is 90 m with a total reservoir capacity of 1.282 × 10 8 m 3 and regulating storage capacity of 1.072 × 10 8 m 3 . The project is located in the south of Yunnan-Guizhou Plateau and its topography is high in the north and low in the south. The highest elevation in the north is about 1600 m. It is gradually descending to the height of 1100-1200 m in the dam site. The main tasks of the project are a comprehensive utilization for water supply and irrigation in urban and rural areas, and power generation. The yearly average temperature is 16.3 °C with a maximum temperature of 36.5 °C and a minimum temperature of −4.7 °C. The annual average rainfall and evaporation are 1437.6 mm and 1512.6 mm, respectively.

Hydrogeological Conditions
There are mainly three layers of medium-thick limestone in the dam site with the stratigraphic symbols of T 2 g 3-2-1 , T 2 g 3-2-2 , and T 2 g 3-2-3 ( Figure 2). The total thickness is 103.4 m. The underlying strata are medium-thick limestone interbedded with thinbedded shale (T 2 g 3-1-2 ). In addition, the overlying strata are medium-thick limestone with mudstone shale (T 2 g 3-3 ) with a thickness of 33.6 m and thick dolomite with dolomitic limestone (T 2 g 3-4 ).
(1) Permeability analysis of rock mass The permeability of rock mass is closely related to its weathering integrity. Rock mass in weathering and unloading zones around the valley has good permeability. It decreases obviously with the increase of rock mass integrity ( Table 1). The Lu values in Table 1 are obtained by water pressure test, which is defined as the inflow rate per meter in unit time, when the maximum pressure of the test section is 1 MPa. The test was proposed by M.L. Ngeon from France in 1933 to estimate the necessity of grouting in rock mass of dam foundation.  (1) Permeability analysis of rock mass The permeability of rock mass is closely related to its weathering integrity. Rock mass in weathering and unloading zones around the valley has good permeability. It decreases obviously with the increase of rock mass integrity ( Table 1). The Lu values in Table 1 are obtained by water pressure test, which is defined as the inflow rate per meter in unit time, when the maximum pressure of the test section is 1 MPa. The test was proposed by M.L. Ngeon from France in 1933 to estimate the necessity of grouting in rock mass of dam foundation. Because of the different lithology in the dam site, the development of karst is different, which constitutes a multi-level aquifer system. For example, the T2g 2 layer belongs to relative aquiclude. The T2g 3-1-1 , T2g 3-1-2 , and T2g 3-3 layers are the fractured weak permeable layers. The T2g 3-2-1 , T2g 3-2-2 , and T2g 3-2-3 are the moderate permeable layers in karst-fractured systems. In addition, T2g 3-4 is the strong karst permeable layer. Therefore, the permeability of rock mass has heterogeneous anisotropy and is controlled by karst development, weathering, and fragmentation degree of rock mass. The zones with high permeability are mainly distributed in three geological bodies: (1) Strong and fractured dissolution weathering zone near the surface; these zones were disclosed in the adits PD-4, PD-12, and PD-13 in the left bank and PD-9 and PD-10 in the right bank (We edited the adits found in the study area as PD-1, PD-2..., PD-N). Serious dripping occurred in these zones. (2) Karst Because of the different lithology in the dam site, the development of karst is different, which constitutes a multi-level aquifer system. For example, the T 2 g 2 layer belongs to relative aquiclude. The T 2 g 3-1-1 , T 2 g 3-1-2 , and T 2 g 3-3 layers are the fractured weak permeable layers. The T 2 g 3-2-1 , T 2 g 3-2-2 , and T 2 g 3-2-3 are the moderate permeable layers in karstfractured systems. In addition, T 2 g 3-4 is the strong karst permeable layer. Therefore, the permeability of rock mass has heterogeneous anisotropy and is controlled by karst development, weathering, and fragmentation degree of rock mass. The zones with high permeability are mainly distributed in three geological bodies: (1) Strong and fractured dissolution weathering zone near the surface; these zones were disclosed in the adits PD-4, PD-12, and PD-13 in the left bank and PD-9 and PD-10 in the right bank (We edited the adits found in the study area as PD-1, PD-2..., PD-N). Serious dripping occurred in these zones.
(2) Karst pipe and corrosion fracture zones; serious dripping was found in the section of 53. Owing to the influences of weathering and unloading, the permeability of rock mass decreases obviously with the increase of burial depth of rock mass (Table 2). At the depth of 0-60 m, the permeability of local rock mass is moderate. About 120 m below the surface, the rock mass is basically impermeable (Lu < 1). (2) Types of aquifers The groundwater types in the study area are mainly pore water in loose accumulation layer, fractured water, and fractured karst groundwater. Pore water mainly distributes in shallow surfaces and gullies. The permeability of fractured and karst groundwater is different due to the lithological difference of each aquifer. Karst is widely distributed in carbonate area. Waterfalls and sinkholes can be seen everywhere. Therefore, there are many springs in reservoir areas, and groundwater resources are abundant. However, in mudstone and shale areas, springs are less exposed and karst is not developed. The permeability of each formation is shown in Table 3.  The groundwater level in the dam site is greatly affected by rainfall. If the rainfall is large, the groundwater level is high, and the level ranges from 3 m to 20 m. According to the relationship between the groundwater level of boreholes and the river level, the groundwater hydraulic gradient is calculated. On the left bank, hydraulic gradient is 0.33 for borehole ZK-l0, 0.16 for borehole ZK-23, and 0.18 for borehole ZK-28. On the right bank, hydraulic gradient is 0.47 for borehole ZK-l6, 0.33 for borehole ZK-22, and 0.27 for borehole ZK-27. So, the gradient is smaller on the valley slopes, but larger on the bank slopes. In the riverbed, the levels are influenced by deep circulation conditions of groundwater, and confined water appears in many places. For example, groundwater tables are 982.75 m for ZK-18, 954.96 m for ZK-19, 956.15 m for ZK-20, and 982.45 m for ZK-24. They are 2.96-30.75 m above the river level. Therefore, the groundwater levels on both banks are higher than that of the river, and groundwater recharges to the river water. Field investigation shows that most springs have larger discharge in rainy season and smaller discharge in dry season. It shows that groundwater is mainly recharged by rainfall. According to the groundwater levels on the both banks, there exists groundwater watershed in the dam site. The outcrop of all springs is higher than the river level. Therefore, the groundwater is discharged to the river in the study area under the natural conditions ( Figure 3).

(4) Hydrogeological characteristics of karst
The karst development is controlled by lithology and structure in the study area. In all pure dolomite and limestone, karst development is relatively strong. For example, sinkholes K l7 , K 18 , K 19 , and Longdang karst pipes on the left bank are mainly developed in thick massive limestone T 2 g [3][4][5][6] . However, the karst development is weak in the heavier argillaceous layers of T 2 g 3-3 and T 2 g 3-2-2 . There is only a small amount of micro-dissolution. Also, the development directions of Nos. 1 and 5 karst pipes on the left bank and Nos. 2 and 4 karst pipes on the right bank in dam site are controlled by strata. The development of karst pipes is shown in Table 4. Furthermore, the elevation of Karst development has obvious zonation. The sinkholes K l7 , K 18 , and K 19 are developed at 1100-1120 m. The springs g1 and g11 are mainly distributed at 960-980 m on the left bank. On the right bank, the karst development is concentrated at three elevations of 1050-1070 m, 1010-1020 m, and 960-980 m.
for borehole ZK-l0, 0.16 for borehole ZK-23, and 0.18 for borehole ZK-28. On the right bank, hydraulic gradient is 0.47 for borehole ZK-l6, 0.33 for borehole ZK-22, and 0.27 for borehole ZK-27. So, the gradient is smaller on the valley slopes, but larger on the bank slopes. In the riverbed, the levels are influenced by deep circulation conditions of groundwater, and confined water appears in many places. For example, groundwater tables are 982.75 m for ZK-18, 954.96 m for ZK-19, 956.15 m for ZK-20, and 982.45 m for ZK-24. They are 2.96-30.75 m above the river level. Therefore, the groundwater levels on both banks are higher than that of the river, and groundwater recharges to the river water. Field investigation shows that most springs have larger discharge in rainy season and smaller discharge in dry season. It shows that groundwater is mainly recharged by rainfall. According to the groundwater levels on the both banks, there exists groundwater watershed in the dam site. The outcrop of all springs is higher than the river level. Therefore, the groundwater is discharged to the river in the study area under the natural conditions (Figure 3).   Because of the particularity of karst development, each layer of karst is connected with each other. The lower karst develops along the upper karst. The outlet K 4 of karst pipe No. 2 has been suspended 60 m above the water level of the present river. With the down cutting of valleys and the decrease of groundwater level, karst also develops to the depth. Karst caves and No. 3 Karst Pipe disclosed in PD-3 Drift have been a part of the lower passage of No. 2 karst pipe. Karst development is weak in the deep riverbed. According to the borehole data of the riverbed in the dam site, there are few corrosion phenomena within 100 m. In addition, the permeability of rock mass is also very small. Therefore, it can be inferred that the relative lower limit of karst development is 30-50 m below the riverbed. In the two banks and watershed areas, the relative lower limit is controlled by lithology and tectonics. These karst pipes are the main channels for groundwater discharge, and their hydrodynamic conditions are complex. The calculation shows that the hydraulic gradients are about 8% for of Nos. 1 and 4 karst pipes, 4% for Nos. 2, and 20-28% for Nos. 3 and 5 (Table 4). Due to the wide distribution of karst strata at both banks of the dam site, karst is relatively developed. When the reservoir is impounded to a normal water level, these parts between the Nizao gully and the dam site have the possibility of leakage around the dam foundation and abutment by the karst fractures.

Methods
In fractured karst systems, a 3D mathematical model is employed to describe the movement of groundwater flow in fractured medium and a conduit model is used to describe it in karst medium. The two models are coupled using the conditions for continuous groundwater flow between the walls of karst conduit and rock mass matrix around the conduit.

Mathematical Model of Groundwater Flow
The transient groundwater flow in heterogeneous anisotropic fractured medium can be simulated using the Boussinesq equation [42]: where x, y, and z are the space coordinates [L]; t is time [T]; µ is the specific yield for unconfined aquifer and specific storage for confined aquifer; h f is the hydraulic head in the fractured system [L]; K fx , K fy , K fz are the hydraulic conductivities along coordinate axes [L/T]; ε is the evaporation and precipitation recharge of phreatic surface [L 3 /T]; h f0 is the initial groundwater level [L]; h f1 is the boundary groundwater level [L]; n is the normal direction of the boundary surface; K fn is the hydraulic conductivity in the normal direction of the boundary surface [L/T]; q f (x, y, z, t) is the unit area flux of the second type boundary [L 3 /T·L 2 ], it is positive for groundwater inflow and negative for outflow; Ω is the domain of the fractured medium, and Γ 1 and Γ 2 are domains of the first and second type boundary. (according to Abbreviations).

Karst Pipes Model and Determination of Its Hydraulic Conductivity
Groundwater flow in karst pipes is generally non-Darcy flow owing to the fast water velocity. The hydraulic head loss, h l , can be given by Darcy-Weisbach equation: where λ is head loss coefficient along the conduit, L is the length of the tube [L], d is the conduit diameter [L], u is the average velocity [L/T], and g is the gravitational acceleration [L/T 2 ]. If the porosity (n) of the conduit is assumed to be 1, the seepage velocity, V, is equal to the average velocity. In addition, the hydraulic gradient (J) can be expressed by head loss: (according to Abbreviations).
Substituting Equation (3) and the seepage velocity (V) into Equation (2), it can be rewritten as: Therefore, Equation (4) can be written in the form of Darcy's Law: A converted hydraulic conductivity, K c , is defined as: In order to calculate K c , the key is to calculate the corresponding λ according to the groundwater flow pattern. The λ values depend on Reynolds number (Re) and relative roughness (∆) of the tube walls. λ values can be given as [41]: For laminar flow (Re < 2320), For turbulent flow (Re ≥ 2320), The seepage velocity can be calculated according to the flux and diameter of the karst conduit. Then Re and λ values are obtained by the seepage velocity and Equations (7) and (8). Finally, the converted hydraulic conductivity from Equation (6) can be determined.
The K c is the main hydraulic conductivity along the direction of the karst pipe. In fact, there is a certain angle between the direction of the pipe and the coordinate axis. Therefore, the K c value must be converted to tensor form in the global coordinate system, that is, geodetic coordinate system, of the computational domain. A local coordinate system, ox y z , is established to parallel the positive direction of the ox -axis with the axial direction of karst pipe. The permeabilities of karst conduit are mainly along the ox'-axis owing to its strong anisotropy. So, under the local coordinate system, ox y z , the hydraulic conductivity tensors of karst pipe, [K] e L , can be expressed as [41]: Rotating the direction of the local coordinate system with the same as that of the global coordinate system, based on the matrix correlation theory and conception, the hydraulic conductivity tensors of the global coordinate system, [K] e G , can be calculated with that of the local coordinate system: where [P] is the orthogonal matrix, and it can be calculated using the local coordinates of karst pipe. Assuming    α 1 = cos α, cos β, cos γ α 2 = cos α , cos β , cos γ α 3 = cos α , cos β , cos γ And assuming that the starting and ending coordinates of the axis of karst pipe are M 1 (x 1 , y 1 , z 1 ), M 2 (x 2 , y 2 , z 2 ), respectively. The directional cosine of vector (α 1 ) or ox'-axis can be given as: Because the ox'-axis, oy'-axis, and oz'-axis are orthogonal to each other, α 1 , α 2 , and α 3 are also orthogonal, that is Therefore, vector α 2 and α 3 satisfy as: We can get the following algebraic equation: Assuming a = cosα, b = cosβ, and c = cosγ, if a is not equal to zero, the orthogonal matrix can be expressed as: The expressions of orthogonal matrices for b = 0 and c = 0 can also be obtained by the same method. Therefore, the hydraulic conductivity tensors in global coordinates can be determined by Equation (10). Finally, the principal hydraulic conductivities along three coordinate axes are obtained in karst conduit system.

Coupled Conditions
For fractured karst systems, the 3D continuum and 1D conduit models are applied to describing the fractured and karst media, respectively. The coupling model is solved using the continue conditions of the hydraulic heads at the common boundaries. The coupling conditions can be expressed as where Γ c is the common boundary of the karst and fracture media, and h k is the hydraulic head in the karst system [L]. Therefore, before the hydraulic heads are calculated in the study area, those in the common boundary must be solved. The iteration method is employed to solve Equations (1) and (17). (according to Abbreviations).

Parameter Inversion Method
In order to obtain the size and direction of the hydraulic conductivities in the study area, the inversion method of discontinuity control is used to determine those of riverbed and left and right bank rocks [21]. The objective function is established by the least square method.
where K i j are the permeability parameters to be calculated, the superscript i is the i-th sub-partition according to permeability of rock mass, and i = I, II, . . . , NNO, NNO is the total number of parametric partitions. The subscript j is the j-th parameter in the i-th sub-partition, and j = 1, 2, . . . , NK, NK is the total number of parameters, ω k is the k-th weight function of observed hydraulic heads, M is the number of observation holes, and h c n and h 0 n are the calculated and observed hydraulic heads.

Calculation Method of Hydraulic Conductivity in Fracture-Karst Media
The values mainly depend on the permeability of fractures in rock mass, spatial distribution characteristics of karst, and fracture discontinuity. Hydraulic conductivity tensor can be used to indicate the size and direction of permeability of anisotropic media. It is determined by the parameter inversion method. The equivalent permeability coefficient required by this method can be determined with the data of water pressure test in the boreholes. The Boussinesq equation was employed to calculate the equivalent permeability coefficient. When the test section for the water pressure test is close to the impermeable layer, the equation can be expressed as: If the test section is far from the impermeable layer, the equation is given as: where ω is the permeable rate, the value is equal to 100 Lu, L is the length of the test section, and r 0 is the radius of the test hole.

Field Tests and Parameter Determination
(1) Permeability coefficient of the unsaturated zone The test is to raise the water head artificially by injecting water into the test pit. It is an in-situ test method for measuring the permeability of loose rock and soil. The permeability coefficient of the unsaturated zone is measured by a water injection test of double ring. Two permeability tests were conducted in this study with the numbers of T1 and T2. The coordinates of T1 are 104 • 54 36.51" E and 25 • 12 12.59" N with the elevation of 1089.90 m. In addition, the coordinates of T2 are 104 • 55 16.01" E and 25 • 11 18.07" N with elevation of 1026.7 m. The data of tests are shown in Figure 4. The final steady infiltration rate at the end of T1 was 4.95 cm 3 /s. The area of the ring was 490.625 cm 2 . The water level in the pit is 10 cm, and infiltration depth of groundwater was 40 cm. Hence, the calculated permeability coefficient of the unsaturated zone was 0.0081 cm/s. The infiltration rate for T2 was 14.1 cm 3 /s and the infiltration depth was 35 cm. The calculated permeability coefficient was 0.022 cm/s. The average value of 0.015 cm/s was regarded as the permeability coefficient of the unsaturated zone.
(2) Division of Permeability coefficient zone According to the characteristics of lithology and discontinuity in the dam site, the media were divided into three zones in the vertical direction. They are strong, weak, and slight weathering zones from top to bottom, respectively. The media were also influenced by river erosion. There existed weathering, unloading, and fracture dissolution in a certain range of valley slopes on both banks. The adit data showed that the integrity of rocks is better in the mountain body than near the valley. Therefore, the mountain bodies of both banks are divided into three zones longitudinally. They are strong, weak, and slightly permeable zones, respectively ( Figure 5). Therefore, there are nine different material zones. The equivalent permeability coefficients are listed in Table 5. ring. Two permeability tests were conducted in this study with the numbers of T1 and T2. The coordinates of T1 are 104°54′36.51″ E and 25°12′12.59″ N with the elevation of 1089.90 m. In addition, the coordinates of T2 are 104°55′16.01″ E and 25°11′18.07″ N with elevation of 1026.7 m. The data of tests are shown in Figure 4. The final steady infiltration rate at the end of T1 was 4.95 cm 3 /s. The area of the ring was 490.625 cm 2 . The water level in the pit is 10 cm, and infiltration depth of groundwater was 40 cm. Hence, the calculated permeability coefficient of the unsaturated zone was 0.0081 cm/s. The infiltration rate for T2 was 14.1 cm 3 /s and the infiltration depth was 35 cm. The calculated permeability coefficient was 0.022 cm/s. The average value of 0.015 cm/s was regarded as the permeability coefficient of the unsaturated zone.  (2) Division of Permeability coefficient zone According to the characteristics of lithology and discontinuity in the dam site, the media were divided into three zones in the vertical direction. They are strong, weak, and slight weathering zones from top to bottom, respectively. The media were also influenced by river erosion. There existed weathering, unloading, and fracture dissolution in a certain range of valley slopes on both banks. The adit data showed that the integrity of rocks is better in the mountain body than near the valley. Therefore, the mountain bodies of both banks are divided into three zones longitudinally. They are strong, weak, and slightly permeable zones, respectively ( Figure 5). Therefore, there are nine different material zones. The equivalent permeability coefficients are listed in Table 5.   (3) Calculation of hydraulic conductivity for karst conduit Five karst pipes were considered in the study area ( Figure 5). According to the measuring data from Table 4 and Equations (6) and (10) Figure 5. Different permeability zones in the study area.  Five karst pipes were considered in the study area ( Figure 5). According to the measuring data from Table 4 and Equations (6) and (10), the calculation results of the main hydraulic conductivities for karst conduits are listed in Table 6.

Hydraulic Conductivities of Fracture-Karst Media by Inversion Analysis
Taking the center of the arch dam as the coordinate origin, the east was the x-axis, the north was the y-axis, and the vertical direction was the z-axis. The model area extended about 600 m upstream, 500 m downstream to Longdang River, and 700 m to watershed on both banks ( Figure 5). It can be seen from Figure 3, groundwater was discharged into rivers under natural conditions. The groundwater flow direction was almost perpendicular to the river, so it could be treated as a streamline boundary. The watersheds in the two banks were also treated as the second boundary conditions. The Mabie River and Longdang River in the left bank was considered as the first boundary condition. The model area was divided into 2,339,812 nodes and 4,334,022 elements. The discrete graph of the model area is shown in Figure 6.

Hydraulic Conductivities of Fracture-Karst Media by Inversion Analysis
Taking the center of the arch dam as the coordinate origin, the east was the x-axis, the north was the y-axis, and the vertical direction was the z-axis. The model area extended about 600 m upstream, 500 m downstream to Longdang River, and 700 m to watershed on both banks ( Figure 5). It can be seen from Figure 3, groundwater was discharged into rivers under natural conditions. The groundwater flow direction was almost perpendicular to the river, so it could be treated as a streamline boundary. The watersheds in the two banks were also treated as the second boundary conditions. The Mabie River and Longdang River in the left bank was considered as the first boundary condition. The model area was divided into 2,339,812 nodes and 4,334,022 elements. The discrete graph of the model area is shown in Figure 6. Assuming that the permeability coefficients of five karst pipes remained unchanged (Table 6), the permeability coefficient of the impervious curtain was set to 0.00864 m/d. The groundwater levels of observation holes ZK1, ZK2, ZK3, ZK5, and ZM10 were selected as the measured values. The calculated groundwater levels of these boreholes can be obtained by solving the coupled model using the iterative method. Permeability coefficients were adjusted continuously. When the errors between calculated and measured groundwater levels were small, the parameters at this time could be considered as accurate parameters (Table 7). Fitting curves between calculated and measured values are shown in Figure 7. It can be seen from Figure 7 that the calculated water levels fitted well with the measured water level. Generally, the calibrated model can better capture the hydrogeological characteristics of karst and fractured medium in the study area (Figure 8). Assuming that the permeability coefficients of five karst pipes remained unchanged (Table 6), the permeability coefficient of the impervious curtain was set to 0.00864 m/d. The groundwater levels of observation holes ZK1, ZK2, ZK3, ZK5, and ZM10 were selected as the measured values. The calculated groundwater levels of these boreholes can be obtained by solving the coupled model using the iterative method. Permeability coefficients were adjusted continuously. When the errors between calculated and measured groundwater levels were small, the parameters at this time could be considered as accurate parameters (Table 7). Fitting curves between calculated and measured values are shown in Figure 7.
It can be seen from Figure 7 that the calculated water levels fitted well with the measured water level. Generally, the calibrated model can better capture the hydrogeological characteristics of karst and fractured medium in the study area ( Figure 8).

Groundwater Level Analysis of Dam Site during the Reservoir Operation
During the reservoir operation, the seepage control and drainage measures for underground buildings were used in the dam site. Drainage holes have been installed in an underground powerhouse. Impervious curtain was arranged in the dam foundation and both banks and extended to the position of karst pipes. It mainly discusses the variable of

Groundwater Level Analysis of Dam Site during the Reservoir Operation
During the reservoir operation, the seepage control and drainage measures for underground buildings were used in the dam site. Drainage holes have been installed in an underground powerhouse. Impervious curtain was arranged in the dam foundation and both banks and extended to the position of karst pipes. It mainly discusses the variable of groundwater flow field of karst pipeline impervious curtain under normal and invalid conditions. The normal impounded level of the reservoir is 1030 m. In addition, the water level behind the dam is 965.47 m. Contour maps of groundwater level for underground powerhouse and karst pipes are shown in Figure 9. It can be seen that contour shape, trend, and distribution of groundwater level accurately reflected the characteristics of seepage control, drainage measures, and boundary conditions in corresponding areas (Figure 9a). Owing to the action of anti-seepage treatment and drainage measures, the seepage control effect was obvious. The general trend of groundwater level was convergence from upstream and downstream to an underground powerhouse area. The groundwater level fell obviously near it. It can be seen from Figure 9b,c that the groundwater level decreased obviously after passing through the curtain and drainage of the underground powerhouse. Where there was a curtain, the water level isoline was relatively close, which reflected the water-blocking effect of the impervious curtain. There occurred regular bending of the groundwater level in the karst pipe area. In addition, the bending place was the direction of the karst pipe. The groundwater level in the pipe decreased obviously after passing through the impervious curtain.
Groundwater in the study area was mainly discharged from the mountain bodies on both banks to the Mabie River. However, after reservoir impoundment, groundwater was mainly discharged from the upper reservoir to the lower reaches of the dam and the powerhouse through the two banks and dam foundation. The groundwater level of the powerhouse has decreased by about 2-5 m. Water conveyance pipeline has risen by 10-20 m and 15-25 m for mountainous bodies on both banks. Therefore, the natural flow field of groundwater has been greatly changed by the water conservancy project.  Groundwater in the study area was mainly discharged from the mountain bodies on both banks to the Mabie River. However, after reservoir impoundment, groundwater was mainly discharged from the upper reservoir to the lower reaches of the dam and the powerhouse through the two banks and dam foundation. The groundwater level of the powerhouse has decreased by about 2-5 m. Water conveyance pipeline has risen by 10-20 m and 15-25 m for mountainous bodies on both banks. Therefore, the natural flow field of groundwater has been greatly changed by the water conservancy project.

Leakage Analysis of Underground Powerhouse and Karst Pipes
There are five karst pipes in the dam site, of which Nos. 2 and 3 are finally merged into a karst pipe. Karst pipes 2, 3, and 4 are distributed on the right bank, and Nos. 1 and 5 on the left bank ( Figure 5). The water flow in the pipes was discharged into the reservoir area under natural conditions. However, the water level of the reservoir was higher than the elevation of pipes after reservoir impoundment. The reservoir water will be discharged into the pipes and flowed out of the reservoir owing to the difference of water levels, which will affect the reservoir impoundment. The leakage of karst pipes and the underground powerhouse was calculated under the normal and failure impervious curtain ( Table 8). The reservoir level of cases 1 and 2 is 1030 m, and that of Case 3 and 4 is 985 m. It can be seen from the Table 8 that the leakage in the case of impervious failure (Case 2) was larger than that of normal seepage control (Case 1), which shows that plugging karst pipes was very necessary. Because the vertical curtain did not reach the elevation of No. 5 karst pipe, there was little difference in leakage. It is suggested that the impervious curtain should be deepened at No. 5 karst pipe. Whether the seepage control of karst pipes is normal or not has little effect on the leakage of the underground powerhouse. During the operation of the reservoir, the leakage of the powerhouse was estimated to be about 1000-3000 m 3 /d.

Conclusions
According to the hydrogeological and karst hydrogeological conditions, a model coupling karst-fracture media and karst pipes has been developed. A 3D mathematical model of heterogeneous anisotropy was used to describe groundwater flow in the rock mass matrix, and the 1D conduit model was built based on non-Darcy's flow for karst pipes. The permeability of karst pipes was larger along the development direction of pipes. Basically, the permeability coefficients in the other two directions can be neglected. However, there was a certain angle between the coordinate systems and these pipes, Therefore, for 1D pipe flow, its permeability coefficients were a second-order tensor. The initial values for hydraulic conductivities of fracture-karst media can be obtained by the data of water pressure tests in the boreholes. The relatively accurate values were determined by the inversion method. The calibrated model was used to predict the variation of groundwater flow field and leakage in the dam site area during the reservoir operation. The results show that the coupled model can capture the hydrogeological characteristics of karst and fractured medium and reflect the features of seepage control and drainage measures in the study area. The groundwater level of powerhouse has decreased by about 2-5 m. That of the water conveyance pipeline has risen by 10-20 m and 15-25 m for both banks during the reservoir operation. The leakage of karst pipes for impervious failure was larger than that for normal seepage control. In addition, the leakage of the powerhouse was estimated to be about 1000-3000 m 3 /d.  is the length of test section; r 0 is the radius of test hole.