# An Estimation Model for Hydraulic Conductivity of Low-Permeability and Unfilled Fractured Granite in Underground Water-Sealed Storage Caverns

Authors: Yangbing Cao; Weiguo Gong; Xiangxiang Zhang; Junxi Chen; Zhenping Huang
Journal: Geofluids
 (2021)

---

## Abstract

The permeability of rock mass is closely related to the stability and safety of underground structure, especially in underground water-sealed storage caverns. With regard to the estimation approaches in predicting the hydraulic conductivity of fractured granite in water-sealed storage caverns, there are some limitations of parameter selection leading to poor applicability. Focusing on the contribution of the water conduction fractures (WCF) to the hydraulic conductivity, we attempted to propose a novel model, the CA model, for estimating its hydraulic conductivity based on the fracture orientation index and the normal stress index by analyzing the borehole wall imaging results and borehole water-pressure test results in the site of underground water-sealed storage caverns. The results indicated that the proposed model is suitable for low-permeability and unfilled fractured granite, exhibiting good effectiveness by clarifying the relation between geomechanical parameters and hydraulic behavior. Further, the parameters upon which the proposed model is based are representative and easy to obtain, which has certain guiding significance and reference value for analyzing the permeability characteristics of similar rock masses.


---

## Body

## 1. Introduction

At present, many underground water-sealed storage caverns are located in the low-permeability and unfilled fractured granite area. A large number of random joints result in distinct heterogeneity and anisotropy in the permeability of fractured granite. With the abovementioned features of granite, most underground water-sealed storage caverns for oil or liquefied petroleum gas make the groundwater flow in complex fracture networks, forming a stable seepage field by constructing a water curtain system, so as to achieve the long-term storage of oil or liquefied petroleum gas. However, groundwater seepage may not only weaken the granite strength but form local high-pressure seepage resulted from water-sealed conditions, resulting in some geological disasters such as block falling or block sliding.Therefore, the groundwater seepage characteristics of fractured granite in these water-sealed storage cavern areas become a critical issue to study the stability of surrounding granite and the evolution of seepage field. Specifically, the construction of the seepage medium model, as well as the determination of hydraulic conductivity, is the most essential question. Currently, the seepage medium model of the rock mass can be divided into three categories: the equivalent continuum medium model [1–3], double-medium model [4, 5], and discrete fracture network (DFN) model [6–8]. However, the equivalent continuum model is only applicable to large-scale loose or broken rock with a representative elementary volume (REV). In addition to a small-scale REV, the input parameters of the dual medium model are complex and possess a wide range of assumptions. Meanwhile, the two types of continuous medium models do not involve the specific spatial distribution of the water conduction fractures (WCF), and they are difficult to characterize the heterogeneity and discontinuity of seepage in the rock mass, making them unsuitable for the underground water-sealed storage caverns. Note that the WCF is a class of rock mass fracture belonging to the primary seepage zone. Generally, primary seepage zones are some fractures with high permeability or distinct systematic distributions as well as a combination of the two in rock mass. The DFN model assumes that flow and solute transport occur primarily within fractures [9]. Therefore, by accurately describing the geometric parameters (orientation, spacing and aperture, etc.) and hydraulic parameters (hydraulic aperture, etc.) of each fracture in the fracture network, the geometric and hydraulic characteristics of the WCF can be described. Thus, the heterogeneity and seepage anisotropy of the fractured rock mass can be determined. This model has been widely used in the stability and seepage analysis of high and steep rock slopes and underground caverns with small research areas and requires high research degree and calculation accuracy [10, 11]. In summary, the fractured granite of the underground water-sealed storage cavern can be examined as a fractured seepage network comprising WCF, enabling seepage calculations to be conducted using the DFN model with a calibrated hydraulic conductivity.Methods for determining the hydraulic conductivity of fractured rock mass mainly include the field hydraulic test method, inverse analysis method, and empirical estimation method. The field hydraulic test method [12, 13] obtains the hydraulic conductivity of rock mass via an analytical calculation of water-pressure test and pumping test data in situ. However, because of its high cost as well as technology and time requirements, its precision and scope are limited, and it is difficult to reflect the hydraulic conductivity characteristics of the rock mass of the entire study area. Meanwhile, the inverse analysis method [14, 15] is based on the monitoring data of groundwater flow rate or velocity. The hydraulic conductivity of a rock mass is calculated via analytical or numerical methods, but it often provides a no-unique solution, which affects the accuracy. Finally, the empirical estimation methods mainly include the fracture network measurement method, fracture network numerical test method, and geological index estimation method. Among these, in the fracture network measurement method [16–18], the geometric parameters of fractures in the measurement and statistic area are substituted into the permeability tensor formula to obtain the initial value of the equivalent permeability tensor, which must then be corrected using the hydraulic test method or inverse analysis method. The fracture network numerical test method [19–21] was proposed by Long et al. [19] to determine the REV and equivalent permeability tensor of rock mass. However, the understanding of the existence and size effect of the REV is not uniform, which makes it difficult to establish the numerical model.The geological index estimation method can be divided into single index and comprehensive index estimation methods. Currently, many scholars use a single index, such as depth [22, 23], P-wave velocity (Vp) [24], geoelectrical parameters [25], or rock quality designation (RQD) [26, 27], to establish an empirical formula to predict the hydraulic conductivity of rock mass. Moreover, existing estimation models based on the comprehensive index have also made great progress. Hsu et al. [28, 29] defined a comprehensive index HC that considers the RQD, depth index (DI), lithology permeability index (LPI), and gouge content designation (GCD) index and estimated the hydraulic conductivity of fractured sedimentary rock by establishing a power function relationship between the hydraulic conductivity (K) and HC index. Song et al. [30] proposed the RMP index, which considers the RQD, rock integrity designation (RID), fracture aperture designation (AD), and LPI and established a power function relationship to estimate the hydraulic conductivity of the underground water-sealed caverns for liquefied petroleum gas (LPG) engineering site. Chen et al. [31] developed the ZRF model, which considers the buried depth (Z), RQD, and GCD, in order to estimate the hydraulic conductivity of the granite in a hydropower station. Thus, geological index estimation models are simple and practical and have a certain application value in stability and permeability research of rock mass engineering. However, as geological indices of rock mass are adopted mainly to characterize rock mass mechanics, the indices may not be entirely applicable to hydraulic behavior [32]. Therefore, the proposed models lack key parameters, such as the fracture orientation and hydraulic aperture, which control the hydraulic conductivity of rock mass and cannot form a connection between the seepage model of the fractured media and permeability of the fractured rock mass.Based on the field test results of borehole water-pressure tests and borehole wall imaging tests of rock mass in a underground water-sealed storage cavern, this paper discusses the existing estimation indices and models and innovatively introduces geometric characteristic parameters such as the fracture orientation, number, and aperture of the WCF. Finally, we proposed a new model for estimating hydraulic conductivity that uses normal stress on the WCF and fracture orientation as its basic parameters. In addition, this model examines the influence of RQD, fracture spacing, and density on the permeability of the rock mass by considering the number of fractures. Based on the estimation results, it is proven that the proposed model can accurately reflect the permeability characteristics of fractured granite in the study area, thereby providing a reference for estimating the hydraulic conductivity and engineering application of other rock masses.


## 2. Hydraulic Conductivity Distribution Characteristics of Rock Mass in the Underground Water-Sealed Storage Caverns

### 2.1. Geological and Hydrogeological Setting

