Seepage of Groundwater in an Underground Fractured Rock Mass and Its Sustainable Engineering Application

Abstract

Due to the existence of tiny cracks in rock, underground engineering has begun to consider how to divert a large amount of groundwater. To divert groundwater more effectively, it is necessary to master the seepage characteristics of fluids in the micropores of rocks. Based on rock samples obtained from an underground engineering site, this paper analyzes the microscopic pore structure of the rock through a combination of laboratory tests and numerical simulations and inputs this information into a computer model. The fluid seepage state in the rock under different conditions is simulated in the computer model, and parameters such as the fluid seepage velocity in the rock are obtained. Afterwards, it has been verified by engineering practice that the smallest remaining water inflow can reach 0.06‰. The results of this paper can effectively guide the discharge of groundwater to better manage water resources, greatly reduce the pollution of groundwater in construction and production environments, and reduce the pollution caused by grouting projects. Furthermore, the cleanliness and safety of underground engineering construction and production could be ensured.

1. Introduction

In many countries, groundwater disasters affect all aspects of human life. These include floods in residential areas [1], inundation of underground construction projects [2], and other disaster impacts. Groundwater disasters also seriously affect economic development [3]. The physical and chemical hazards of groundwater seriously affect human life. How to solve this problem and make rational use of groundwater is an urgent problem to be solved [4,5,6].

Regarding the multiple hazards of a large amount of groundwater under complex geological conditions and the problem of sustainable development, many scholars have carried out research from different aspects [7,8,9,10].

To solve this problem, the method of injecting cement slurry and grouting to plug the water is generally used to plug rock pores [11,12,13]. This is due to the fact that, in underground engineering, the main circulation path of groundwater is the pores in the rock.

Regarding grouting materials, Gao found that although cement slurry can seal groundwater to a certain extent, it also causes problems such as high energy consumption and large carbon dioxide emissions. Therefore, he developed cement slurries with different admixtures and different proportions, and the research results indicate high performance, low cost, and environmental protection [14]. Gao considered the interface transition zone between the aggregate and cementitious material, and the results showed that materials such as carbon nanotubes can be used as excellent additives to reduce the propagation of cracks [15].

The abovementioned scholars have conducted in-depth research on cement grouting and water shutoff, and some have also considered certain environmental protection measures and have made positive research progress.

In the grouting process, Han et al. [16], studied the effect of viscosity on cement grouting and the results show that the reduction of grouting pressure can improve the grouting efficiency. Han’s work focused on the site conditions of the Pearl River Project and, considering the influence of time-varying viscosity and other factors, the results show that the proposed new grouting algorithm effectively helps grouting to plug water [17].

The above studies have two drawbacks: Firstly, this method can only block part of the groundwater and reduce groundwater damage, but it does not substantially solve the problem of groundwater outflow. This is because there are still tiny cracks in the rock that ordinary cement slurry cannot enter, and groundwater still flows out of the tiny cracks [18,19]. Secondly, the mixing of groundwater and cement slurry during grouting affects the properties of the cement slurry, and the final results are often worse than the research results [20]. This is because a large amount of groundwater will dilute the cement slurry, thereby reducing the setting time of the cement slurry. The cement slurry does not condense in the pores of the rock for a long time and cannot achieve the effect of blocking the groundwater, which affects the entire grouting and water-blocking project, making the groundwater hazard situation still severe [21].

Therefore, an increasing number of scholars have begun to study groundwater drainage (draining groundwater to other places or collecting it for resource utilization). Hong [22], derived the coupling equation of stress and seepage, and applied it to high-speed railway tunnels to reduce the impact of groundwater on the project. Zhou et al. [23], numerically simulated the influence mechanism of tunnel water pressure on lining, and the proposed new method effectively guided the tunnel construction. Zheng et al. [24], carried out experiments on the project site, and proposed improvement measures to effectively solve the problem of air leakage in water diversion, which greatly improved the construction efficiency.

In a grouting water plugging project, Ren [18,25], installed part of the drainage pipe utilizing the field experiment method to draw the groundwater out of the grouting project, which effectively reduced a large amount of groundwater during the grouting process. The way in which to address the amount of grouting engineering groundwater presents new feasible ideas. Based on engineering examples, Jin et al. [26], considered the problem of slurry solidification and verified the drainage design with experiments which effectively reduced the impact of groundwater on grouting. To determine the influence of groundwater on grouting, Pan et al. [27], considered more factors of grouting pipes and the characteristics of underground water inflow, and proved that drainage is also an important means to ensure project safety. Pan et al. [28], considered diverting grouting and optimizing the construction process. The results show that drainage is the best method for sustainable development when the water inflow is large.

Scholars have found that these methods have two main advantages. For instance, large amounts of groundwater can be diverted, which can reduce pollution and hazards to engineering construction and avoid the problem of groundwater flowing out of tiny cracks in the grouting method of water blocking, which is conducive to the engineering environment and safety assurance. In addition, the drained groundwater can be stored for domestic purposes, such as hydroelectric power generation. This not only saves water resources and increases resource utilization but also enables better water resource management [29]. However, drainage is not easy, and drainage tubes and other equipment should not be arranged blindly. If the main part of the drainage pipe and other equipment are placed in a location with small groundwater flow, the location with a large flow of groundwater cannot effectively divert the groundwater; therefore, a large amount of groundwater will still gather there and the result will still cause groundwater pollution and harm. Therefore, it is necessary to rationally arrange drainage pipes and other equipment to efficiently divert groundwater [30]. The key to rationally arranging drainage pipes and other equipment is to clarify the pore structure characteristics of the rock and then grasp the flow characteristics of fluid seepage. In addition, most drainage engineering examples studied by scholars did not exceed a depth of 1000 m. The deep underground geological conditions below 1000 m (hereafter referred to as the deep) are complex, and the geological conditions of the shallow part are quite different. In other words, the fluid seepage characteristic at low depths cannot be directly used for reference in deep underground engineering. This further increases the difficulty of drainage operations in deep underground projects.

