Estimating dewatering in an underground mine by using a 3D finite element model

Groundwater inflow to an underground mine will seriously affect its mining plan and engineering geology safety. Groundwater models are powerful tools commonly used in the mines to develop dewatering strategies. Many mines in the Kolwezi area have been present since the 1950s, and groundwater flow patterns have been significantly influenced by mining activities. A mining plan is developed for an underground mine with overturned syncline strata in Kolwezi, Congo. Previous groundwater models using layered homogeneous media lowered model accuracies. A new three-dimensional groundwater model using FEFLOW, consisting of a combined regionally and locally geology models integrating 16 hydrogeological cross-sections and borehole logging data, are formulated to predict the underground dewatering in the study area. A 31-days pumping tests with 3 pumping wells and 28 observation wells are carried out to estimate the hydrogeological properties. The simulated water level data match the observed data rather well. Under 8 scenarios of possible well designs, the model predicts a possible dewatering capacity greater 23,900 m3/d at the initial stage of mining. The concept of the model and its application can be a reference for other mines with complex geology for mining safety in the region of interest.


Introduction
Mineral exploitation and mining activities play a vital role in the economic development of the Democratic Republic of Congo (DRC). According to the USGS 2015 Minerals Yearbook, the DRC accounted for 51% of Africa's copper mine production, and the mining and mineral processing sector in the DRC accounted for an estimated 22.1% of its gross domestic product (GDP). Many of the most important mining operations in the DRC consist of copper and cobalt production in the copper-belt region of the southern Katanga province. The belt stretches for 250 km between Kolwezi and Lubumbashi. Industrial copper production started in 1911 with the Union Minière du Haut Katanga (which became Gécamines), and the Musonoi mine has been known for many years in the western part of the Katanga Copperbelt. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 During mineral deposit exploitation, dewatering of subsurface strata is one of the most important tasks to ensure an engineering geologically safe and efficient mining operation [1]. Under complex geological condition, analytical solution [2,3] and analytical element models cannot be available, and thus numerical models are efficient tools to provide a dewatering plan in porousmedium areas [4][5][6][7] and fractured and karst regions [8][9][10][11]. The zonation of hydrogeology parameters is usually made based on borehole data and data interpolation [12]. Unlike layered soils [13], zonation is hard to achieve by data interpolation for complex geological setting, such as overturned synclines. 3D geology models, which integrate many kinds of field investigation data, are necessary and usually used tools to demonstrate the distribution of formations in mines activities [14,15]. Importantly, geological heterogeneity is associated with open mine voids, tectonic faults and fractures, and traditionally structured grids have great difficulty representing high heterogeneity. Geology models, from existing geological modeling software such as Petrel, GeoModeller and Leapfrog, provide a visual demonstration and knowledge of the distribution of formations, especially controlling aquifers or faults. However, few researches are concerned with combining groundwater flow models with 3D geology models in complicated aquifers. MOD-FLOW-USG [16] and FEFLOW 7 [17] have been developed to describe complex geology conditions with unstructured grid systems. Meanwhile, uncertainties are still present in aquifer geometry and its characteristics due to limited investigations and complexity of geological conditions. Hassen and Gibson et al. [18] construct a 3D geology model using GeoModeller software and evaluate the influence of aquifer connectivity and the role of faults and geology in groundwater flow within the aquifer system; however, to our knowledge, few groundwater flow models with coupled geology model are reported. The coupling of flow model and geology model are effective to improve the accuracy of numerical model because it provides both the visual display of geology conditions and groundwater flow pattern, and thus accurate estimation of dewatering will minimize the risks during the deep excavation of mining.
The Kolwezi megabreccia contains Cu-Co deposits hosted in folded and brittle-fractured structures of the Mines Subgroup. As indicated by the surface map and cross-sections [19][20][21], these structures are aligned more or less parallel to the E-W/60-SW-NE orientation of the thrust nappe [22,23]. The Musonoi mine is a newly planned underground mine in the Kolwezi Copper Deposit, DRC. The syncline strata in the Musonoi mine are overturned, with complex geologic and hydrogeological conditions. Several mines are present in the Kolwezi nappe area, including the KOV (Kamoto, Oliveira and Virgule) open pit [24], the Kamoto underground mine, Kamoto East, the old Musonoi open pit, the Dikuluwe open pit, the Mashamaba West and East open pits, the Sicomines mine, and the Kolwezi copper mine. Intense mining activities alter the natural path of groundwater flow and strengthen the groundwater connections among different mines in the Kolwezi copper mine. Intense mining activities alter the natural path of groundwater flow and strengthen the groundwater connections among different mines in the Kolwezi nappe area. Groundwater flow models have been developed to explore the dewatering strategies in the Kolwezi nappe area. The first public computer modeling of groundwater flow [25] is performed for the Dikuluwe and Mashamaba West pits and indicates that the progress of dewatering in 1985 would not be satisfactory for the existing mining plans. Africa Geo-Environmental Services (Pty) Ltd. develops groundwater flow models that consider the dewatering of KOV, Kamoto East, Kamoto Underground, Dikuluwe, Mashamaba East and Mashamaba West in 2006. SLR Consulting (South Africa) (Pty) Ltd. takes the Kolwezi nappe area as the model area to estimate the dewatering in the newly planned Musonoi mine in 2013 and 2017 using FEFLOW software. Such models simplify the strata into layered properties, and the simplifications affect model accuracies to some extent.
Although a lot of field investigations are carried out in the local mine area, few are present at the regional scale. However, the model boundary is required to extend beyond the mine area so as to avoid the influences of adjacent mining activities and man-made boundary on the accuracy of groundwater flow model. So the main objectives of this research are to combine groundwater flow model with regional and local geology model, and then explore the dewatering in underground mine using the established model. First, a hydrogeological conceptual model is designed for the underground mining activities. Second, the 3D geology model is constructed at regional and local scales and thus the parameterization will be carried out. Then, data will be prepared and then the groundwater flow model will be calibrated. Last, the established model is used to predict the dewatering in the Musonoi mine. The results of this paper can be effectively used for references in developing an efficient dewatering program.