In this study, the underground water-sealed storage cavern area is situated in the eastern coastal area of Shandong Province, China. The research object of this study is a water-sealed storage cavern for propane on the east side of the project area. As shown in Figure1, the study area is located in a transition zone between hills and alluvial plains. The altitude of the ground surface ranges from 15 m to 40 m, the topographic slope gradient is less than 5°, and the overall dip direction of the slope is westward. Because of the construction of ancillary facilities on the ground, the surface morphology is required to create a flat with the final elevation of approximately 30 m. In the study area, fine-grained monzonitic granite and marble occurred in the northeast of it, and the surface is covered with thin Quaternary overburden. The grained monzonitic granite was formed in the early stages of the Yanshan period which has a blocky structure and is generally medium to coarse grained. It is relatively hard rock accounting for more than 80% of the cavern rock mass. The rock mass of the slightly weathering layer has high strength, and the fractures range from extremely low developed to low developed. This is the main rock constructing the cavern, and the basic quality grade of the rock mass is II-III. The cavern construction is greatly influenced by secondary faults and joints distributed in the engineering area. According to the survey results of faults and discontinuities, the faults within the cavern site mainly include fault F9 and fracture zone P9. The remaining faults or fracture zones are far from the caverns and are small scale, barely influencing the cavern construction.Figure 1
Geological and hydrogeological setting: (a) the regional area and (b) the study area.The hydrogeological survey revealed that the groundwater is mainly recharged by the vertical infiltration of atmospheric precipitation and lateral groundwater in mountainous areas, whereas it is mostly discharged by surface runoff with small infiltration. As illustrated in Figure1, groundwater flows from higher altitudes areas of the eastern mountains to the western and northern regions and finally into the adjacent river and ocean. According to the laboratory and in situ test results, the granite matrix has very low permeability and porosity. Because of the development of joints, the low-permeability value of the rock matrix, and the large hydraulic gradient in the slightly weathered granite, groundwater mainly flows through the complex fracture network. Moreover, it appears to be the primary pathway in some broken areas and fracture concentrated belts. Therefore, the hydrogeological conditions of the study area highly obey the assumption of the fracture network model.In the fresh or slightly weathered granite, the geostress was measured in two vertical boreholes using the hydraulic fracturing method. The classical theory of hydraulic fracturing was proposed by Hubbert and Willis [33]. To date, it has been deeply developed in theory and widely used in the geostress measurement and the stimulation technique of oil and gas reservoirs in the petroleum industry [34–36]. The conditions of this study area obey the assumptions of the hydraulic fracturing classical theory, including the test section of the rock which is homogeneous, isotropic, linear elastic, intact, and impermeable, and there is one principal stress direction parallel to the borehole axis. Because the final altitude of the ground surface will be flat, it is assumed to be a plane with a 30 m altitude in the following regression analysis. Therefore, the calculation of the geostress can be transformed to a planar stress issue with a circular hole in an infinite plate [34]. In addition, it was predetermined that the tensile stress is positive, and the compressive stress is negative. The maximum and minimum horizontal principal stresses are, respectively, 13.71 MPa and 8.05 MPa in the corresponding design altitude range of the cavern floor. Thus, it can be judged that the site is an area of medium and low in situ stress in accordance with the Chinese standard for engineering classifications of rock mass [37]. After tensor transformation, tensor averaging, and a linear regression analysis, the direction of the maximum principal stress was determined to be NE78° along the horizontal plane. Assuming that the weight density of the rock mass was 27 kN/m3, the fitting relationships between the principal stresses and altitude are as follows:
(1)σH=0.0361h−6.0651,σh=0.0206h−3.9489,σv=0.027h−0.81,where σH, σh, and σv are the maximum principal stress (MPa), intermediate principal stress (MPa), and minimum principal stress (MPa), respectively, and h is the altitude (m).Borehole wall imaging tests were conducted in 28 boreholes using borehole televiewer technology and the latest digital panoramic borehole camera system (DPBCS) [38]. The DPBCS were employed as follows: a digital high-definition camera device in front of the probe can record a continuous, magnetically orientated, digital, 360° color image of the borehole wall. The depth and plane orientation are measured by the wheel and magnetic electronic compass of the DPBCS, respectively, and the annular borehole wall image is obtained. Then, in the flattened pattern, the annular borehole wall image is converted into a two-dimensional panoramic image in the order of N-E-S-W-N. Using the established coordinate system, it is reduced to a real borehole wall image. Finally, image interpretation software can be used to process the image to obtain the orientation, location, geometrical aperture, spacing, and filling content. After completing abovementioned steps, it was found that the fractures of the slightly weathered rock mass were basically closed, and there was no filled medium inside. However, not all the data from DPBCS and interpretation software are available, and geometrical aperture is also difficult to accurately process as the use of rigorous mathematical formulas and interpretation algorithms produce errors. Therefore, this study mainly conducted permeability analyses and hydraulic conductivity estimations based on the fracture orientation, number, spacing, and other geometric parameters.


### 2.2. Analysis of Borehole Water-Pressure Test Data

In accordance with the Chinese code of water-pressure tests in boreholes for water resources and hydropower engineering (CCWPT) [39], a total of 609 test sections were then subjected to double-packer water-pressure tests adopted intelligent water-pressure testers. Each test section is divided into three levels of pressure and three or five stages owing to the depth of the test section. Figure 2 shows the basic curve between the test pressure (P) and input flow (Q) of the test section when five stages existed. According to the CCWPT, if the absolute differences in the flow rate between stages 4 and 2 (or between stages 3 and 1 if there are only three stages) and between stages 5 and 1 are less than 1 L/min, or the relative errors are less than 5%, the boost and the buck phases in the P-Q curve can be considered to be basically coincident. In this study, the flow rate of every stage was invariably less than 1 L/min; so, the absolute difference was less than 1 L/min. Thus, the P-Q curves conform to the type of laminar flow, meaning that the fracture water flowing in the low-permeability and unfilled fractured granite was considered to be laminar flow.Figure 2
Diagram of stable flow rate and testing pressure variation at the five test stages.Owing to the fact that the groundwater flow through the fracture was laminar, the Lugeon value of the rock mass in the test sections can be calculated using formula (2), except for the test sections with high permeability around the fracture zone or a joint dense development zone.
(2)q=QmaxL∙Pmax,where q (Lu) is the Lugeon value, L (m) is the length of the test section, Qmax (L/min) is the maximum flow rate for each test stage (the third stage if five stages exist and the second stage if three stages exist), and Pmax (MPa) is the maximum testing pressure for each test stage (the third stage if five stages exist and the second stage if three stages exist).The distribution of the Lugeon values in the fresh or slightly weathered granite is illustrated in Figure3. A three-parameter Weibull distribution was used to fit the results, revealing a shape parameter of 0.80, a scale parameter of 0.061, and a location parameter of 0.001. In addition, the results show that the permeability of the cavern site is mainly distributed in 0.001–0.20 Lu, with an average of 0.082 Lu and a maximum of 1.098 Lu. Therefore, according to the CCWPT, the rock mass in the study area is mainly low permeability or extremely low permeability.Figure 3
Distribution characteristics of Lugeon values in fresh or slightly weathered granite.According to the method suggested in the CCWPT, when the test section is below the groundwater level,q<10 Lu and the fracture flow are laminar. Thus, the equivalent hydraulic conductivity of the rock mass in the test section can be calculated using formula (3), which was proposed by Hvorslev [40]. Figure 4 illustrates the distribution characteristics of the hydraulic conductivity of 599 test sections in the fresh or slightly weathered granite area (without considering the test sections with larger K values around the fracture zone or joint dense development zone). A three-parameter Weibull distribution was used to fit the data. Based on the fit, the shape parameter was found to be 0.851, the scale parameter was 7.2×10−4, and the location parameter was 3.6×10−6. After obtaining the distribution model and statistical characteristics of the random variables, the Monte Carlo method can be used for generation in computer software [34].
(3)K=Q2πLHlnL2r+1+L2r2,where K is the hydraulic conductivity, Q is the flow rate, r is the radius of the borehole, H is the testing pressure head, and L is the same as above.Figure 4
Distribution characteristics of hydraulic conductivity in fresh or slightly weathered granite.


### 2.3. Applicability of the Single Index Estimation Method of Hydraulic Conductivity

At present, many empirical models use a single index, such as depth (Z or DI), P-wave velocity of rock mass (Vp), or rock quality designation (RQD), to estimate the hydraulic conductivity of rock mass. These indices are mostly based on the relationship between the degree of fracture development and the corresponding index. Then, we explored whether the above indices are suitable for estimating the hydraulic conductivity of fractured granite in the study area.According to the results of the water-pressure tests, Figure5 demonstrated the distribution of hydraulic conductivity of rock mass with depth. It can be seen that the hydraulic conductivity is mainly concentrated in the range of 10−5–10−2 m/d, and it is irregular at each depth, reaching even three orders of magnitude differences at the same depths. Further, it was found that with increasing depth, the hydraulic conductivity of each borehole exhibited the following changes:
(1)
In a small portion of the boreholes, with increasing depth, the hydraulic conductivity of the rock mass remained basically unchanged, only experiencing a very slight fluctuation(2)
In a small portion of the boreholes, with increasing depth, the hydraulic conductivity of the rock mass remained basically unchanged, and only in dense local fractures or large fractures did the hydraulic conductivity increase sharply(3)
In most boreholes, the hydraulic conductivity of the rock mass changes irregularly with increasing depth, and the relationship between the hydraulic conductivity and depth is disorderedFigure 5
Scatter diagram of hydraulic conductivity distribution in study area.The permeability of the granite at the site is closely related to the distribution characteristics of the fractures and the integrity of the rock mass, which can be characterized by the RQD and RID, in which the RQD refers to the cumulative length of core pieces longer than 100 mm in a run (RS) divided by the total length of the core run (RT) under standard drilling processes, and RID is obtained using the following:
(4)RID=VpVr2,where Vp is the elastic P-wave velocity of the rock mass, and Vr is the elastic P-wave velocity of the rock block. As RQD and RID are both, to a certain extent, influenced by development density, size, and aperture of fractures, the relationship between rock permeability and the fracture distribution and rock integrity can be described using RQD and Vp.Based on the core logging results of the boreholes and the test results of the P-wave velocity, the average P-wave velocity of the rock mass in each borehole was found to be between 4.8 km/s and 5.2 km/s. Further, it was found that most boreholes are similar to theK-RQD and K-Vp curves of ZK55 shown in Figures 6 and 7, respectively, which both exhibit no obvious regularity and large dispersion.Figure 6
Relationship curve between hydraulic conductivity (K) and rock quality designation (RQD).Figure 7
Relationship curve between hydraulic conductivity (K) and P-wave velocity (Vp).In summary, the single index method of estimating hydraulic conductivity, including the index of depth,RQD, and P-wave velocity, is not suitable to directly estimate the hydraulic conductivity of the rock mass in the study area. Therefore, based on the paradigm of the comprehensive index estimation method, a new estimation model to estimate the hydraulic conductivity of low-permeability and unfilled fractured granite in underground water-sealed storage caverns must be proposed considering the comprehensive effect of multiple indices.


