Abstract
A coal mine underground reservoir, composed of a large number of goaf groups, is effective in achieving the protective utilization of coal mine water. Mastering the hydraulic connection between the goaf groups and the law of seepage is the key to the design and safe operation of coal mine underground reservoirs. Coal/rock sample seepage tests were conducted in Wanli No. 1 Mine in the Shendong mining area according to the evolution law of the lateral supporting pressure of the coal pillars. From the tests, the function of vertical stress in relation to the permeability of coal/rock samples was obtained and the distribution law of permeability of porous media in the goaf was revealed. On this basis, the computational fluid dynamics (CFD) numerical calculation of seepage properties of the coal/rock pillars was conducted. The results show that with the decrease of water level in the goaf and the increase of coal pillar width, the flow range and flow velocity keep decreasing. At the water level of 18.65, 28.65, and 38.65 m, the critical coal pillar widths for the goaf water to penetrate into the adjacent goaf are 30, 50, and 70 m, respectively. When the coal pillar width is less than 10 m, the water can bypass the elastic core area at the bottom coal seam to the adjacent goaf area; when the width exceeds 20 m, the water cannot enter the adjacent goaf area through the bottom coal pillar. On the basis of CFD simulation study on seepage properties of coal/rock pillars between goafs, this paper conducts calculations for all the goafs in the 4-2 coal-panel block of Wanli No. 1 Mine and reveals the seepage path and seepage law in the underground reservoir of this mine.
1. Introduction
An underground reservoir, formed by underground pores in sand and gravel, rock cracks or karst caves in a certain range, has the advantages of large storage capacity, small evaporation loss, and not occupying land surface. Japan’s Nagasaki has built the world’s first underground reservoir, though the storage capacity is only 9,000 m3 [1]. The aquifer storage and recovery (ASR) program implemented in the United States stores water by preinjecting it into the aquifer and extracts it when needed [2, 3]. Since the 1950s, Amsterdam, the Netherlands, has been diverting water from the Rhine River to canals and ditches on natural sand dunes to replenish groundwater [4]. These practices can shed light on water resource protection in western China where the ecological environment is fragile, especially areas with large-scale coal mining activities. The traditional water-preserving mining technology adopts the “blocking method” to protect the aquiclude, but this technology greatly limits the mining efficiency and recycling rate of coal mines. For this reason, Shenhua Group put forward a new concept of “diverting, storing, and utilizing water” with coal mine underground reservoirs to achieve water resource protection and utilization [5].
Wanli No. 1 Mine, located in Dongsheng District, Ordos City, Inner Mongolia, is in a typical arid and semiarid area with a fragile ecological environment. In view of this, water conservation is introduced to build underground reservoirs that have been adopted in the Shendong mining area. An underground reservoir water conservation project was implemented in the goaf of south and north wings of coal panel No. 4-2, where some coal pillars were left between adjacent goafs after the mining of working face. These structures had a significant impact on the hydraulic connection of this area. Scholars from China and abroad have, at present, studied the seepage properties of fractured rock from both theoretical and practical perspectives. Ma systematically studied the properties of broken rockmass and saturated broken rockmass, including the swelling and permeability properties, porosity compaction, and the deformation modulus. He further studied seepage property of broken rock in the goaf with numerical simulation [6]. Yang et al. suggested that water flow should experience three physical processes during the water inrush after the mining rockmass failure, i.e., Darcy laminar flow in the aquifer, non-Darcy high-speed flow in the rockmass, and the Navier-Stokes turbulent flow in the roadway [7, 8]. Based on the key layer theory of rock strata control studies, Cao et al. examined the properties of water flow in the fracture and its main channel distribution model [9] from the perspective of fracture pattern and fracture distribution in the overlaying rockmass. In searching for methods preventing water damage in the abandoned goaf of Yangquan No. 3 Mine, Xu et al. focused on predicting the roof water-conducting crack in the lateral boundary of the goaf and proposed a method to judge the lateral boundary of the roof water-conducting crack in relation to the scale of horizontal deformation [10]. Fink et al. studied the permeability properties of shale under the influence of both slip flow and pore elastic effects in the range of large pore pressure. They established a permeability model taking slip flow and pore elastic effects into account [11]. Ghanizadeh et al. first revealed that the permeability of early Jurassic shale in northern Germany decreased nonlinearly with the increase of effective stress [12].
The existing studies on seepage in the goaf mainly focus on the prevention and control of water hazards, but quite little has been done on the hydraulic connection, the seepage law, and the seepage path between the goaf of coal mine underground reservoirs. Therefore, taking the Wanli No. 1 Mine in the Shendong mining area as the engineering background, this study designed tests according to the evolution law of the lateral supporting pressure of the coal pillars in the area and conducted seepage tests of the coal/rock samples to obtain the function of the coal/rock permeability with the change of vertical stress. Following that, the CFD numerical calculation of the seepage properties in the goaf area was conducted, which revealed the hydraulic connection, the seepage law, and the seepage path between the adjacent goaf. The results could provide guidance for optimizing the design of coal mine underground reservoirs and hence ensuring the safe operation of them.
2. Coal/Rock Pillar Seepage Test
2.1. Test Equipment
In order to study the law of permeability change of coal pillars and its overlaying rock with coal/rock stress in Wanli No. 1 Mine, a three-axis servo seepage testing machine (Figure 1) was used for detecting coal/rock deformation under seepage action. It consists of five parts, namely, confining pressure loading system, axial loading system, gas injection system, flow measurement system, and computer control system.

