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

Abstract

A coupled model has been developed to simulate groundwater flow 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 flow 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 flux 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 flow field 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³/d, and the seepage control of karst pipes had little influence on the leakage of underground powerhouse.

1. 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 karstification for different boundary conditions. Huang et al. [21,22] developed a coupling model based on the non-overlapping domain decomposition method to simulate the groundwater flow in fractured bedrock areas. Kaufmann et al. [23] used the finite element method to study the flow evolution in a fractured, porous karst aquifer on a two-dimensional mesh of irregularly spaced nodal points. Chang et al. [24] developed a coupled groundwater model CE (CFPv2 + ERCH), consisting of a simple reservoir model generating the recharge source function, assuming the concept of an epikarst horizon (ERCH) with CFPv2. It was a discrete conduit-continuum groundwater flow model of the phreatic karst aquifer. Among these numerical modeling methods, the coupling model is the most commonly used [25]. The coupled-continuum pipe-flow (CCPF) model is a dual flow system consisting of a matrix representing the bulk mass of permeable limestone and a conduit system representing the karst conduit network. The flow exchange between the two systems is controlled by differences of hydraulic heads, the hydraulic conductivity, and the geometric setting [26,27]. For the karst pipe, its hydraulic conductivity plays an important role in simulating groundwater flow. How to accurately describe the diversion capacity of groundwater flow is the focus and difficulty of many experts and scholars. Warren and Root [28]; Narasimhan [29]; Lei [30]; and de Rooij [31] coupled the fractured medium and the karst pipes by the exchange flow rate between two flow systems. The exchange flow rate was proportional to the difference of hydraulic heads in the fractured and pipe systems. Chen and Teutsch et al. [32,33] investigated the validity of the CCPF model for flow in a karst aquifer. Giese et al. [34,35,36,37] applied large-scale tests to solve the problem of hydraulic parameter evaluation scale caused by significant differences in hydraulic parameters of karst aquifers. However, most of these research methods underestimated the permeability of karst pipelines and did not consider the hydrodynamic process in the pipeline clearly. The description of water flow characteristics was relatively vague. Aiming at this phenomenon, Chen et al. [38] proposed the concept of equivalent hydraulic conductivity. In addition, the expression of equivalent hydraulic conductivity (EHC) was derived for the different flow conditions of non-Darcy flow. The non-Darcy flow was coupled into the model of karstic triple void media.

It is worth noting that these coupling models mainly regard karst pipelines as one-dimensional models. Actually, the strikes of these pipes may not align with the global coordinate system [39,40,41]. Therefore, the main purpose of this paper is to derive the expression of the hydraulic conductivity tensors for karst pipes based on the equivalent hydraulic conductivity from Chen [38]. Then, the model of karst pipe with the hydraulic conductivity tensors was coupled into the continuous media model that consisted of pores, micro-fissures, and mesoscale fissures. The calibrated coupling model was applied to simulate the groundwater flow in fractured karst systems.

2. Studied Site

2.1. 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². 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⁸ m³ and regulating storage capacity of 1.072 × 10⁸ m³. 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.

2.2. Hydrogeological Conditions

There are mainly three layers of medium-thick limestone in the dam site with the stratigraphic symbols of T₂g^3-2-1, T₂g^3-2-2, and T₂g^3-2-3 (Figure 2). The total thickness is 103.4 m. The underlying strata are medium-thick limestone interbedded with thin-bedded shale (T₂g^3-1-2). In addition, the overlying strata are medium-thick limestone with mudstone shale (T₂g^3-3) with a thickness of 33.6 m and thick dolomite with dolomitic limestone (T₂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. Percentage of Lu values for different weathering zones (%).

Weathered Zones	5–10 Lu	1–5 Lu	<1 Lu
Weak weathered zone	28	4	6
Slightly weathered zone	44	13	6
Intact rock mass	28	83	88
Weak weathered zone	28	4	6

Lu is equivalent to 10⁻⁵ cm/s.

). 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 T₂g² layer belongs to relative aquiclude. The T₂g^3-1-1, T₂g^3-1-2, and T₂g^3-3 layers are the fractured weak permeable layers. The T₂g^3-2-1, T₂g^3-2-2, and T₂g^3-2-3 are the moderate permeable layers in karst-fractured systems. In addition, T₂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.2–55.3 m in the adit PD-10 and 53.3–57.5 m in the adit PD-11. (3) Fractured and bedding dissolution zones; there are the dripping phenomenon in the section of 74.3–74.7 m in the adit PD-9 and 43–49 m in the adit PD-10.

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. Percentage of Lu values for different burial depths (%).