## 2.1. Geological and Hydrogeological Setting

In this study, the underground water-sealed storage cavern area is situated in the eastern coastal area of Shandong Province, China. The research object of this study is a water-sealed storage cavern for propane on the east side of the project area. As shown in Figure1, the study area is located in a transition zone between hills and alluvial plains. The altitude of the ground surface ranges from 15 m to 40 m, the topographic slope gradient is less than 5°, and the overall dip direction of the slope is westward. Because of the construction of ancillary facilities on the ground, the surface morphology is required to create a flat with the final elevation of approximately 30 m. In the study area, fine-grained monzonitic granite and marble occurred in the northeast of it, and the surface is covered with thin Quaternary overburden. The grained monzonitic granite was formed in the early stages of the Yanshan period which has a blocky structure and is generally medium to coarse grained. It is relatively hard rock accounting for more than 80% of the cavern rock mass. The rock mass of the slightly weathering layer has high strength, and the fractures range from extremely low developed to low developed. This is the main rock constructing the cavern, and the basic quality grade of the rock mass is II-III. The cavern construction is greatly influenced by secondary faults and joints distributed in the engineering area. According to the survey results of faults and discontinuities, the faults within the cavern site mainly include fault F9 and fracture zone P9. The remaining faults or fracture zones are far from the caverns and are small scale, barely influencing the cavern construction.Figure 1
Geological and hydrogeological setting: (a) the regional area and (b) the study area.The hydrogeological survey revealed that the groundwater is mainly recharged by the vertical infiltration of atmospheric precipitation and lateral groundwater in mountainous areas, whereas it is mostly discharged by surface runoff with small infiltration. As illustrated in Figure1, groundwater flows from higher altitudes areas of the eastern mountains to the western and northern regions and finally into the adjacent river and ocean. According to the laboratory and in situ test results, the granite matrix has very low permeability and porosity. Because of the development of joints, the low-permeability value of the rock matrix, and the large hydraulic gradient in the slightly weathered granite, groundwater mainly flows through the complex fracture network. Moreover, it appears to be the primary pathway in some broken areas and fracture concentrated belts. Therefore, the hydrogeological conditions of the study area highly obey the assumption of the fracture network model.In the fresh or slightly weathered granite, the geostress was measured in two vertical boreholes using the hydraulic fracturing method. The classical theory of hydraulic fracturing was proposed by Hubbert and Willis [33]. To date, it has been deeply developed in theory and widely used in the geostress measurement and the stimulation technique of oil and gas reservoirs in the petroleum industry [34–36]. The conditions of this study area obey the assumptions of the hydraulic fracturing classical theory, including the test section of the rock which is homogeneous, isotropic, linear elastic, intact, and impermeable, and there is one principal stress direction parallel to the borehole axis. Because the final altitude of the ground surface will be flat, it is assumed to be a plane with a 30 m altitude in the following regression analysis. Therefore, the calculation of the geostress can be transformed to a planar stress issue with a circular hole in an infinite plate [34]. In addition, it was predetermined that the tensile stress is positive, and the compressive stress is negative. The maximum and minimum horizontal principal stresses are, respectively, 13.71 MPa and 8.05 MPa in the corresponding design altitude range of the cavern floor. Thus, it can be judged that the site is an area of medium and low in situ stress in accordance with the Chinese standard for engineering classifications of rock mass [37]. After tensor transformation, tensor averaging, and a linear regression analysis, the direction of the maximum principal stress was determined to be NE78° along the horizontal plane. Assuming that the weight density of the rock mass was 27 kN/m3, the fitting relationships between the principal stresses and altitude are as follows:
(1)σH=0.0361h−6.0651,σh=0.0206h−3.9489,σv=0.027h−0.81,where σH, σh, and σv are the maximum principal stress (MPa), intermediate principal stress (MPa), and minimum principal stress (MPa), respectively, and h is the altitude (m).Borehole wall imaging tests were conducted in 28 boreholes using borehole televiewer technology and the latest digital panoramic borehole camera system (DPBCS) [38]. The DPBCS were employed as follows: a digital high-definition camera device in front of the probe can record a continuous, magnetically orientated, digital, 360° color image of the borehole wall. The depth and plane orientation are measured by the wheel and magnetic electronic compass of the DPBCS, respectively, and the annular borehole wall image is obtained. Then, in the flattened pattern, the annular borehole wall image is converted into a two-dimensional panoramic image in the order of N-E-S-W-N. Using the established coordinate system, it is reduced to a real borehole wall image. Finally, image interpretation software can be used to process the image to obtain the orientation, location, geometrical aperture, spacing, and filling content. After completing abovementioned steps, it was found that the fractures of the slightly weathered rock mass were basically closed, and there was no filled medium inside. However, not all the data from DPBCS and interpretation software are available, and geometrical aperture is also difficult to accurately process as the use of rigorous mathematical formulas and interpretation algorithms produce errors. Therefore, this study mainly conducted permeability analyses and hydraulic conductivity estimations based on the fracture orientation, number, spacing, and other geometric parameters.


## 2.2. Analysis of Borehole Water-Pressure Test Data

In accordance with the Chinese code of water-pressure tests in boreholes for water resources and hydropower engineering (CCWPT) [39], a total of 609 test sections were then subjected to double-packer water-pressure tests adopted intelligent water-pressure testers. Each test section is divided into three levels of pressure and three or five stages owing to the depth of the test section. Figure 2 shows the basic curve between the test pressure (P) and input flow (Q) of the test section when five stages existed. According to the CCWPT, if the absolute differences in the flow rate between stages 4 and 2 (or between stages 3 and 1 if there are only three stages) and between stages 5 and 1 are less than 1 L/min, or the relative errors are less than 5%, the boost and the buck phases in the P-Q curve can be considered to be basically coincident. In this study, the flow rate of every stage was invariably less than 1 L/min; so, the absolute difference was less than 1 L/min. Thus, the P-Q curves conform to the type of laminar flow, meaning that the fracture water flowing in the low-permeability and unfilled fractured granite was considered to be laminar flow.Figure 2
Diagram of stable flow rate and testing pressure variation at the five test stages.Owing to the fact that the groundwater flow through the fracture was laminar, the Lugeon value of the rock mass in the test sections can be calculated using formula (2), except for the test sections with high permeability around the fracture zone or a joint dense development zone.
(2)q=QmaxL∙Pmax,where q (Lu) is the Lugeon value, L (m) is the length of the test section, Qmax (L/min) is the maximum flow rate for each test stage (the third stage if five stages exist and the second stage if three stages exist), and Pmax (MPa) is the maximum testing pressure for each test stage (the third stage if five stages exist and the second stage if three stages exist).The distribution of the Lugeon values in the fresh or slightly weathered granite is illustrated in Figure3. A three-parameter Weibull distribution was used to fit the results, revealing a shape parameter of 0.80, a scale parameter of 0.061, and a location parameter of 0.001. In addition, the results show that the permeability of the cavern site is mainly distributed in 0.001–0.20 Lu, with an average of 0.082 Lu and a maximum of 1.098 Lu. Therefore, according to the CCWPT, the rock mass in the study area is mainly low permeability or extremely low permeability.Figure 3
Distribution characteristics of Lugeon values in fresh or slightly weathered granite.According to the method suggested in the CCWPT, when the test section is below the groundwater level,q<10 Lu and the fracture flow are laminar. Thus, the equivalent hydraulic conductivity of the rock mass in the test section can be calculated using formula (3), which was proposed by Hvorslev [40]. Figure 4 illustrates the distribution characteristics of the hydraulic conductivity of 599 test sections in the fresh or slightly weathered granite area (without considering the test sections with larger K values around the fracture zone or joint dense development zone). A three-parameter Weibull distribution was used to fit the data. Based on the fit, the shape parameter was found to be 0.851, the scale parameter was 7.2×10−4, and the location parameter was 3.6×10−6. After obtaining the distribution model and statistical characteristics of the random variables, the Monte Carlo method can be used for generation in computer software [34].
(3)K=Q2πLHlnL2r+1+L2r2,where K is the hydraulic conductivity, Q is the flow rate, r is the radius of the borehole, H is the testing pressure head, and L is the same as above.Figure 4
Distribution characteristics of hydraulic conductivity in fresh or slightly weathered granite.


## 2.3. Applicability of the Single Index Estimation Method of Hydraulic Conductivity