2.2. Test Design
During mining of a coal seam, the variation in stress in a coal/rock mass significantly influences the deformation and failure characteristics of the coal/rock mass itself. Figure 2 shows how the variation of the stress field influences the distribution of mining-induced fractures after mining of the coal seam. After the coal seam is mined, the stresses in the surrounding rock are redistributed. From the center of the coal pillar to both sides of the coal pillar, the stress of coal and rock gradually changes from the triaxial stress state to the uniaxial stress state. The bearing capacity of the coal/rock mass near the coal pillar boundary decreases, and the stress is transferred to the interior of the coal/rock pillar, forming lateral abutment pressure. Under the effect of lateral abutment pressure, plastic failure occurs in the coal/rock masses, significantly raising the permeability of the coal masses. As shown in Figure 2, the coal pillar in the No. 4-2 Coal Panel of Wanli No. 1 Mine (henceforth referred to as “No. 4-2”) is about 20 m in width, with an elastic core area in the middle and a certain length of plastic zones on both sides. The permeability of the plastic failure zone is high, meaning that water in the goaf can easily penetrate it. Similarly, the permeability of the elastic core area is comparatively low, which means water on one side can only pass through the coal/rock pillar to the other side when the pressure discrepancy is large enough. Due to the narrow width of the upper rock pillar, the whole rock pillar is in a plastic state, and the elastic core area disappears.

The obtained coal/rock sample from No. 4-2 was cut and processed into specimens with a diameter of 50 mm and a height of 100 mm. The change of vertical stress was simulated by loading and unloading of axial pressure, and the change of horizontal stress was simulated by loading and unloading of confining pressure. Air tightness between the airway system, and the sample was checked before the test. A specimen coated with silicone rubber and enclosed within a heat shrinkable tube was put into the seepage experiment machine. Axial pressure was applied slightly first to press the specimen, then oil was injected to remove the air in the pressure chamber, and the oil return valve was closed. In the first stage of the test, the axial pressure and confining pressure of the coal/rock sample were loaded at a rate of 0.01 MPa/s to 6 MPa; in the second stage, the confining pressure (6 MPa) was kept constant, and the axial pressure was loaded to 30 MPa at a rate of 0.01 MPa/s; in the third stage, the axial pressure and confining pressure are both unloaded at the speed of 0.01 MPa/s, until the coal/rock sample is destroyed. In the second stage of the test process, that is, when the axial pressure and the confining pressure are both at 6 MPa, the pressure relief valve was adjusted to stabilize the cylinder pressure at 1 MPa, and then, the pipeline valve was opened. After the flow rate was stable, the rate was recorded. During the test, the flow rate was recorded for every 0.5 MPa change in axial pressure.
2.3. Test Results
Permeability is an important parameter in the law of solid-liquid coupling seepage. According to Darcy’s law, the permeability of the sample can be expressed as [13] where is permeability, m2; is helium flow rate, m3/s; is height of the coal sample, m; is helium dynamic viscosity (18.9 μPa·s); is cross-sectional area of the coal sample, m2; is inlet gas pressure, Pa; is outlet gas pressure (0.1 MPa).
By taking the specimen parameters in this test and the helium gas flow data into equation (1), the relationship was obtained between coal/rock permeability and axial pressure in pressure-loading and pressure-unloading stages, as shown in Figure 3.