Therefore, the characteristics of how the fluid seeps into the deep underground rock, where the flow is large, where the flow velocity is fast, and how these locations are distributed are all important to the drainage of groundwater. This manuscript aims to effectively divert groundwater, reduce groundwater hazards, reduce environmental pollution caused by groundwater, solve the pollution and safety hazards of groundwater influx to underground resource mining and other projects, and help clean the production of underground resource mining and other projects. At the same time, the use of water resources should be increased to contribute to the sustainable development of resources. We mainly explore the fluid seepage characteristic in deep underground rock and perform the following work. In this study, rock samples were obtained from the site of a deep underground project. And this research does not rely solely on engineering experience, but comprehensively uses various methods such as experiments, numerical simulations, and field verifications.

A fluid seepage simulation under different conditions was carried out on the computer, and the fluid flow characteristic (this does not refer to specific research laws, etc., but refers to the description of the fluid flow state that can be applied to engineering and is sufficient to guide the construction content such as the layout of engineering drainage pipes.) in the deep underground engineering rock was explored. Since this study not only examined the macropores but also included the micropores in the rock, the results are more accurate and credible. Through a test in the actual project, the reliability of the conclusion was verified. The results obtained can effectively guide the drainage of groundwater to carry out better water resource management and can greatly reduce the pollution and harm from groundwater to the environment during project construction and production, as well as solving the problems of pollution and safety hazards caused by overflowing groundwater to underground resource mining and other projects, and helping the clean production of underground resource mining and other projects. At the same time, this study contributes to the sustainable development of water resources.

The main idea of this research is shown in Figure 1.

2. Methodology

2.1. Brief Introduction of the Pore Structure Experimental Method

The microscopic pore structure characteristics of the rock were mainly studied by scanning X-ray diffraction, casting thin sections and high-pressure mercury intrusion experiments.

The sandstone sample was obtained from the construction site of an underground project (depth exceeding 1300 m). According to different experimental requirements, the rock samples were cut into shapes that met the needs of the experiment. The depth of the engineering site for this study was constant and samples were taken at 5 different locations during the study. Since the results of the five studies were largely consistent, not all of them are listed.

X-ray diffraction: through X-ray diffraction experiments, the content and characteristics of various minerals can be determined.

Casting thin sections: through the casting thin section experiment, the pore throat type, particle size and contact method can be identified under a polarized light microscope.

High-pressure mercury intrusion: High-pressure mercury intrusion is one of the most important methods for studying the characteristics of rock pore throat structures. It can obtain the pore throat characteristics and porosity of rock samples.

2.2. Rock Lithological Analysis

According to X-ray diffraction and casting thin section experiments, the mineral composition and relative content of rock samples are shown in Figure 2.

The mineral composition of the rock not only affects the physical and chemical effects during the formation of the rock but also affects seepage in the pores by affecting the pore structure. Quartz and feldspar edges are more likely to generate inter-granular pores with a larger radius, seepage easily passes through such pores, and various impurities and cements fill the pores and throats to hinder seepage.

The rock particle-size classification standard and rock particle-size structure parameters were derived through X-ray diffraction, casting thin sections and observational statistics of rock samples. It was found that for rock sample sorting, the particles were mainly medium sand, some coarse sand, and some fine sand. The percentages were 67.74%, 30.92% and 1.32%, respectively. The weathering degree of feldspar was relatively low in structure. The largest particle size was 1.45 mm, and the main particle-size range was 0.35 mm to 0.8 mm. The contact mode of the particles was mainly line contact.

2.3. Rock Pore Characteristics Analysis

2.3.1. Pore Types

Based on the testing and analysis of the rock samples based on casting thin section technology, the pore types mainly included residual intergranular pores, intergranular dissolved pores, intragranular dissolved pores, mold pores, and microcracks. Among them, the most important pore types were residual intergranular pores and intergranular dissolved pores.

Intergranular pores refer to pores surrounded by particles in the rock. During diagenesis, the intergranular pores were squeezed under the action of in situ stress and filled with interstitials to form residual intergranular pores. The current research also mostly refers to residual intergranular pores. The shapes were mostly triangles, polygons, and irregular shapes of varying sizes (a₁, a₂). The average pore diameter was 10~150 μm.

Intergranular dissolution pores showed that the edges of the particles were corroded in a bay-like or irregular shape, mainly from feldspar dissolution (b₁, b₂), followed by interstitials (c₁, c₂), (d₁, d₂), rock debris (e₁, e₂), and the edge dissolution of quartz particles (f₁, f₂).

The intragranular dissolved pores were mainly formed by partial dissolution of rock particles (g₁, g₂). The mold pores were formed after the soluble minerals were completely dissolved (h₁, h₂), retaining the original shape and contour of the particles.

There were two main types of microcracks: one was distributed along the inside of the debris or the edges of the particles, and the other was distributed between the particles [31]. The width was generally from a few tenths of a micrometer to a few tens of micrometers (i₁, i₂).

2.3.2. Throat Types

The throat is a connecting channel between two pores which has a constricted, sheet, curved, or tube bundle shape. A pore can connect to one throat or multiple throats. The pores can pass through the throats and connect with multiple throats. The other pores are connected.

Necked throats were more common in rock supported by grains. Compaction reduced the pores and throats, and the sizes of the pores and throats were significantly different, as shown in Figure 3a.

Flaky throats and curved flaky throats were more common in linear contact and concave–convex contact rocks. After compaction, the particles appeared to have inlaid contacts, and the throat was thin and curved, as shown in Figure 3b.

2.4. Porosity Analysis of Rock

In this research, high-pressure mercury intrusion analysis was performed on 5 samples, and the porosity parameters of the 5 samples are listed in

Table 1. Measurement of high-pressure mercury in 5 rock samples.

	1	2	3	4	5
Porosity (%)	22.31	19.64	20.47	18.56	16.35
Permeability (10−3 μm2)	302.02	235.18	204.35	141.52	98.35

in conjunction with the subsequent research in this study.

According to the results of mercury intrusion experiments, the porosity range of the rock samples is 16.35~22.31%, which is the pre-value segmentation of the pores in the subsequent CT experiments as the division standard.