At present, many empirical models use a single index, such as depth (Z or DI), P-wave velocity of rock mass (Vp), or rock quality designation (RQD), to estimate the hydraulic conductivity of rock mass. These indices are mostly based on the relationship between the degree of fracture development and the corresponding index. Then, we explored whether the above indices are suitable for estimating the hydraulic conductivity of fractured granite in the study area.According to the results of the water-pressure tests, Figure5 demonstrated the distribution of hydraulic conductivity of rock mass with depth. It can be seen that the hydraulic conductivity is mainly concentrated in the range of 10−5–10−2 m/d, and it is irregular at each depth, reaching even three orders of magnitude differences at the same depths. Further, it was found that with increasing depth, the hydraulic conductivity of each borehole exhibited the following changes:
(1)
In a small portion of the boreholes, with increasing depth, the hydraulic conductivity of the rock mass remained basically unchanged, only experiencing a very slight fluctuation(2)
In a small portion of the boreholes, with increasing depth, the hydraulic conductivity of the rock mass remained basically unchanged, and only in dense local fractures or large fractures did the hydraulic conductivity increase sharply(3)
In most boreholes, the hydraulic conductivity of the rock mass changes irregularly with increasing depth, and the relationship between the hydraulic conductivity and depth is disorderedFigure 5
Scatter diagram of hydraulic conductivity distribution in study area.The permeability of the granite at the site is closely related to the distribution characteristics of the fractures and the integrity of the rock mass, which can be characterized by the RQD and RID, in which the RQD refers to the cumulative length of core pieces longer than 100 mm in a run (RS) divided by the total length of the core run (RT) under standard drilling processes, and RID is obtained using the following:
(4)RID=VpVr2,where Vp is the elastic P-wave velocity of the rock mass, and Vr is the elastic P-wave velocity of the rock block. As RQD and RID are both, to a certain extent, influenced by development density, size, and aperture of fractures, the relationship between rock permeability and the fracture distribution and rock integrity can be described using RQD and Vp.Based on the core logging results of the boreholes and the test results of the P-wave velocity, the average P-wave velocity of the rock mass in each borehole was found to be between 4.8 km/s and 5.2 km/s. Further, it was found that most boreholes are similar to theK-RQD and K-Vp curves of ZK55 shown in Figures 6 and 7, respectively, which both exhibit no obvious regularity and large dispersion.Figure 6
Relationship curve between hydraulic conductivity (K) and rock quality designation (RQD).Figure 7
Relationship curve between hydraulic conductivity (K) and P-wave velocity (Vp).In summary, the single index method of estimating hydraulic conductivity, including the index of depth,RQD, and P-wave velocity, is not suitable to directly estimate the hydraulic conductivity of the rock mass in the study area. Therefore, based on the paradigm of the comprehensive index estimation method, a new estimation model to estimate the hydraulic conductivity of low-permeability and unfilled fractured granite in underground water-sealed storage caverns must be proposed considering the comprehensive effect of multiple indices.


## 3. Hydraulic Conductivity Estimation Model of Low-Permeability and Unfilled Fractured Granite

### 3.1. Model Construction and Parameter Selection

The permeability of fracture media in the DFN depends almost entirely on the spatial structure of the fracture system; meaning, the relationship between the fracture properties and permeability can be clarified by establishing a simple medium model of fractured rock mass. Snow [16] assumed that there was a water-conducting fracture network composed of n groups of directional fractures in a medium of pure fractured rock mass, and that the fracture water took a laminar flow. Therefore, the permeability tensor of each group of the WCF was superimposed based on the cubic law of the parallel plate model, and considering geometric factors, such as the fracture orientation, spacing, and hydraulic aperture, the following permeability tensor calculation formula was proposed, which is referred to as the snow model:
(5)KSNOW=g12v∑i=1nbi3Si1−nxi2−nxinyi−nxinzi−nyinxi1−nyi2−nyinzi−nzinxi−nzinyi1−nzi2,where KSNOW (m/s) is the hydraulic conductivity tensor of the rock mass, g (m/s2) is the gravitational acceleration, ν (m2/s) is the kinematic viscosity, bi (m) is the average hydraulic aperture of group i, Si (m) is the space of group i, nxi, nyi, and nzi are components in the x, y, and z directions, respectively, of the unit normal vector of group i fracture, and n is the number of fractures.Snow [16] considered the number of fractures to reflect the spacing and density of the fractures by superposing the hydraulic conductivity tensor of all the WCF, which is actually the contribution of rock integrity parameters to the hydraulic conductivity. Further, the snow model is suitable for laminar flow and clean fractures, which are consistent with the characteristics of the granite in the study area. Therefore, the snow model can be simplified into two parts: the orientation of all fractures and the hydraulic aperture. According to the superposition principle of the hydraulic conductivity tensor and considering the number of fractures, the fracture orientation index and hydraulic aperture index were tentatively selected to establish a new model for predicting the hydraulic conductivity of the low-permeability and unfilled fractured granite in underground water-sealed storage caverns. The hydraulic aperture is mainly controlled by normal stress, contact area, geometrical aperture, and roughness of the fracture, and it is difficult to determine the hydraulic aperture directly. Previous relevant studies [34, 41, 42] show that the quantitative relationship exists distinctly between normal stress and fracture hydraulic aperture. Fortunately, the normal stress of fractures is relatively easy to obtain in engineering; so, it was selected to replace the index of hydraulic aperture for the construction of the estimation model. Therefore, the normal stress, as a comprehensive index, combined with the orientation index of the discontinuity can comprehensively characterize the hydraulic behavior of the fracture.The hydraulic conductivity tensor of the rock mass is based on the permeability characteristics of the WCF. However, according to the results of the borehole water-pressure tests and borehole wall imaging tests in the study area, it was found that not all fractures are the WCF, and most have no seepage significance. As the orientation of the WCF determines the primary seepage direction of the rock mass, the fracture network model needs to be characterized emphatically. Therefore, in the case of existing fracture orientation index, the proposed model introduces orientation index of the WCF to characterize the contribution of the hydraulic characteristics of the WCF to the permeability of the rock mass.


### 3.2. Novel Model for Estimating Hydraulic Conductivity of Fractured Granite in Water-Sealed Storage Caverns

Based on the snow model and fracture network model elements, this study used the fracture orientation (denoted asA or Aj) and the normal stress (σn) of the WCF as the basic index to characterize the geometric characteristics of the fracture and the hydraulic properties of the WCF, respectively, and proposed an estimation model for hydraulic conductivity of the low-permeability and unfilled fractured granite in underground water-sealed storage caverns, called the CA model:
(6)K=α∙A+β∙Cwcf+γ,where α, β, and γ are empirical constants, and K (cm/s) is the hydraulic conductivity of each test section. The index A represents the contribution of the orientation of all fractures to the hydraulic conductivity, and Cwcf is an integrated hydraulic index representing the contribution of all the WCF to hydraulic conductivity of rock mass. Therefore, the name of the CA model is taken from the main indices Cwcf and A. In terms of the unit of two parameters, the index A proposed in this study is dimensionless, and the unit of the index Cwcf is MPa-3. Although the unit is different, the model only uses its numerical value.
(1)
The method for obtaining the indexA. It is assumed that the contribution of the orientation of n fractures to the hydraulic conductivity of the rock mass at a certain depth range is An. Based on the principle of permeability tensor superposition, the following formula can be obtained:(7)An=∑i=1nAi=∑i=1n1−nxi2−nxinyi−nxinzi−nyinxi1−nyi2−nyinzi−nzinxi−nzinyi1−nzi2.If the dip direction of fracture i is αi, the dip angle is βi, the x-axis points to the east, the y-axis points to the north, the z-axis points upwards, and a space rectangular coordinate system ENZ(xyz) is constructed. As shown in Figure 8, the x, y, and z components of the normal vector of the ID i of the fracture are expressed as follows:
(8)nxi=sinαi∙sinβinyi=cosαi∙sinβinzi=cosβi.Figure 8
Schematic diagram of unit normal vector of fracture orientation in geographic coordinate system.The eigenvalues ofAn was calculated and denoted as A1, A2, and A3, of which the three eigenvectors are the principal directions of the eigenvalues. The fracture orientation index A is defined as the geometric mean value of the eigenvalue An, as follows:
(9)A=A1∙A2∙A33,where the physical meaning of A is the contribution of the orientation of all fractures to the hydraulic conductivity in the depth range of the test section. Therefore, the orientation index of fractures can be obtained by simultaneous formulas (7), (8), and (9).
(2)
The method for obtaining the indexCwcf. To obtain Cwcf, the identification and selection of WCF should be conducted first to determine the mathematical relationship between the hydraulic parameters of the WCF and the hydraulic conductivity of the rock mass. Firstly, based on borehole wall image tests and water-pressure tests data, the hydraulic conductivity of each test section and internal fractures orientation was obtained. Secondly, the orientation ranges of fractures with a small hydraulic conductivity in the test section need to be screened out, and it was considered that it does not belong to the primary seepage zone. Finally, it was removed from the orientation range of the test sections with large hydraulic conductivity, and what were retained were the WCF and its orientation in each test section. Based on the DFN model of a real discrete fracture network generated using a Monte Carlo simulation by Min et al. [6], Figure 9 represents the WCF selected from fracture network and internal groundwater seepageFigure 9
The WCF and the direction of groundwater flow in fracture system (the blue line represents the WCF, and the black line represents ordinary discontinuities).Although we have selected the normal stress index (σn) instead of the hydraulic aperture (eh) of the WCF to build the new model, it still plays a role as a bond between the hydraulic conductivity of rock mass and normal stress on the fracture. Based on the research results of Cao et al. [34] on the fracture hydraulic aperture of slightly weathered granite which located nearby the study area, it was assumed that a rock mass test section with upper and lower surfaces parallel to its contained fracture and with a distance of L between the two interfaces. Other parameters were kept the same as previously mentioned. The hydraulic conductivity Kr of the fractured rock mass is the superposition of the hydraulic conductivity of each fracture in the test section as follows:
(10)Kr=∑i=1ngeh312vL·Ai.Thus, the basic relationship between hydraulic aperture and normal stress was obtained, in which the quantitative relationship is given by fitting the average variation characteristics [34].
(11)eh=292.1σn+1.78.According to formulas (10) and (11), the relationship between hydraulic conductivity and normal stress on the WCF of the fractured granite in the study area can be summarized as follows:
(12)Kr≈∑i=1ng12vL·aAiσni3,where a is constant. Therefore, the contribution index of the WCF to the hydraulic conductivity is
(13)Cwcf=∑j=1NAjσnj3.To facilitate the calculation in formula (13) and relate the physical meaning of Aj to A, Aj represents the contribution value of the orientation of all the WCF in the test section to the hydraulic conductivity of the rock mass, and the calculation method is the same as that used for the index A. When there is no WCF in the test section, that is, the number N of the WCF in the test section is 0, Cwcf=0. As the index A includes the contribution of some WCF orientation to the hydraulic conductivity of the rock mass, and because the relationship between the normal stress of the fracture and the hydraulic conductivity of the rock mass obtained by formula (12), it is appropriate to divide Aj into σnj3, which not only avoids the issue that the geometric parameter of the WCF repeatedly lead to excessive weights in the estimation of the hydraulic conductivity but also possesses some theoretical and logical basis. The introduction of correction parameters α, β, and γ does not make the three terms of formula (6) have a large difference, especially when Cwcf=0.The normal stressσnj of the WCF is calculated using the following steps. First, the fracture depth can be obtained from the location from the borehole wall image. Combined with the altitude of the borehole orifice, the fracture altitude can be calculated. Then, the principal stresses σH, σh, and σv at the fracture location can be obtained by substituting the fracture altitude into formula (1). Next, a unified spatial coordinate system of the principal stresses is constructed, as shown in Figure 10, wherein the H-axis positively points to the direction of the maximum principal stress, the h-axis positively points to the direction of the intermediate principal stress, and the v-axis points positively upward. According to second law of Cauchy stress, using the geometric relationship, the dip angle β and the adjusted dip direction α of the fracture are transformed into the rectangular coordinate system of the principal stress space to obtain the normal stress of the WCF:
(14)σn=nH2σH+nh2σh+nv2σv,where nH, nh, and nv are the unit normal vectors of the fracture in the principal stress space coordinate system.Figure 10
Schematic diagram of the orientation of discontinuity orientation in principle stress space.In summary, the hydraulic parametersAj of the WCF are calculated by formulas (7), (8), and (9), just like index A. Substituting σnj calculated by formula (14) and Aj into formula (13), Cwcf can be easily acquired. According to the CA model, the empirical constants α, β, and γ can be gained from linear regression. Thereby, the estimated hydraulic conductivity K of the CA model is finally obtained.


