Numerical Modeling of Groundwater Flow in Karst Aquifer, Makeng Mining Area

Problem statement: Making iron ore is one of the largest mining proje cts in south east of china. Due to the development of mining activities in that site, it has become necessary to increase t he depth of exploration. Increase the exploration dept h makes the mining tunnels subjected to the karst water inrush. Approach: A hydrological and a hydrogeological model for the Makeng area have been developed, which yield information on relevant para meters such as groundwater recharge and margins’ lateral inflow, to estimate aquifer yield. USGS flo w code, MODFLOW 2000, was used to produce the numerical model. Collected GIS based information wa s synthesized in a finite difference numerical model. The regional steady and transient-state flow was calibrated under pre-development conditions assuming an equivalent porous medium approach. Results: Water budget calculations show that the total groundwater flow in regional aquifers amounts to 2.98 mm year. Infiltration from precipitation provides 61.7% of the groundwater supply, while 21% comes from lateral inflow and the remaining 17.3% is induced recharge from surface waters. Disc harge from regional aquifers occurs through springs outflow 88.5% and flow to streams 11.5%. Conclusion/Recommendations: Although the karstic nature of the limestone aquifer the equival ent porous-medium flow model is appropriate to represent hydraulic heads and recharge/discharge re lationships on a regional scale. The results of thi s study can be used to predict the required amounts o f pumping and the possible locations to dewater the groundwater around the mining tunnels.


INTRODUCTION
China is rich in mineral resources. Up to now, geologists have confirmed reserves more than 158 kinds of world's known minerals can be found in china, putting China as the third in the world in total reserves and the second largest mineral producer (Rani and Chen, 2009). Because of the current large-scale development of mining, it becomes necessary to increase the depth of exploration to extract the minerals. Increase the exploration depth makes the mining tunnels subjected to the karst water inrush (Wu et al., 2004). Therefore, the clear identification of hydrogeological conditions of karst aquifers is required as an important task in dewatering design.
Makeng iron ore is one of the largest mining projects in china, has three mining levels +420, 300 and +200 m above the sea level. The average groundwater level around the mining area is about +450 m. at present, abundance karst water seeps into the first level mining tunnels. The groundwater seeps through the tunnels connections with the faults and the karst caves.
At present, the mining activities reduced gradually affected with that water inrush. To continue the exploration a dewatering system should be implemented, taking in consideration the surface-water interaction. The dewatering process will depend on pumping and a drainage system around the mining tunnels to reduce the groundwater level for about 300 m. That required a complex study for all the hydrogeological conditions. Then applying a numerical simulation to determine the groundwater flow directions, forecasting the amount of inflow to the drainage system and calculate the amounts and locations for the required pumping.
Application of numerical models in karst aquifers is more problematic, because of the high heterogeneity (Singhal and Gupta, 1999;Bridget et al., 2003;Nico and David, 2007). The first stage of this study is the description of the hydrogeological framework. The study proceeded with the development of the conceptual model of regional groundwater flow. MODFLOW (McDonald and Harbaugh, 1988) was used to simulate the 3-dimensional regional groundwater flow. The calibration of the model parameters was conducted under steady and transientstate flow conditions. The numerical model was then used to estimate the regional water budget. This manuscript is the second study describing the hydrogeological components of Makeng mining project. The first study presents the estimation of aquifer recharge from springs' flow records (Rani and Chen, 2009).