Study area
The study area is located at the south of the equator in the Katanga plateau of the DRC (Fig 1). The study area has a savanna climate. The annual mean temperature is approximately 21.2˚C, the highest temperature is approximately 31.9˚C, and the lowest temperature is approximately 9.4˚C. The average annual precipitation from 1979 to 2017 is approximately 1144.90 mm and the average annual evaporation is approximately 1860.00 mm. The maximum annual precipitation is 1826.00 mm (in 1986), and the minimum annual precipitation is 803.50 mm (in 1993). Precipitation mainly happens from November to March of the following year and accounts for more than 85% of the annual precipitation. The dry season is from May to September, with little monthly precipitation of less than 5 mm. The overall terrain is high in the south and low in the north, with varying elevations from 1250 m to 1550 m. There are three main rivers in the study area, the Luilu River, the Musonoi River and the Dilala River. The Luilu River (including the Ptopoto River) and the Musonoi River flow towards the north. Catchment areas for the Luilu and Musonoi Rivers are 236 km 2 and 190 km 2 , respectively. The Dilala River surrounds the east and north sides of the mine area and finally joins the Musonoi River in the northwest of the mine area. During the dry season, the measured runoff of the Dilala River on the east side of the mine area is about 0.085 m 3 /s. The strata in the study area are of the Late Proterozoic Katanga supergroup, which can be subdivided into the upper Kundelungu group and the lower Roan group (host strata). Regional geology is mainly controlled by the Katanga Fold-Nappe structure. The secondary structure in the nappe body is complex with folded fractures. The nappe structure has an important influence on the water content and spatial distribution of the aquifer. The strata in this area are mainly the Katanga series and the quaternary. The Katanga series mainly includes the Roan group (R), the Nguba group (Ng) and the Kundelungu group (Ku). The series of geology from the young to the old can be seen in Table 1.