## 3.1. Model Construction and Parameter Selection

The permeability of fracture media in the DFN depends almost entirely on the spatial structure of the fracture system; meaning, the relationship between the fracture properties and permeability can be clarified by establishing a simple medium model of fractured rock mass. Snow [16] assumed that there was a water-conducting fracture network composed of n groups of directional fractures in a medium of pure fractured rock mass, and that the fracture water took a laminar flow. Therefore, the permeability tensor of each group of the WCF was superimposed based on the cubic law of the parallel plate model, and considering geometric factors, such as the fracture orientation, spacing, and hydraulic aperture, the following permeability tensor calculation formula was proposed, which is referred to as the snow model:
(5)KSNOW=g12v∑i=1nbi3Si1−nxi2−nxinyi−nxinzi−nyinxi1−nyi2−nyinzi−nzinxi−nzinyi1−nzi2,where KSNOW (m/s) is the hydraulic conductivity tensor of the rock mass, g (m/s2) is the gravitational acceleration, ν (m2/s) is the kinematic viscosity, bi (m) is the average hydraulic aperture of group i, Si (m) is the space of group i, nxi, nyi, and nzi are components in the x, y, and z directions, respectively, of the unit normal vector of group i fracture, and n is the number of fractures.Snow [16] considered the number of fractures to reflect the spacing and density of the fractures by superposing the hydraulic conductivity tensor of all the WCF, which is actually the contribution of rock integrity parameters to the hydraulic conductivity. Further, the snow model is suitable for laminar flow and clean fractures, which are consistent with the characteristics of the granite in the study area. Therefore, the snow model can be simplified into two parts: the orientation of all fractures and the hydraulic aperture. According to the superposition principle of the hydraulic conductivity tensor and considering the number of fractures, the fracture orientation index and hydraulic aperture index were tentatively selected to establish a new model for predicting the hydraulic conductivity of the low-permeability and unfilled fractured granite in underground water-sealed storage caverns. The hydraulic aperture is mainly controlled by normal stress, contact area, geometrical aperture, and roughness of the fracture, and it is difficult to determine the hydraulic aperture directly. Previous relevant studies [34, 41, 42] show that the quantitative relationship exists distinctly between normal stress and fracture hydraulic aperture. Fortunately, the normal stress of fractures is relatively easy to obtain in engineering; so, it was selected to replace the index of hydraulic aperture for the construction of the estimation model. Therefore, the normal stress, as a comprehensive index, combined with the orientation index of the discontinuity can comprehensively characterize the hydraulic behavior of the fracture.The hydraulic conductivity tensor of the rock mass is based on the permeability characteristics of the WCF. However, according to the results of the borehole water-pressure tests and borehole wall imaging tests in the study area, it was found that not all fractures are the WCF, and most have no seepage significance. As the orientation of the WCF determines the primary seepage direction of the rock mass, the fracture network model needs to be characterized emphatically. Therefore, in the case of existing fracture orientation index, the proposed model introduces orientation index of the WCF to characterize the contribution of the hydraulic characteristics of the WCF to the permeability of the rock mass.


## 3.2. Novel Model for Estimating Hydraulic Conductivity of Fractured Granite in Water-Sealed Storage Caverns