3. Rock 3D Structure Creation

A previous study clarified the two-dimensional characteristics of the mineral composition and micropore structure of the rock samples studied in this research. However, the two-dimensional characteristics of the pore structure alone cannot completely and clearly describe the characteristics of the microscopic pore structure of the rock, nor can they be used to simulate the microscopic seepage under the three-dimensional pore structure. Therefore, this study combined the CT-scanning experiment to build a 3D digital core based on this and then created a 3D rock structure. This lays the foundation for subsequent microscopic seepage simulation research under the three-dimensional pore structure.

To avoid ambiguity, a simple explanation of digital core technology is as follows: digital core technology is based on two-dimensional scanning electron microscopy or three-dimensional CT scan images, using computer image processing technology to complete the three-dimensional reconstruction of the rock through certain algorithms; that is, to create a real three-dimensional rock structure.

3.1. CT Experiment

The CT equipment used in this experiment is the SkyScan2211 multirange nano-CT system. The equipment is manufactured by a precision instrument equipment factory in Shanghai, China.

Since the resolution of the scanned two-dimensional image has a negative correlation with the size of the rock sample, the smaller the rock sample is, the higher the scanned resolution. Therefore, this study combined previous research and references to cut the rock into a cube sample with a side length of 10 mm. In this way, pores with a radius of more than 1 μm could be scanned, which conforms to previous research and guarantees the scanning resolution.

Before the formal CT-scanning experiment, a pre-experimental scan of the rock sample was needed. In the pre-experiment, various instrument parameters were calibrated to detect in advance whether the resolution of the CT scan image met the requirements, which also saved time and reduced the scrap rate in the formal experiment. After the CT scan was completed, 1200 continuous two-dimensional grayscale images and image information were obtained. Figure 4 shows the two-dimensional grayscale image displayed after processing.

The image in Figure 4 consists of gray areas of varying degrees. The image has a black background. The irregular black parts in the circular rock sample indicate pores, and the blacker the image is, the lower the density. Gray indicates the rock matrix. The white part represents rock minerals. The figure shows that the sample had developed a certain amount of larger pores and contained more minerals.

3.2. 3D Digital Core Model Creation Process

Avizo software [32], was used to process a series of two-dimensional grayscale images scanned by CT, synthesize a three-dimensional model, and reconstruct a three-dimensional digital core model of rock pores. That is, the continuous two-dimensional CT data scanned by Avizo software (it comes with processing software for the device after CT scan) were superimposed, the pixels were connected to form a three-dimensional data volume, and then the three-dimensional data volume was visualized. The reconstruction of the three-dimensional digital core model from the two-dimensional gray image data is shown in Figure 5.

The process of constructing a three-dimensional digital core was to use CT scanning to obtain a rock sample and obtain a series of two-dimensional grayscale images that could reflect the microscopic pore characteristics. All two-dimensional grayscale images were processed by Avizo large-scale visualization software, and a series of two-dimensional grayscale images were superimposed. After multiple restoration and reconstruction processes, a three-dimensional digital core was created, as shown in Figure 5.

3.3. CT Image Processing

3.3.1. Filter Processing of CT Images

The scanned image had system noise, which not only reduced the quality of the image but also affected the subsequent research and analysis and required noise reduction processing. Commonly used methods include Gaussian filtering, median filtering, and mean filtering. Among them, median filtering is the best processing method and the most widely used processing method. The median filter can not only smooth the noise efficiently but can also keep the image boundary clear. Therefore, this research selected the median filter to reduce the noise of the original gray image. Figure 6 shows the effect diagrams after the three filtering processes.

3.3.2. CT Image Segmentation Processing

Although larger pores could be directly seen, smaller pores could not be seen directly. The boundary between the pores and rock matrix was blurred, and it was impossible to visually distinguish the edges of the pores and the rock matrix. The different structures and material compositions of the rock were reflected in different gray levels in the two-dimensional gray image. At the boundaries of the pores, rock matrix, and high-density mineral components (high-density nodules), the grayscale changed suddenly. Although the naked eye could not find the location of the change, computer technology could identify the change. Therefore, a special algorithm was used to process and segment the pores and rock matrix in the two-dimensional grayscale image. This research used the morphological watershed method [33], to segment CT images based on the measured porosity in a previous study.

3.3.3. Opening Operation and Closing Operation of CT Images

The two iterative algorithms of the opening operation and closing operation could soften the sudden change in brightness of the image, weaken the brightness of the bright spots in the image, and connect the discontinuous edges on the image. In addition, the effect of noise could be eliminated, and the edges of the image could be made clearer. In Figure 7, the rock matrix is dark blue and purple, the pores are white, and the high-density nodules are red. The processing results are shown in Figure 7b.

3.4. Extraction and Characterization of the Representative Elementary Volume

3.4.1. Representative Elementary Volume Analysis

Scholars have studied the effect of fractured rock series and different pore volume scales through various means [34], and the representative elementary volume (REV) was used in this study. REV can represent the overall characteristics of the rock and can meet the requirements of the study of seepage in the pore structure of the rock mass. Therefore, this research used Avizo three-dimensional visualization software to statistically extract the porosity of different volumes of cubes at different positions on the three-dimensional digital core for REV analysis.

With the aid of Avizo software, the coordinate axis was calibrated for the three-dimensional digital core model, and three different positions were selected in the three-dimensional digital core model, denoted as O₁, O₂, and O₃ points. Taking these three points as the center point of the cube unit, its size was increased along the coordinate axis, and the porosity was calculated. In this way, the porosity curves of three cube elements with different positions and different volumes could be obtained, and the characteristic unit REV could be found by analyzing the change law of porosity with cube elements. The cutting diagram of the characterization unit REV is shown in Figure 8, and the porosity value of the rock sample is shown in

Table 2. Unit cell porosity with O₁, O₂ and O₃ as the central voxel point (%).

Voxel Point	Unit Body Side Length Size/μm
	100	200	300	400	500	600	700	800
