Mine water inrush is one of the important factors threatening safe production in mines. The accurate understanding of the mine groundwater flow field can effectively reduce the hazards of mine water inrush. Numerical simulation is an important method to study the groundwater flow field. This paper numerically simulates the groundwater seepage field in the GaoSong ore field. In order to ensure the accuracy of the numerical model, the research team completed 3,724 field fissure measurements in the study area. The fracture measurement results were analyzed using the GEOFRAC method and the whole-area fracture network data were generated. On this basis, the rock mass permeability coefficient tensor of the aquifer in the study area was calculated. The tensor calculation results are used in the numerical model of groundwater flow. After calculation, the obtained numerical model can better represent the groundwater seepage field in the study area. In addition, we designed three different numerical models for calculation, mainly to explore the influence of the tensor assignment of permeability coefficient on the calculation results of water yield of the mine. The results showed that irrational fathom tensor assignment would cause a significant deviation in calculation results.

Prediction of water inrush in mines not only involves complex geological problems, but must also be combined with mine production (He & Yang 2008). Groundwater numerical simulation technology is widely used in the evaluation and research of groundwater resources. Under the promotion of the comprehensive management concept of water resources system, the coupling with surface water simulation is becoming more increasingly close (Labolle et al. 2003; Liu & Mou 2016) and provides an effective basis for the unified management of water resources. According to the conceptual model of mine hydrogeology, numerical simulation technology is used to establish the numerical model of groundwater seepage field in the study area, and the calculation and prediction of water yield of the mine can effectively prevent mine flooding accidents.

The rock mass medium in the natural state is discontinuous, and the groundwater is mainly seepage in the rock mass fracture network. In cases where dimensions of the works or scale of application are large compared to the heterogeneities of the medium, flow properties are computed with the aid of an homogeneous model of permeability. The notion of effective permeability is used for a medium that is statistically homogeneous on a large scale beyond a minimum volume called representative elementary volume (REV) (Pouya & Fouche 2009). The permeability coefficient tensor is a second-order symmetric tensor. The permeability coefficient tensor is only the equivalent of the seepage flow. The rock is equivalent to a continuous osmotic medium, which can be borrowed from the analysis method of soil seepage, which is the main reason why the equivalent continuous medium is widely used in rock mass seepage (Zhang 2008; Wu et al. 2014). How to use effective penetration tensors in groundwater numerical models is also a problem to be discussed. Use of error will cause significant deviation calculation results.

The research object is the GaoSong ore field in Gejiu City, Yunnan Province, China. It is located in the south-central part of the Yunnan-Guizhou Plateau (Figures 1 and 2). The study area has the following specific hydrogeological characteristics. (1) Its groundwater type is dominated by karst water and bedrock fissure water. Due to the obvious directionality of the distribution of bedrock karst and fissures in space, the permeability coefficient of the fracture medium exhibits strong anisotropy. (2) The mine tunnel drains the groundwater, forms a groundwater depression cone, and lacks groundwater level observation data. When verifying the numerical model, only the observation data of the 1,360 mid-water yield can be used. (3) There is no aquifuge in the study area, and the groundwater is all generalized to the unconfined aquifer.

Figure 1

Regional hydrogeological map. 1. Strongly water-bearing karst aquifer 2. Medium water-bearing karst aquifer 3. Medium water-bearing porous aquifer 4. Weakly water-bearing igneous fissure water aquifer 5. Weakly water-bearing clastic rock fissure water aquifer 6. Strongly water-bearing clastic rock fissure water aquifer 7. Underground river 8. Direction of groundwater flow 9. Groundwater divide 10. Fault 11. Stratum boundary 12. River 13. Lake 14. City 15. Water-rich block.

Figure 1

Regional hydrogeological map. 1. Strongly water-bearing karst aquifer 2. Medium water-bearing karst aquifer 3. Medium water-bearing porous aquifer 4. Weakly water-bearing igneous fissure water aquifer 5. Weakly water-bearing clastic rock fissure water aquifer 6. Strongly water-bearing clastic rock fissure water aquifer 7. Underground river 8. Direction of groundwater flow 9. Groundwater divide 10. Fault 11. Stratum boundary 12. River 13. Lake 14. City 15. Water-rich block.

Close modal
Figure 2

Terrain of the research area.

Figure 2

Terrain of the research area.

Close modal