Based on the snow model and fracture network model elements, this study used the fracture orientation (denoted asA or Aj) and the normal stress (σn) of the WCF as the basic index to characterize the geometric characteristics of the fracture and the hydraulic properties of the WCF, respectively, and proposed an estimation model for hydraulic conductivity of the low-permeability and unfilled fractured granite in underground water-sealed storage caverns, called the CA model:
(6)K=α∙A+β∙Cwcf+γ,where α, β, and γ are empirical constants, and K (cm/s) is the hydraulic conductivity of each test section. The index A represents the contribution of the orientation of all fractures to the hydraulic conductivity, and Cwcf is an integrated hydraulic index representing the contribution of all the WCF to hydraulic conductivity of rock mass. Therefore, the name of the CA model is taken from the main indices Cwcf and A. In terms of the unit of two parameters, the index A proposed in this study is dimensionless, and the unit of the index Cwcf is MPa-3. Although the unit is different, the model only uses its numerical value.
(1)
The method for obtaining the indexA. It is assumed that the contribution of the orientation of n fractures to the hydraulic conductivity of the rock mass at a certain depth range is An. Based on the principle of permeability tensor superposition, the following formula can be obtained:(7)An=∑i=1nAi=∑i=1n1−nxi2−nxinyi−nxinzi−nyinxi1−nyi2−nyinzi−nzinxi−nzinyi1−nzi2.If the dip direction of fracture i is αi, the dip angle is βi, the x-axis points to the east, the y-axis points to the north, the z-axis points upwards, and a space rectangular coordinate system ENZ(xyz) is constructed. As shown in Figure 8, the x, y, and z components of the normal vector of the ID i of the fracture are expressed as follows:
(8)nxi=sinαi∙sinβinyi=cosαi∙sinβinzi=cosβi.Figure 8
Schematic diagram of unit normal vector of fracture orientation in geographic coordinate system.The eigenvalues ofAn was calculated and denoted as A1, A2, and A3, of which the three eigenvectors are the principal directions of the eigenvalues. The fracture orientation index A is defined as the geometric mean value of the eigenvalue An, as follows:
(9)A=A1∙A2∙A33,where the physical meaning of A is the contribution of the orientation of all fractures to the hydraulic conductivity in the depth range of the test section. Therefore, the orientation index of fractures can be obtained by simultaneous formulas (7), (8), and (9).
(2)
The method for obtaining the indexCwcf. To obtain Cwcf, the identification and selection of WCF should be conducted first to determine the mathematical relationship between the hydraulic parameters of the WCF and the hydraulic conductivity of the rock mass. Firstly, based on borehole wall image tests and water-pressure tests data, the hydraulic conductivity of each test section and internal fractures orientation was obtained. Secondly, the orientation ranges of fractures with a small hydraulic conductivity in the test section need to be screened out, and it was considered that it does not belong to the primary seepage zone. Finally, it was removed from the orientation range of the test sections with large hydraulic conductivity, and what were retained were the WCF and its orientation in each test section. Based on the DFN model of a real discrete fracture network generated using a Monte Carlo simulation by Min et al. [6], Figure 9 represents the WCF selected from fracture network and internal groundwater seepageFigure 9
The WCF and the direction of groundwater flow in fracture system (the blue line represents the WCF, and the black line represents ordinary discontinuities).Although we have selected the normal stress index (σn) instead of the hydraulic aperture (eh) of the WCF to build the new model, it still plays a role as a bond between the hydraulic conductivity of rock mass and normal stress on the fracture. Based on the research results of Cao et al. [34] on the fracture hydraulic aperture of slightly weathered granite which located nearby the study area, it was assumed that a rock mass test section with upper and lower surfaces parallel to its contained fracture and with a distance of L between the two interfaces. Other parameters were kept the same as previously mentioned. The hydraulic conductivity Kr of the fractured rock mass is the superposition of the hydraulic conductivity of each fracture in the test section as follows:
(10)Kr=∑i=1ngeh312vL·Ai.Thus, the basic relationship between hydraulic aperture and normal stress was obtained, in which the quantitative relationship is given by fitting the average variation characteristics [34].
(11)eh=292.1σn+1.78.According to formulas (10) and (11), the relationship between hydraulic conductivity and normal stress on the WCF of the fractured granite in the study area can be summarized as follows:
(12)Kr≈∑i=1ng12vL·aAiσni3,where a is constant. Therefore, the contribution index of the WCF to the hydraulic conductivity is
(13)Cwcf=∑j=1NAjσnj3.To facilitate the calculation in formula (13) and relate the physical meaning of Aj to A, Aj represents the contribution value of the orientation of all the WCF in the test section to the hydraulic conductivity of the rock mass, and the calculation method is the same as that used for the index A. When there is no WCF in the test section, that is, the number N of the WCF in the test section is 0, Cwcf=0. As the index A includes the contribution of some WCF orientation to the hydraulic conductivity of the rock mass, and because the relationship between the normal stress of the fracture and the hydraulic conductivity of the rock mass obtained by formula (12), it is appropriate to divide Aj into σnj3, which not only avoids the issue that the geometric parameter of the WCF repeatedly lead to excessive weights in the estimation of the hydraulic conductivity but also possesses some theoretical and logical basis. The introduction of correction parameters α, β, and γ does not make the three terms of formula (6) have a large difference, especially when Cwcf=0.The normal stressσnj of the WCF is calculated using the following steps. First, the fracture depth can be obtained from the location from the borehole wall image. Combined with the altitude of the borehole orifice, the fracture altitude can be calculated. Then, the principal stresses σH, σh, and σv at the fracture location can be obtained by substituting the fracture altitude into formula (1). Next, a unified spatial coordinate system of the principal stresses is constructed, as shown in Figure 10, wherein the H-axis positively points to the direction of the maximum principal stress, the h-axis positively points to the direction of the intermediate principal stress, and the v-axis points positively upward. According to second law of Cauchy stress, using the geometric relationship, the dip angle β and the adjusted dip direction α of the fracture are transformed into the rectangular coordinate system of the principal stress space to obtain the normal stress of the WCF:
(14)σn=nH2σH+nh2σh+nv2σv,where nH, nh, and nv are the unit normal vectors of the fracture in the principal stress space coordinate system.Figure 10
Schematic diagram of the orientation of discontinuity orientation in principle stress space.In summary, the hydraulic parametersAj of the WCF are calculated by formulas (7), (8), and (9), just like index A. Substituting σnj calculated by formula (14) and Aj into formula (13), Cwcf can be easily acquired. According to the CA model, the empirical constants α, β, and γ can be gained from linear regression. Thereby, the estimated hydraulic conductivity K of the CA model is finally obtained.


## 4. Engineering Application and Reliability Analysis of CA Model

### 4.1. Parameter Fitting of the CA Model in the Study Area

Based on the interpretation results of borehole wall imaging and water-pressure test data of low-permeability and unfilled fractured granite in situ, we calculated the values ofA, Cwcf, and Kin‐situ (the in situ test values of hydraulic conductivity from double-packer borehole water-pressure test) of each test section in ZK61 and ZK69, as listed in Table 1. The fracture strike range with low hydraulic conductivity (K<5×10−5m/d) was removed from the strike range with high hydraulic conductivity (K>5×10−3m/d), and that of the WCF was obtained as 0°–20°, 30°–50°, 330°–340°, and 350°–360°. In addition, according to the results of the borehole wall image, it was found that the dip direction of the dominant fractures in the rock mass is mainly distributed in the intervals 80°–90°, 110°–140°, and 270°–300°.Table 1
Test results of boreholes in the study area.
BoreholeDepth [m]IndexAIndexCwcfKin−situ [m/d]ZK6131-416.89907.91×10−541-5112.02409.75×10−351-617.58403.60×10−461-718.0872.53×10−34.29×10−471-816.5631.34×10−25.67×10−481-913.7735.85×10−34.01×10−591-1016.5435.44×10−33.91×10−4101-1112.3867.30×10−31.74×10−4111-1212.9727.41×10-34.25×10−4121-1311.79101.46×10−4131-1415.6853.55×10−26.85×10−4151-1613.3947.01×10−34.12×10−4161-1713.1871.16×10−31.55×10−4171-1816.4294.24×10−31.90×10−4181-1913.6964.12×10−33.99×10−4191-2012.8513.48×10−32.97×10−4201-2118.0355.02×10−32.09×10−4ZK6933-435.5922.90×10−21.23×10−442-526.7611.30×10−21.22×10−451-615.2541.86×10−29.54×10−560-705.3674.30×10−32.36×10−469-792.8268.72×10−32.01×10−478-886.1471.80×10−34.48×10−487-978.9333.27×10−31.94×10−496-1065.0946.54×10−36.34×10−4105-1157.7383.07×10−38.57×10−4114-12410.4731.96×10−23.63×10−4123-1335.8953.42×10−34.13×10−4132-1423.95804.37×10−4141-1510.63301.11×10−4150-1604.93101.04×10−4159-1691.1851.24×10−21.48×10−4168-1782.7138.79×10−41.30×10−4177-1874.5744.13×10−32.19×10−4186-1968.0451.34×10−12.65×10−3195-2058.2791.33×10−12.08×10−3204-2141.4049.13×10−21.15×10−3213-2234.0966.33×10−38.40×10−4The CA model was used to fit the test data of ZK61 and ZK69, and the model parameters are listed in Table2. The spatial variation characteristics of the regional geological conditions and permeability may be the reason for the different fitting parameters of the CA model in two boreholes. Given that ZK61 and ZK69 are located in the same engineering site, they have similar geological and hydrogeological conditions, and the hydraulic conductivity is mainly concentrated in the magnitude of 10-4–10-3 m/d. Moreover, the CA model was used to fit the two boreholes tests data simultaneously, and the estimation results are in accordance with the test results in situ (R2=0.85). Therefore, the model proposed in this study is suitable for projects that have similar engineering geological conditions to that of the water-sealed storage caverns. According to the aforementioned research, the CA model may be more appropriate to estimate the hydraulic conductivity of rock mass with low permeability and internal unfilled fractures.Table 2
Fitting results of CA model parameters.
BoreholesαβγR2ZK612.89×10−61.30×10−21.95×10−40.78ZK691.84×10−51.40×10−29.88×10−50.74All2.04×10−51.38×10−27.40×10−50.85


### 4.2. Verification and Reliability Analysis of the CA Model

