Research on Seepage Control of Jurong Pumped Storage Hydroelectric Power Station

: Based on the geological and hydrogeological conditions of the Jurong Pumped Storage Hydroelectric Power Station (JPSHP), a 3D groundwater ﬂow model was developed in the power station area, which took into account the heterogeneity and anisotropy of fractured rocks. A control inversion method for fractured rock structural planes was proposed, where larger-scale fractures were used as water-conducting media and the relatively intact rock matrix was used as water-storage media. A statistical method was used to obtain the geometric parameter values of the structural planes, so as to obtain the hydraulic conductivity tensor of the fractured rocks. Combining the impermeable drainage systems of the upper storage reservoir, underground powerhouse and lower storage reservoir, the 3D groundwater seepage ﬁeld in the study area was predicted using the calibrated model. The leakage amounts of the upper storage reservoir, powerhouse and lower storage reservoir were 710.48 m 3 /d, 969.95 m 3 /d and 1657.55 m 3 /d, respectively. The leakage changes of the upper storage reservoir, powerhouse and lower storage reservoir were discussed under the partial and full failure of the anti-seepage system. The research results provide a scientiﬁc basis for the seepage control of the power station, and it is recommended to strengthen the seepage control of the upper and lower storage reservoirs and the underground powerhouse to avoid excessive leakage and affect the efﬁciency of the reservoir operation. and hydrogeological


Introduction
In recent years, with the increase in water conservancy and hydropower projects, seepage problems have attracted widespread attention. When there are problems with seepage control, the efficiency of the project can be affected and even cause safety and stability problems such as landslides, submersion and permeability failure. An example of this is the famous French Malpasset dam failure. Some researchers have studied the cause of the incident. They found that the presence of fractures and the permeability of the rock mass eventually led to the collapse of the Malpasset dam, which also reflects the necessity of paying attention to fractures when conducting seepage control in water conservancy and hydropower projects [1,2]. The hydraulic conductivity or hydraulic conductivity tensor of a fractured rock is an important hydrogeological parameter when studying the movement of groundwater in a rock mass [3].
The field hydrogeological test is a quick method used to understand the permeability characteristics of fractured rock in the field, which mainly includes a drainage test [4], piezometric test [5][6][7], pumping test [8] and micro-water test [9]. The above methods mainly include test principles and methods, equipment development, test procedure improvement, test data processing and interpretation, error analysis, influence factors, numerical reservoir adopts reinforced concrete panels for seepage control, and the bottom of the re voir adopts a geomembrane for seepage control. A reservoir road is set up along the per eter of the upper storage reservoir with an elevation of 271.80 m. The water convey system is located in the main peak of Lunshan mountain wit total length of 1327.89-1363.72 m (from the inlet and outlet of the upper storage reserv to the inlet and outlet of the lower storage reservoir) ( Figure 2). It includes the inlet a outlet of the upper storage reservoir, the upper flat section of the diversion, the divers pressure chamber, the diversion shaft, the lower flat section of the diversion, the divers steel fork pipe and the branch pipe before the underground powerhouse, the tailwa tunnel and the inlet and outlet of the lower storage reservoir. The inlet and outlet of upper storage reservoir is located on the left bank with a floor elevation of 220 m. T diversion system adopts a three-hole and six-machine shaft arrangement. The inlet a outlet of the lower storage reservoir are located at a small ridge on the right bank wit floor elevation of 50 m. The size of the underground powerhouse is 246. 5  The climate of Jurong is a north subtropical monsoon climate. There are high temp atures and rain in summer, warm and humid winters, high rainfall and abundant wa resources. The annual average temperature is 17.1 °C, and the annual average precip tion is 1222.3 mm.

Terrain and Landforms
The reservoir area is located in the south of the Yangtze River and the north of Taihu Lake. The site is located in the Guanyin and Lunshan mountain ranges on the south of the middle section of the Ningzhen mountain range, which runs from northeast to southwest, and is a low mountainous terrain with high terrain in the west and low terrain in the east. The south of the site is the Jurong Basin, and its ground elevation is mostly 35-50 m. The top elevation is 50-80 m and has the typical characteristics of riverbed phase and river floodplain phase, which is situated in Class I and Class II base terraces. The south-east is the Gao Li Mountain range with a summit elevation of 425 m, with the Lun Shui reservoir sandwiched between it and the Lun Shan Mountain range. The north-east and north of the Lun Shan Mountain consist of low hills, the summits of which are mostly between 115 and 200 m in elevation. The study area spans the main peak of Lunshan from west to east. The summit of Lunshan is 400.4 m in elevation. The ridge generally extends from west to east and is about 1.8 km long with an elevation of 60-400 m. The northwest of Lunshan is the recently exploited Jurong Forestry Quarry.
The upper storage reservoir is located on the southwest of Lunsan Mountain in a large sentinel ditch. The ditch flows to the southeast with a short extension and a narrow ditch bottom. The ground elevation is 90-110 m. The ditch at the dam site is an asymmetrical V shape, and the reservoir is basin-shaped. The ditch is surrounded by ridges to the east, north and southwest in a near circle, consisting of peaks and passes. The highest point around the reservoir is at the northeast corner of the reservoir, which is the main peak of Lunshan. The peak on the west is 375.1 m, and the ground elevation of the pass is 288. .60 m. The mouth of the ditch is located in the south-east. The overall slope of the bank is 25-40°. There are a number of small shallow erosion gullies and an uneven terrain in the reservoir basin with an elevation of 100-140 m.
The lower storage reservoir is located on the northeast of Lunshan Mountain, and the natural reservoir basin is composed of the Sister Bridge Ditch and a small washout ditch on its north side. The Sister Bridge Ditch flows in a south-easterly direction, and the bottom of the ditch is flat with a ground elevation of 56-82 m. The south bank, or right bank, is Lunshan Mountain. The north bank, or left bank, has wide and gentle hills with a small topographic slope and a summit elevation of 150-226 m. These small ridges are mostly in the direction of SN to NWW. The largest branch ditch is the Tiger Dam Ditch in front of the dam, which meanders from north to south to the dam site and joins the Sister Bridge Ditch at the dam site into Lunshan Reservoir. The left bank of the Laohu Dam ditch is the head of the left dam with a ridge elevation of 134 m. The gully outside the reservoir is well developed.
The water transmission and power generation system is buried within the Lunshan Mountain, crossing the main peak of Lunshan Mountain from south-west to north-east with the ground elevation ranging from 70 to 400 m. The highest peak is the main peak of Lunshan Mountain. The ridge of Lunshan Mountain is gentle, and the branch ditches along both sides are undeveloped and shallowly cut. The slope of inlet and outlet for the