Conceptual model
Due to the long-term mining activities, the natural groundwater flow paths in each mine area are significantly altered. The hydraulic connections of the groundwater system in each mine area are much closer to the mining process due to the intense mining. Therefore, it is unreasonable to establish a numerical model of the groundwater flow in a certain mine area, and, based on the study of the digital elevation model of the study area, hydrology toolboxes in QGIS software are used to divide the watershed boundary. The model boundary, shown in Fig  1, considers the structure of the nappe body to minimize the influence caused by the setting of

Flowchart of the study
The flowchart is shown in Fig 2. Firstly, the model area is defined and conceptual model is established. Secondly, model area is discretized using FEFLOW [26] software into grid systems, and a layered groundwater flow model is used. Thirdly, all data, including borehole logging, pumping tests data, local geology model and hydrogeological reports, are prepared, and then a 3D geology model will be developed regionally. Next, other model data is collected and used to carry out model calibration by pumping tests in the Musonoi area. Finally, the model is used to predict the dewatering flux for the underground Musonoi mine.

Hydrostratigraphy and aquifer parameters
The geological conditions are very complex in the study area, especially the Kolwezi nappe structure. According to previous geological research reports [20], 16 geological profiles ( Fig  3A and 3B) are determined to divide the nappe bodies. The formation is reversed due to the presence of the overturned syncline, and the thickness differs in different locations. The geological modeling software is used to display the plane and section geological maps in the threedimensional spatial coordinate system for better presentation ( Fig 3B). Then, a three-dimensional geological model can be established to determine the parameters of the formation after raster map digitalization. The hydrogeological parameters are closely related to the structure and the degree of fragmentation caused by the stratigraphic structure. According to pumping test data, the hydraulic conductivity of the medium corresponding to the layer of CMN and RSF is approximately 0.65 m/d where the breccia zones are developed. The average hydraulic conductivity of the medium corresponding to the MWashya and RGS formation is 0.12 m/d, and the average hydraulic conductivity of media corresponding to the Kando and Kundelungu formation is approximately 0.026 m/d.

Interaction between the rivers and the aquifer
According to the field investigation, the flow rates at each gauging station are basically the same from the start and end station along the flow of the Musonoi River. Considering the possible existence of drainage from the KOV pit to the mainstream of the river, it can be speculated that the hydraulic exchange between the Musonoi River and groundwater is present but not very strong. The interaction between the Musonoi River and groundwater is calculated by the groundwater level, the vertical hydraulic conductivity of the semipermeable underlying layer of the river, and the river stage [27], which is represented as the third boundary type in FEFLOW software. Little information about the seasonal Dilala River and the Luilu River is collected, and it is assumed that groundwater gains the recharge from the rivers and the infiltration rate is estimated by 5 to 10% of the runoff in this study.

Sources and sinks
At the upper boundary of the system, the aquifer mainly receives recharge from precipitation, rivers, and the wastewater from tailing reservoirs. Groundwater discharge items mainly include mining drainage, drainage to abandoned pits, drainage to the river, pit leakage and groundwater evaporation. Groundwater evaporation occurs in the shallow depth areas of the groundwater. Recharge to the aquifer is mainly from precipitation during the rainy season. The precipitation infiltration coefficient is closely related to rainfall, the lithological structure of the aeration zone and buried depth of the groundwater level. According to previous studies, the rainfall infiltration recharge coefficient is generally between 10% and 25% but is only 3% in the RAT rock zone, 75% in some open pits, and up to 25% in other sandstone, waste craters and tailing ponds. In the model, the infiltration recharge rate is calculated by multiplying the annual average precipitation infiltration by the infiltration coefficient.

Water levels and pumping tests
Observation wells are distributed in the model area (Fig 1).