In this paper, a reasonable numerical model of groundwater seepage field will be established according to the specific conditions of the mine. The use of the model focuses on the following aspects: the assignment of the groundwater permeability coefficient tensor, the generalization of the fracture as a hydrogeological boundary, and the generalization of complex tunnels.

Hydrogeological structure

The hydrogeological unit strata are dominated by Triassic carbonate strata, and the geological structures are mainly northeastward and east-west (Li 2011). The groundwater in the study area is mainly discharged to the Honghe River by the underground dark river (Figure 1).

According to the regional hydrogeological survey data, the domain of the study area is determined by the four faults of the Gesong fault, the Beiyingshan fault, the Gejiu fault, and the Jiajieshan fault. The aquifer in the study area is an unconfined aquifer, and the aquifer medium is carbonate rock with a total thickness greater than 1,000 m. The carbonate rock is exposed directly on the surface and there is no aquifuge in the upper part. The lower part of the carbonate rock is granite, which is an aquitard in the study area (Figure 3). The shape of the granite roof has a significant impact on the groundwater seepage field.

Figure 3

The granite top tray in the research area.

Figure 3

The granite top tray in the research area.

Close modal

Groundwater recharge and discharge

Groundwater recharge sources in the study area include rainfall and underground runoff. The karst depression, funnel, falling water cave, and dissolution fissures in the study area are mostly developed, and the rainfall infiltration coefficient is 65–85%. The atmospheric rainfall vertically supplies groundwater, and the seasonal variation of the groundwater level varies greatly, which can be greater than 150 m. The impact of the groundwater level in the central part of the study area was greatly affected by the mine draining. Under normal circumstances, groundwater in the study area is all collected in the 1,360 tunnels. The surrounding groundwater is supplied to the groundwater aquifer in the study area through the water-conducting fault (Kang & Yang 2017). When the groundwater level is high during the rainy season, the aquifer in the study area can be drained to the lower terrain.

Boundary conditions

The boundaries of the study area are the Gesong fault, the Beiyingshan fault, the Gejiu fault, and the Jiajieshan fault. According to regional hydrogeological data, the Gesong fault in the north of the study area and the Beiyingshan fault in the south are large-scale water-conducting faults, but the water level in the fault zone varies greatly and is not suitable as a fixed water head boundary. Therefore, the Gesong fault and Beiyingshan fault are the boundaries of fixed flow in the study area. The water-carrying section of the two boundaries can be estimated based on the shape of the underlying granite. In addition, the western section of the Gesong fault is blocked by an intrusive rock mass and partially serves as the confining boundary.

The east side of the Gejiu fault is carbonate rock, and the west side is a sandstone and a shale aquifuge. Therefore, the Gejiu fault is used as a confining boundary (Figure 4).

Figure 4

Cross section of Gejiu fault. 1. Clay 2. Silty mudstone and shale 3. Limestone 4. Karst cave 5. Bore and water level 6. Gejiu fault.

Figure 4

Cross section of Gejiu fault. 1. Clay 2. Silty mudstone and shale 3. Limestone 4. Karst cave 5. Bore and water level 6. Gejiu fault.

Close modal

There is a large area of granite bulge near the Jiajieshan fault in the eastern part of the mining area, and the granite is an aquitard. In addition, the elevation of the surface water body in the Mengzi Basin in the eastern part of the mining area is lower than the 1,360 m, so the Jiajieshan fault in the eastern part of the study area serves as the confining boundary.

Aquifer property

The permeability coefficient of the fracture medium exhibits strong anisotropy and the vertical karst development in the study area is strong. It is obviously more reasonable to use the permeability coefficient tensor to indicate the permeability of the aquifer in the study area.

In order to obtain the permeability tensor of the aquifer in the study area, the research team carried out a large number of structural fracture measurements. The fracture measurement items include fracture's dip, the angle of dip, the trace length, the fracture development density, and the aperture. In order to ensure the uniformity of data measurement, the study area is divided into a grid of 200 m × 200 m. At least one set of valid data within each grid is guaranteed during the fracture measurement. The total area measured is about 133 km2, and the data volume covers 80% of the study area. A total of 3,724 surface fissure measurement data were obtained, and the research team performed statistical analysis on these measurements indoors.

According to the statistical law of fractures, using the GEOFRAC (Liu & Koike 2007; Liu et al. 2013, 2018) method, the research team simulated the discrete fracture network in the study area (Figure 5).

Figure 5

All the fractures simulated by GEOFRAC.

Figure 5

All the fractures simulated by GEOFRAC.