Terrain and Landforms
The reservoir area is located in the south of the Yangtze River and the north of Taihu Lake. The site is located in the Guanyin and Lunshan mountain ranges on the south of the middle section of the Ningzhen mountain range, which runs from northeast to southwest, and is a low mountainous terrain with high terrain in the west and low terrain in the east. The south of the site is the Jurong Basin, and its ground elevation is mostly 35-50 m. The top elevation is 50-80 m and has the typical characteristics of riverbed phase and river floodplain phase, which is situated in Class I and Class II base terraces. The south-east is the Gao Li Mountain range with a summit elevation of 425 m, with the Lun Shui reservoir sandwiched between it and the Lun Shan Mountain range. The north-east and north of the Lun Shan Mountain consist of low hills, the summits of which are mostly between 115 and 200 m in elevation. The study area spans the main peak of Lunshan from west to east. The summit of Lunshan is 400.4 m in elevation. The ridge generally extends from west to east and is about 1.8 km long with an elevation of 60-400 m. The northwest of Lunshan is the recently exploited Jurong Forestry Quarry.
The upper storage reservoir is located on the southwest of Lunsan Mountain in a large sentinel ditch. The ditch flows to the southeast with a short extension and a narrow ditch bottom. The ground elevation is 90-110 m. The ditch at the dam site is an asymmetrical V shape, and the reservoir is basin-shaped. The ditch is surrounded by ridges to the east, north and southwest in a near circle, consisting of peaks and passes. The highest point around the reservoir is at the northeast corner of the reservoir, which is the main peak of Lunshan. The peak on the west is 375.1 m, and the ground elevation of the pass is 288.30-313.60 m. The mouth of the ditch is located in the south-east. The overall slope of the bank is 25-40 • . There are a number of small shallow erosion gullies and an uneven terrain in the reservoir basin with an elevation of 100-140 m.
The lower storage reservoir is located on the northeast of Lunshan Mountain, and the natural reservoir basin is composed of the Sister Bridge Ditch and a small washout ditch on its north side. The Sister Bridge Ditch flows in a south-easterly direction, and the bottom of the ditch is flat with a ground elevation of 56-82 m. The south bank, or right bank, is Lunshan Mountain. The north bank, or left bank, has wide and gentle hills with a small topographic slope and a summit elevation of 150-226 m. These small ridges are mostly in the direction of SN to NWW. The largest branch ditch is the Tiger Dam Ditch in front of the dam, which meanders from north to south to the dam site and joins the Sister Bridge Ditch at the dam site into Lunshan Reservoir. The left bank of the Laohu Dam ditch is the head of the left dam with a ridge elevation of 134 m. The gully outside the reservoir is well developed.
The water transmission and power generation system is buried within the Lunshan Mountain, crossing the main peak of Lunshan Mountain from south-west to north-east with the ground elevation ranging from 70 to 400 m. The highest peak is the main peak of Lunshan Mountain. The ridge of Lunshan Mountain is gentle, and the branch ditches along both sides are undeveloped and shallowly cut. The slope of inlet and outlet for the upper storage reservoir is slightly steeper, with a slope of 30-35 • . The inlet and outlet of the lower storage reservoir is gentler. The overlying rock thickness of the upper flat section is 65-160 m. The overlying rock thickness of the lower flat section is 180-270 m. It is 95-180 m thickness in the underground powerhouse.

Stratigraphic Lithology
The stratigraphy exposed in the study area is classified from old to new as: (1) The stratum of the Sinian Dengying Formation (Z 2 dn) is mainly gray-white, medium-thicklayered fine-crystalline dolomite, which is mainly distributed in the middle of the main underground powerhouse, the secondary underground powerhouse and the main transformer cave. At an elevation of 51. dolomite and fragmented dolomite. The argillaceous dolomite is thin layered, and the surface is strongly weathered. Muddy filling is more common in contact zones and layers, which is a soft interbedded layer. The contact surface is dominated by cuttings with mud. It is distributed in strips. According to the analysis of the horizon-slice map, at an elevation of 51.5 m, this rock mass is mainly distributed on the right end wall of the powerhouse, accounting for about 2% of the powerhouse rock mass, and the main transformer caves and busbar cave are not distributed. At an elevation of 37 m, this rock mass is distributed outside the power station area.

Geological Structure
The geological structure of the study area is complex. Folding is not developed. The stratigraphy is monoclinic. The main tectonic traces are mainly faults, joints and fractures. According to the attitude of structural planes, the first type is in the NW direction, such as F 6 and F 8 . The width is 0.5-2 m. The lithology is tectonic conglomerate and fractured rock. The contact zone is mostly seen as dissolution fractures. The second type is in the NWW direction with the width of 0.5-2 m. The lithology is tectonic conglomerate and fractured rock. The contact zone is mostly seen as dissolution fractures. The third type is in the NNW direction, such as F 7 and F 11 . The width is 1-3 m. The lithology is tectonic conglomerate and fractured rock. The fractures are mostly filled by the intrusion of rock veins. The fourth type is in the NE direction with the width of 2-3 m, such as F 2 and F 3 . The lithology is tectonic conglomerate and fractured rock with minor fragmented chalk and faulted mud.