O1	14.49	16.36	17.55	19.88	19.85	19.7	19.83	19.71
O2	14.25	15.88	18.15	19.85	19.76	19.82	19.75	19.85
O3	14.87	15.46	17.95	19.75	19.68	19.98	19.7	19.75

. The obtained porosity and the cell side length size of the relationship curve are shown in Figure 9.

Figure 9 shows that as the cube unit increased, the corresponding porosity also changed. When the cube unit increased to a certain value, the porosity of the rock sample gradually stabilized. The cube unit of this value was the desired REV. The porosity of the rock sample studied in this research tended to be stable when the side length was 400 μm, this value was the side length of the REV. At this time, the average porosity of the rock sample was 19.79% and the error with the measured porosity was 0.76%, which met the research requirements. Therefore, this research selected the REV with a size of 400 μm × 400 μm × 400 μm from the three-dimensional digital core model as the research object, built the microscopic pore finite element model, and conducted the seepage simulation study.

3.4.2. Micropore Model Creation of Representative Elementary Volume Analysis

The two-dimensional grayscale image and distribution information data scanned by the CT software system were imported into Avizo software and the three-dimensional visualization of the grayscale image of the rock sample was realized through the Volume Rendering module, as shown in Figure 10a. The Extract Subvolume module was used to select the size and position of the REV. The REV was 400 μm × 400 μm × 400 μm voxels. For the center of the voxel, the O₂ point was selected to segment the REV, as shown in Figure 10b. Since the resolution of the CT-scanning device was 3 μm/voxel, the true three-dimensional size of the rock mass was 1200 μm × 1200 μm × 1200 μm.

According to previous research, REV was extracted from the 3D digital core model, and the whole REV was divided into pore structure and skeleton structure (rock matrix). Blue represents the rock skeleton (rock matrix), and yellow represents the pore structure of the rock. Figure 11 shows the rock micropore structure model based on REV.

3.4.3. Representative Elementary Volume Pore Structure Connectivity Analysis

For the quantitative characterization of pore structure and seepage simulation to proceed smoothly, pore connectivity analysis was required to distinguish between connected pores and closed pores. The connected pores were connected through Axis Connectivity and then the Separate Objects module was used to separate the skeleton (matrix) of the connected pores from the pore structure. If one wants to show closed pores, one needs to use the Lable Analysis module to adjust the parameters.

As shown in Figure 12, blue pores represent connected pores and pores in other colors represent isolated pores. The total pores are the sum of the isolated pores and the connected pores. The total pores accounted for 19.79% (porosity) of the target study area. Among them, the connected pores accounted for the largest proportion, approximately 75.63% of the total pores, and the isolated pores accounted for approximately 24.37% of the total pores, as shown in

Table 3. Statistical table of model pore connectivity.

Pore Type	Total Pores	Isolated Pores	Connected Pores
Target area	19.79%	4.82%	14.97%
Total pores	1	24.37%	75.63%

. Figure 12 and Table 3 show that the micron-sized pores had better connectivity in three-dimensional space.

4. Results (Numerical Simulation under Different Conditions)

In this study, the following three common conditions were used as examples to study the multidirectional entrance and multidirectional exit seepage.

• The positive direction of “X + Y” was set as the entrance, the negative direction of “X + Y” was set as the exit, and seepage was set in both directions.
• The positive direction of “X + Y + Z” was set as the entrance, the negative direction of “X + Y + Z” was set as the exit, and seepage was set in the three directions simultaneously.
• The positive directions of “X, −X, Y, −Y, Z,−Z” were set as the entrance and the other five wall faces were set as the exit faces for seepage.

In this paper, the ANSYS software Fluent module is used to perform boundary conditions and solver settings for the finite element model of the microscopic pore structure. Combined with previous studies and references, it can be seen that 3D microscopic seepage studies generally do not consider the influence of pore shape and seepage effects on confining pressure changes. Based on the previous research and the research needs of this paper, the pressure inlet and the pressure outlet are set, and the pore wall is a non-slip boundary without the influence of confining pressure. The results of previous studies show that when the Reynolds number is lower than a certain critical value, the fluid flow is laminar; otherwise, the fluid flow is turbulent, and the turbulent flow model is more suitable for seepage in microscopic pores. In this paper, the common k-ε turbulence model is used to describe the flow in the pores. To simplify the simulation and analysis process, the following assumptions are made:

(1)

Water only flows in the pores of the rock mass and will not penetrate into the matrix;

(2)

Water is an incompressible fluid that flows continuously, and the temperature is constant during the flow of pores;

(3)

Water is only affected by gravity and pressure.

4.1. Simulation Study of Two-Way Entrance and Two-Way Exit Seepage

Under the condition that the entrance pressure was 5 MPa and the exit pressure was 0, the model was subjected to multidirectional seepage simulation, and the results are shown in Figure 13.

Figure 13 shows that, compared with the one-way channel seepage, the seepage velocity in the case of the “X + Y” two-way channel was significantly higher. The average velocity in the X direction was 31.02 m/s, which was 72.8% higher than the average velocity of the X forward single-channel seepage. The average velocity in the Y direction was 30.69 m/s, which was 33.6% higher than the average velocity of the single-channel seepage, characteristic in the Y direction. The water flows in the X and Y directions affected each other, and part of the water flow originally flowing out from the negative X direction flowed out from the negative Y direction. Part of the water flow originally flowing out from the negative direction of Y flowed out from the negative direction of X to form a stream. This created the pores through which a large amount of water flowed, but now only a small part of the water flowed through.

The velocity contours of the X-direction section and Y-direction section of the “X + Y” bidirectional channel seepage is shown in Figure 14, and the velocity changes on the corresponding section were obtained. The velocity of the cross-section in the X direction and Y direction of the “X + Y” bidirectional channel seepage and the cross-sectional velocity of the unidirectional channel are shown in

Table 4. The average velocity distribution of each section of two-way channel seepage and one-way channel seepage.

		X = 0 Section (m/s)	X = 0.2 Section (m/s)	X = 0.4 Section (m/s)	X = 0.6 Section (m/s)	X = 0.8 Section (m/s)	X = 1.0 Section (m/s)	X = 1.2 Section (m/s)