Since this study has proposed a new model, it is necessary to test it, which is an inevitable process of scientific research. Moreover, the superiority of the new model needs to be demonstrated from the comparison with other models. Therefore, existing indices or models should be summarized and used to predict the hydraulic conductivity of low-permeability and unfilled fractured granite. The prediction results of CA model and previous models are compared with theKin‐situ to explicate its superiority.
(1)
The HC model and the RMP model. Existing comprehensive indices or estimation models mainly consider rock integrity (characterized by the RID or RQD), depth index (DI orZ), lithology permeability index (LPI), gouge content designation (GCD), and fracture aperture designation (AD) as basic parameters. To be specific, the HC index and the RMP index are the most representative, as follows:(15)HC=1−RQD1−GCDDILPI,(16)RMP=1−RQD1−RIDADLPI,where LPI=0.15 is for granite through the classification of rock permeability [28]. Among these, the GCD is calculated using the following equation:
(17)GCD=RGRT−Rs,where RG is the total length of filling content and GCD=0 for the unfilled granite in the study area. In addition, the HC index defined DI as the buried depth index:
(18)DI=1−LCLT,where LC is a depth that is located at the middle of a double-packer test interval in the borehole, and LT is total borehole length. However, as the hydraulic conductivity of the rock mass is independent of the length of the borehole, the selected index DI may be unreasonable. The fracture aperture parameter (AD) in the RMP index is defined as the ratio of the sum of fracture aperture dA and the length L of the test section. Since the two models are based on the HC index and RMP index, respectively, they can be called the HC model and RMP model. Based on two comprehensive indices, the two estimation models were used in estimating the hydraulic conductivity of sedimentary rocks with gouge filled and water-sealed storage caverns with unfilled granite, respectively.
(19)K=bHCc,(20)K=mRMPn,where b, c, m, and n are all fitting parameters. Using the HC model, the hydraulic conductivity of ZK76 at different depths in the study area was estimated, and the results are listed in Table 3.Table 3
HC model estimation results of ZK76.
Z [m]1-RQDDI1-GCDLPIKHCmodel [m/d]350.360.8401.00.1503.55×10−5450.330.7951.00.1502.91×10−5550.30.7491.00.1502.35×10−5650.180.7031.00.1501.07×10−5750.230.6581.00.1501.36×10−5850.310.6121.00.1501.86×10−5950.190.5661.00.1508.51×10−61050.180.5211.00.1507.04×10−61250.20.4291.00.1506.24×10−61350.140.3881.00.1503.32×10−61550.130.2921.00.1502.02×10−61650.110.2471.00.1501.27×10−61750.150.2011.00.1501.47×10−61850.240.1551.00.1501.97×10−61950.140.1101.00.1505.79×10−72050.230.0641.00.1505.46×10−7Based on the empirical RMP model, the hydraulic conductivity of the ZK76 at different depths was estimated, and the estimation results are listed in Table4.
(2)
The ZRF model. In addition to considering RQD, the ZRF model introducedZ and GCD to predict the hydraulic conductivity of fractured rock mass at the site of Yagen-II Hydropower Station in China. The expression of the ZRF model is(21)lgK=κ+λ·PD+ω·lgZ,(22)PD=1−RQD1−GCD,where κ, λ, and ω are empirical constants. The hydraulic conductivity of the ZK76 at different depths was estimated by the ZRF model, and the results are listed in Table 5.
(3)
Analysis of the CA model and applicability comparison with other models. Based on the water-pressure test of ZK76, the hydraulic conductivity of rock mass in each depth section was calculated according to formula (3). According to the interpretation results of wall imaging of ZK76, the location and orientation of all fractures in each depth section were obtained. Based on the data of orientation and hydraulic conductivity, the location and orientation of WCF were then picked out. Equations (7)–(9) were used to calculate A and Aj; the normal stress σnj of the WCF was calculated by combining formulas (1) and (14). According to Table 2, the fitting parameters of the CA model are as follows: α=2.04×10−5, β=1.38×10−2, and γ=7.40×10−5. Substituting the above parameters and basic indices into the CA model(23)K=α∙A+β∙∑j=1NAjσnj3+γ.Table 4
RMP model estimation results of ZK76.
Z [m]1-RQD1-RIDADLPIKRMPmodel [m/d]350.360.256.72×10−40.154.63×10−4450.330.181.04×10−30.154.72×10−4550.300.141.23×10−30.154.04×10−4650.180.119.75×10−40.151.72×10−4750.230.117.91×10−40.151.78×10−4850.310.111.12×10−30.153.11×10−4950.190.081.31×10−40.152.41×10−51050.180.088.97×10−40.151.22×10−41250.200.085.92×10−40.159.29×10−51350.140.083.35×10−40.154.17×10−51550.130.086.73×10−40.157.15×10−51650.110.082.48×10−40.152.60×10−51750.150.083.07×10−40.154.10×10−51850.240.086.92×10−40.151.25×10−41950.140.108.38×10−40.151.12×10−42050.230.101.25×10−30.152.43×10−4Table 5
ZRF model estimation results of ZK76.
Z [m]1-RQD1-GCDPDKZRFmodel [m/d]350.361.00.363.01×10−4450.331.00.332.06×10−4550.31.00.301.47×10−4650.181.00.186.53×10−5750.231.00.237.63×10−5850.311.00.311.07×10−4950.191.00.195.00×10−51050.181.00.184.34×10−51250.21.00.204.18×10−51350.141.00.142.81×10−51550.131.00.132.36×10−51650.111.00.112.00×10−51750.151.00.152.38×10−51850.241.00.243.75×10−51950.141.00.142.05×10−52050.231.00.233.25×10−5The prediction result of hydraulic conductivity of ZK76 has been listed in Figure11, which also involved the estimation results of the HC model, RMP model, ZRF model, and in situ measured values of ZK76. It can be seen that the new CA model provides significantly better results than the existing models. Specifically, the root mean square error of the CA model is 3.76×10−5m/d, which is far less than that of the HC, RMP, and ZRF models (respectively, 1.46×10−4m/d, 1.02×10−4m/d, and 1.24×10−4m/d), and the relative errors are 45.0%, 95.8%, 64.1%, and 70.6%, respectively. Thus, the CA model is more suitable and accurate in estimating the hydraulic conductivity of low-permeability and unfilled fractured granite in underground water-sealed storage caverns.Figure 11
Comparison between the estimated values of hydraulic conductivity.


## 4.1. Parameter Fitting of the CA Model in the Study Area

Based on the interpretation results of borehole wall imaging and water-pressure test data of low-permeability and unfilled fractured granite in situ, we calculated the values ofA, Cwcf, and Kin‐situ (the in situ test values of hydraulic conductivity from double-packer borehole water-pressure test) of each test section in ZK61 and ZK69, as listed in Table 1. The fracture strike range with low hydraulic conductivity (K<5×10−5m/d) was removed from the strike range with high hydraulic conductivity (K>5×10−3m/d), and that of the WCF was obtained as 0°–20°, 30°–50°, 330°–340°, and 350°–360°. In addition, according to the results of the borehole wall image, it was found that the dip direction of the dominant fractures in the rock mass is mainly distributed in the intervals 80°–90°, 110°–140°, and 270°–300°.Table 1
Test results of boreholes in the study area.
BoreholeDepth [m]IndexAIndexCwcfKin−situ [m/d]ZK6131-416.89907.91×10−541-5112.02409.75×10−351-617.58403.60×10−461-718.0872.53×10−34.29×10−471-816.5631.34×10−25.67×10−481-913.7735.85×10−34.01×10−591-1016.5435.44×10−33.91×10−4101-1112.3867.30×10−31.74×10−4111-1212.9727.41×10-34.25×10−4121-1311.79101.46×10−4131-1415.6853.55×10−26.85×10−4151-1613.3947.01×10−34.12×10−4161-1713.1871.16×10−31.55×10−4171-1816.4294.24×10−31.90×10−4181-1913.6964.12×10−33.99×10−4191-2012.8513.48×10−32.97×10−4201-2118.0355.02×10−32.09×10−4ZK6933-435.5922.90×10−21.23×10−442-526.7611.30×10−21.22×10−451-615.2541.86×10−29.54×10−560-705.3674.30×10−32.36×10−469-792.8268.72×10−32.01×10−478-886.1471.80×10−34.48×10−487-978.9333.27×10−31.94×10−496-1065.0946.54×10−36.34×10−4105-1157.7383.07×10−38.57×10−4114-12410.4731.96×10−23.63×10−4123-1335.8953.42×10−34.13×10−4132-1423.95804.37×10−4141-1510.63301.11×10−4150-1604.93101.04×10−4159-1691.1851.24×10−21.48×10−4168-1782.7138.79×10−41.30×10−4177-1874.5744.13×10−32.19×10−4186-1968.0451.34×10−12.65×10−3195-2058.2791.33×10−12.08×10−3204-2141.4049.13×10−21.15×10−3213-2234.0966.33×10−38.40×10−4The CA model was used to fit the test data of ZK61 and ZK69, and the model parameters are listed in Table2. The spatial variation characteristics of the regional geological conditions and permeability may be the reason for the different fitting parameters of the CA model in two boreholes. Given that ZK61 and ZK69 are located in the same engineering site, they have similar geological and hydrogeological conditions, and the hydraulic conductivity is mainly concentrated in the magnitude of 10-4–10-3 m/d. Moreover, the CA model was used to fit the two boreholes tests data simultaneously, and the estimation results are in accordance with the test results in situ (R2=0.85). Therefore, the model proposed in this study is suitable for projects that have similar engineering geological conditions to that of the water-sealed storage caverns. According to the aforementioned research, the CA model may be more appropriate to estimate the hydraulic conductivity of rock mass with low permeability and internal unfilled fractures.Table 2
Fitting results of CA model parameters.
BoreholesαβγR2ZK612.89×10−61.30×10−21.95×10−40.78ZK691.84×10−51.40×10−29.88×10−50.74All2.04×10−51.38×10−27.40×10−50.85


## 4.2. Verification and Reliability Analysis of the CA Model