Hydrogeological Conditions
(1) Permeability of rock mass The statistical results of conventional water pressure tests show that 54.6% of the rock mass is slightly permeable, 29.5% is weakly permeable rock mass, 10.2% is weakly permeable rock mass and 5.6% is moderately permeable rock mass. The permeability of the various sets of strata is basically similar. Among them, the Qixia Formation, the Qinglong Formation, the Zhouchong Village, the Longtan Formation, the Shangdang Formation, the Mufushan Formation and the Paotaishan Formation are relatively weak in water permeability, with slightly permeable rock mass accounting for more than 60%. The weakly permeable formations of the Lunshan Formation, the Honghuayuan Formation, the Dengying Formation and the Guanyintai Group reach 46-49%, and the slightly permeable rocks make up less than 50%. The permeability of non-soluble rock in Yangchong Formation is mainly weak and slightly permeable, very small and moderately permeable, without strong permeable rock mass. The impure carbonate rock mass is mainly weak to slightly permeable, and there is no rock mass with medium or above permeability. The medium and above permeable rock mass in the carbonate rock contains caves or wide fractures, dissolution fissures, etc. According to the relationship between flow rate and pressure, the

Hydrogeological Conditions
(1) Permeability of rock mass The statistical results of conventional water pressure tests show that 54.6% of the rock mass is slightly permeable, 29.5% is weakly permeable rock mass, 10.2% is weakly permeable rock mass and 5.6% is moderately permeable rock mass. The permeability of the various sets of strata is basically similar. Among them, the Qixia Formation, the Qinglong Formation, the Zhouchong Village, the Longtan Formation, the Shangdang Formation, the Mufushan Formation and the Paotaishan Formation are relatively weak in water permeability, with slightly permeable rock mass accounting for more than 60%. The weakly permeable formations of the Lunshan Formation, the Honghuayuan Formation, the Dengying Formation and the Guanyintai Group reach 46-49%, and the slightly permeable rocks make up less than 50%. The permeability of non-soluble rock in Yangchong Formation is mainly weak and slightly permeable, very small and moderately permeable, without strong permeable rock mass. The impure carbonate rock mass is mainly weak to slightly permeable, and there is no rock mass with medium or above permeability. The medium and above permeable rock mass in the carbonate rock contains caves or wide fractures, dissolution fissures, etc. According to the relationship between flow rate and pressure, the fracture state did not change during the water pressure test. At the same time, two sets of high water pressure tests were carried out in the underground powerhouse area, and the maximum water permeability was 2.95 Lu, which indicated that the rock mass of the underground powerhouse is weakly to slightly permeable.
(2) Groundwater types There are three main groundwater types in the study area: pore water, fracture water and karst water, with fracture-karst water being the main type.

Fracture water
Fracture water is mainly found in bedrock fractures and fault zones. It is divided into two categories: tectonic fracture water and weathered fracture water. The thickness of the aquifer is large, recharged by atmospheric precipitation, drained in the form of springs, controlled by faults and fractures and affected by karsts. It is mainly distributed in veins and bands, especially at the low groundwater level in the Lunshan Mountains. Dripping and seepage along faults, fractures and karsts are more common. The degree of development of tectonic fracture water is determined by the mechanical properties and shale content of the rock. Strata dominated by quartz sandstones are brittle, and the fractures are developed under the action of structure, and the water richness is good. Strata dominated by siltstone are soft, with a high argillaceous content, undeveloped fractures, which are easily filled and closed, and poor water content. Weathered fracture water generally decreases with the increase in depth. It is widely distributed on the plane. The shallow fractures are connected in a network. Due to the influence of topography, stratum lithology, geological structure, surface water and groundwater, it varies from place to place. The groundwater in the weathered zone mainly recharges by the fracture water in the hills, and has weak pressure bearing capacity, and the water inflow of a single well is less than 10 m 3 /d.