(a) Loading stage

(b) Unloading stage
The coal/rock permeability was closely related to the coal/rock deformation and damage. As shown in Figure 3(a), when axial pressure increases while confining pressure is fixed, the gas flow channel is closed, and the coal/rock permeability and the gas flow rate decreases gradually. Therefore, the permeability is the lowest when the lateral support pressure of the coal pillar is at its peak position. Compared with rock sample, the initial permeability of coal sample is higher, but when axial pressure increases, the drop of coal permeability is also greater due to its softness and the compaction of its internal pores. Figure 3(b) shows that with the decrease of axial pressure and confining pressure, the permeability of coal/rock samples continues to increase. In the early stage of the unloading stage, due to the high confining pressure, the coal/rock sample still maintains its strength, with a slight permeability increase. In the later stage of the unloading stage, featured with the low confining pressure, the coal/rock sample becomes unstable. After reaching its strength limit, its permeability increases sharply. After the unstable state, the scale of increase in permeability of the coal sample is the smallest, only 263 times higher than its minimum permeability, while that of the siltstone is the largest, about 139,519 times higher than its own minimum value. In the pressure loading and unloading stages, the relationship between coal/rock sample permeability and axial pressure can be fitted with an exponential function, and the fitting results are shown in Table 1.
3. CFD Numerical Calculation of Seepage Properties of the Coal/Rock Pillar
3.1. Establishment of Numerical Calculation Model
3.1.1. Geometric Model and Boundary Conditions
Strata of different hardness and thickness may play two controlling roles by affecting part of the strata above or the whole strata above, where the former is called the secondary key stratum and the latter the primary key stratum [14]. The sample in this simulation is a coal pillar from a section of 20 m in width between goaf No. 42111 and 42101 of Wanli No. 1 Mine. According to the ZK1416 borehole in the above area, the Key Strata Distinguishing software (KSPB) is used to examine the key strata of the overlaying rock in this area. The result is shown in Figure 4.

According to Figure 4, we established a goaf model (Figure 5) with a geometric structure of 115.65 m in height and 550 m in length. In terms of solid mechanics calculations, the top interface of the model is a free boundary, the side boundaries are roller bearing boundaries, and the bottom interface a fixed boundary. In terms of fluid calculation, the target areas of this model are and , where is a coal/rock pillar with its left boundary pressure at Pa and is the goaf area with its water level being . So the pressure on the left is . Assuming that the water ratio of the initial state of area is 0 with an initial pressure of 0 Pa, and the initial water level of area is , and the initial water pressure of the area is . The permeability is determined according to Table 1.