Close modal
The two-dimensional distributions of the simulated borehole fractures and the linear features of lineaments and microcracks can be linked by permeability tensor analysis to estimate permeability of geological media and its change with scales. Following Oda's equation, the permeability tensor, kij of degree three is expressed by
formula
where δij is a delta function, λ a constant expressing the continuity of fractures.
formula
where V is the volume of the study region, t the hydraulic aperture of the fracture, n a normal vector of the fracture plane. For two-dimensional analysis, V and r are replaced by the area and trace length of the fracture, and the effect is only considered for the n. Assuming that the fracture shape is a plate, the λ was defined as 1/12 (Oda et al. 1987).

Using the simulated discrete fracture network, the research team calculated the permeability coefficient tensor of the aquifer rock mass in the study area (Oda et al. 1987; Yang et al. 2013). During the calculation, the REV (representative elementary volume) of the discrete fracture network of the rock mass in the study area was analyzed, indicating that REV is present and the value is 36 m. Finally, the permeability tensor of the discrete fracture network of the aquifer in the GaoSong ore field is , and the main infiltration direction is 135°.

The karst development in the vertical direction of the study area shows obvious zoning characteristics. The surface karst is strongly developed in the area, within 300 m from the surface the karst rate is >1.5%, and the karst rate is >1% within 500 m. The degree of karst development gradually decreases with increasing depth (Li 2011). According to the law of vertical karst development, we divided the aquifer in the study area into two layers vertically.

Mathematical model

There is an REV value in the fracture network of the rock mass in the study area, so we can use the rock mass permeability coefficient tensor to equalize the fracture medium of the rock mass to an anisotropic continuous medium. In this medium, the law of groundwater movement is in accordance with Darcy's law. Therefore, the mathematical model of the groundwater seepage field in the study area can use the ordinary three-dimensional unsteady flow groundwater model.
formula
where ω is groundwater recharge or discharge; Ss is specific storage of the aquifer; is initial water level; Ω is groundwater runoff area; q1(x,y,z,t) unsteady-flow boundary flow; Kxx, Kyy, Kzz represent the permeability coefficient in the x,y,z direction; n is normal direction outside the boundary of the runoff area. Specific to the model, boundary of fixed flow: AB, EF; confining boundary: AF, BC, CD, DE (Figure 6(a)).
Figure 6

Conceptual model coverage.

Figure 6

Conceptual model coverage.

Close modal

Establishment and calibration of groundwater flow model

The geostrome of the study area are generalized into three aquifers. The lithology of the first and second aquifers is carbonate rock, which has different permeability; the third aquifer is granite, which acts as an aquitard. The first aquifer is affected by topography and karst development, and its horizontal permeability is anisotropic, divided into three types of zones (Figure 6(a), Table 1); and the recharge coefficient of the first aquifer is also inconsistent. Therefore, the precipitation recharge in the study area is divided into two zones (Figure 6(c)), and the average monthly infiltration in each zone is calculated based on the precipitation (Figure 7). Other aquifer property parameters are shown in Table 1.

Table 1

Material properties of the aquifer system

PropertyAquifer1_1Aquifer1_2Aquifer1_3Aquifer2Aquifer3
Horizontal K(m/d) 0.64 0.45 0.40 0.38 0.038 
Horizontal anis. (Kx/Ky) 1.49 1.49 1.49 1.49 1.00 
Vertical anis. (Kx/Kz) 0.67 1.00 1.25 1.25 1.00 
Specific yield 0.07 0.06 0.047 0.047 0.002 
Porosity 0.07 0.06 0.047 0.047 0.002 
PropertyAquifer1_1Aquifer1_2Aquifer1_3Aquifer2Aquifer3
Horizontal K(m/d) 0.64 0.45 0.40 0.38 0.038 
Horizontal anis. (Kx/Ky) 1.49 1.49 1.49 1.49 1.00 
Vertical anis. (Kx/Kz) 0.67 1.00 1.25 1.25 1.00 
Specific yield 0.07 0.06 0.047 0.047 0.002 
Porosity 0.07 0.06 0.047 0.047 0.002 
Figure 7

Observed vs. simulated water yield of mine and precipitation.

Figure 7

Observed vs. simulated water yield of mine and precipitation.

Close modal

The mine tunnels are mainly horizontal with an elevation of 1,360 m. The number of tunnels is large and the structure is complicated. According to the size of the tunnel and the waterproof measures, the conductivity of the tunnel is different. We divide the tunels in the model into three types and we assign different parameter values to the three types of tunnels (Figure 6(d) and Table 2).