Since this study has proposed a new model, it is necessary to test it, which is an inevitable process of scientific research. Moreover, the superiority of the new model needs to be demonstrated from the comparison with other models. Therefore, existing indices or models should be summarized and used to predict the hydraulic conductivity of low-permeability and unfilled fractured granite. The prediction results of CA model and previous models are compared with theKin‐situ to explicate its superiority.
(1)
The HC model and the RMP model. Existing comprehensive indices or estimation models mainly consider rock integrity (characterized by the RID or RQD), depth index (DI orZ), lithology permeability index (LPI), gouge content designation (GCD), and fracture aperture designation (AD) as basic parameters. To be specific, the HC index and the RMP index are the most representative, as follows:(15)HC=1−RQD1−GCDDILPI,(16)RMP=1−RQD1−RIDADLPI,where LPI=0.15 is for granite through the classification of rock permeability [28]. Among these, the GCD is calculated using the following equation:
(17)GCD=RGRT−Rs,where RG is the total length of filling content and GCD=0 for the unfilled granite in the study area. In addition, the HC index defined DI as the buried depth index:
(18)DI=1−LCLT,where LC is a depth that is located at the middle of a double-packer test interval in the borehole, and LT is total borehole length. However, as the hydraulic conductivity of the rock mass is independent of the length of the borehole, the selected index DI may be unreasonable. The fracture aperture parameter (AD) in the RMP index is defined as the ratio of the sum of fracture aperture dA and the length L of the test section. Since the two models are based on the HC index and RMP index, respectively, they can be called the HC model and RMP model. Based on two comprehensive indices, the two estimation models were used in estimating the hydraulic conductivity of sedimentary rocks with gouge filled and water-sealed storage caverns with unfilled granite, respectively.
(19)K=bHCc,(20)K=mRMPn,where b, c, m, and n are all fitting parameters. Using the HC model, the hydraulic conductivity of ZK76 at different depths in the study area was estimated, and the results are listed in Table 3.Table 3
HC model estimation results of ZK76.
Z [m]1-RQDDI1-GCDLPIKHCmodel [m/d]350.360.8401.00.1503.55×10−5450.330.7951.00.1502.91×10−5550.30.7491.00.1502.35×10−5650.180.7031.00.1501.07×10−5750.230.6581.00.1501.36×10−5850.310.6121.00.1501.86×10−5950.190.5661.00.1508.51×10−61050.180.5211.00.1507.04×10−61250.20.4291.00.1506.24×10−61350.140.3881.00.1503.32×10−61550.130.2921.00.1502.02×10−61650.110.2471.00.1501.27×10−61750.150.2011.00.1501.47×10−61850.240.1551.00.1501.97×10−61950.140.1101.00.1505.79×10−72050.230.0641.00.1505.46×10−7Based on the empirical RMP model, the hydraulic conductivity of the ZK76 at different depths was estimated, and the estimation results are listed in Table4.
(2)
The ZRF model. In addition to considering RQD, the ZRF model introducedZ and GCD to predict the hydraulic conductivity of fractured rock mass at the site of Yagen-II Hydropower Station in China. The expression of the ZRF model is(21)lgK=κ+λ·PD+ω·lgZ,(22)PD=1−RQD1−GCD,where κ, λ, and ω are empirical constants. The hydraulic conductivity of the ZK76 at different depths was estimated by the ZRF model, and the results are listed in Table 5.
(3)
Analysis of the CA model and applicability comparison with other models. Based on the water-pressure test of ZK76, the hydraulic conductivity of rock mass in each depth section was calculated according to formula (3). According to the interpretation results of wall imaging of ZK76, the location and orientation of all fractures in each depth section were obtained. Based on the data of orientation and hydraulic conductivity, the location and orientation of WCF were then picked out. Equations (7)–(9) were used to calculate A and Aj; the normal stress σnj of the WCF was calculated by combining formulas (1) and (14). According to Table 2, the fitting parameters of the CA model are as follows: α=2.04×10−5, β=1.38×10−2, and γ=7.40×10−5. Substituting the above parameters and basic indices into the CA model(23)K=α∙A+β∙∑j=1NAjσnj3+γ.Table 4
RMP model estimation results of ZK76.
Z [m]1-RQD1-RIDADLPIKRMPmodel [m/d]350.360.256.72×10−40.154.63×10−4450.330.181.04×10−30.154.72×10−4550.300.141.23×10−30.154.04×10−4650.180.119.75×10−40.151.72×10−4750.230.117.91×10−40.151.78×10−4850.310.111.12×10−30.153.11×10−4950.190.081.31×10−40.152.41×10−51050.180.088.97×10−40.151.22×10−41250.200.085.92×10−40.159.29×10−51350.140.083.35×10−40.154.17×10−51550.130.086.73×10−40.157.15×10−51650.110.082.48×10−40.152.60×10−51750.150.083.07×10−40.154.10×10−51850.240.086.92×10−40.151.25×10−41950.140.108.38×10−40.151.12×10−42050.230.101.25×10−30.152.43×10−4Table 5
ZRF model estimation results of ZK76.
Z [m]1-RQD1-GCDPDKZRFmodel [m/d]350.361.00.363.01×10−4450.331.00.332.06×10−4550.31.00.301.47×10−4650.181.00.186.53×10−5750.231.00.237.63×10−5850.311.00.311.07×10−4950.191.00.195.00×10−51050.181.00.184.34×10−51250.21.00.204.18×10−51350.141.00.142.81×10−51550.131.00.132.36×10−51650.111.00.112.00×10−51750.151.00.152.38×10−51850.241.00.243.75×10−51950.141.00.142.05×10−52050.231.00.233.25×10−5The prediction result of hydraulic conductivity of ZK76 has been listed in Figure11, which also involved the estimation results of the HC model, RMP model, ZRF model, and in situ measured values of ZK76. It can be seen that the new CA model provides significantly better results than the existing models. Specifically, the root mean square error of the CA model is 3.76×10−5m/d, which is far less than that of the HC, RMP, and ZRF models (respectively, 1.46×10−4m/d, 1.02×10−4m/d, and 1.24×10−4m/d), and the relative errors are 45.0%, 95.8%, 64.1%, and 70.6%, respectively. Thus, the CA model is more suitable and accurate in estimating the hydraulic conductivity of low-permeability and unfilled fractured granite in underground water-sealed storage caverns.Figure 11
Comparison between the estimated values of hydraulic conductivity.


## 5. Discussion

Existing comprehensive indices or models are not suitable for the low-permeability and unfilled fractured granite in underground water-sealed storage caverns as the aforementioned LPI and GCD lead to a weak representation and poor correlations between the geological parameters and the hydraulic conductivity, resulting in large estimation errors. Although the RMP model has been proposed for fractured granite of LPG caverns in adjacent areas, it commonly overestimated the contribution of rock integrity to the hydraulic conductivity by simultaneously adding the RQD and RID parameters. Therefore, the randomness of discontinuities causes large fluctuations in the hydraulic conductivity of the rock mass at different depths by affecting the integrity of rock mass. In addition, the fracture aperture designation is difficult to measure, which may be another cause of errors.Compared with existing models, the CA model considers the characteristics of the seepage fracture network model of the rock mass. Based on the snow model, the number of fractures is introduced to reflect the integrity of the rock mass, and the orientation of fractures is considered to reflect the permeability of the rock mass. Moreover, the proposed model considers the characterization of mechanical and hydraulic properties of the rock mass and has a certain theoretical basis. Simply put, the CA model has only two basic parameters: fracture orientation index (A) and normal stress index (σn), which can be easily obtained using test data in situ. Thus, it is reliable to estimate the hydraulic conductivity of the cavern granite. However, this model is dependent on the number and accuracy of the in situ data, and it may be difficult to identify the WCF. Therefore, an accurate screening method for the WCF must be explored.A sensitivity analysis of the fracture orientation index and normal stress index of the model indicated that the change in the hydraulic conductivity caused byσn is much larger than that caused by A when the two parameters change by the same multiple. This shows that the influence of the normal stress of the fracture on the rock mass permeability is greater than that on the fracture orientation. This is because the normal stress controls the magnitude of the permeability of the rock mass, whereas the fracture orientation determines the seepage direction. Further, it shows that the normal stress index selected in this study is reasonable to construct the new model.


## 6. Conclusions

(1)
Existing models have difficulty accurately reflecting the hydraulic conductivity characteristics of low-permeability and unfilled fractured granite in underground water-sealed storage caverns. The rock mass was regarded as a fractured network comprising WCF; using the test data of borehole wall imaging and water-pressure tests of fractured granite of an underground water-sealed storage caverns, the fracture orientation index and normal stress index were selected as basic parameters, and the CA model for estimating hydraulic conductivity was established by introducing the WCF contribution indexCwcf(2)
The hydraulic conductivity of ZK76 in the study area was estimated using the CA model to determine its reliability. Compared with the existing models, the estimation results of the CA model were closest to the in situ measured values which have a high accuracy. Further, it is extremely suitable for estimating the hydraulic conductivity of low-permeability and unfilled fractured granite in underground water-sealed storage caverns which can also be applied for similar rock mass(3)
The paper proposed a screening method of the water conduction fractures (WCF) that mainly consists of the following steps. Firstly, the appropriate method should be used to determine the hydraulic conductivity of each test sections. Secondly, the strike and orientation of internal fractures were obtained based on survey results of borehole wall image tests and water-pressure tests. Thirdly, the fracture strike range in the test section with low hydraulic conductivity was removed from the strike range with high hydraulic conductivity, and the WCF in rock mass was then identified. In addition, the study indicated that the WCF orientation can be determined only by the dominant fractures orientation in the study area of the underground water-sealed storage caverns(4)
A sensitivity analysis of the parameters of the CA model revealed that the influence of the normal stress of the fracture on the permeability of rock mass is greater than the fracture orientation



---

*Source: 6423048-2021-09-27.xml*