The parameter values are shown in Table 2.
3.1.2. Simulation Schemes
The model plans to calculate the width d of the coal pillar in this section and the water level in area. According to the abovementioned plan, the orthogonal simulation schemes are designed by setting the water level in area at 18.65, 28.65, and 38.65 m and the coal pillar width at 10, 20, 30, 40, 50, 60, and 70 m. The simulation schemes are shown in Table 3.
3.2. Analysis of Simulation Results
3.2.1. Distribution Law of Permeability in Goaf
By arranging survey lines at 2.55, 7.55, 12.55, 17.55, 22.55, 27.55, 32.55, and 37.55 m from the coal seam floor along the coal seam inclination, the permeability of each point was obtained and hence the changing pattern (Figure 6) of permeability of the overlaying rock coal pillars at different widths. It shows that stress on the coal/rock in the goaf has a significant impact on permeability, which is mainly manifested in the following five aspects: (1)When the coal pillar width is 10 m, the permeability below 12.55 m of the survey line drops sharply because the original fractures of the rockmass are compacted, which indicates that the area below 12.55 m is an elastic core area with extremely low permeability. As the coal pillar width increases, the permeability at higher survey lines drops significantly too, indicating that the height of the elastic core area is increasing. The permeability of the elastic core area of the rock pillar is between and m2 and that of the coal pillar is between and m2(2)The permeability of the rock above the elastic core area is relatively high, and the coal/rock body in this range is destroyed due to the reduction of horizontal stress. When the coal pillar width is 10 m, the permeability below 12.55 m is low, but the permeability above 17.55 m is higher, which indicates that the lateral boundaries of the water-conducting fractures on both sides of the coal pillar meet at the height between 12.55 m and 17.55 m from the coal seam floor. As the coal pillar width increases, the permeability in the plastic failure zone increases slightly, and the permeability in this range is between and m2(3)Under the support of the coal pillar, the coal/rock near the mining boundary is not fully compacted, and the permeability in this area is relatively high. Because the coal pillar width is 10 m which is smaller than boundary coal pillar width (40 m) in the model, the permeability above the goaf is distributed asymmetrically. As the coal pillar width increases, the above phenomenon gradually disappears with the permeability starting to distribute symmetrically. The peak permeability near the mining boundary is between and m2(4)In contrast with the coal/rock naturally deposited near the mining boundary of goaf, the permeability of the fractured coal/rock in the middle and lower parts of the goaf is relatively low, because it is recompacted. But such permeability (~ m2) is still greater than that of the coal pillar in the plastic failure zone(5)Affected by the mining of the lower coal seam, the stress of the upper coal/rock is relieved, but permeability of the coal/rock is increased. However, as the degree of stress unloading on the coal in this area is less than that on the coal near the mining boundary, its permeability (~ m2) is lower than that of the latter area, but significantly higher than the fractured coal/rock in the middle and lower parts of the goaf

(a) Coal pillar width is 10 m

(b) Coal pillar width is 30 m

(c) Coal pillar width is 50 m

(d) Coal pillar width is 70 m
In summary, the order of coal/rock permeability in the goaf is: natural accumulation area near the mining boundary > low stress area in the middle and upper part of the goaf > recompaction area in the middle and lower part of the goaf > plastic failure zone above the coal pillar > the elastic core area in the middle of the coal/rock pillar. From the perspective of the overall distribution of coal/rock permeability in the goaf, the boundary area is connected to the high permeability area in the middle and upper part of the goaf, forming a round hat-shaped structure in space, which corresponds to the “O-shape” theory and the “ellipse parabolic zone” theory in mining cracks studies.
3.2.2. The Law of Fluid Movement in the Goaf
The pressure gradient is the ratio of the pressure difference between two points in the flow area to the distance between them. According to the theory of fluid seepage in porous media, water seepage in the coal happens at a certain starting . The starting of seepage in coal/rock can be expressed as [16]:
In this formula, is the starting , Pa/m; is permeability, 10-3 μm2, which can be obtained through tests.
When , the water body is in a static state; when , the water body is in a flowing state. The ratio of to in the fluid calculation range is derived, and the flowing area, where , is colored and labeled, while the static area, where , is not colored and labeled. So we obtained the cloud map of r and fluid streamline distribution of different water levels and different coal pillar widths (Figure 7). The red lines in Figure 7 represent the water flow in the goaf, and the thickness of the line represents the seepage velocity. (1)Figure 7(a) indicates that when the water level is at 18.65 m and the coal pillar width is 10 m, both the goaf and the coal pillar are flowing areas. However, the rock above the coal pillar forms an elastic core area, also called static zone. In addition, when the coal pillar width is 20 m, in the upper part of the elastic core area where , the water in the goaf overflows the elastic core area. When the coal pillar width exceeds 30 m, in most fluid calculation areas where , the water in the goaf is in a static state(2)According to Figure 7(b), when the water level is at 28.65 m and the coal pillar width is 20 m, the coal pillar becomes a static area, and the water in the goaf can only move above the elastic core area of the rock pillar to the adjacent goaf. When the coal pillar width increases to 50 m, the coal/rock pillar and the goaf almost entirely become static areas, and there is only downslope flow on the left side of the elastic core area(3)It can be seen from Figure 7(c) that when the water level is at 38.65 m and the coal pillar width is 20 m, a flowing zone appears in the coal pillar, but it is blocked by the static zone, so continuous seepage is not actually formed. When the coal pillar width increases to 70 m, the static zone cuts off the entire seepage area, and there is almost no water flow in the goaf