Table 2

The conductance of tunnels

ClassesConductance ((m2/d)/m)
Exits of the main tunnels 5.16 
Main tunnels 25.80 
Other tunnels 10.32 
ClassesConductance ((m2/d)/m)
Exits of the main tunnels 5.16 
Main tunnels 25.80 
Other tunnels 10.32 

The permeability tensor of the aquifer must be treated when the model is meshed. According to the main direction value of the permeability coefficient tensor, we rotate the mesh when we build the model. The specific way of rotation is that the mesh is rotated 135° around the Z axis. 135° is the main direction of groundwater infiltration (Figure 8).

Figure 8

Grid of model. Model A grid rotation 135°, Model B and C grid rotation 0°.

Figure 8

Grid of model. Model A grid rotation 135°, Model B and C grid rotation 0°.

Close modal

The precipitation and water yield of mine data in the numerical model use actual data. The data were recorded from January 1, 2014 to December 31, 2014. The precipitation is from the data of the Gejiu City Meteorological Bureau. The observation point of the water yield of the mine is located at the exit of the 1,360 tunnel, and its water inflow represents the entire mine displacement (Figure 7).

The model mesh size is 100 m × 100 m, and the number of cells is 21,885; simulation package uses MODFLOW 2005 version with GMS. Model type is unsteady flow. The simulation model stresses periods from January 1, 2014 to December 31, 2014. The numbers of time steps amounts to 60.

The groundwater level of the second aquifer is significantly affected by the drainage effect of the tuunels (Figure 9). Centered on the tunnels, the groundwater aquifer formed a groundwater depression cone with a maximum depth of 150 m. Groundwater level depletion also occurred in the middle of the second aquifer. We believe that this phenomenon is caused by granite protrusions.

Figure 9

Simulation results of aquifer 2.

Figure 9

Simulation results of aquifer 2.

Close modal

Comparing the simulated and the observed water yield, we find that the two curves are basically the same (Figure 7). This phenomenon shows that the numerical model of groundwater we have established is basically in line with the actual situation. However, there are still large deviations. The main reasons for the analysis are as follows:

  • (1)

    The fractures of various scales in the study area develop strongly, so the permeability coefficient tensor of the aquifer in the study area will not have only one value; we only carry out the calculation of the permeability coefficient tensor of the two-dimensional plane. Actually, it is a three-dimensional infiltration tensor value.

  • (2)

    The study area is a karst aquifer with karst pipes and caves. These karst pipes and caves cause strong inhomogeneities in groundwater infiltration.

  • (3)

    The fixed flow boundary is an estimated value and it is not very accurate.

This is the first time we have carried out numerical simulation of groundwater in the study area based on the permeability coefficient tensor. The main purpose of this simulation work is to find out the various problems in the numerical simulation of groundwater. After a thorough analysis of this problem, it will assist in more accurate simulation calculations. Therefore, we believe that this calculation of deviation is an acceptable range. To this end, the impact of other aspects will continue to be discussed below.

In order to observe the influence of the permeability coefficient tensor on the water yield of the mine in the numerical model, we designed three grid models for analysis and comparison from the perspective of grid rotation angle and horizontal permeability coefficient anisotropy (Kx/Ky).

  • (1)

    Simulation Model A: the model grid is rotated 135° about Z_axis (Figure 8), and the horizontal permeability coefficient anisotropy (Kx/Ky) is 1.49.

  • (2)

    Simulation Model B: the model grid is rotated 0° about Z_axis (Figure 8), and the horizontal permeability coefficient anisotropy (Kx/Ky) is 1.0.

  • (3)

    Simulation Model C: the model grid is rotated 0° about Z_axis (Figure 8), and the horizontal permeability coefficient anisotropy (Kx/Ky) is 1.49.

The three models use the unsteady flow calculation mode and then compare the calculation results.

It can be obtained from the comparison curve of the calculation results that the grid rotation angle has the most obvious influence on the calculation results (Figure 10, Model A vs. Model C). The permeability coefficient anisotropy also has an effect on the calculation results, but the effect is smaller.

Figure 10

Simulated water yield of mine.

Figure 10

Simulated water yield of mine.

Close modal