Model domain, mesh and model layers
The model domain has an area of approximately 623.63 km 2 . The study area is divided into 12,740 nodes and 24,988 elements, with a total of 18 model layers (the total number of nodes is 242,060 and the number of elements is 449,784). In the study area, the region of the Kolwezi nappe is refined, and the Musonoi mine area is refined again. Observation wells and pumping wells are considered as fixed points, and the mesh is divided along the shape of the river network. The dry season in the study area is from May to September, and the time period for model calibration is set from September to December. The initial distribution of water level is thus calculated based on groundwater flow simulation under the steady state and then adjusts after the input of model data [27]. The model calibration is carried out based on the estimated hydrogeological parameters, the trial-and-error method and FePEST tool. The multilayered wells were represented as the multilayered well module in FEFLOW 7. The model adopts automatic time-step control via the forward Adams-Bashforth/backward trapezoid (AB/TR) time integration scheme.

3D geology model
The ground elevation data used in this model are obtained by highly resolution elevation data from BIGEMAP software (http://www.bigemap.com/) with a horizontal resolution of 7 meters, which are publicly available. For the regional geology model, the depth range is taken from the ground surface down to 2000 m. The method is established based on grid systems as follows: 1. Establish a regional geological model in the study area. First, the plane geological map and 16 cross-sections are digitized. Then, the geological body is generated by using 3D data integration and Datamine RM software from all geological maps. The horizon clips of profiles at the 18 model layer depths are derived.
2. Develop a local geological model in the area of interest. The area of interest is approximately 6 km 2 where the borehole logging data is used. The refined 3D geological model in the area of interest is produced using borehole data; the breccia zone is separately described in the geology model.
3. Couple the regional and local geology models and generate the zonation of hydrogeological parameters. The regional and local 3D geologies are coupled by the same grid system, and the partition of lithology in the study area is then formulated. Fig 5 shows the detailed zonation of hydrogeological parameters at the regional (Fig 5A and 5B) and local (Fig 5C and  5D) scales.

Anisotropy and estimation of hydraulic parameters
High heterogeneity is present because of the overturned syncline in the area of interest. As indicated by the pumping tests, the maximum drawdown is approximately 61 m in the PW01 well, 18 m in the GM11 well and 18 m in the TB01 well; however, the distance between the pumping well of PW01 and the observation well is 130 m from the PW01 to the GM11 well and 44 m from the PW01 to the TB01 well, which suggests that anisotropy is still present in the study area and should be considered. According to the distribution of the formation, the angle between the direction of the anisotropy principal axis of hydraulic conductivity and the axis of the Cartesian coordinate system is generalized as the same value in the study area, and the angles of the principal axes are 60, 20 and 20 degrees from the x, y and z axes, respectively. The hydraulic conductivity of the long anisotropy principal axis in the plain is estimated by the pumping test data from engineering boreholes in the neighboring mine areas, and the value for each model layer is shown in Table 2. The hydraulic conductivity of the short anisotropy principal axis in the plain is set as the same as that in the long axis. The hydraulic conductivity of the vertical anisotropy principal axis is set at 1/10 of that in the long axis. Specific yield varies from 0.01 to 0.05 and is set to 0.02, 0.02, 0.04, and 0.05 for the CMN, SD, Kamoto, and Breccia zones, respectively. Storativity varied from 0.00001 to 0.0001 m -1 and is set to 0.0002, 0.0002, 0.0002, and 0.0001 m -1 for the CMN, SD, Kamoto, and Breccia zones (fractured zone), respectively.

Sensitivity analysis of hydraulic conductivities and boundary conditions
Hydraulic conductivity (K) and recharge boundaries such as precipitation and surface water infiltration are critical factors for model calibration, and their uncertainties are required to assess. In this study, local sensitivity is calculated based on the single-factor variation method. The main factors include hydraulic conductivities in the Ku zone, the RAT zone, the Kamoto zone, the SD zone, the CMN zone and the Dipeta zone, and infiltration from rivers and reservoirs. The base scenario (S0) is set as initial estimated parameters. Scenarios are set by a factor with 25%, 75%, 125% and 150% of the initial parameters (shown in Table 3). The simulation  (1)) is defined as the root mean square error of the simulated water level difference from between the scenario and the base scenario.
where SA si is the sensitivity index for the scenario Si (i = 1,2,3,4), m; N is the total number of simulated values for 28 wells; h si and h s0 are simulated water level from the scenario Si and the base scenario S0, respectively.
Sensitivities of hydraulic conductivities for six zones, and surface water infiltration on groundwater level changes are shown in Fig 6. When hydraulic conductivities in the CMN, SD and RAT zones change, significant groundwater level changes are found, whereas hydraulic conductivities in the Kamoto, Dipeta, and Ku zones have smaller sensitivity. The scenario Cm1 has the biggest sensitivity (4.99 m). However, the sensitivity of river and reservoir infiltration is low (less than 0.10 m). So the hydraulic conductivities in the CMN, SD and RAT zones are the controlling factors which should be considered in the calibration period.