“X + Y” bidirectional	X	15.58	19.16	26.59	34.02	47.21	43.66	30.92
	Y	23.20	38.65	37.89	28.52	31.94	27.13	27.48
Unidirectional	X	9.03	16.43	21.73	22.14	23.09	18.86	14.36
	Y	16.08	22.80	25.99	27.71	25.87	24.20	18.16

.

The seepage velocities of different cross-sections in the X and Y directions are compared in Figure 14 with the seepage velocities of the various cross-sections of the X and Y forward single channels. Not only the size of the flow velocity in each section in the X and Y directions is shown but also the change in the velocity gradient in each section. This is the same as the velocity distribution streamline and velocity distribution vector trend in Figure 13. The flow lines in the pore channels through which the fluids in both directions pass were denser, and the flow velocity was greater. Flowing through the same pore channel, compared with unidirectional flow, the flow velocity of bidirectional flow changed greatly. Table 4 shows that the seepage rate of the “X + Y” two-way channel was higher than that of the one-way seepage in each section.

To further quantitatively study the seepage characteristic of the pore model, the permeability (according to existing research [35], it is considered and calculated in combination with Darcy’s law) and velocities of the X and Y exit faces under the conditions of the “X + Y” two-way channel seepage were extracted and calculated and plotted into tables, as shown in

Table 5. The permeability of the X and Y exit surfaces of the “X + Y” bidirectional channel seepage.

Entrance Surface	Permeability (10−3 μm2)	Percentage (%)	Exit Surface	Permeability (10−3 μm2)	Percentage (%)
X	337	43.31	X	391	50.33
Y	441	56.69	Y	386	49.67

and

Table 6. The average velocity distribution of the X and Y exit surfaces of the “X + Y” bidirectional channel seepage.

Entrance Surface	Entrance Velocity (m/s)	Exit Surface	Exit Velocity (m/s)
X + Y	15.58	X	23.20
	30.92	Y	27.48

.

Table 5 shows that during the seepage process of the “X + Y” bidirectional channel, the seepage characteristic in the X direction and Y direction crossed flows. During the stream flow, the permeability of the X exit accounted for a larger proportion than that of the Y exit. This shows that there were pore channels with larger radii and shorter lengths from the Y entrance surface to the X exit surface, which were easier for water to pass through than the pore channels from the X entrance surface to the Y exit surface.

In addition, the permeability data in Table 5 shows that its value is relatively close to the permeability of sample 1 in Table 1, but both are greater than the values in Table 1. This may be due to the assumption that the numerical simulation is an ideal state, and the considered Darcy’s law usually assumes that the flow rate of the seepage channel is the same, while the seepage channel has different subtle changes in the actual seepage process.

Table 6 shows that during the seepage process of the “X + Y” bidirectional channel, the seepage velocities in the X direction and Y direction were also different. In addition, the seepage velocity in the Y direction was greater, indicating that the pore channels in the Y direction were more tortuous, and the radius of the pore channels at different positions changed more.

4.2. Three-Way Entrance and Three-Way Exit Seepage Simulation Study

Figure 15 shows that the seepage velocity of the “X + Y + Z” multidirectional channel had a significant increase compared with the single channel seepage. The average velocity in the X direction in the “X + Y + Z” multidirectional channel seepage was 34.30 m/s, which was 91.14% higher than the average velocity of the X forward single channel seepage. The average velocity in the Y direction of the “X + Y + Z” multidirectional channel seepage was 34.87 m/s, which was 51.81% higher than the average velocity of the Y positive single channel seepage. The average velocity in the Z direction of the “X + Y + Z” multidirectional channel seepage was 33.89 m/s, which was 66.29% higher than the average velocity of the Z positive single channel seepage. The water flows in the X, Y, and Z directions affected each other, and part of the water entering from the positive direction of X flowed out from the negative direction of Y and the negative direction of Z. The water entering in the positive direction of Y and Z also had the same characteristics as the direction of X, and the water in the three directions, X, Y, and Z, flowed with each other.

Compared with the seepage characteristic in the “X + Y” bidirectional channel, the seepage velocity and path of the “X + Y + Z” multidirectional channel underwent greater change. Affected by the seepage pressure in the Z direction, the seepage paths in the positive X and Y directions became shorter, and the seepage velocity increased significantly. For example, the “X + Y” bidirectional channel seepage had an average seepage velocity of 30.86 m/s, and the “X + Y + Z” multichannel seepage had an average seepage velocity of 34.35 m/s.

The characteristics of the pressure field and velocity field of the “X + Y” and “X + Y + Z” multidirectional channel seepage simulations not only confirmed the characteristics of simulation research by [36,37], but also the seepage characteristic of the research was more comprehensive.

The velocity cloud diagram of the X-, Y-, and Z-direction cross-sections of the “X + Y + Z” multidirectional channel seepage is shown in Figure 16.

Table 7. The average velocity distribution of each section of multidirectional channel seepage and unidirectional channel seepage.

		X = 0 Section (m/s)	X = 0.2 Section (m/s)	X = 0.4 Section (m/s)	X = 0.6 Section (m/s)	X = 0.8 Section (m/s)	X = 1.0 Section (m/s)	X = 1.2 Section (m/s)
“X + Y + Z” multi-directional	X	21.93	24.15	32.17	34.75	48.04	47.98	31.10
	Y	31.22	40.20	36.69	40.88	31.92	29.45	33.72
	Z	12.31	26.76	36.14	45.89	30.09	37.76	48.31
Unidirectional	X	9.03	16.43	21.73	22.14	23.09	18.86	14.36
	Y	16.08	22.80	25.99	27.71	25.87	24.20	18.16
	Z	12.92	24.01	23.24	22.47	23.75	21.57	14.71

shows the velocity of each section of the “X + Y + Z” multidirectional channel seepage in the X, Y, and Z directions and the velocity of each section of the unidirectional channel.