Karst water
Karst is the main water-bearing space where groundwater occurs in carbonate strata, and it occupies a major position in the study area. The study area is characterized by weak to moderately strong karst strata, surface erosion recesses, surface erosion fractures, solution channels and small solution holes. Drill holes and exploration caves revealed the presence of karst fractures and caves at certain depths. Intermittent water gushing in the fault zone and tracer test results in the borehole show that there are crevices or local pipeline flows in the Lunshan rock mass. Therefore, the karst water in the study area is diverse, including karst-fracture water, pore-fracture water and crevice-fissure-pipe water, etc. It is the main type of groundwater migration in the study area, and atmospheric precipitation is the main recharge source. Water and ditch water are also sources of low-elevation karst water, which are discharged in the form of karst springs or crevice water. According to the outcropping conditions, lithological characteristics, geological structure and topography of carbonate rock formations, combined with regional hydrogeological data and spring water outcrops, water richness is divided into two levels: formations with strong water richness and formations with weak water richness.  The topography of the study area is large undulations, and rainfall is the main source for groundwater recharge. In the exposed bedrock mountains, groundwater recharge sources are dominated by depressions, caves and vertically dispersed infiltration from surface solution gaps, fractures and gullies. The variation of groundwater level is closely related to the amount of rainfall. Groundwater runoff is mainly controlled by topography and tectonics, especially the control of bedding fissures, water-resistant rock veins and fractures, forming spring areas. Groundwater in the gullies is constrained by topographic and hydrological conditions and drains to the Gangwash and gully runoff. From the analysis of the topography, geological conditions and the distribution of spring outcrops, the short runoff pathway and local recharge and discharge is one of the main characteristics of bedrock groundwater transport in the area.
Discharge conditions are dominated by evaporation, springs and inter-stratigraphic drainage, which are the main modes of groundwater discharge in the study area. There are deep wells for groundwater extraction in Judong farm and the army at the western end of Guanyintai Mountain, and shallow wells for mining shallow pore water on farms such as Zimeiqiao and Shangmeng. The terrain in the mountains is steep. The rock outcrops receive infiltration recharge from rainfall, and groundwater runoff circulation conditions are well established, with short pathways and fast flow rates, so that groundwater discharge in the form of springs from bedrock is extremely common, especially for karst water. In the contact zone between soluble and non-soluble rocks, controlled by tectonic and topographic cutting conditions, springs are frequently exposed, with a maximum flow of over 1000 m 3 /d for Xianren Cave Spring. These springs are controlled by NNW-, NNE-and NWW-oriented fractures and veins. This makes the bedrock groundwater show a certain vein-like flow or strip-like flow of runoff and discharge pattern. Evaporation from low-lying gully and valley areas and the extraction of folk wells are also groundwater drainage routes. In addition, the infiltration of groundwater along the water barrier in mountainous areas is Water 2022, 14, 141 9 of 24 also a means of drainage. There is a weak runoff discharge to the deeper and eastern parts along the Sister Bridge fracture zone.

Mathematical Model of Rock Seepage and Its Solution
There is some karst in the JPSHP, but its connectivity is generally poor and no obvious karst pipes are formed, so the influence of rock fractures is mainly considered. When studying the groundwater movement of rock in the area of a hydraulic building, the waterbearing medium can be studied as a continuous medium if the ratio of the average spacing of fractures in the reservoir area to the minimum size of the building is less than 1:20. The ratio of the average spacing of the fractures in the JPSHP to the minimum size of the building is much less than 1:20, thus the mathematical model of a continuous medium can be used to simulate and calculate the seepage field. Consider the following assumptions: (i) the water is incompressible; (ii) the aquifer material is incompressible; (iii) external loads (e.g., atmospheric pressure, overburden) on the aquifer are constant; (iv) the groundwater is flowing slowly (Reynolds number is less than unity). The 3D movement of groundwater satisfies the following control equation: where ∇ is the Hamilton operator, K is the hydraulic conductivity tensor, H is the hydraulic head at any point in the seepage field The governing Equation (1) for groundwater movement in a rock mass is solved by the following initial and boundary conditions: Initial conditions can be given as First boundary conditions or Dirichlet conditions are Second boundary conditions or Neumann conditions can be expressed as where H(x, y, z, t) is the hydraulic head at any point (x, y, z) in the region at some time t [L], φ 0 (x, y, z), φ 1 (x, y, z, t) and q 0 are known functions and n is the direction of the outer normal to the boundary Γ 2 . Using the Garlerkin finite unit method, the final algebraic equation for the entire seepage problem is solved by where {H} is the unknown nodal hydraulic head array, {F} is the known right end term and [G] is the overall penetration matrix or conduction matrix. dH/dt can be approximated with the finite difference formula based on the Forward Euler method [32]: The equation can be rewritten as using the implicit differential format: That is,

Drainage Hole Simulation and Cross-Sectional Flow Rate Determination
The JPSHP has four drainage corridors and holes, including the top, upper, middle and bottom levels. A rod unit is used to deal with these drainage holes, and the method is illustrated with a typical unit. For the hexahedral unit (Figure 4), it is assumed that the node 1485 is the profile where the drainage holes are located. The rod unit penetration matrix is assigned as that is, add or subtract p s to the diagonal elements a 11 , a 44 , a 88 and a 55 of the permeability matrix of a typical unit, respectively.
where b is the well spacing, γ ω is the well radius, K is the hydraulic conductivity, K 1 , K 2 and K 3 are the three main hydraulic conductivities of the rock mass and ∆F is the area of the unit profile where the hole is located.
Water 2022, 13, x FOR PEER REVIEW 10 of 24 The equation can be rewritten as using the implicit differential format: That is,

Drainage Hole Simulation and Cross-Sectional Flow Rate Determination
The JPSHP has four drainage corridors and holes, including the top, upper, middle and bottom levels. A rod unit is used to deal with these drainage holes, and the method is illustrated with a typical unit. For the hexahedral unit (Figure 4), it is assumed that the node 1485 is the profile where the drainage holes are located. The rod unit penetration matrix is assigned as that is, add or subtract ps to the diagonal elements a11, a44, a88 and a55 of the permeability matrix of a typical unit, respectively.
where b is the well spacing, γω is the well radius, K is the hydraulic conductivity, K1, K2 and K3 are the three main hydraulic conductivities of the rock mass and ΔF is the area of the unit profile where the hole is located. By solving Equation (8), the seepage field over the entire seepage area can be obtained. As the partial derivative of the hydraulic head on the unit nodes with respect to the coordinates cannot be obtained directly when solving the seepage field using the finite element method, the intermediate section method is used when calculating the section flow ( Figure 5). Using the node 8 flow element as an example, the equation for the seepage through a section S is given as By solving Equation (8), the seepage field over the entire seepage area can be obtained. As the partial derivative of the hydraulic head on the unit nodes with respect to the coordinates cannot be obtained directly when solving the seepage field using the finite element method, the intermediate section method is used when calculating the section flow ( Figure 5). Using the node 8 flow element as an example, the equation for the seepage through a section S is given as where S is the overflow section and n is the unit vector in the direction of the positive normal to the overflow section.  (10) where S is the overflow section and n is the unit vector in the direction of the positive normal to the overflow section. For any eight-node hexahedral isoparametric element, the midsection abcd is chosen as the overflow section S and S is projected onto the yoz, zox, and xoy planes, noted as Sx, Sy, and Sz, respectively; then, the seepage through the midsection of the element is given as For any eight-node hexahedral isoparametric element, the midsection abcd is chosen as the overflow section S and S is projected onto the yoz, zox, and xoy planes, noted as S x , S y , and S z , respectively; then, the seepage through the midsection of the element is given as