Model calibration method
The calibration period is set from 8  from 28 observation wells in the area of interest. The 28 observation wells mentioned above are selected as the fitting object to estimate the hydrogeological parameters. The objective function is set for this purpose as in Eq (2), and drawdown is used for metrics because the accurate head distribution in the Musonoi mine is not available. Under the constraints of the upper and lower limits of each parameter, the objective function G is the minimum. If the difference between simulated and observed water level during the calibration period is not sufficiently small, the values of some parameters are adjusted and the objective function is gradually minimized. First, the trial-and-error method was used to manually adjust the parameters. FePEST in FEFLOW was then used to estimate the parameters.
where p 1 , p 2 , . . ., p n are the hydrogeological parameters to be estimated; N g and N tc are the number of observation wells and the time with maximum drawdown, respectively; ω h is the weight of the groundwater level, here is set to 1 for 28 observation wells in the area of interest; S(i,j) is the simulated drawdown at the j th time steps for observation well i (m), and S g (i,j) is the observed drawdown at the j th time steps for observation well i (m).

Comparison of observed and simulated drawdown
The simulated drawdowns are compared with the observed measurements, and the statistics are shown in Table 4. In the area of interest, approximately 41% of the absolute differences between the observed and simulated drawdown are within 0.5 m, approximately 56% of the absolute differences are less than 1.0 m, and approximately 22% of the absolute differences are bigger than 5.0 m. A total of 12 observed wells are selected to demonstrate the comparison of drawdown between the simulated and observed measurements over the calibration period (shown in Fig 7). During the pumping test period, the maximum drawdowns are located around wells PW01, PW02, PW03, TB01, TB02 and TB03. Fig 8 shows the comparison between the observed and simulated maximum drawdowns among six wells (PW01, PW02, PW03, TB01, TB02 and TB03) and the relative errors. The relative error of simulated maximum drawdowns is within 10% for wells PW01, PW02, TB01 and TB02 and approximately 15% for wells PW03 and TB03.

Calibrated hydraulic parameters
The hydraulic conductivity along the long anisotropy principal axis in the plain for each parameter zonation of each model layer in the study area after the calibration is shown in Table 5. The hydraulic conductivity along the long anisotropy principal axis is from approximately 0.05 to 0.43 m/d for the breccia zone, which is lower than the estimated value. Additionally, hydraulic conductivity from layers 5 to 8 is sensitive to the objection function in Eq (2) and is thus refined after calibration. The hydraulic conductivity of the short anisotropy principal axis in the plain is approximately 90% of the value in the long axis. The hydraulic conductivity of the vertical anisotropy principal axis is set at 1/15 of that in the long axis.

