A three-dimensional mathematical model is used to simulate the groundwater flow in the aquifer system in the study area. The mathematical equation is as follows:
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
Where kxx, kyy and kzz are the permeability coefficients (m/d) along the x, y and z axes respectively; H is the water level elevation of shallow groundwater system and middle-deep groundwater system (m); S is the water storage rate of point (x, y, z) (1/m); W is a source-sink term (1/d) and a function of time and space W=W(x, y, z, t); T is time (d); N is the normal vector within the boundary of the second kind; Q(x, y, z, t) is the variation function of flow in time and space; F2(x, y, z, t) is a function of water level change in time and space; D is an infiltration basin jointly determined by shallow groundwater system and middle-deep groundwater system; γ 1 is the known water level boundary of seepage calculation domain d; γ 2 is the known flow boundary of seepage calculation domain d.
(2) Numerical simulation model of groundwater
In order to evaluate the groundwater resources in the study area more accurately and predict the state of groundwater system and its changing trend under the condition of coal mine development and utilization, a numerical simulation model of groundwater flow in the study area is established by using Feflow software based on finite element method. After identification and inspection, the groundwater seepage law in the study area is simulated, and the future groundwater dynamics are analyzed and evaluated.
1. Simulation software selection and grid generation
The software used in this work is Feflow software based on finite element method developed by WASY company in Germany. To establish the numerical simulation model of groundwater flow in the study area, the simulation area should be triangulated first. In addition to following the general principles of division (for example, the internal angles of triangular units should not be obtuse as far as possible, and the area difference between adjacent units should not be too large), the following actual conditions should also be considered: ① Fully consider the boundaries of the study area, lithologic zoning boundaries, faults, etc. ; (2) Observation wells and water sources should be placed on the nodes of subdivision units as far as possible; (3) In areas where the hydraulic gradient changes greatly and the research is focused on, the cutting should be properly densified. The divided simulation area * * * has 4473 nodes and 6996 units (Figure 5-5). The simulated aquifer includes coal14-k3 (Ⅱ) aquifer and Ordovician limestone aquifer. The three-dimensional structure of the study area is shown in Figure 5-6.
Figure 5-5 Schematic Diagram of Triangulation Network in the Study Area
Figure 5-6 Three-dimensional structure of the study area
Because the boundary of the calculation area is the fault boundary, not the boundary of regional stratigraphic change, there are some errors in boundary treatment. By processing and analyzing the observation data of groundwater quantity and water level in the study area and its periphery for many years, and referring to relevant research results, the boundary of model calculation is preliminarily determined.
(1) Vertical boundary: The upper boundary of the calculation simulation area is the Quaternary confined aquifer, which is a water exchange boundary with changing position, and there is a close hydraulic relationship among atmospheric precipitation, surface water and groundwater. The lower boundary of the calculation simulation area is the floor of the confined aquifer, which is composed of bedrock (mudstone, shale, etc.). ) the permeability is extremely poor, and it is generalized as a water-proof boundary.
(2) Lateral boundary: the groundwater watershed of confined aquifer is treated as zero flux boundary; Groundwater movement can be generalized as three-dimensional spatial flow; The input and output of groundwater system change with time and space, so groundwater is an unstable flow; This parameter varies with space, reflecting the heterogeneity of water-bearing medium, but it has no obvious directionality, so it can be generalized as isotropic. To sum up, the groundwater flow system in the study area can be summarized as heterogeneous, isotropic, three-dimensional spatial structure and unstable groundwater flow system.
2. Finite element numerical solution of groundwater system simulation.
The finite element method was first put forward by China mathematician Feng Kang and others in, and was mainly used in mechanics at first. At the end of the 20th century, the finite element method was introduced to solve the problem of groundwater flow. With the rapid development of computer technology, finite element method has become an effective method to solve complex hydrogeological seepage problems. The finite element discretizes the differential equations solved continuously in the region into linear algebraic equations by subdivision interpolation, and replaces the exact solutions of the differential equations with approximate numerical solutions. According to different mathematical equations, the finite element method can be divided into Ritz finite element, equilibrium finite element and Galerkin finite element.
Ritz finite element is based on variational principle and subdivision interpolation. Variational principle is to transform the solution of partial differential equation describing groundwater movement into a functional extreme value problem; Subdivision interpolation is to divide the studied seepage area into point, line and area units in geometry, assume the head value of nodes, then use some form of interpolation method to interpolate the units according to the actual situation, and finally form the interpolation of the whole unit geometry by constructing the head expression of each unit. Based on the variational principle, this method transforms the problem of finding the extreme value of a functional into the problem of solving a group of multivariate linear algebraic equations by using the interpolation of the whole unit set, so as to get the solution of the seepage problem. Balanced finite element method is to triangulate the seepage area from the perspective of small balance, and then establish the equilibrium equation of the area composed of the connecting line between the center of gravity of each triangle and the midpoint of the corresponding side with any node as the common vertex, and then get a linear algebraic equation. For the whole research area, a system of linear equations will be obtained. Although the balanced finite element method is simple and intuitive, it is considered to be inaccurate in mathematics. Galerkin finite element method discretizes the continuous differential equation by residual weighting method, so as to find the numerical solution of the control equation satisfying the hydrogeological constraints. Its mathematical principle is as follows.
(1) principle of residual weighting method: residual weighting method is an effective method to solve approximate solutions of differential equations, and its mathematical model is as follows:
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
Its analytical solution must satisfy the differential equation and the corresponding constraints of the equation in this region, that is, at any point in the above formula. For complex problems, analytical solution is not easy to solve, so to find an analytical solution that can meet a certain accuracy, it is generally an approximate solution of analytical solution with a preset error range.
Currently using the trial function:
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
Where: φi is a known function selected according to specific requirements; φ 1, φ2, …, φn are linearly independent basis functions; αi is the undetermined coefficient.
Because it is an approximate solution of the equation, it will be substituted into the above formula, so that the right end of the equation is not 0, resulting in a remainder called R, that is, L (V) = R.
If several weight functions Ω1,Ω 2, …, Ω n corresponding to the undetermined coefficient are selected by some method, and the weighted integer value of the remainder r is required to be equal to 0 with the weight function, then an equation group composed of several equations can be obtained, namely
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
Where: ω is the integral interval, which can be any one of line domain, area domain or volume domain. This formula is a residual weighted integral formula. By solving this formula, several undetermined coefficients α 1, α2, ..., αn can be determined, and the approximate solution of the mathematical model can be obtained.
(2) Galerkin method: In practical application, there are many ways to select the weight function ω i. If the selected weight function is the same as the basic function φi, it can be written as i= 1, 2,3, …, n, and this formula is called Galerkin equation.
Galerkin finite element method can solve the three-dimensional flow problem of groundwater. It is assumed that the calculation area ω is a heterogeneous anisotropic water-bearing medium, and the boundary is composed of γ = γ 1+γ 2, where γ 1 is the first boundary of the water-bearing system and γ 2 is the second boundary. If the selected spatial coordinate system is consistent with the direction of the main permeability coefficient of the aquifer, the mathematical model of three-dimensional groundwater flow can use the following equation:
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
In the above formula: H(x, y, z) is the distribution of groundwater head in aquifer; Kxx, kyy and kzz are the main permeability coefficients of aquifer respectively; T is a time variable; ω is the seepage space of the study area; γ 1 and γ 2 are the boundaries of the first and second study areas; W is the algebraic sum of source term and sink term; Ss is the water storage rate of aquifer; H0 is the head distribution of the initial flow field; H 1 is the first kind of boundary head distribution; Q is the second boundary flow.
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
According to the residual weighting theory, the trial function is taken, and the formula satisfies the definite solution conditions of the equations. Where {φ I (x, y, z), I = 1, 2, ..., n} is the constructed basis function group, and {Ci(t), I = 1, 2, ..., n} is the time correlation coefficient. For a fixed t, these coefficients can be solved by the following equation:
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
Variables y and x have similar expressions, namely:
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
The above equation is written as a matrix of the following form
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
The elements of matrix [H] and [P] in the formula are respectively
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
The elements of the column vector f are
Numerical simulation of groundwater movement and pollution in typical coal mines Application of Feflow and Modflow
As long as the basis function {φ I} is determined, all the above variables can be calculated, so that the approximate solution of the groundwater flow motion control equation can be solved.
In practical application, the main steps of Galerkin finite element method to solve the groundwater motion equation include: ① dividing the research area into finite seepage units in space to form a grid of seepage units; (2) Calculating the plane area, spatial volume and dielectric constant of seepage units according to the divided node coordinates, wherein the node coordinates can be absolute coordinates or relative coordinates, and establishing a water conducting matrix of each unit; (3) according to the flow on the head expression of each node, the known flow of each unit is moved to the node according to the principle of equal sharing, that is, the equivalent flow; (4) establishing a continuous equation for each unknown node head and solving the equation.
All the above steps are automatically completed by Feflow.
3. Treatment of simulation period
In this study, the flow field in June 5438 +2005 10 is taken as the initial flow field for the calculation and evaluation of groundwater resources, and the flow field in June 5438 +2007 10 is taken as the fitting flow field for the calculation and evaluation of groundwater resources. The simulation period is from 65438+ 10 in 2005 to 65438+ 10 in 2007, and one month is a time period, each time period includes several time steps, which is automatically controlled by the model.
4. Treatment of definite solution conditions
Initial conditions: According to the groundwater level observed in June, 5438+October, 20051October, the initial water level of the aquifer is obtained by interpolation and extrapolation.
Boundary conditions: the parameters of each flow boundary mainly consider the flow field at the initial stage and the end of the simulation, and the boundary with detailed data matches the inflow and outflow of the boundary. The time step is controlled by the program, and the error is strictly controlled in each operation.
Overview of water inrush points:
According to the ratio of water inflow at each water inrush point, the statistical data of water inflow of roadway is distributed to the water inrush point, which is generalized as a pumping well (Figure 5-7) to simulate the drainage effect of roadway on aquifer.
Figure 5-7 is a schematic diagram of the distribution of water inrush points in roadway after pumping well.
5. Identification and test of the model
The process of model identification and verification is an extremely important step in the whole simulation, and it usually needs to modify parameters and adjust some source and sink items repeatedly to achieve the ideal fitting effect. This method of model identification and test is also called trial estimate correction method, which belongs to one of the indirect methods of inverse parameter calculation.
By running the calculation program, the temporal and spatial distribution of groundwater level under given hydrogeological parameters and various equilibrium terms can be obtained. Through the flow field in quasi-contraction period and the diachronic curve of long observation well, the equilibrium items such as hydrogeological parameters and boundary values are identified, which makes the established model more in line with the hydrogeological conditions in the study area, thus quantitatively studying the recharge and discharge in the simulation area more accurately and predicting the groundwater level.
The identification and verification of the model mainly follow the following principles: ① The simulated groundwater flow field should be basically consistent with the actual groundwater flow field, that is, the simulated groundwater level isoline should be similar to the measured groundwater level isoline; ② The dynamic process of simulated groundwater should be basically similar to the measured dynamic process, that is, the shape of simulated groundwater level hydrograph should be similar to the actual groundwater level hydrograph; (3) From the point of view of balance, the simulated groundwater balance changes are basically consistent with the actual situation; ④ The determined hydrogeological parameters should conform to the actual hydrogeological conditions.
Because the hydrogeological conditions in Linnancang Mine are complex, the boundary conditions are short of quantitative data, and the degree of hydrogeological exploration in the mining area is low, the observation wells with systematic long-term data are few and unevenly distributed, and there are many uncertain factors in the model, which brings great difficulties to the parameter identification of the model. In order to overcome the fuzziness and instability of parameter identification caused by insufficient data, the observation data of hydrogeological conditions, hydrodynamic field and tectonic stress field in this area are studied in detail to help analyze the spatial distribution law of parameters and determine the value space of parameters. In the process of forward and backward modeling, various means are adopted to make the solution and prediction of parameters, water level and water quantity conform to the objective reality as far as possible. See fig. 5-8 ~ fig. 5- 13 for the change of head height of each observation well during the simulation.
Of the 7 observation wells, the observation data of No.6 observation well is only 4 months, so only 6 observation wells are actually used (the observation data of No.7 observation well is 15 months). The simulated value of 139 times is compared with the observed value. According to the calculation results, the error of fitting results is shown in Table 5- 15.
Figure 5-8 1 # Fitting Diagram of Water Level of Observation Well
Figure 5-9 Fitting Diagram of Water Level of 2 # Observation Well
Figure 5- 10 3# observation well water level fitting diagram
Figure 5- 1 1 4# observation well water level fitting diagram
Figure 5- 12 5# observation well water level fitting diagram
Figure 5- 13 7# observation well water level fitting diagram
Table 5- 15 Error of Simulation Results
Table 5- 16 Statistical Table of Fitting Errors
Hydrogeological parameters include the permeability coefficient and water supply of aquifers in all directions. Through fitting, the fitting diagram of seepage field between Ordovician limestone karst aquifer and coal aquifer 14-K3 is obtained (Figure 5- 14 ~ Figure 5- 15). As can be seen from the figure and table 5- 16, the fitting effect is very good.
Figure 5- 14 Fitting Diagram of Seepage Field of Ordovician Limestone Aquifer
Figure 5- 15 Fitting Diagram of Seepage Field of Coal Seam Aquifer 14-K3
According to the long observation well fitting curve and the flow field fitting curve, the established simulation model basically meets the accuracy requirements of the model, conforms to the hydrogeological conditions of the study area, basically reflects the dynamic characteristics of the groundwater system, and can be used for groundwater level prediction.
Finally, the determined hydrogeological parameters are shown in Figure 5- 16 ~ Figure 5-2 1 respectively.
Figure 5- 16 X-direction permeability coefficient of Ordovician limestone aquifer
Fig. 5- 17 permeability coefficient of Ordovician limestone aquifer in y direction
Figure 5- 18 Z-direction permeability coefficient Figure 5. 18 Z-direction permeability coefficient in Ordovician limestone aquifer
Figure 5-/KOOC-0/9 Coal Seam/KOOC-0/4-K3 (Ⅱ) X-direction Permeability Coefficient Figure 5-/KOOC-0/9 Aquifer/KOOC-0/4-K3 (Ⅱ) X-direction Permeability Coefficient
Figure 5-2014k3 (Ⅱ) Permeability coefficient of aquifer in Y direction Figure 5.2014k3 (Ⅱ) Permeability coefficient of aquifer in Y direction
Figure 5-2 1 coal seam14-k3 (Ⅱ) Z-direction permeability coefficient of aquifer Figure 5-2 1Z-direction permeability coefficient14-k3 (Ⅱ)
6. Seepage field prediction
Using the previously determined model, the simulation period is adjusted, and the seepage field of 14-K3 coal seam aquifer and Ordovician limestone aquifer after 10 is predicted under the condition that the current mining conditions and drainage conditions are unchanged (Figure 5-22 ~ Figure 5-23).
Fig. 5-22/KOOC-0/4-K3 Coal Seam Aquifer/KOOC-0/0 Seepage Field Prediction Diagram of Tu Tu 5.22/KOOC-0/4-K3 Coal Seam Aquifer Ten Years Later
Fig. 5-23 10 Seepage Field Prediction of Ordovician Limestone Aquifer.
7. Water inrush analysis
According to the water inrush coefficient method, the water inrush risk is determined by the water pressure and the thickness of the aquiclude. The higher the water head, the thinner the effective water-resisting layer and the greater the risk of water inrush. According to the prediction results of seepage field, the risk of water inrush is analyzed below.
Based on the identified seepage model, the seepage field changes of 14 sandstone aquifer and 20 19 Ordovician limestone aquifer are predicted (see Figure 5-24 and Figure 520). The water level of 14-K3 aquifer is basically flat, showing the characteristics of high in the north and low in the south. In 20 19, the groundwater gathered due south, and the water head was still at a high level, especially in the aquifer, which was about -290.67 m (see Figure 5- 19). Under other conditions unchanged, the groundwater in the 14-K3 aquifer can easily enter the roadway through the vertical water filling channel under high pressure, resulting in water inrush or water inrush, so the water inrush risk is greater. In the northern region, due to the high hydraulic gradient of underground water in the coal 14-K3 aquifer and the dense isograms of underground water level, once the water in the coal 14-K3 aquifer is affected by mining cracks, it is likely that the water in the coal 14-K3 aquifer will first enter the mining cracks, which may lead to the increase of pore water pressure and water inrush. Therefore, attention should be paid to the observation and control of mining cracks in this area. The water level variation law of Ordovician limestone aquifer is basically similar to that of 14-K3 aquifer, showing the characteristics of high in the north and low in the south. 20 19 groundwater is collected to the south, and the collection area is less than 14-K3 aquifer. However, the water head is still high, especially in the aquifer, which is about -3.67 m (see Figure 5-23). Under other conditions unchanged, the groundwater in Ordovician limestone karst aquifer can easily enter the roadway through the vertical water filling channel under high pressure, resulting in water inrush or water inrush, so the risk of water inrush is greater. In the northern area, because of the large hydraulic gradient, dense isograms of underground water level and high water pressure in the Ordovician limestone karst aquifer, the water in the Ordovician limestone karst aquifer is likely to enter the mining fracture preferentially through the water-conducting structure, resulting in the increase of pore water pressure and water inrush, so it is necessary to make long-term dynamic observation on the water level in the Ordovician limestone aquifer.
8. Mine water inflow prediction
Using the previously identified model, the groundwater level flow field in June 2005 (1) is brought into the model as the initial water level. Four imaginary wells are evenly distributed in the mining area, and the pumping capacity of each well is adjusted to make the stope water level drop to a reasonable position. Then the total pumping capacity of the four wells is the future mine water inflow, and the running time of the model is the future drainage time.
For the specific drainage water level, besides the water pressure on the floor of 14 coal seam, it is also necessary to consider the characteristics of large Ordovician limestone water volume and difficult drainage, and the drainage of aquifer in 14-K3 coal seam is the key factor to be considered. Four pumping wells are evenly distributed in the mining area (Figure 5-24), and the pumping capacity of each well is 600m3/d, that is, the current drainage capacity is increased by 10%, 30% and 50% respectively. It is judged that the water level of the coal 14-K3 aquifer has obviously decreased, and the drainage economy (Figure 5-25 ~ Figure 5-) shows that the water level of the 14-K3 aquifer has obviously decreased under the condition of strong drainage.
Figure 5-24 Location of New Pumping Hole
Fig. 5-25 14k 3 Aquifer Seepage Field 10% after the pumping quantity is increased; Fig. 5.25 14-K3 Coal Seam Aquifer Seepage Field 10% after the pumping quantity is increased.
Fig. 5-26 Seepage Field of Aquifer 14K3 after 30% pumping increase Fig. 5.26 Seepage Field of Aquifer 14-k3 after 30% pumping increase.
Fig. 5-27 Seepage field of 14K3 aquifer with pumping capacity increased by 50%.
As can be seen from the figure, the effect of dewatering and depressurization is obvious, but relatively speaking, improving the drainage capacity by 30% is similar to reducing the groundwater level by 50%. Therefore, in the future mining scheme, the scheme of improving drainage capacity by 30% should be chosen. According to this plan, when the drainage capacity is increased by 30%, the economic pressure is relatively small and the effect is quite obvious. It is suggested to increase the drainage pressure by 30%.
9. Analysis of water inflow prediction results
Due to the powerful function of Feflow, its fluid flux analysis module can carry out Flu xinside and Flux outside analysis on a given block in a given time period in the simulation period, and predict the water inflow in areas A, B, C, D, E and F6 of coal seam aquifer 14-K3 in different years in the simulation period (Table 5- 17, Figure 5-28).
Table 5- 17 Prediction Results of Water Inflow in Coal 14-K3 Aquifer
Figure 5-28 zoning map of water inflow
According to the prediction results of water inflow from coal 14-K3 aquifer provided in the above table, the linear trend in the above table is analyzed (Figure 5-29 ~ Figure 5-34):
Figure 5-29 A Area Water Inflow Prediction Result Diagram
Figure 5-Prediction Results of Water Inflow in Area 5-30 B
Figure 5-3 1 C Area Water Inflow Prediction Result Diagram
Figure 5-32D Prediction Results of Water Inflow in Area
Fig. 5-33 forecast results of water inflow in e area.
Figure 5-34 F Area Water Inflow Prediction Result Diagram
It can be seen from the forecast diagrams of water inflow in Area A, Area B, Area C, Area D, Area E and Area F that the water inflow in Area B, Area C, Area D, Area E and Area F fluctuates, but it can be seen from the linear diagram that the trend of water inflow is increasing. Therefore, it is necessary to dewater and depressurize the coal 14K3 aquifer.
See Figure 535 for the specific water inrush path, from which it can be seen that:
Figure 5-35 Plan of Water Inrush Path
Figure 5-36 perspective view of water inrush path of existing water inrush point