Hydrogeological characteristics of the study area:
The making iron-mining project is located in mountain region at the southern part of Longyan City, west of Fujian province. The landscape is dominated by isolated mountains with altitudes range from +350 m to +1170 above the sea level. The mining area extended from the southern part of the study area, through a transition zone with an elevation averages from +800 at the southeastern part, to +420 m in the eastern part beside Xima River. The annual maximum rainfall is 2348 mm; the minimum rainfall is 1188 mm, the daily maximum rainfall 322 mm, rainfall concentrated in summer and autumn. The annual average temperature is 20°C. The potential evaporation also varies between 1930.6-1166.7 mm year −1 . The study area is igneous rocks groundwater basin with a total thickness of around 150-700 m, comprising Permian sandstone and a Paleozoic limestone system, separated by mudstone aquitard. The aquifers system is bounded with a faults zone, these faults are assumed as a model boundary. Other faults inside the study area are vertical and acts as flow barriers. The flow system of the study area represents a sequence of two aquifers and one aquitard. All three aquifers contribute to the regional groundwater flow system extending from the south eastern mountains toward the Qilai spring at the north western part. The upper unconfined sandstone aquifer is 0-400 m thick and is composed of sands and gravels lenses. The second, aquifer, which is 0-800 m thick, with common thickness in the range of 300 m, its lithology is mainly siliceous limestone, shale limestone and pure limestone and dolomitic limestone. The limestone aquifer outcrop the surface with a weathering zone in the south eastern part of the study area and plays an important role in receiving and passing the downward precipitation recharge over the entire area. It is separated from the upper unconfined sandstone aquifer by a 0-50 m thick mudstone aquitard.

Conceptual model:
A total of about 346 borehole data points, are available from previous investigation for making area, allowed full use of the advantage of GMS to build the conceptual model within its software environment. Analysis of the data followed by semiautomatic preparation of cross-sections with 3 hydrostratigraphic units consisting of 2 aquifers and 1 aquitard. The maps of the top and bottom of each of the hydrostratigraphic units were created by kriging of borehole data and outcrop boundaries. The three potentiometric surfaces of the three aquifers were constructed in the same way and afterwards imported into the model as initial heads. It can be seen that the regional groundwater flow is from the SW boundary, across the west fault zone towards the NE outflow boundary. The head pattern of the upper aquifer is influenced by the surface topography and frequently receiving rainfall recharge. The hydrogeological setting of the study area is characterized by large horizontal and vertical heterogeneity. The recharge conditions of the flow system were characterized by 3 main external sources of water: (1) lateral inflow from the SW boundary fault; (2) rain recharge; (3) rivers bed infiltration in some sections. A substantial contribution of rain recharge and lateral inflow from the boundary fault along the SW margin of the study area to the overall water budget of the study area was documented by (Rani and Chen, 2009).