(a) Water level at 18.65 m

(b) Water level at 28.65 m

(c) Water level at 38.65 m
The coal pillar width in No. 4-2 was 20 m. During the mining process, water once flowed in from adjacent goafs. When the water level of the adjacent goaf became lower than 18 m, the water influx decreased and gradually stopped, indicating that when the coal pillar width is 20 m, the critical water level is 18 m, which is consistent with the numerical calculation results in Figure 7(a).
4. CFD Numerical Calculation of Seepage Law in the Goafs of Wanli No. 1 Mine
In order to rationally design coal mine underground reservoirs and ensure the safe operation of them, it is necessary to understand clearly the seepage law and seepage path of the water in goaf. So, we conducted CFD numerical calculations for the whole goaf groups of No. 4-2. After assigning values to the permeability of the overlaying rock in the goaf with coupling function of coal/rock permeability and vertical stress obtained from the seepage test, the seepage law and seepage path in the underground reservoir was finally obtained.
4.1. Establishment of Numerical Calculation Model
4.1.1. Geometric Model and Boundary Conditions
First, a three-dimensional surface was generated according to the elevation points from the floor of No. 4-2, and then, the three-dimensional surface was scanned vertically to produce a geometric structure model (Figure 8). This model was 2400 m in width, 4150 m in length, and 115.65 m in thickness, including 14 goafs with 13 coal pillars between them. Next, the elastoplastic calculation module was chosen for the simulation of solid mechanics calculation, in which the top interface of this model was considered a free boundary, the sides around the model as roller bearing boundaries, the bottom interface of the model as a fixed boundary, and the internal boundaries of the model as free boundaries. In terms of fluid calculation, since the water level at the outer boundary of the goaf in goaf No. 42109 was 35 m, the pressure at this inlet boundary was assigned Pa. As lateral boundary of goaf No. 42301 was open to the atmosphere, it was set as pressure outlet with a pressure of 0 Pa. In addition, it was assumed that except for working face No. 42109, the initial water ratio of other areas of the model was 0, and the initial pressure was 0 Pa.

This model can investigate the spatial distribution of seepage law and seepage path of reservoir water (the values for the parameters were shown in Table 2). Assuming that the seepage in the overlaying area of the goaf complies with Darcy’s law, the permeability of overlaying rock can be determined according to the coal/rock permeability and axial stress fitting function in Table 1. The coal pillar width (in Figure 9) between the south and north wings of No. 4-2 is about 210 m, much longer than the critical width of 70 m, which means that the water in the goaf will not flow into the adjacent goaf. That is to say, the south and north goaf of No. 4-2 are separated by coal pillars. However, there are roadways connecting the south and north wings, as shown in Figure 9. The red dots in Figure 9 are waterproof walls with great strength and good airtightness. Therefore, goaf No. 42301 and 42103 in the south wing cannot be connected to goaf No. 42102 in the north wing. The blue dots in Figure 9 are firewalls used as temporary walls with fairly poor airtightness. Therefore, goafs No. 42104, 42105, 42106, 42107, and 42108 in the south wing and No. 42102 in the north wing are connected to each other, and goafs No. 42110, 42109, and 42111 are connected as well.

4.2. Analysis of Simulation Results
4.2.1. The Distribution Law of Water Pressure on the Floor of the Goaf
A cloud diagram of the pressure distribution on the floor of the goaf (Figure 10) shows that the coal pillar width between the south and north wings is very large. The northern floor of the coal pillars is under high water pressure of 194 kPa with a water level at 19.4 m, due to the blocking of boundary coal/rock pillars. Similarly, the floor at the northern boundary of No. 42301 is also under high water pressure of 163 kPa with a water level at 16.3 m. In addition, the ground drilling ZK-1 on No. 42301 working face reveals that the water level here is 15 m, which is consistent with the calculation results of this model, indicating that the model parameters are selected reasonably and the simulation results are accurate and reliable.