Parameter Inversion Method
(1) Inversion method based on structural surface control Based on the error between the observed and calculated hydraulic heads at a number of known coordinate positions in the study seepage area, the least squares method is used to establish an objective function: where K i j is the coefficient, the superscript i indicates the i-th subzone according to the rock permeability, i = I, II, . . . , NK, NK is the total number of subzones, the subscript j denotes the j-th coefficient in the i-th subzone, j = 1, 2, . . . , MK, MK is the total number of parameters for a subzone and ω k is the weight function for the k-th head observation Obviously, a set of parameters reflecting the objective fractured rock seepage pattern, K i j , needs to be configured such that the objective function E tends to a minimum. Considering the actual field work, the structural plane attitudes are represented by the strike, dip direction and dip angle of the structural planes. For a 3D seepage problem, if the x, y and z axes are chosen to be east, north and vertical, respectively, the expression for the permeability tensor in a Cartesian coordinate system can be rewritten as − sin β j sin 2 γ j cos β j − cos β j sin γ j cos γ j − sin β j cos j sin 2 γ j 1 − sin 2 β j sin 2 γ j − sin β j sin γ j cos γ j − cos β j sin γ j cos γ j − sin β j sin γ j cos γ j 1 − cos 2 γ j   (13) where β j and γ j are the dip direction and dip angle of the structural planes in group j of the rock mass, respectively, K ej is the equivalent hydraulic conductivity [L/T] for the j-th group of fractures and M is the number of structural plane groups in the rock mass that control groundwater seepage. From tectonic geology, the development and distribution of fracture planes in a given area is regular and is developed in groups. Therefore, statistical methods can be used to firstly count the development and distribution patterns of structural planes in the study area, so that the distribution function of the structure planes and the statistical values of the geometric parameters are obtained. In this way, the matrix in the tensor expression can be uniquely determined and the unknowns are only the coefficients, K ej . If the hydraulic conductivity can be back-calculated for each group of fractures, then the permeability tensor can be uniquely determined from the tensor expressions. Among the subzones with a total number of NK, the hydraulic conductivity of the j-th group of fractures in the i-th subzone can be expressed as K i ej , and K i ej > 0, i = I, I I, . . . , NK; j = 1, 2, . . . , MK As the correlation between the objective function E and the independent variable K i ej cannot be described by an explicit expression, it is easier to solve the inverse problem using an optimal method. Each step of this method can be solved in a similar way to the positive problem, initially using the 3D finite element method to find the approximate distribution of the hydraulic head based on a set of trial solutions K i ej . The parameters are then substituted with the most recently obtained parameters and continuously tuned so that the calculated groundwater levels gradually approximate the observed values, that is, E tends to a minimum.

(2) Treatment of groundwater levels in observation holes
If an observation hole is drilled in a water-bearing rock, water can enter the hole across the thickness of the rock through which the observation hole is drilled. Therefore, the observed groundwater level, H o , is not the true groundwater level at a particular point, but rather the average value at each point on the section of the water-bearing rock uncovered by the borehole in the observation hole: where L is the thickness of the water-bearing rock uncovered by the observation hole [L]; H is the water level at different depths (locations) in the observation hole [L]. When calculating with the finite element method, the local coordinates of the observation hole in a typical element are first generated from the overall coordinates. Then, according to the interpolation formula, the water level at a point in the borehole in the vertical direction is calculated. In this way, the accumulative value of the water level of the borehole in the typical element can be approximated as The calculated water level, H c , in the entire observation hole is where H + and H − are the groundwater levels at the top and bottom for the i-th typical element disclosed by the observation hole, respectively. ∆L i is the vertical length [L] of the i-th element uncovered by the observation hole, n is the total number of elements and L can be expressed as (

3) Determination of the weight function ω
There is a certain amount of error between the calculated and observed water levels at the position corresponding to the observation hole. Additionally, there are various reasons for this error, including errors in the numerical calculation itself, errors in the simplification of the model, errors in the formation of the observation hole, errors in observation, etc. Theoretically, if the mathematical model established to describe groundwater movement is identical to the actual groundwater movement and the observed groundwater level is not disturbed in any way, the calculated and observed water levels should be equal at the same point, excluding the numerical errors in the numerical calculation itself. However, due to the complexity of the geological body and the randomness of the water control structure, although the overall view of groundwater movement basically conforms to the established mathematical model, but specifically to a certain borehole and a certain local point, the actual groundwater level is not necessarily equal to the calculated value, so that there is an inherent error between the calculated and observed water levels, which cannot be eliminated by improving the calculation method and adjusting the model. Therefore, in order to correctly simulate the spatial distribution of the actual groundwater movement from a macroscopic perspective and the flow field as a whole, and to eliminate the accuracy of the entire calculation result caused by the error of individual points, the difference between the water level of each observation hole and the calculated water level is multiplied by the weight function, so that the calculation as a whole meets the control accuracy requirements.
Suppose there are n observation holes in the study area, the observed water level of the i-th observation hole is H o i and the corresponding calculated water level is H c i . Define the weight function: The n weight functions are arranged in order from largest to smallest, while the water level error is arranged in order from smallest to largest, and the renumbered weight functions and water level error are brought into the objective function Equation (12), which can be solved to obtain the objective function value. Because the size of this weight function ω k (k = 1,2,···,n) and its arrangement is with the parameter inversion process and the water level error H o k − H c k 2 dynamic self-adjustment, the difference between the observation and calculated water levels multiplied by the weight function can automatically and effectively eliminate the error of individual points, so that the calculation in the overall control accuracy can achieve the requirements.