The seepage velocities of different cross-sections in the X, Y, and Z directions are compared in Figure 16 with the seepage velocities of various cross-sections in the X, Y, and Z forward single channels. The seepage velocity in the three directions of the “X + Y + Z” multidirectional channel was significantly higher, and the velocity change range of each section was also changed. This was the same as the velocity distribution streamline and velocity distribution vector trend in Figure 15.

The fluid originally flowing out in the negative X direction was divided into three parts: flowing out from the X negative direction, the Y negative direction and the Z negative direction. The fluid originally flowing out in the X negative direction was divided into three parts, which flowed out from the Y negative direction, X negative direction, and Z negative direction. The fluid originally flowing out in the Z negative direction was divided into three parts, which flowed out from the Z negative direction, X negative direction, and Y negative direction.

The flow lines in the pore channels where the fluids in the three directions converge were denser, and the speed was also greater. Table 7 shows that the “X + Y + Z” multidirectional channel seepage had a significant increase in the speed of each cross-section compared to the one-way seepage and the two-way channel seepage. The velocity sequence of single-way channel seepage, “X + Y” two-way channel seepage, and “X + Y + Z” multiway channel seepage in the X, Y, and Z directions was ordered as “X + Y + Z” multiway channel seepage > “X + Y” two-way channel seepage >X, Y, Z one-way seepage.

To study the seepage characteristic of the pore model, the permeability, and velocities of the X, Y, and Z exit surfaces under the condition of “X + Y + Z” multidirectional channel seepage were extracted and calculated, as shown in

Table 8. The permeability of the X, Y and Z exit surfaces of the “X + Y + Z” bidirectional channel seepage.

Entrance Surface	Permeability (10−3 μm2)	Percentage (%)	Exit Surface	Permeability (10−3 μm2)	Percentage (%)
X	445	37.88	X	394	33.58
Y	534	45.49	Y	421	35.88
Z	195	16.63	Z	358	30.54

and

Table 9. The average velocity distribution of the X and Y exit surfaces of the “X + Y” two-way channel seepage.

Entrance Surface	Entrance Velocity (m/s)	Exit Surface	Exit Velocity (m/s)
X + Y + Z	21.93	X	31.10
	31.22	Y	33.72
	12.31	Z	48.31

.

Table 8 shows that during the seepage process of the “X + Y + Z” multidirectional channel, the seepage in the three directions, X, Y, and Z crossed flows. In the seepage process, the proportions of the permeability of the three exits, X, Y, and Z were ranked as Y > X > Z. This indicates that the fluid entering from the entrance surface flowed out more from the Y exit. That is, under the three-way seepage pressure, there was a pore channel “X (Z) positive → Y negative”, which made it easier for fluid to pass through. The Z exit had the smallest permeability, indicating that the pore channels from the X, Y, and Z positive directions to the Z negative direction were more tortuous, and there were pore throats with a smaller radius, which made it difficult for fluid to pass. In addition, the permeability data in Table 8 shows that its value is relatively close to the permeability of sample 1 in Table 1, but most of the values are greater than those in Table 1. This may be due to the assumption that the numerical simulation is an ideal state, and the considered Darcy’s law usually assumes that the flow rate of the seepage channel is the same, while the seepage channel has different subtle changes in the actual seepage process. Table 9 shows that the average exit velocity in the three directions was sorted as Z > Y > X.

4.3. Simulation Study of One-Way Entrance without Wall Seepage

“X, −X, Y, −Y, Z, −Z” were set as the entrance and the other five wall surfaces were set as the exit surfaces. The flow conditions and velocity distribution were obtained by simulation, and the simulation results are shown in Figure 17 and Figure 18.

Figure 17 and Figure 18 show that due to the different development statuses of the pore structure at different entrance surfaces, the development status was also different in different directions. Therefore, when multidirectional channel seepage was performed, different entrance faces were selected, and the seepage channels through which the fluid passed were also different. Compared with unidirectional seepage, there were more seepage paths for water flow due to the presence of five exit faces. However, there were more pore channels without fluid passing through them. The seepage characteristic of multidirectional channels was the same as that of unidirectional, two-way, and three-way channels. That is, in the process of fluid seepage, the water flow chose pore channels with larger radii and shorter seepage paths to pass through. In a single pore, the flow velocity was largest closer to the center of the channel, the flow velocity was smaller closer to the pore wall, and the seepage velocity was zero at the position close to the pore wall.

The X positive direction was taken as the entrance and the other five wall surfaces were taken as the exits. The entrance face pressure was set to 5 MPa and the exit face pressure was set to 0 to perform the seepage simulation. The results are shown in Figure 19, Figure 20 and Figure 21.

Figure 20 shows that the pressure distribution of each section in the X direction shows that the pressure distribution on the same section is not uniform. The pressure cloud diagram of each section shows a trend in which the pressure in the middle area of the pores is larger and diffuses to the surroundings and gradually decreases. Among them, X = 0 mm (pressure entrance), and most of the area of the cross-section color is dark red. However, it appears blue in two areas. These two areas are the corners of the pore model. No fluid passes through the corners. As the distance from the entrance end increases, the area of the dark red area of the X = 1.2 mm cross-section is greatly reduced, and the yellow and green areas increase. From the entrance face to the exit face, the pressure value gradually decreases along the seepage path. Figure 21 shows that the color of the cross-section near the entrance surface has yellow, green, and blue gradient distributions, indicating that the entrance flow rate is higher. The color of the cross-section near the exit surface is mostly blue, indicating that the exit velocity is low. After the fluid enters the model from the entrance surface, it preferentially flows through the pore channels with a shorter distance and a larger radius. Most fluid flows out at the Y and Z exits, and only a small part of the fluid flows out at the negative X exit.

To study the seepage characteristic of the pore model, the permeability, and velocities of the exit surfaces in the X, Y, and Z directions under multidirectional channel seepage conditions are listed in

Table 10. Permeability distribution of each exit surface under different seepage entrance surface conditions.

Entrance Surface	Permeability (10−3 μm2)	Exit Surface	Permeability (10−3 μm2)	Percentage (%)
X	191	−X	50.8	10.65
		Y	50.2	10.53
		−Y	138.4	29.02
		Z	120.9	25.35
		−Z	116.5	24.43