Numerical model:
The objectives of the numerical model were to evaluate and test the validity of various interpretations about flow system. These include: (1) the location and type of flow-system boundaries, (2) the definition of recharge areas and (3) variations in interpretation of hydrogeological frame work. The numerical groundwater flow model MODFLOW 2000 was used. The aquifer system domain was horizontally divided into 80 rows and 80 columns resulting in a uniform size of cells of 83×83 m. The model consists of 4 layers simulating the principal hydrostratigraphic units' presents the two regional aquifers and the aquitard. The top layer (layer #1) represents at the same time quartz-sandstone aquifer, the mudstone and limestone when present. The following layer (Layer #2) the mudstone aquifer and the highly fractured rock from the limestone aquifer which present. The last two layers, (Layer #3) continue the highly fractured rock and (Layer #4) the medium fracture rocks of the limestone aquifer. Surface and bottom elevations for all layers are interpolated from the boreholes data.
The model boundaries were selected to correspond as much as possible with the natural geomorphological and hydrological features (Fig. 1). A specified flux was assigned along the SW boundary simulating the groundwater inflow from the adjacent sandstone aquifer. This flux was constrained with the estimated potentiometric surface along the limit. Elsewhere, lateral boundaries coincident with topographic divides are assumed at the same time as groundwater divides and are simulated as no-flow boundaries. The main Faults inside the study area were simulated by the barrier package for all the four aquifers.
The river package simulates the effects of flow between surface water features and groundwater systems. Water levels in the rivers and streams' stages were approximated from the monitoring gauges and the available topographic map. The drain package of MODFLOW was used to simulate Qilai spring at NW corner of the model area.
The recharge rates were assigned as in (Rani and Chen, 2009) for the upper aquifer and the pinched part of the limestone aquifer using the MODFLOW recharge package. The initial hydraulic conductivity distribution for the limestone aquifer (layer # 3 and 4) were assigned as shown in (Fig. 2) (Rani and Chen, 2009). For the upper unconfined quartz-sandstone aquifer there was a lack of hydraulic head and pumping test data therefore heads and K were assumed on the basis of lithological knowledge and fragmented information from shallow wells.

Steady-state simulation:
The steady state simulation based on the predevelopment data in the study area. Investigations of the Making area started as early as 1977.   [1977][1978][1979][1980], was no pumping stresses and the water level data are extensive during the period. From the water level measurements in the study area, it appears the pre-development did not impose significant changes to the groundwater flow system, the aquifer system is assumed to have been in state of dynamic equilibrium. Because no data were available for the north area, a few "fictitious" wells were introduces along the profile to ensure that a realistic interpolation was obtained from the contouring software. Water levels at these "fictitious" wells were set proportional to the elevation of Qilai spring.

Calibration of the unstressed steady-state groundwater flow model:
The observation data used in the calibration process was the water table elevations from observation wells. Calibration is achieved when the error is within the estimated error interval (±0.15 m) of the observed value. For the initial simulation, the values of hydraulic conductivity assigned were as shown in (Fig. 2). Initial rainfall recharge amounts have been assigned based on the values documented by (Rani and Chen, 2009). These initial values were adjusted during calibration where the degree of fit between model simulations and field measurements was quantified by statistical means. For this model, the estimated RMS of error was 1.7 which is very small relative to the total head loss (1.8 m). Once these criteria were satisfied, the model was considered calibrated. It is important to note that these criteria of calibration were met with the residual mean being close to zero and the ratio of residual standard deviation to the overall range in head is adequately less than 10%.

The transient-state simulation:
The transient simulation was run from the initial un-pumped steady state to examine the effect of the pumping well discharging 21 L/S from layer 3 and 4 from December 5, 1977 (08:00 AM) until December 10, 1977 (10:00 AM). The discharging of water through the pumping well was simulated as specified flow boundary using the specified pumping rate. The pumping well is a sink and is represented in the model by a node. The stress periods selected to match the dates/times whenever input data in the map module changed as well as whenever groundwater flow collected. Five stress periods were simulated with 1 time step of simulation assigned to each stress period. Results of transient flow simulation were presented in the form of head contour map at the end of each time step. There are 11 observation wells in addition to the pumping well trapping the deep aquifer units (model layer 3 and 4).

Calibrating the transient-state groundwater flow simulation:
For each transient model run, an analysis of the observed versus computed water levels was conducted to determine the accuracy of the simulation. Three methods were used for determining the level of accuracy. The first method involved calculating the mean of the residuals. This provides a measure of the bias of the distribution, indicating whether the simulation was over or under-estimating the water table as a whole. The second method involved the calculation of the root-mean-square or standard deviation of the residuals, which provides a measure of the squared differences in measured and computed water levels, a sensitive test of the range of differences. The third method involved calculating the standard error between the observed and computed values. The standard error calculates the mean of the absolute values of the residuals. This provides a more realistic measure of the average difference between the observed and computed values. The calibration procedure used with the transient simulation resembled the one used with calibrating the steady-state simulation described above. The degree of fit between model simulations and field measurements was evaluated based on the principal calibration statistics: The acceptable residual should be a small fraction of the difference between the highest and lowest heads across the site and be based on: (1) the ratio of the Root Mean Squared (RMS) of error to the total head loss should be small; (2) head differential of <5% for the residual mean and standard deviation; and <10% for the ratio of the standard deviation to total head change. Once these criteria were satisfied, the model was considered calibrated. For this model, the estimated RMS of error was 2.1 m which is very small relative to the total transient head loss (total drawdown) which was 4.61 m Sensitivity analysis: A sensitivity analysis was performed for the following calibrated parameters: specified flux at the SE boundary, recharge flow rates and hydraulic conductivities. The applied flux at the SE was varied between 500 m 3 day −1 less and higher the calibrated value. The performance of the model in terms of the calibration targets was stable up to 1000 m 3 day −1 . For values greater than 2000 m 3 day −1 , the groundwater head increased in such a way that exceeds the surface elevations, which is unrealistic.
For the recharge flow rates, an arbitrary range of variation was selected (R ± 0.5R). NRMS, RMS and ARM were rather sensitive to changes in recharge rates. Each of the 29 calibrated hydraulic conductivity (K) zones was multiplied by the following factors: 0.01, 0.1, 0.5, 5, 10 and 100. Model results were particularly sensitive to changes in the hydraulic conductivity in the middle area. In terms of flow components, a variation from 0.5 up to 10 K resulted in ground water inflows through the SE boundary within the defined range (1000-1560 m 3 day −1 ).

RESULTS AND DISCUSSION
The model was considered a calibrated base model by adjusting the hydraulic parameters of the limestone aquifer until the computed water levels were reasonably well matched to observed potentiometric surface elevation contours. Flow out of the model area by rivers and Qilai spring was also adjusted until it was reasonably well matched to observed data.
Hydraulic connection between groundwater of the first, quartz-sandstone aquifer and surface water of the rivers was simulated as head dependent boundary condition using the MODFLOW River Package. A good hydraulic connection between rivers and the limestone aquifer was modeled by enlargement of vertical leakance, the segments which pass through the pinched part of the limestone was also enlarged.
The calibrated values for hydraulic conductivities (29 zones) vary between 2×10-4 and 6.7 m 3 day −1 . Whereas the calibrated values for the storage coefficient vary from 5×10-3 up to 3×10-1. The observed Vs simulated groundwater heads for the steady-state and transient-state models are shown in (Fig. 3). For the case of the steady-state model (Fig. 3a), NRMS, ARM and RMS values show an acceptable calibration when compared with the previously defined calibration targets. In the case of the transient model, (Fig. 3b) shows the results for the complete period of the pumping test 5 days. From (Fig. 3b) it is possible to see that the calibration targets are also met, but with slightly poorer quality than the steady case. Most of the residuals are contained in the range between −2 and +2 m for a large part of the modeled area. For the northern area, which corresponds to the thinnest part of the aquifer, the range of the error terms correspond to -1 to +1 m (2 m) and this represents 1% of the thickness of the aquifer. For the rest of the modeled domain, the range of the error terms is less than 2% of the average aquifer thickness (thickness 250-300 m).
Simulated groundwater heads clearly overestimate the observed values in the middle of the domain where head residuals are on the order of 2-4 m. In the west of the domain, observed values are underestimated on the order of 1-2 m. In the north of the model area, the residuals are on the order of 1 m.

CONCLUSION
A regional groundwater flow modeling study was conducted in Making mining area, Fujian province, china. This study represents a significant advance in the understanding of the groundwater flow in the region using equivalent porous medium approach, particularly when considering that no quantitative analysis has been performed to date. That model represents the first step towards design a dewatering system to continue mining activities in the region.
The general finding was that the limestone aquifer received lateral influx from the outside sandstone rockfracturing pattern and the interconnection along the south east margin. Finally, the numerical model was calibrated against the measured potentiometric surface under the assumed steady and transient-state conditions. The quantitative estimates of the dynamic groundwater parameters show that the total flow in the region amounts to 2.98 Mm 3 year −1 . The regional areal recharge is 1.67 Mm 3 year −1 . On the other hand, 0.34 Mm 3 year −1 are discharged naturally to surface waters and 0.57 Mm 3 year −1 are recharged to the limestone aquifer as lateral inflow.