Groundwater flow simulation and its application in GaoSong ore field, China

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.

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.

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.

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.

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).

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.

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 km², 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).

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.

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.

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).

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).

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.

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 (K_x/K_y).

• (1)

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

• (2)

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

• (3)

  Simulation Model C: the model grid is rotated 0° about Z_axis (Figure 8), and the horizontal permeability coefficient anisotropy (K_x/K_y) 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.

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°.

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 (K_x/K_y) 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.