−X	188	X	53.7	11.43
		Y	218.4	46.47
		−Y	43.2	9.18
		Z	41.2	8.77
		−Z	113.3	24.11
Y	199	X	50.8	10.21
		−X	201.8	40.60
		−Y	85.5	17.21
		Z	68.8	13.84
		−Z	89.9	18.10
−Y	186	X	119.6	25.69
		−X	73.2	15.72
		Y	87.1	18.70
		Z	107.1	23.00
		−Z	78.6	16.88
Z	149	X	99.1	26.63
		−X	35.6	9.56
		Y	78.3	21.04
		−Y	110.2	29.62
		−Z	48.9	13.15
−Z	217	X	112.3	20.72
		−X	128.1	23.64
		Y	149.3	27.56
		−Y	89.8	16.57
		Z	62.3	11.49

and

Table 11. The distribution of seepage velocity on each exit surface under different seepage entrance surface conditions.

Entrance Surface	Average Entrance Velocity (m/s)	Exit Surface	Exit Velocity (m/s)	Average Exit Velocity (m/s)
X	34.583	−X	7.026	20.6018
		Y	12.702
		−Y	15.977
		Z	35.866
		−Z	31.438
−X	42.748	X	5.269	17.6374
		Y	26.690
		−Y	24.719
		Z	7.797
		−Z	23.712
Y	45.422	X	12.682	20.865
		−X	29.116
		−Y	10.168
		Z	30.282
		−Z	22.077
−Y	43.740	X	28.140	20.1734
		−X	17.460
		Y	9.063
		Z	27.752
		−Z	18.452
Z	40.490	X	23.620	13.4746
		−X	5.464
		Y	14.505
		−Y	17.860
		−Z	5.924
−Z	46.331	X	18.280	22.8756
		−X	37.607
		Y	31.675
		−Y	18.366
		Z	8.450

.

Table 10 shows that due to the differences in the development of the pore structure in various directions, the fluid flow direction tended to the exit surface with a larger pore channel radius and a shorter seepage path. By calculating the percentage, it can be found that there was always an optimal seepage exit surface among all exit surfaces corresponding to each seepage entrance surface.

The permeability data in Table 10 have a large change, which may be due to the change of the calculated value due to the change of the size of the seepage channel in the complex seepage process. The data changes in Table 10 are generally about two to three times (except for very few data), which is similar to the data in Table 1 (about three times changes). In addition, if the cell size in this study is 500 μm or other sizes, the permeability data may change accordingly, which will be further studied in the future.

A total of 29.02% of the fluid entering from the +X surface exited from the exit of the −Y surface. The order of the proportions of the fluid on each exit surface was ranked as −Y > +Z > −Z > −X > +Y. This shows that the −Y exit surface was the best percolation exit surface of the +X entrance.

A total of 46.47% of the fluid entering from the −X surface flowed out from the +Y surface exit, and the proportion of the fluid at each exit surface was in the order of Y > −Z > X > −Y > +Z. This shows that the +Y exit surface was the best percolation exit surface of the −X entrance.

A total of 40.60% of the fluid entering from the +Y surface flowed out from the exit of the −X surface, and the proportions of the fluid at each exit surface were in the order of −X > −Z > −Y > +Z > X, which indicates that the −X exit surface was +, which was the best percolation exit surface for the Y entrance.

A total of 25.69% of the fluid entering from the −Y surface flowed out from the +X exit surface, and the proportions of the fluid on each exit surface were sorted as +X > +Z > Y > −Z > −X. This shows that the +X exit surface was the best percolation exit surface for the −Y entrance.

A total of 29.62% of the fluid entering from the +Z plane flowed out from the exit of the −Y plane, and the proportion of the fluid at each exit plane was in the order of −Y > +X > +Y > −Z > −X. This shows that the −Y exit surface was the best percolation exit surface of the +Z entrance.

A total of 27.56% of the fluid that entered from the −Z plane flowed out from the exit of the +Y plane, and the proportions of the fluid at each exit surface were in the order of +Y > −X > +X > −Y > Z. This shows that the +Y exit surface was the best percolation exit surface of the −Z entrance.

The pore channel from the entrance surface to the optimal percolation exit surface is called the optimal percolation surface channel. The best seepage surface channels corresponding to different entrance surfaces are +X → −Y, −X → +Y, +Y → −X, −Y → +X, +Z → −Y, and −Z → +Y. Among the six groups of best seepage surface channels, up to 46.47% of the fluid passed through the −X → + Y seepage surface channel. This value was the highest proportion of fluid passing through the seepage channels corresponding to all entrance surfaces. This pore channel was the easiest channel for fluid to pass through.

From Table 11, it was found that due to the difference in pore structure development, the corresponding seepage velocities of different entrance faces were quite different. The velocity of the seepage entrance was between 34.58 m/s and 46.33 m/s, and the velocity of the seepage exit was between 13.47 m/s and 22.88 m/s. The order of the speed of each entrance surface was −Z > +Y > −Y > −X > +Z > +X. The factor that affected the velocity of the entrance surface was the area of the pores participating in the seepage. The pore area of the entrance surface participating in the seepage characteristic was negatively correlated with the seepage velocity; that is, the +X pore area participating in the seepage characteristic in each seepage entrance surface was larger, and the pore development was better.

5. Discussion (Engineering Applications)

To verify the reliability of the research results, actual tests were carried out in this project.

In combination with the actual work and research conditions, the technicians at the engineering site came to the conclusion through discussion and research that: if the seepage condition of the large-scale rock mass at the engineering site is to be obtained, the minimum research unit should be controlled in the range similar to the particle size of the rock sample and choose the smallest size possible within the range. Therefore, a cell size of 0.4 mm was ultimately chosen. This both meets the range requirements (0.35–0.8 mm) and ensures a smaller element model size.