Results of Parameter Inversion
The computational domain is bounded by the upper and lower storage reservoirs.  Table 1. According to the weathering characteristics and structural plane development in the power station area, the study area is divided into three different permeability zones in the profile, (I) strong-moderate weathering zone, (II) weak weathering zone and (III) slightly weathered zone.
of −30 m. The upper storage reservoir treats the water level of the observation hole as the first type of boundary condition. The lower storage reservoir treats the ditch water as the first type of boundary condition, and both sides and the bottom of the model are treated as no flow boundary conditions. The permeability of the rock in the JPSHP area is mainly controlled by the structural planes, and several groups of faults and fractures are developed. In the inversion, the faults and four groups of joints and fractures, which have a controlling effect on the permeability of the underground numbers, are mainly considered. The fracture statistics are shown in Table 1. According to the weathering characteristics and structural plane development in the power station area, the study area is divided into three different permeability zones in the profile, (I) strong-moderate weathering zone, (II) weak weathering zone and (III) slightly weathered zone.  There are 25 long-term observation holes for groundwater level in the study area, and the locations of some observation holes are shown in Figure 6. The groundwater levels of eight observation holes were selected to invert the model parameters, including four in the upper reservoir and its surroundings, two in the vicinity of the diversion tunnel  There are 25 long-term observation holes for groundwater level in the study area, and the locations of some observation holes are shown in Figure 6. The groundwater levels of eight observation holes were selected to invert the model parameters, including four in the upper reservoir and its surroundings, two in the vicinity of the diversion tunnel and underground powerhouse and two in the lower reservoir, basically covering the entire study area.
According to the established groundwater flow model, the finite element method was used for calculation. The calculated values in the observed borehole were obtained corresponding to 12 January 2004 ( Table 2). The objective function value, E, was 0.48, corresponding to the optimized solution of the hydrogeological parameters. The calculated water level of each borehole and the corresponding weight function value, measured value and the error between the calculated and measured values are shown in Table 2. It can be seen from Table 2 that the fitting errors between the observed and calculated water levels are generally small, most of the boreholes have high calculation accuracy and only a few points have large errors. Additionally, the largest errors are found in ZK7 and ZK14, where the calculated values are almost 2 m lower than the observed values. Due to the strongly heterogeneous anisotropy of the fractured rock, there are some errors between the calculated and observed values. The main values of the hydraulic conductivity tensor, the main directions and the inversion values of the average permeability coefficient in each permeability zone are shown in Table 3. The calculation results of permeability tensor, the main value of the tensor and the average value show that the trend of the parameter magnitudes in each permeability partition is basically consistent with the weathering of the fractured rock mass and the hydrogeological characteristics. Overall, the permeability parameters of the rock masses gradually decrease with increasing depth. The average hydraulic conductivity ranges from 0.308 m/d in the strongly and moderately weathered zone to 0.062 m/d in the slightly weathered rock mass.  In order to verify the correctness of the inversion model, the long-observed water levels of the borehole for a total of 12 time periods from December 2003 to April 2004 were selected and fitted to the calculated water levels. The permeability tensor obtained from the inversion was used to calculate the groundwater levels for each time period. Using time t as the horizontal coordinate and hydraulic head as the vertical coordinate, the observed water level H o and the calculated water level H c for the same observation hole were plotted on the graph to obtain a fitted graph (Figure 7). From the fitted graphs, it can be seen that the calculated water levels in each borehole fit well with the observed water levels. It shows that the inversion model can well reflect the actual conditions in the study area and that the values of the hydrogeological parameters obtained from the inversion are reasonable.

Simulation of 3D Groundwate Seepage in the Powerstation Area
(1) Calculation range The study area mainly comprises the upper storage reservoir, underground powerhouse and lower storage reservoir, extending to the watershed on both sides and to a depth of −30 m elevation. Based on the permeability zoning of the rock in the study area and the setting of impermeable curtains during the reservoir operation period, the study area was divided into four zones according to permeability, namely a strong-moderate weathering zone, weak weathering zone, slightly weathered zone, the upper and lower storage reservoir and underground powerhouse impermeable curtains. The hydraulic conductivity of the geotextile at the bottom of the upper storage reservoir is 1 × 10 −11 cm/s. The hydraulic conductivity of the reinforced concrete panel around the reservoir is 1 × 10 −7 cm/s. The hydraulic conductivity of the clay cover of the lower storage reservoir is taken as 1 × 10 −7 cm/s. The hydraulic conductivity values of the other three divisions are taken as shown in Table 3. Considering the topography, the arrangement of hydraulic buildings and the arrangement of drainage holes, etc., the study area was divided into a total of 1603 super units. Due to the different calculation accuracies required for different engineering parts, the super units were dissected with unequal density, with a total of 19,409 nodes and 17,100 units, and the mesh map is shown in Figure 8. Taking the south side of the area as the origin of the coordinates, the coordinate system with the x-axis pointing to true east, the y-axis pointing to true north and the z-axis pointing vertically upward was used as the calculation coordinate system. level Ho and the calculated water level Hc for the same observation hole were plotted on the graph to obtain a fitted graph (Figure 7). From the fitted graphs, it can be seen that the calculated water levels in each borehole fit well with the observed water levels. It shows that the inversion model can well reflect the actual conditions in the study area and that the values of the hydrogeological parameters obtained from the inversion are reasonable.