4.2.2. Seepage Law in Goaf Groups
In order to obtain the seepage law in goaf groups, the ratio of the in the area to the at the starting point is calculated. If , it means that the fluid is greater than the starting , and the water in the goaf is in a flowing state; if , it indicates that the fluid is less than the starting , and the water in the goaf is in a static state. (1)The value cloud diagram and the flow direction within 0 ~ 5.1 m from the coal seam floor (see Figure 11(a)) reveal that the coal pillars permeability both in the section and at the boundary of the goaf is low, and the fluid is in a static state as . Such phenomenon appears only in the compaction zone in the middle of part of the goaf, while within the goaf, the flow range at the height of 0 ~ 5.1 m is relatively large. As to fluid movement, the water may flow in completely different directions, such as either flowing from the surroundings to the center point or vice versa. The main reason for this is the force of gravity, which means, within 0 ~ 5.1 m from the coal seam floor, gravity plays a dominant role in the direction of fluid movement, and the water normally moves from a high position to a lower one. In addition, due to high permeability of the roadways, the flow velocity in the lanes is quite high(2)Figure 11(b) shows the value cloud map and the flow direction within 15-20 m from the coal seam floor. The result shows that the water level in the goaf begins to exceed the height of the elastic core, and the fluid in it flows from the high-pressure area to the low-pressure area under (3)Figure 11(c) shows the value cloud diagram and the flow direction within 25-30 m from the coal seam floor. The result shows that, compared with that of 15-20 m, the range of water flow is enlarged, which shows that the permeability at the 25-30 m layer is large, and the flow velocity in the goaf is high(4)The seepage path in Wanli No. 1 Mine is: 42109-42110-42112-42111-42101-42102- roadway connecting the north and south wings-“O” ring fissure zone of 42108, 42107, 42106, 42105, 42104, and 42103-42301-42113

(a) The range of 0 ~ 5.1 m from the coal seam floor

(b) The range of 15~20 m from the coal seam floor

(c) The range of 25~30 m from the coal seam floor
5. Conclusions
(1)A coal/rock sample was collected from No. 4-2 of Wanli No. 1 Mine, and tests were designed based on the evolution laws of the lateral supporting pressure on coal pillars. Through tests on a three-axis servo seepage testing machine, the permeability of the coal/rock pillar was obtained as a function of the vertical stress, which provides a basis for studying the seepage properties of the coal/rock pillars between adjacent goafs of the underground reservoir(2)CFD simulation on seepage properties of coal/rock pillar showed that as the water level in the goaf decreased and the coal pillar width increased, the water flow range and flow velocity continued to decrease. At the water level of 18.65, 28.65, and 38.65 m, water could enter adjacent goaf when coal pillar widths were of 30, 50, and 70 m, respectively. In addition, when the coal pillar width was 10 m, the coal pillar was penetrated by high permeability area, and water in the goaf could bypass the elastic core area through coal pillars at the bottom to the adjacent goaf. As the coal pillar width increases, water could no longer enter the adjacent goaf through the coal pillar at the bottom(3)A seepage model for the goaf groups of underground reservoirs was established and numerical calculations performed. The results were verified by the in situ water level data. Numerical calculation results showed that at the height of 0 ~ 5.1 m from the coal seam floor, water in the goaf flowed from a higher place to a low one under gravity. At the height of 15~20 m from the floor, water in the goaf flowed from high-pressure area to low-pressure area under the pressure gradient. At the height of 25~30 m from the floor, water in the goaf flowed at a high velocity in a wide range, as this area’s permeability was high. Based on the seepage area and seepage direction on different layers of the goaf, the seepage path of the underground reservoir in Wanli No. 1 Mine was obtained
Data Availability
All relevant data are available from FigShare, under the DOI: (https://figshare.com/s/e72c1f8e7c8ed0bcf963).
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
We thank X. Chen for assistance with the experiments and Y. Dong for valuable discussion. Financial support for this work was provided by the Independent Research Project of State Key Laboratory of Coal Resources and Safe Mining (No. SKLCRSM19X009) and the National Natural Science Foundation of China (No. 52174212).