Establishment of 1. numerical model
1. Hydrogeological conceptual model
As mentioned above, the karst water system in Niangziguan spring area is surrounded by a relatively closed "spontaneous flow" type full drainage system except for a small amount of undercurrent near the spring mouth. The medium structure of the aquifer is dominated by dissolved pores and fissures, and there are anisotropic water-bearing rock masses with karst caves in some areas. According to the groundwater flow field, Niangziguan Spring generally sinks fan-shaped from north, west and south to east, which can be further divided into two recharge areas: the confluence area from Yangquan to Chengju and the discharge area at the spring mouth. Corresponding to the first three areas, there are three strong runoff zones, one is the northern karst groundwater level strong runoff zone from Yuxian along Wenhe to Cheng Ju, the other is the southern karst groundwater level strong runoff zone from Xiyang along the contact zone between Middle Ordovician and Carboniferous and Permian to Pingding Minmetals to Yangquan, and the third is the central karst groundwater level strong runoff zone from Yangquan to Chengju to Quankou confluence area. According to the characteristics of karst media in the north, it can be considered that two-dimensional plane flow obeys Darcy's law. The attenuation coefficient of spring flow in Niangziguan is between N 10-3 ~ N 10-8, which also shows that the karst development is mainly fractures and caves, and the existence of non-Darcy flow is not excluded in some large karst cave development areas, and these special areas are not treated specially in the calculation. In addition, the model ignores the movement of groundwater basin and does not consider the vertical movement in groundwater.
2. Mathematical model and its solution
According to the hydrogeological conceptual model of karst water in Niangziguan spring area, the mathematical model expression of water migration is as follows
Environmental problems and protection of karst groundwater in northern China
Considering that the variation of groundwater level relative to aquifer thickness is very small and can be ignored, the aquifer thickness can be approximately regarded as a certain value, and its migration model is as follows.
Environmental problems and protection of karst groundwater in northern China
Where: h is the head of the aquifer; Tx and Ty are hydraulic conductivities in x direction and y direction; Q is the vertical water exchange capacity of aquifer; H0(x, y) is the initial flow field of karst groundwater; H 1(x, y, t) is the first boundary water point; μ is the water storage coefficient; Q(x, y, t) is the second kind of boundary single-width replenishment.
3. Discrete mathematical model
The irregular finite difference method is used to divide the definite solution of the above-mentioned karst groundwater calculation area, and many small triangular bins are used to approximate the groundwater level surface. The whole spring area is divided into 597 triangular bins and 343 nodes, and the discrete definite solution of the above groundwater seepage is:
Environmental problems and protection of karst groundwater in northern China
In which: β is a β-binding composed of I-node equilibrium domain; P is the number of bins composed of I-node equilibrium domain; , which is the water level at the vertex t and t- 1 of β bin J node; Δ t is the time step; D is a water-conducting substrate;
=(×ai×Al+×bi×Bl)/(8×δβ)
E is the water storage matrix; e0i =- 1/( 16×△β)[(ai×Aj+bi×Bj)×(+)+(ai×Ak+bi×Bk)×()]
δ β is the area of β bin;
Ai=Yj-YK, Aj = YK- yi, Ak = yi -Yj.
Bi=Xj-Xk，Bj=Xk-Xi，Bk=Xi-Xj
X and y are the vertical and horizontal coordinates of the three vertices of β bin.
4. Revision and solution of mathematical model
Limited by observation data, the simulation calculation period is from June 1982 to February 1983, and the calculation time unit is months.
In the process of simulation calculation, for the river infiltration with linear seepage characteristics, firstly, the river is naturally divided according to the surface element that the river crosses, and the midpoint of the cutting line segment is taken as the infiltration point, which approximately replaces the linear infiltration; Then, according to the position of the equivalent seepage point in the panel and the distribution principle of the concentrated load of the triangular panel to each support point (the vertex of the triangular panel), it is added to the three vertices to calculate the seepage quantity. The exchange volume between production wells and other vertical points is treated in this way.
According to the difference of rainfall in different places, four rainfall zones are divided in the calculation, which are roughly divided into the middle and upper reaches of Wenhe River, the middle and upper reaches of Taohe River (including Nanchuan River), Songxi River and the lower reaches of Wenhe River. Due to the lag effect of vadose zone thickness in different areas, according to the comparison of the dynamic relationship between rainfall and groundwater level, the rainfall lag from drainage area, confluence area to recharge area is 1 ~ 3 months respectively.
Through the parameter adjustment of karst hydrogeology zoning, the groundwater flow field, 1 1 long groundwater level observation well and spring flow are simulated and calculated, and finally 33 parameter zoning (Table 8-5 and Figure 8- 13) and fitting curves as shown in Figure 8- 14 are obtained.
Table 8-5 Zoning Table of Hydrogeological Parameters of Niangziguan Spring Area Unit: m2/d
Second, the calculation of response matrix
In 1970s, from the point of view of groundwater management, in order to pursue the goal of the optimal development state of groundwater with unified flow field, people combined the groundwater numerical seepage model with optimization technology to form a groundwater system management model. In order to ensure that the relationship between water level constraint and output in the management model obeys the seepage equation of groundwater flow, predecessors proposed "nested method" (aguado, Remsen, 1974) and "response matrix method" (Haidari, 1982) to establish the relationship between them. The essence of water vulnerability assessment is the expression of the influence degree of water exploitation amount at each point in the system on the water quantity (or water level) at a specific target point, and it is a problem of calculating the relationship between water level (or flow reduction) at different points and water exploitation amount, so the response matrix method is used to link the two in this calculation. The seepage model of karst groundwater in spring area described by Formula (8-2) can be decomposed into the following two definite solutions:
Figure 8- 13 Zoning Diagram of Numerical Simulation Parameters of Karst Water System in Niangziguan Spring Area
Environmental problems and protection of karst groundwater in northern China
Figure 8- 14 Summary of numerical simulation curves of karst water system in Niangziguan spring area
and
Environmental problems and protection of karst groundwater in northern China
Where: q = ε+p; ε is rainfall infiltration, river leakage and vertical recharge of reservoir; P is the exploitation amount of karst groundwater in the system.
Equation (8-4) represents the natural water level field without groundwater exploitation in the management area under the given initial and boundary conditions; Equation (8-5) describes the variation field of water head under mining conditions, regardless of initial value and boundary value. Easy to prove:
Environmental problems and protection of karst groundwater in northern China
When the output P is controlled in a certain range that is not enough to cause the system parameters to change greatly, the hydrogeological characteristics described in Equation (8-5) are a stationary linear time-invariant system with homogeneous boundary and initial equation, and the water level drop can be characterized by the unit impulse response function of output, and the mathematical expression of the input-output relationship of the system can be obtained. According to the properties of linear time-invariant system, the groundwater level decline can be expressed by the convolution of its exploitation amount and the unit impulse response of the system:
Environmental problems and protection of karst groundwater in northern China
Where: Q(t) is the exploitation amount; B(t) is the unit impulse response function of the system; 5(t) is the decrease of depth caused by mining.
For (8-7) discretization, there are
Environmental problems and protection of karst groundwater in northern China
According to the superposition principle, for the linear system of formula (8-8), at the end of n period, the depth drop caused by I water sources at node K is
Environmental problems and protection of karst groundwater in northern China
Where: S(k, n) is the depth drop of K node at the end of n period; Q(i, j) is the exploitation amount of water source in J phase I; β(k, i, n-j+ 1) is the response matrix.
The physical meaning represented by the response matrix can be understood as follows: I water source pumps water at the unit pumping capacity in the first time period, and when pumping stops in the later time period, the remaining water level of K node drops at the end of each time period.
Given the unit pulse output, the water level response matrix of any point in the system at each time period can be obtained by solving Equation (8-5).
Three, the division of karst water system aquifer water source protection areas
Niangziguan spring is currently the water source of Yangquan tap water lifting project and Pingding county water lifting project. The planned maximum pumping capacity (excluding the closed Niangziguan Power Plant) is 5.09m3/s, and the local water consumption in Niangziguan Town is 0.5m3/s, so the total water consumption will reach 5.59m3/s (Table 8-6). According to the measured average flow of spring water after 2000, it is 6.452m3/s, and the residual flow after reaching the maximum pumping capacity is less than1.0m3/s. There is ecological water in Mianhe irrigation area downstream of the spring water and along the Yehe River in Hebei Province. The continuous attenuation of spring water flow will inevitably threaten the sustainable utilization of spring water, so maintaining the minimum spring water flow of 6.0 ~ 6.5m3/s should be an important goal of karst water protection. In order to achieve this goal, it is necessary to take Niangziguan spring flow as the protection core, divide the reduction of spring flow caused by pumping of karst aquifer in the system, that is, the difference in time and space of Niangziguan spring's sensitive response to pumping of karst aquifer in the system, and evaluate the vulnerability of spring flow to aquifer in the system.
Table 8-6 Investigation Statistics of the Present Situation of Niangziguan Spring Pumping Project
1. Water vulnerability assessment of aquifer
Through the response matrix based on the numerical seepage model of karst water in the system, the water level response values of other arbitrary points at any time can be obtained when the aquifer of the system is pumped at any point or multiple points. In plane space, it can be obtained by changing the number of nodes of numerical method (or the number of nodes that make up a bin), and in time, it can be obtained by adjusting the time step. For a given time step, it needs to be given in combination with the dynamic characteristics of karst water system.
The dynamic characteristics of karst water system in Niangziguan spring area are mainly restricted by atmospheric precipitation recharge, so we adopt the following conditions in the selection of vulnerability assessment period:
1) spring water discharge area. Compared with the direct extraction of spring water, it instantly decreases 100% in the extremely sensitive area of water quantity.
2) 1 year, where the spring water decreases by more than 80%. Due to the uneven distribution of precipitation in the system during the year, the spring flow and karst groundwater level show regular changes in a year.
3) Spring water in moderately sensitive waters is reduced to more than 80% in the longest continuous dry season. As far as water quantity is concerned, the influence of continuous drought years on the dynamic characteristics of karst water and related hydrogeological environmental problems is the most important. Therefore, the longest period when the continuous precipitation is lower than the multi-year average is regarded as the vulnerability assessment period. According to the dynamic analysis of precipitation series in the system, there have been two periods in history when the average precipitation for five consecutive years is lower than that for many years. Therefore, the five-year period is chosen as the evaluation period.
4) Remaining low sensitive area of water in the whole system. In the delineation of five protected areas of karst water system in Soignies River Management Area, Florida, USA, the karst aquifer area reaching the basin area in 20 years is determined as the protected area. There are some human factors in this definition, and it is more reasonable to extend it to the whole spring area.
2. Division of systematic karst aquifer water source protection areas
Maintaining the minimum flow of Niangziguan Spring above 6.0m3/s should be one of the important objectives of aquifer water protection in this system. In order to achieve this goal, it is necessary to make a general development plan of systematic karst water according to the results of vulnerability assessment. At present, there is almost no remaining available flow, so the exploitation of karst groundwater in the first and second protected areas, which have great influence on spring flow, should be strictly restricted. According to the fourth water-saving protection zone division scheme in Chapter VII (Figure 7-3), * * * divides four-level karst aquifer water-saving protection zones (Figure 8- 15), namely:
Figure 8- 15 Zoning Map of Karst Aquifer Water Vulnerability (Protected Area) of Karst Water System in Niangziguan Spring Area
1) Yolanda nature reserve area13.87km2.
2) The first-class protection zone of karst aquifer water source covers an area of 20.73km2
3) The area of the karst aquifer water secondary protection zone is 137.6438+02km2.
4) Quasi-protected area of system water, with an area of 7024.47km2