It should be noted here that the engineering site already has geological survey reports on various aspects of rock heterogeneity analysis, analysis of different strata, and hydrogeological conditions. If the specific conditions of large-scale in-situ rocks are to be obtained according to the small unit model, the size of the unit model should be within the range of the particle size of the rock sample that has been measured (0.35–0.8 mm).

To ensure the accuracy and reliability, the size of the element model should be selected as small as possible, so that it can be finally calculated into the large-scale model (when the size is too large, accidental phenomena are prone to occur, and if a large number of accidental phenomena are calculated in the large-scale model, there will be a large error).

After selecting the unit model size of 0.4 mm, the engineering field technicians discussed the research results of this manuscript in combination with the existing multi-faceted geological and hydrological survey reports. Due to the need to consider the heterogeneity of different rock layers, the regionality of water influx, and the effect of in-situ stress under this geological condition, the engineering site technicians classified the different rock conditions in the geological report, and based on the results obtained in this study (for example: if it is similar to the unit model taken in the study, the results of the manuscript will be applied; if it is different from the unit model taken in the study, the inverse calculation will be carried out through the study results of the manuscript), and the final comprehensive study has obtained the distribution of the water inflow distribution of the rock mass in the project. According to the distribution of water inflow, etc., the installation of drainage pipes and other equipment can be further guided.

In the project, three test points were selected and monitored to verify the feasibility of the research conclusions. For the three areas with small, medium, and large water inflows in a state of no drainage, three measuring points, 1, 2, and 3, were set, respectively.

According to the characteristics of the rock distribution and hydrological conditions in the area, combined with the aforementioned analysis of the seepage characteristic, the site layout of drainage equipment was carried out for the three test points (because drainage pipes and other equipment needed to be adjusted at any given time according to the actual situation on site, the layout of the drainage equipment plan is not provided). Then, the groundwater drainage operation began and the total water influx was recorded at different times.

The water output at different time periods of the three measuring points before and after drainage was compared, as shown in

Table 12. Water inflow at different measuring points and times when drainage equipment was not arranged (m³).

	Test Point 1	Test Point 2	Test Point 3
1 day	886.6	2031.5	4178.2
7 days	251.3	13,894.1	29,206.4
15 days	13,245.4	29,987.3	62,759.6
30 days	26,467.9	60,026.3	125,347.5

and

Table 13. The water inflow (m³) at different measuring points and times after the drainage equipment was arranged (m³) (the value in brackets is the ratio per thousand of the inflow of water under this condition to the inflow of water without drainage equipment under the same conditions).

	Test Point 1	Test Point 2	Test Point 3
1 day	3.1 (3.50‰)	6.8 (3.35‰)	8.2 (1.96‰)
7 days	2.4 (9.55‰)	8.5 (0.61‰)	7.9 (0.27‰)
15 days	1.9 (0.14‰)	7.1 (0.24‰)	6.9 (0.11‰)
30 days	2.5 (0.09‰)	5.9 (0.10‰)	7.5 (0.06‰)

.

Comparing Table 12 and Table 13, it can be found that after the drainage equipment was arranged according to this method the amount of overflowing groundwater was greatly reduced. The maximum point only reached 9.55‰ (less than one percent), and the minimum point only reached 0.06‰. Under different test points and different time conditions, the remaining water overflowing volume was minimal. This shows that the method effectively alleviated the problem of overflowing groundwater. This shows the feasibility of the conclusions of this paper.

It is worth mentioning that rock fissures have a certain influence on seepage [38]. Moreover, aggregates also have certain effects on rock or stone bodies [39,40], but these issues were not considered in this study. While this study has yielded good results, the authors believe that further research into issues such as rock fracturing can lead to better results in the future.

6. Conclusions

This research starts from the measurement and characterization of the real three-dimensional pore structure of rocks and studies the microscopic pore structure characteristics through scanning X-ray diffraction, casting thin sections, high-pressure mercury intrusion, and CT rock-scanning experiments. Using a computer numerical simulation method, the rock microscopic pore structure is restored in the computer model, and the pore structure model seepage simulation study at the microscale is carried out. Exploring the characteristic of fluid seepage in rocks effectively resolves the groundwater drainage problem in the project and further contributes to the sustainable development of water resources.

• The sandstone sample particles are mainly medium sand. The mineral composition of the rock sample contains mostly quartz, followed by feldspar, and the contents of cuttings, argillaceous base, and siliceous cement are lower.
• Rock pores are mainly residual intergranular pores, intergranular dissolved pores, intragranular dissolved pores, mold pores, and a small number of microcracks. The main types of throats are flaky and curved throats, followed by constricted throats, and there is a small amount of bundle-shaped throats. The rock developed multiscale pore throats and the pore radius is large. The porosity of the sample is 16.35–22.31%.
• The main distribution range of the pore radius is 30~80 μm. The main distribution range of the throat radius is 4~14 μm. The main distribution range of the pore coordination number is 6~17. The length of the throats is mainly distributed between 10 and ~300 μm. The pore shape factor is mainly distributed between 0.02 and ~0.07.
• For multidirectional seepage simulation studies, the seepage velocity is higher than that of the unidirectional flow. In the case of multidirectional seepage, fluids in multiple directions flow in series. Through multidirectional seepage simulation, one can also find the best seepage channel.
• Through verification in an actual project, groundwater diversion based on the research results of this paper can greatly reduce the amount of groundwater influx. This shows that the previous research conclusions can effectively guide groundwater drainage. In other words, it is feasible to study the characteristic of rock seepage in the area first and then to divert groundwater. This method effectively solves the pollution and safety hazards caused by overflowing groundwater in underground resource mining and other projects, reduces the pollution caused by grouting projects, and helps clean the production of underground resource mining and other projects. Stored groundwater can increase resource utilization and contribute to the sustainable development of water resources. At the same time, the implementation of this method ensures the safety of the working space of underground resource mining and other projects. After underground resource mining and other projects are completed, the operation space can still be fully utilized instead of directly discarded, which is in line with the concept of sustainable development. This method also has important significance for the clean production and safe construction of other underground projects. In addition, the study can also be further studied from the permeability and Darcy’s law to obtain more accurate research results.