Burial Depth (m)	5–10 Lu	1–5 Lu	<1 Lu
20	13	45	42
40	2	66	32
60	1	66	33
80	0	75	25
100	0	59	41
120	0	29	71
>120	0	0	100

1 Lu is equivalent to 10⁻⁵ cm/s.

). 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. Permeability statistics of different strata.

Stratigraphic Age	Stratigraphic Symbol	Stratigraphic Thickness (m)	Lithology	Permeability of Rock Mass
Paleogene	E2l	>220	Conglomerate and sandstone	Moderate
Triassic	T2g4	>150	Dolomitic limestone	Strong
	T2g3	570	Marl and argillaceous dolomite	Moderate
	T2g2	120	Mudstone and dolomite interbedding	Relative impermeability
	T2g1-2	80	Argillaceous dolomite	Weak karst permeable layer
	T2g1-1	240	Limestone and dolomitic limestone	Strong karst permeable layer
	T1yn2	50	Mudstone with limestone	Relative impermeability
	T1yn1	284	Muddy limestone	Strong
	T1f	625	Marlite	Impermeable layer
Permian	P2c	2.5–20	Limestone	Strong
	P2l2	234	Sandstone with clay	Impermeable layer
	P2l1	38	Bioclastic limestone	Strong
	P1m	318–544	Bioclastic limestone	Strong

.

(3) Recharge and discharge of groundwater

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₁₈, K₁₉, and Longdang karst pipes on the left bank are mainly developed in thick massive limestone T₂g^3-6. However, the karst development is weak in the heavier argillaceous layers of T₂g^3-3 and T₂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. Characteristics of karst pipe at dam site.

Number of Karst Pipeline	Location	Export Elevation of Karst Pipeline (m)	Length of Karst Pipe (m)	Export Flow of Karst Pipeline (m3/s)	Characteristics of Karst Pipeline
1	Left bank	983.70	2000	0.005	Development along the strike of strata T2g3-1-1
2	Right bank	1014.37	600	0.030	Development along the strike of strata T2g3-2-1
3	Right bank	980.31	150	0.020	Development along the layer surface of strata T2g3-6
4	Right bank	973.63	600	0.00002	Development along the layer surface of strata T2g3-1-1
5	Left bank	958.50	1000	0.015	Development along the layer surface of strata T2g3-6

. Furthermore, the elevation of Karst development has obvious zonation. The sinkholes K_l7, K₁₈, and K₁₉ 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.

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

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

3.1. Mathematical Model of Groundwater Flow

The transient groundwater flow in heterogeneous anisotropic fractured medium can be simulated using the Boussinesq equation [42]:

$$\{\begin{array}{c}\mu \frac{\partial h_f}{\partial t}=\frac{\partial }{\partial x}\left(K_{\mathit{fx}}\frac{\partial h_f}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{\mathit{fy}}\frac{\partial h_f}{\partial y}\right)\frac{\partial }{\partial z}\left(K_{\mathit{fz}}\frac{\partial h_f}{\partial z}\right)+\epsilon ,x,y,z\in \mathsf{\Omega },t\ge 0\\ \begin{array}{c}{h_f\left(x,y,z,t\right)|}_{t=0}{=h}_{f0},x,y,z\in \mathsf{\Omega },t\ge 0,\\ {h_{f1}\left(x,y,z,t\right)|}_{{\Gamma }_1}{=h}_{f1},x,y,z\in {\Gamma }_1,t\ge 0\\ {K_{fn}\frac{\partial h_f}{\underset{\partial n}{\to }}|}_{{\Gamma }_2}{=q}_f\left(x,y,z,t\right),x,y,z\in {\Gamma }_2,t\ge 0\end{array}\end{array}$$

(1)

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³/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³/T·L²], it is positive for groundwater inflow and negative for outflow; Ω is the domain of the fractured medium, and Γ₁ and Γ₂ are domains of the first and second type boundary. (according to Abbreviations).

3.2. 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:

$$h_1=\lambda \frac{L}{d}\frac{u^2}{2g}$$

(2)

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

$$J=\frac{h_1}{L}$$

(3)

Substituting Equation (3) and the seepage velocity (V) into Equation (2), it can be rewritten as:

$$\mathit{JL}=\mathsf{\lambda }\frac{L}{d}\frac{\mathit{VV}}{2g}$$

(4)

Therefore, Equation (4) can be written in the form of Darcy’s Law:

$$V=\frac{2\mathit{gd}}{\lambda V}J$$

(5)

A converted hydraulic conductivity, K_c, is defined as:

$$K_c=\frac{2\mathit{gd}}{\lambda V}$$

(6)

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

$$\mathsf{\lambda }=\frac{64}{\mathrm{Re}}$$

(7)

For turbulent flow (Re ≥ 2320),

$$\mathsf{\lambda }={\left[1.74-21\mathrm{g}\left(\frac{\Delta }{d}+\frac{21.25}{{\mathrm{Re}}^{0.9}}\right)\right]}^{-2}$$

(8)

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, ${\left[K\right]}_L^e$, can be expressed as [41]:

$${\left[K\right]}_L^e=\left[\begin{array}{c}\begin{array}{ccc}{K}_c& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\end{array}\right]=diag\left(K_c,0,0\right)$$

(9)

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, ${\left[K\right]}_G^e$, can be calculated with that of the local coordinate system:

$${\left[K\right]}_G^e=\left[P\right]{\left[K\right]}_L^e{\left[P\right]}^T$$

(10)

where [P] is the orthogonal matrix, and it can be calculated using the local coordinates of karst pipe.

Assuming

$$\{\begin{array}{c}{\alpha }_1=\left\{\begin{array}{ccc}\cos \alpha ,& \cos \beta ,& \cos \gamma \end{array}\right\}\\ {\alpha }_2=\left\{\begin{array}{ccc}{\cos \alpha }^{\prime },& {\cos \beta }^{\prime },& {\cos \gamma }^{\prime }\end{array}\right\}\\ {\alpha }_3=\left\{\begin{array}{ccc}{\cos \alpha }^{″},& {\cos \beta }^{″},& \cos \gamma \end{array}\right\}\end{array}$$

(11)

And assuming that the starting and ending coordinates of the axis of karst pipe are M₁ (x₁, y₁, z₁), M₂ (x₂, y₂, z₂), respectively. The directional cosine of vector (α₁) or ox’-axis can be given as:

$$\{\begin{array}{c}\cos \alpha =\frac{x_2-x_1}{\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_1-z_2\right)}^2}}\\ \cos \beta =\frac{y_2-y_1}{\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}}\\ \cos \gamma =\frac{z_2-z_1}{\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}}\end{array}$$

(12)

Because the ox’-axis, oy’-axis, and oz’-axis are orthogonal to each other, α₁, α₂, and α₃ are also orthogonal, that is