Further analysis of the penetration ellipse of the three models can better illustrate this calculation (Pouya & Fouche 2009; Wu et al. 2014) (Figure 11). The north and south boundaries of the numerical model are boundaries of fixed flow, and some of the tunnels are roughly parallel to the boundaries. When the main infiltration direction of the aquifer is consistent with the direction of the hydraulic gradient, according to Darcy's law, water yield of themine is the largest. Among the three models, only the main infiltration direction of Model A has an angle of 135° with the direction of the hydraulic gradient, so water yield of the mine of Model A is the smallest. The horizontal permeability coefficient anisotropy (Kx/Ky) of Model B is 1.0, and all directions are the main infiltration directions, so water yield of the mine of Model B is the largest. In addition, we compare the values of water yield of the mine in each period of Model A and Model C. The ratio between the two ranges from 0.568 to 0.767. The value of cos135° is 0.707 after ignoring the negative sign. Obviously, the calculation results of Model A and Model B are different, the main reason being that the grid is rotated by 135°.

Figure 11

Polar diagram of penetration ellipse.

Figure 11

Polar diagram of penetration ellipse.

Close modal

We can assume the following understanding: When using GMS for numerical simulation of groundwater seepage field, if you want to consider the influence of the permeability coefficient tensor value of the rock mass, make the grid rotation direction consistent with the main direction of the permeability tensor, and the simulation result will be closer to real groundwater seepage. Only modifying the horizontal permeability coefficient anisotropy (Kx/Ky) does not reflect the influence of the permeability coefficient tensor of the aquifer.

The effect of the main directionality of the plane permeability coefficient tensor can be reflected by the rotation of the model grid when using GMS. The assignment method of the permeability coefficient tensor in the numerical model of groundwater is very critical. When building the model, not only the non-uniformity of the permeability coefficient is considered, but also the relationship between the main direction and the mesh needs to be observed. Tunnels' generalization as close as possible to the actual situation is equally important for the accuracy of the water yield calculation.

We will improve the work in the following aspects: (1) the permeability coefficient tensor of the fracture network should consider the karst pipeline; (2) the permeability coefficient tensor of the study area cannot have only one value, but should be assigned by partition; and (3) rock permeability test and pumping are required tests, etc., used to correct the calculation results of the rock mass permeability coefficient tensor.

This work was supported by the National Natural Science Foundation of China Projects (NSFC) numbers 40902058 and 41562017.

He
X.-W.
&
Yang
J.
2008
Problem and countermeasure of mine water resource regeneration in China
.
Journal of China Coal Society
33
(
1
),
63
66
.
Kang
D.
&
Yang
Z.
2017
The hydrogeological investigation in GaoSong field of Gejiu tin mine
.
China Mining Magazine
26
(
S1
),
191
194
.
Labolle
E.-M.
,
Ahmed
A. A.
&
Fogg
G. E.
2003
Review of the integrated groundwater and surface-water model (IGSM)
.
Ground Water
41
(
2
),
238
246
.
Li
R.-W.
2011
Characteristics and control factors of groundwater of Gejiu orefield
.
Yunnan Geology
30
(
1
),
64
66
.
Liu
C.
,
Ni
C.
,
Yan
Y.
&
Tan
L.
2013
Three-dimensional simulation of rock fractures by geostatistical method: a case study of Gaosong field in Yunnan province
.
Acta Geoscientica Sinica
34
(
3
),
347
353
.
Liu
C.
,
Kubo
T.
,
Lu
L.
,
Koike
K.
&
Zhu
W
, .
2018
Spatial simulation and characterization of three-dimensional fractures in Gejiu tin district, southwest China, using GEOFRAC
.
Natural Resources Research.
https://doi.org/10.1007/s11053-018-9381-8.
Oda
M.
,
Hatsuyama
Y.
&
Ohnishi
Y.
1987
Numerical experiments on permeability tensor and its application to jointed granite at Stripa mine, Sweden
.
Journal of Geophysical Research
92
(
B8
),
8037
8048
.
Wu
Y.
,
Liu
Q.
,
Li
J.
&
He
J.
2014
Effect of correlation between joint trace and aperture on hydro-mechanical behavior of fractured rock mass
.
Chinese Journal of Rock Mechanics and Engineering
33
(
S2
),
3555
3562
.
Yang
J.-P.
,
Chen
W.-Z.
,
Wu
Y.-X.
&
Tan
X.-J.
2013
Numerical study on equivalent permeability tensor of fractured rock masses
.
Chinese Journal of Geotechnical Engineering
35
(
6
),
1183
1188
.
Zhang
Y.
2008
Rock Hydraulics and Engineering
.
China Water & Power Press
,
pp
.
112
124
.