Simulation of 3D GroundwateS Seepage in the Powerstation Area
(1) Calculation range The study area mainly comprises the upper storage reservoir, underground powerhouse and lower storage reservoir, extending to the watershed on both sides and to a depth of −30 m elevation. Based on the permeability zoning of the rock in the study area and the setting of impermeable curtains during the reservoir operation period, the study area was divided into four zones according to permeability, namely a strong-moderate   and the arrangement of drainage holes, etc., the study area was divided into a total of 1603 super units. Due to the different calculation accuracies required for different engineering parts, the super units were dissected with unequal density, with a total of 19,409 nodes and 17,100 units, and the mesh map is shown in Figure 8. Taking the south side of the area as the origin of the coordinates, the coordinate system with the x-axis pointing to true east, the y-axis pointing to true north and the z-axis pointing vertically upward was used as the calculation coordinate system.  (2) Calculation schemes In order to know the groundwater movement pattern in the power station area during the operation period, the groundwater seepage field during the operation period was simulated. Considering the failure of the curtain due to chemical and physical effects during the operation period of the reservoir and the drainage of the underground powerhouse, simulations were made to calculate the failure of the impermeable curtain in the upper storage reservoir, the underground powerhouse and the lower storage reservoir during the operation period, and the specific calculation conditions and schemes are shown in Table 4.

Seepage Field Calculations and Evaluation of Results
(1) Comparison of seepage fields in normal and natural conditions Scheme f01 is the normal reservoir operation and effective anti-seepage of the upper and lower storage reservoirs and underground powerhouse as well as drainage. Scheme f02 does not consider any form of anti-seepage and drainage. It can be seen in Figure 9 that the groundwater level in front of the underground powerhouse under scheme f02 is about 130-160 m, about 80 m higher than under the drainage condition of scheme f01. In the underground powerhouse, the free surface of groundwater decreases globally due to the drainage of the underground powerhouse. The groundwater level is basically the same as the elevation of the bottom in the underground powerhouse ( Figure 10) and shows that the drainage system of the underground powerhouse is effective. The hydraulic gradients for the upper and lower storage reservoirs and the underground powerhouse are calculated separately in Table 5. For scheme f02, the average hydraulic gradients are 0.0015, 0.1476 and 0.0036 for the upper storage reservoir, underground powerhouse and lower storage reservoir, respectively. For scheme f01, the hydraulic gradients are 0.2032, 0.3308 and 0.0019, respectively, with an increase in the hydraulic gradient due to the presence of the underground powerhouse drainage holes. It can be seen in Table 5 that the underground powerhouse has the largest hydraulic gradient, followed by the upper storage reservoir, and the lower storage reservoir has the smallest hydraulic gradient. Generally, the allowable hydraulic gradient of RCC dams is 1.5-2.5, and the allowable hydraulic gradient of the anti-seepage curtain can reach 3-4. It can be seen from Table 5 that the calculated maximum hydraulic gradient does not exceed 1, which indicates that during the operation of the reservoir, the seepage prevention curtain will not be damaged. These values are acceptable.        (2) Leakage analysis The leakage was calculated for the upper storage reservoir, underground powerhouse and lower storage reservoir under various schemes. The results are shown in Table 6. During normal operation of the reservoir, scheme f01, the leakage from the upper, underground powerhouse and lower storage reservoirs are 710.48 m 3 /d, 969.95 m 3 /d and 1657.55 m 3 /d, respectively. The underground powerhouse anti-seepage fails, that is, when the hydraulic conductivity of the impermeable curtain is taken to be the hydraulic conductivity of the underground powerhouse surrounding rock, the underground powerhouse leakage is 8571.72 m 3 /d, approximately 8.84 times that under normal anti-seepage conditions, indicating that the underground powerhouse anti-seepage system can effectively reduce the drainage of the underground powerhouse. If the leakage is too large, it will seriously affect the benefits of reservoir power generation, especially pumped storage power stations. The water in the upper reservoir is extremely precious, and an appropriate leakage is acceptable. For example, the leakage in the upper reservoir is 710.48 m 3 /d for f01; the leakage accounts for 0.4 ten thousandth of the total storage capacity in the upper reservoir, which is less than one ten thousandth, indicating that the leakage value of the upper reservoir is acceptable.
The failure of the underground powerhouse anti-seepage has no effect on the upper storage reservoir. The leakage of the lower storage reservoir slightly increases, mainly due to the drainage system of the lower reservoir being very close to the underground powerhouse. The failure of the upper storage reservoir anti-seepage, scheme f04, results in a larger increase in leakage from the upper storage reservoir, greater than 35.65 times that of the normal condition. Leakage in the underground powerhouse increases by 2 times, and leakage in the lower storage reservoir has basically no effect. When the lower storage reservoir anti-seepage fails, scheme f05, the leakage in the lower storage reservoir increases by 11.5 times, and the leakage in the underground powerhouse increases by approximately three times. The upper storage reservoir has basically no effect. Obviously, the anti-seepage of the lower storage reservoir has a greater impact on the leakage of the underground powerhouse than the upper storage reservoir. When the upper storage reservoir is at a dead level of 239 m, the leakage in the upper storage reservoir is 241.03 m 3 /d, which is about three times less than when the storage level is 267 m. The leakage from the reservoir is positively correlated with the reservoir level, but the change in the upper storage reservoir level has little effect on the leakage in the underground powerhouse and lower storage reservoir. Therefore, even taking into account the failure of the anti-seepage, the change in water level in the upper storage reservoir mainly affects the leakage in the upper storage reservoir and has very little effect on the leakage in the underground powerhouse and lower storage reservoir. When the water level of the lower storage reservoir is the dead level of 65 m, the leakage of the lower storage reservoir is 1023.18 m 3 /d, reducing the leakage by about 640 m 3 /d, which is about 66.67% of the normal storage level of 81 m of scheme f01. Similarly, changes in the level of the lower storage reservoir have little effect on the leakage from the upper storage reservoir and the underground powerhouse, mainly affecting the leakage in the lower storage reservoir. Therefore, the impact of leakage in the underground powerhouse mainly comes from the upper storage reservoir, the lower storage reservoir and the underground powerhouse impermeable system, while changes in water levels in the upper and lower storage reservoirs have a minimal impact on the leakage in the underground powerhouse.
(3) Curtain failure and parameter sensitivity analysis The calculation schemes consider the failure of the impermeable curtain for the upper storage reservoir, the underground powerhouse and the lower storage reservoir separately, as well as simultaneous failures, mainly at 100%. It is impossible that the curtain will fail 100% during actual reservoir operation, so we also calculated a partial failure of the curtain, with the failure of hydraulic conductivity for the curtain calculated by the following equation: where w% is the percentage of curtain failure. Curtain failures of 100%, 80%, 50%, 30%, 10% and 0% were considered. It can be seen from the phreatic water contour diagram that the seepage field remains largely unchanged compared to normal operating conditions after the partial or complete failure of the curtain in the underground powerhouse area, but the hydraulic gradient before and after the curtain decreases with curtain failure. Since the hydraulic head of water behind the curtain is mainly controlled by the drainage holes, even if the curtain fails, the seepage field will not change much, but the amount of water discharged through the drainage hole will increase greatly (Table 6). When the curtain of the upper reservoir fails by 30%, the flow rate increases by 94.5%, and when the curtain fails by 50%, the flow rate increases by 212.9%. When the underground powerhouse curtain fails by 30%, the flow rate increases by 79.7%, and when the curtain fails by 50%, the flow rate increases by 163.7%. When the curtain of the lower reservoir fails by 30%, the flow rate increases by 84%, and when the curtain fails by 50%, the flow rate increases by 177.1%. It can be seen that the curtain plays a very important role in reducing the flow into the drainage hole.
In addition, the leakage at the percentage of curtain failure was calculated ( Figure 11). It can be seen from the figure that as the percentage of curtain failure increases, the leakage increases gradually, but not linearly, because the size of the curtain changes the size of the water level and thus affects the hydraulic gradient. For example, when the curtain of the underground powerhouse fails, the hydraulic gradient increases due to the presence of drainage. When the gradient of water flow in the upper and lower storage reservoirs decreases, the multiple of decrease of the hydraulic gradient is less than the multiple of increase in the curtain permeability coefficient. Therefore, the leakage increased. When the curtain fails by 10%, 30%, 50%, 80% and 100%, for the upper reservoir, the leakage increase is 1.24, 1.95, 3.13, 7.98 and 35.65 times, respectively. The leakage increase from the powerhouse is 1.22, 1.80, 2.64, 4.95 and 8.84 times, respectively. The increase in leakage of the lower reservoir is 1.22, 1.84, 2.77, 5.60 and 11.50 times, respectively. Owing to the important role of the anti-seepage curtain, during the construction of the hydropower station, it is recommended to strengthen the anti-seepage of each engineering site and conduct regular inspections to reinforce the failed anti-seepage curtain in time.