Water table change and distribution of maximum drawdown
Usually, the initial head of each node is interpolated from the observed heads in observation wells. There are 18 model layers in the study area, and it is thus hard to obtain the initial head for all model layers because the wells are not well distributed in each layer and over the study area. The initial head is simulated by the coupled head and parameter method [27]. The cone of depression of the water table is clearly simulated and seen in the model area. In the area of interest, groundwater flows from the eastern and southern boundaries. Because the distance between the area of interest and the adjacent mines is short, the groundwater flow pattern will probably be altered after the intense mining activities in the area of interest, and it is uncertain whether the groundwater divide is present between the two mines. The distribution of maximum drawdown during pumping tests can provide insightful understanding of dominant pathways for groundwater flow. The development of cone of depression in response to the pumping is mainly subjected to the area of the fault-controlled fracture zones (about the depth of about 150-350 m. Although the multilayered observation wells are not of the same depths as the boreholes, the maximum drawdowns of all wells are demonstrated in the plain (Fig 9). The distribution of maximum drawdown is highly un-uniform. The long axis of maximum drawdown is approximately 45˚North of East, and the length of influence is approximately 1.50 km. The short axis of maximum drawdown is approximately 45˚North of West, and the length of influence is approximately 1.0 km. Observation wells at the vertical depth of about 150-350 m have the biggest drawdown and then the wells in the shallow depth have the lowest drawdown.
In the area of interest, the drawdown induced by dewatering is so large that drying of the aquifer usually happens. During the calibration, the phenomenon of drying and rewetting of the model layer (Layer 1) is observed (Fig 10). Under the influence of pumping well PW02, the change of groundwater level with time is normal in site P1 but abnormal in site P2 (here the locations of P1 and P2 are very close to PW02). The layer interface between model layer 1 and

PLOS ONE
2 is marked in Fig 10; two sudden water level jumps are present at the beginning and end of the pumping tests, which usually indicates convergence problems and low calculation efficiency.

PLOS ONE
Estimating dewatering using a 3D finite element model

Dewatering flux prediction
According to the draft mining plan, the first phase during the mine service period of the area of interest is to the depth of 140 mL and fully dewatered by pumping wells. Three pumping wells have been drilled during the period of pumping test, namely PW01, PW02 and PW03. The area of cone of depression under current pumping wells are relatively small, and additional wells should be added to extend the area. 140 mL dewatering plan will be achieved within 1 year, and 8 scenarios are thus considered (listed in Table 6).  Fig 11. The simulated results under the eight scenarios are shown in Fig 12. The cone of depression with the drawdown of 140 m is required to cover the area of the mine area (as indicated by the red color in Fig 12). The areas of cone of depression under the eight scenarios are basically within the area of ore body, but the areas of the cone are different for scenarios. The contour map of drawdown pattern is an elliptical shape with a major axis from northeast to southwest. The scenario S5 has the largest area of the cone, while the scenario S1 has the smallest area. As can be seen from the figure, scenario S1, S2, S3, S4, S5, S6, S7 and S8 can form continuous cone area of drawdown with 80 m, 80 m, 90 m, 80 m, 150 m, 140 m, 120 m and 110 m, respectively. The total dewatering flux for scenario S1, S2, S3, S4, S5, S6, S7 and S8, is 16,400 m 3 /d, 17,100 m 3 /d, 18,900 m 3 /d, 15,900 m 3 /d, 28,400 m 3 /d, 22,800 m 3 /d, 23,900 m 3 /d and 18,900 m 3 /d, respectively. The scenario S5 have the highest dewatering flux, whereas the scenario S4 have the lowest flux. Although the area of cone of depression formed by the scenario S5 is large, the invest for drilling will be higher. So the scenario S7 is recommended when considering the area of dewatering and cost for the drilling and operation. It can be seen from the figure that the scenario S7 (Fig 12) can form a cone of depression with the drawdown of 120 m in the middle of 140 mL shaft after one year with maximum drawdown of 320 m, and the drawdown of 200 m at local regions.
Based on the simulated patterns of drawdown changes for the scenario S7, it can be found from the results that the cone of depression at the sixth model layer is the largest, and the maximum drawdown is nearly 320 m. The maximum drawdown decreases from the sixth model layer to the first model layer, and also from the sixth model layer to the eighth. From the existing pumping tests, the maximum production rate of PW03, PW01 and PW02 wells is 2700 m 3 / d, 1200 m 3 /d and 3500 m 3 /d, and the amount of production rate for additional wells in Table 6 can be reasonable. The total water volume of all wells in the middle section of 140 mL shaft is about 23,900 m 3 /d for the scenario S7. So the dewatering capacity with 23,900 m 3 /d will be required at least for the dewatering plan.
Pumping tests are carried out during the dry season, where the infiltration is about 2.19 mm/day. Different seasons such as normal season and wet season will also have an important influence on the dewatering. The model is also used to evaluate the influence of wet and normal season on the change of groundwater level. The precipitation is 1952 mm/year and 1200 mm/year for wet season and normal season, respectively, and then used as model input. More infiltration will increase groundwater level. For the scenario S7, relative to the dry season, the difference of drawdown is about 2-35 m under the wet season, and 1-30 m under the normal season in the area of interest. The dewatering activities will alter the water pressure conditions, and then engineering geologically safety such as collapse may occur. These geologically safety issue should be seriously paid attention during the underground mining. For the safety of  NW101  PW01  NW102  PW02  PB02   S1  2000  2700  2500  1200  2500  3500  2000   S2  2000  2700  2700  1200  3000  3500  2000   PB01  PW03  NW201  NW202  PW01  NW203  PW02  PB02   S3  2000  2700  2500  2500  1200  2500  3500  2000   S4  2000  2700  1500  1500  1200  1500  3500  2000   PB01  PW03  NA01 NW201 NA02  NW202  NA03  PW01  NA04  NW203  NA05  PW02  PB02   S5  2000  2700  2500  2500  2500  2500  1500  1200  1500  2500  1500  3500  dewatering, the change of groundwater levels should be monitored during the implement of dewatering project.

Uncertainty of model results
Due to the high heterogeneity and anisotropy of geology formations, the simulated results have uncertainties. The main effects of the uncertainties are as follows: 1. Uncertainty of geological formations. Although this paper establishes the 3D geology model by coupling the regional 3D geological model based on geological cross-sections with a local refined geology model based on borehole data, local fractures may exist and the distributions of fractures remain still unknown due to limited field investigations.
2. Uncertainty of hydrogeological parameters. Pumping test data in the area of interest are used to calibrate the model. However, the three pumping wells are located at depths of 250-320 m below the ground surface, and observation wells are not well distributed in the area and in the vertical direction. The limited duration of a pumping test (less than 3 months in this study) is usually not sufficient to track the seasonality and changes in groundwater heads. From the calibration on pumping test data in the area of interest, obtained aquifer properties may not be representative over the whole extent of the model layers and model area.
3. Uncertainty of groundwater recharge. The model mainly uses the infiltration from precipitation in the dry season to estimate groundwater recharge, and the precipitation in wet years has different effects on the change of groundwater flow. No complete pumping test was conducted throughout the whole hydrological year, including the wet season, normal season and dry season. Therefore, the estimated dewatering may be underestimated.

Conclusions
This study collects much data, such as precipitation, river runoff, groundwater level and the results of pumping tests, and constructs a three-dimensional groundwater flow model in an area of interest with a 3D geology model used to aid in the prediction of the dewatering of underground mines. Model calibration is conducted using pumping test data, and the characteristics of groundwater flow are discussed. The following conclusions are obtained.
1. Development of a 3D geology model for complicated formations can visualize the distribution of formations and facilitate the zonation of parameters in the groundwater models. The description of complex geological structure in the study area is challenging work. The geological information about sixteen hydrogeological cross-sections and the geology distribution in the plain are first digitized and then adopted to construct a regional geology model by using mining software. The local geology model in the mine area is refined based on borehole logging and fractures. Geology information from the two models is combined to formulate the 3D aquifer model for the study area. The description of the regional complex geological structure provides better conditions for improving the accuracy of the model.

2.
A groundwater flow model was conceptualized, including the scope, the boundary conditions, the representation of multilayered wells, and sources and sinks. Although the area of Writing -review & editing: Zhengqiu Yang.