$$\{\begin{array}{c}{\alpha }_1^T{\alpha }_2=0\\ {\alpha }_1^T{\alpha }_3=0\end{array}$$

(13)

Therefore, vector α₂ and α₃ satisfy as:

$${\alpha }_1^{\mathrm{T}}\mathrm{X}=0$$

(14)

We can get the following algebraic equation:

$${\cos \alpha x}_1{+\cos \beta x}_2{+\cos \gamma x}_3=0$$

(15)

Assuming a = cosα, b = cosβ, and c = cosγ, if a is not equal to zero, the orthogonal matrix can be expressed as:

$$\left[\mathrm{P}\right]=\left[\begin{array}{ccc}a& -\frac{b}{a}\sqrt{\frac{a^2}{a^2{+b}^2}}& -\frac{\mathrm{ac}}{\sqrt{a^2{+b}^2}}\\ b& \sqrt{\frac{a^2}{a^2{+b}^2}}& -\frac{\mathrm{bc}}{\sqrt{a^2{+b}^2}}\\ c& 0& \sqrt{a+b}\end{array}\right]$$

(16)

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.

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

$${h_f\left(x,y,z\right)|}_{{\Gamma }_c}=h_k\left(x,y,z\right)\hspace{1em}\hspace{1em}\left(x,y,z\right)\in {\Gamma }_c,t>0$$

(17)

where ${\Gamma }_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).

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

$$E\left(K_j^i\right)={\displaystyle \sum }_{n=1}^M{\mathsf{\omega }}_n\sqrt{{\left(h_n^c-h_n^0\right)}^2}$$

(18)

where $K_j^i$ 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_n^c$ and $h_n^0$ are the calculated and observed hydraulic heads.

3.5. 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:

$$K =0.528\mathsf{\omega }\lg \left(\frac{0.66L}{r_0}\right)$$

(19)

If the test section is far from the impermeable layer, the equation is given as:

$$K =0.528\mathsf{\omega }\lg \left(\frac{1.32L}{r_0}\right)$$

(20)

where ω is the permeable rate, the value is equal to 100 Lu, L is the length of the test section, and r₀ is the radius of the test hole.

4. Results and Discussion

4.1. 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³/s. The area of the ring was 490.625 cm². 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³/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. Equivalent permeability coefficients for different zones (m/d).

Zones of Vertical Direction	Zones of Horizontal Direction
	Strong Permeable Zone	Weak Permeable Zone	Slight Permeable Zone
Strong weathering zone	0.2–5.0	0.1–0.5	0.1–0.3
Weak weathering zone	0.1–0.5	0.1–0.3	0.05–0.2
Slight weathering zone	0.1–0.3	0.05–0.2	0.005–0.05

.

(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), the calculation results of the main hydraulic conductivities for karst conduits are listed in

Table 6. Main hydraulic conductivities for karst conduits.

Number of Karst Pipes	Main Hydraulic Conductivities (m/d)
	Kx	Ky	Kz
1	4745	582	363
2	5795	1602	241
3	1650	442	485
4	4871	1124	410
5	3368	1226	537

.

4.2. 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. Permeability coefficients of fracture-karst media by inversion method (m/d).

Zones of Vertical Direction	Zones of Horizontal Direction
	Strong Permeable Zone			Weak Permeable Zone			Slight Permeable Zone
	Kx	Ky	Kz	Kx	Ky	Kz	Kx	Ky	Kz
Strong weathering zone	0.49	0.49	0.05	0.39	0.38	0.04	0.28	0.29	0.03
Weak weathering zone	0.39	0.38	0.04	0.29	0.30	0.03	0.15	0.15	0.02
Slight weathering zone	0.28	0.28	0.03	0.11	0.10	0.01	0.03	0.02	0.003

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

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

4.4. 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. Prediction of leakage for karst pipes and powerhouse during the operation period of the reservoir (m³/d).

Cases	Karst Pipes				Underground Powerhouse
	1	2 and 3	4	5
Case 1	0.43	0.36	1.99	313.39	2272.51
Case 2	224.34	162.86	465.95	306.89	2307.62
Case 3	0.16	0.12	0.82	117.92	1205.22
Case 4	80.86	68.30	200.97	113.90	1232.58

). 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³/d.

5. 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³/d.