Conclusions
Through the analysis of hydrogeological conditions in the JPSHP and the results of 3D finite element simulation calculations of the seepage field, the following conclusions can be drawn: (1) The geological and hydrogeological conditions of the JPSHP determine that the permeability of the rock mass is macroscopically heterogeneous and anisotropic. The existence of carbonate and non-carbonate constitutes the heterogeneity of the medium. The coexistence of karst and fracture leads to the obvious anisotropy of the medium, and the weathering effect exacerbates the degree and complexity of the heterogeneous anisotropy of the medium. The rock is cut by faults and fractures to form an uneven fracture-karst aquifer. The permeability of the rock decreases with the

Conclusions
Through the analysis of hydrogeological conditions in the JPSHP and the results of 3D finite element simulation calculations of the seepage field, the following conclusions can be drawn: (1) The geological and hydrogeological conditions of the JPSHP determine that the permeability of the rock mass is macroscopically heterogeneous and anisotropic. The existence of carbonate and non-carbonate constitutes the heterogeneity of the medium. The coexistence of karst and fracture leads to the obvious anisotropy of the medium, and the weathering effect exacerbates the degree and complexity of the heterogeneous anisotropy of the medium. The rock is cut by faults and fractures to form an uneven fracture-karst aquifer. The permeability of the rock decreases with the increasing depth of burial and is stronger at the weathered and unloaded parts of the rock. The permeability of the rock is mainly controlled by the structural plane, and the vertical water permeability generally tends to be stronger than the horizontal water permeability. The fractures play the role of water conduction, and the joints in the lightly weathered and fresh rock mass mainly play the role of water storage. (2) A control inversion method for the structural plane of fractured rock mass is proposed.
For the fractured rock mass, the larger scale faults and fractures are used as waterconducting media, and a relatively complete rock matrix is used as water-storage media. The statistical method is used to obtain the distribution function of the structural plane of the fractured rock mass in the study area and the statistical values of the spatial geometric parameters, so as to obtain the hydraulic conductivity tensor of the fractured rock mass. (3) The inversion results of the hydrogeological parameters of the JPSHP area reflect the characteristics of anisotropic permeability of the rock mass. The main permeability coefficient K zz is the largest, followed by K yy , and K xx is the smallest, that is, the rock mass is permeable in the vertical direction. The permeability is obviously greater than the horizontal permeability, mainly because the dip angle of the structural plane is relatively steep, which is basically greater than the steep dip angle of 60