In order to simplify the well completion process, increase the water yield and reduce the cost, mixed wells are mostly used for groundwater exploitation in this area (Figure 5-2).
At present, the description of mixed wells in groundwater flow model is mostly handled by three-dimensional finite difference groundwater flow model MODFLOW software, which was introduced by the US Geological Survey in 1988 and is widely used in the world.
Figure 5-2 Mixed Pumping Well
MODFLOW suggests that the flow of multi-layer wells must be artificially distributed to each layer in some way ... Well flow should be distributed according to the hydraulic conductivity of each layer, that is, Qi/Qw=Ti/∑T( 1988,1996,2000). Where Qi and Qw are the flow rate and total (wellhead) flow rate of layer I respecTively, and ti and σT are the hydraulic conductivity and total hydraulic conductivity of layer I respectively.
In order to facilitate the discussion without losing its generality, we take a mixed well hole (Figure 5-2) passing through two aquifers as an example for discussion. In this way, the above formula can be expressed as Q 1/Q2=T 1/T2.
For this mixed well treatment method, MODFLOW has not given any theoretical analysis and explanation, and lacks theoretical basis, and its application is also inconsistent with reality (Chen Chongxi et al.,1998; Chen Chong and Jiao J J J, 1999). This is because:
1) The influence of aquifer drainage system number on borehole flow is not so simple. For example, if the lithology (permeability coefficient) changes near the mixed well, or even if the mixed well meets the lithologic lens (Figure 5-3), what impact will it have on the flow distribution of each layer? When the aquifer thickness changes, how to change the distribution of flow? MODFLOW can't answer these common practical questions.
Figure 5-3 The mixed well tube passes through the lithologic lens.
2) The parameters of aquifer affect the flow distribution of mixed wells, and the permeability coefficient is only one of the factors, and the elastic water supply (water storage coefficient) of aquifer should also play a role.
3) The distribution of well pipe flow is not only related to the distribution of hydrogeological parameters of aquifer, but also related to external boundary conditions, well diameter (effective well diameter), the position of suction pipe of pump and the interference of other pumping wells. MODFLOW only pre-allocates the flow of mixed wells according to the number of drainage systems, without considering the role of other factors.
4) The hydraulic conductivity of aquifer generally does not change with time, that is, the hydraulic conductivity ratio T 1/T2 will not change, but the flow ratio Q 1/Q2 is not constant.
5) According to the MODFLOW method, the flow rate Q 1/Q2 of the mixed pumping well is always constant. However, in the simulation process, especially in the prediction, mixed or non-mixed pumping wells can be added or shut down at any time. Under the interference of such well groups, will the original mixed well flow ratio remain unchanged? Obviously impossible.
6) From another perspective, the hybrid observation well is a special case of hybrid pumping well (Qw=0). For a two-layer mixed observation well, the water level in the hole (mixed water level) must be between two aquifers, that is, the mixed observation well can pump water from one aquifer (such as 1 > 0 aquifer) and the other aquifer (2 aquifer). Therefore, Q 1/Q20/T2 > 0. So, how can these two ratios be equal? Some people think that MODFLOW does not say that the above formula can be used to mix observations well. In this regard, it is easy to prove that the above argument is also true when the pumping capacity of the mixed pumping well is small enough to keep the water level in the mixed well between the heads of the upper and lower aquifers.
7) According to the ratio of hydraulic conductivity t of each layer, the flow of each layer is divided artificially in advance. This is lack of theoretical basis, because this method requires that the effective borehole diameter of each layer is equal and the hydraulic gradient at the borehole wall is equal everywhere. These two conditions cannot be controlled artificially and are unknown in advance.
Obviously, it is not appropriate to use the ratio of hydraulic conductivity to give the flow rate of each layer in advance. In essence, this method is not simulation, but "treatment", a kind of "treatment" without considering the mechanism.
"Preventing the simulation distortion of groundwater and improving the simulation effect" is the core task of hydrogeologists, and the mixed well, which is ubiquitous in groundwater flow system, is one of the main problems of simulation distortion at home and abroad, which should be paid enough attention to.
The mixed wells in the study area include pumping wells and observation wells. In the three-dimensional flow field, even if it belongs to homogeneous aquifer (no multi-layer aquifer system), conventional (theoretically non-point filter tube) pumping wells and observation wells belong to mixed well holes. Because the water heads in different depths in the filter tube are not equal, the water in the filter tube will flow vertically, even in the observation well without pumping water. This is the essence of the response of groundwater flow to mixed wellbore.
In this study, the "seepage-pipe flow coupling model" (Chen Chongxi et al., 1992, 1996) is used to describe the mixed drilling, the "seepage" is used to describe the movement of groundwater, and the "pipe flow" is used to describe the water flow in the drilling, which solves the simulation problem of mixed drilling and greatly improves the simulation effect of the model. The basic idea of "seepage-pipe flow coupling model" is as follows (Figure 5-4):
Figure 5-4 Calculation Block Diagram
1) The pipe-flow mixed well is regarded as a "cylindrical lens" with high permeability, and the aquifer exchanges water strongly through this cylinder (wellbore) with high hydraulic conductivity, so the mixed pumping problem can be regarded as a special "overflow system".
2) Find the "permeability coefficient" of the "cylindrical lens" of the overflow system. It should be emphasized that the "permeability coefficient" of "cylindrical lens" is not constant, but changes with the change of flow pattern.
The specific calculation method is as follows:
According to the knowledge of fluid mechanics, when the pipe flow is laminar, the head loss can be calculated by Darcy-Weisbach equation:
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: Δ h is the head loss; F is the coefficient of friction; L is the pipe length, m; D is the inner diameter of the pipe, m; U is the average velocity in the pipeline, m/d; G is the acceleration of gravity, m/d2.
When the pipe flow is laminar
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: re is Reynolds number.
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: ν is the kinematic viscosity of the fluid; μ is the dynamic viscosity of the fluid; ρ is the density of the fluid.
Substitute the formulas (5-2) to (5-4) into the formula (5- 1) to obtain.
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: j is the hydraulic gradient; Gamma is that severity of the fluid.
Write the above formula in the form of seepage velocity. For pipe flow, the void ratio n= 1 and the seepage velocity v=nu=u, so there is
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Comparing this formula with Darcy's law v=KJ, the expression of equivalent permeability coefficient KL of pipe flow in laminar state can be obtained.
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
This formula was obtained by Chen Chongxi (1966) and J.Bear( 1972) respectively.
According to the knowledge of fluid mechanics, when Reynolds number is greater than 100000, the friction coefficient f has nothing to do with Reynolds number Re, but depends on the relative roughness of the inner wall of the pipeline (), and the head loss δ h is proportional to the square u2 of the flow. In the range of 3000 < Re < 10000, there is a section of f=, that is, δH∝u 1.75. Under turbulent conditions, in addition to the above two sections, there are two transition zones, which are divided into four sections. So turbulence is a complicated problem. In order to solve this problem, the concept of equivalent permeability coefficient KN is put forward here.
When the pipe flow is in a turbulent state, the equation (5- 1) can be rewritten as
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
If the equivalent permeability coefficient of pipe flow in turbulent state is defined
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
rule
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
This formula also has the form of Darcy's law.
Originally, the motion law of turbulence is different from Darcy's law, and different flow areas have different forms. After introducing the equivalent permeability coefficient, the motion laws of five flow zones (1 laminar flow zone and four turbulent flow zones) are unified into Darcy's law, which is consistent with the groundwater seepage law. that is
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where Ke is the equivalent permeability coefficient of the coupling model of seepage and pipe flow.
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
It should be noted that the equivalent permeability coefficient KN under turbulent conditions is related to the flow velocity v and the friction coefficient f, which is a variable varying with Reynolds number. In addition, the determination of Reynolds number Re is related to flow velocity V, which in turn depends on friction coefficient F, which in turn depends on Reynolds number Re. Therefore, these three factors must be determined by iterative method, and then the equivalent permeability coefficient Ke must be determined. The specific calculation process is shown in Figure 5-4.
Secondly, the simulation of mixed observation well water level.
For various reasons, there are almost no areas without mixed observation wells. What is a mixed observation well, how to form the water level (mixed water level) in the mixed observation well, how to determine the distribution of initial water head according to the mixed water level, and how to calculate the hydrogeological parameters of each layer by using the mixed water level. , is one of the important problems to prevent simulation distortion and improve simulation (Chen Chongxi, 2003). However, MODFLOW(McDonald et al., 1988) has not analyzed and discussed the mixed observation well, and there is no model of the mixed observation well. It is wrong to simply give up the information of the mixed observation well, because the "mixed observation well" has only zero orifice flow, but there are "pumping" and "water injection" in the mixed observation well.
D sokol (1963) used Thiem's "radius of influence" model to prove that the water level of the mixed observation well is equal to the algebraic average value of the water head of each layer weighted by hydraulic conductivity. There are several problems in this paper: ①Thiem's "influence radius" model can't form a stable flow (Chen Chongxi, 1966, 1975,1983); (2) (2) Sokol implicitly assumes that, although the two aquifers have different permeability coefficients and different flow rates, the "influence radius" is equal; (3) (3) Sokol also implies the assumption that the well pipe resistance can be ignored, that is, the head of the well pipe is equal everywhere. The conclusions drawn from these assumptions are not credible.
Hantushi (196 1) and бочевер et al. (former Soviet Union) (196 1) all think that the water head in the observation well is in the confined aquifer under study.
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where, l' and d' are the elevations of the top and bottom points of the observation well filter tube. This book is called hantush-бочевер equation.
For the study of three-dimensional phreatic flow, Neuman( 1972) also thinks that the head drop depth in the observation well can be regarded as the average value of the drop depth of each point in the filter tube, that is,
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: r is the distance from the observation well to the pumping well; Z 1 and z2 are the elevations of the top and bottom points of the filter tube in the observation well.
The opinions of the above scholars and the determination method of water level in observation wells have been used in hydrogeology now, and there has been no objection for more than 40 years. However, the above methods for determining the water level in observation wells lack a basic analysis of the formation mechanism. Simply analyze, the vertical observation well in three-dimensional flow will have vertical flow due to the unequal water head in the observation well filter tube. The flow of water will lead to the change of water head, and it will also cause the movement of groundwater around the borehole and the redistribution of water head. Therefore, different from the manometer observation well, generally speaking, the observation well is essentially a part of water supply (pumping) filter tube and another part of water drainage (water injection) filter tube, that is, the observation well is not simply reflecting the head of aquifer groundwater, but a well hole with the functions of "pumping" and "water injection", but the flow at the orifice is zero, and the absolute value of "pumping" in the well tube is equal to (if the tiny wellbore storage is ignored)
Obviously, the above-mentioned "comprehensive average water level" of Hantushi, бочевер and Newman does not involve the mechanism of water flow, and it is a pure mathematical method lacking physical foundation. The author uses "seepage-pipe flow coupling model" to simulate the water level (mixed water level) in mixed observation wells. The conditions of this example are:
A homogeneous confined aquifer with equal thickness and infinite horizontal extension is pumped by a hemispherical well at the top at a constant flow rate. The pumping duration is t= 10.02d, and there are complete observation wells (i.e. mixed observation wells) with different diameters at different radial distances r, and the head distribution and flow change in these observation wells are simulated. Aquifer thickness M= 100m, permeability coefficient K= 100m/d, unit elastic specific yield μ s = 0.00001m-0/,pumping flow q = 36000m3/d ... The simulation results are shown in Figure 5-5 and Figure 5-6. As can be seen from Figure 5-5, in the observation well with aperture ≤0. 1m, the head values at different depths are obviously different, which is also very different from the integral average of hantush-бочевер.
Fig. 5-5 Depth distribution of head drop of observation wells with different radial distances r and values sga calculated by hantush-бочевер formula.
Figure 5-6 Flow Distribution of Observation Wells with Different Diameter Spacing
Because the pumping well (point sink) is located at the top of the confined aquifer, the head drop depth S at the top is greater than that at the bottom at all radial distances R, that is, the head H at the bottom is greater than that at the top, so the water in the observation well flows from bottom to top. Based on the principle of water flow continuity, the observation well at the bottom should absorb water (pump water) and the observation well at the top should drain water (inject water). Figure 5-6 shows the flow distribution in the observation well. Figure 5-6a shows that in an observation well with a diameter of 0.20m and a diameter of 2 1. 19m, it is "pumping" below z=35m and "water injection" above 35m; The maximum vertical displacement of observation well exceeds 500m3/d, which is what we say, observation well is not an observation hole!
Based on the above analysis, this study uses the "seepage-pipe flow coupling model" to simulate the water level in the mixed observation well, so that the mixed water level becomes useful information, that is, it is used for initial head distribution and fitting parameters.
3. Dynamic simulation of spring flow
Spring is one of the main forms of converting groundwater into surface water. The emergence of different types of springs and their flow dynamics are key issues for hydrogeologists to analyze, because these data provide extremely important information for understanding hydrogeological conditions. Spring is the main simulation element in the numerical model; Spring flow prediction is an indispensable part of groundwater resources evaluation and management.
The U.S. Geological Survey (200 1) suggested using the Drain module of MODFLOW software to calculate the spring flow iteratively. In essence, the spring water was described as the second boundary, that is, the water head distribution was first calculated without considering the existence of the spring water, and then the spring flow was calculated by multiplying the difference between the groundwater head at the grid point where the spring is located and the elevation of the spring mouth by a certain proportional coefficient (lacking physical significance), and then the water head distribution at the grid point (the second boundary) was recalculated. When the groundwater level is lower than the elevation of the spring mouth, the flow is zero; When the groundwater level is higher than the elevation of the spring mouth, the spring flow is directly proportional to the head difference, and its proportional coefficient must be calculated by flow fitting.
In the previous numerical models in China, the description of spring flow is similar to the input (given) flow of pumping wells. How to deal with this and how to predict it? Spring flow is unknown. In recent years, there are also iterative algorithms similar to those proposed by MODFLOW in China. For example, some research reports describe Shanxi Liulinquan (actually a monoclinal self-flowing inclined rising spring) on the east side of Ordos syncline as a plane two-dimensional flow model, and then calculate the flow of the spring with the formula of complete well under pressure, that is, Q=qs. There are several problems here: ① Liulinquan is an ascending spring, and the deeper the hole near the spring, the higher the groundwater head, which is the characteristic of monoclinic gravity flow inclined ascending spring, which has obvious three-dimensional flow characteristics; (2) It is not appropriate to calculate the spring flow with the formula of complete well under pressure, especially Liu Linquan's flow; ③ The proportional coefficient q is not constant at different times (under different conditions).
"Spring water must be placed in hydrogeological body to study its genetic type, in order to simulate it correctly". Different from MODFLOW, the method proposed in this model takes the elevation of the spring mouth as the first boundary condition (when the spring flow exists), and the grid points of the spring mouth, the surrounding grid points and the lower grid points of the same layer conform to Darcy's law respectively (if the groundwater flow belongs to linear flow, for example, Chen Chongxi et al. (1995), Chen Chongxi et al. (1996), Chen Chongxi et al. Or use non-Darcy's law, such as Chen Chongxi (1995), Cheng Jianmei (1998) and Cheng (1998), and establish the relationship with the principle of water balance. The whole calculation method is based on the flow mechanism, so the spring flow is generated directly through the operation of the model, and there is no need for iterative solution (in some cases, iteration and union are needed). Because of this, it is possible to take the spring flow as the fitting object in the model identification stage. This is very important.
Specifically (Figure 5-7), according to Darcy's law, spring flow Q is the sum of horizontal flow and vertical flow. If the horizontal flow is dominant, it is a downward spring, otherwise it is an upward spring. Specific can be calculated by the following formula:
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: q is the spring water flow; Kzsp is the vertical permeability coefficient of spring flow area; Hb is the head value of the node under the spring mouth; Zsp is the exposed height of the spring; Dzsp is the seepage length (aquifer thickness) in the process of spring flow; S is the control area of the spring mouth node; CA is the correction coefficient of the spring neck; Cs is the spring diameter correction coefficient; Tij and Tik are the average hydraulic conductivity of flow sections ij and ik; ,,, is the length of the line segment; Hi, hj and hk are the head values of nodes I, J and K; E is the equalization unit around node I, such as ipbq quadrilateral.
Figure 5-7 Schematic Diagram of Spring Flow Calculation
In the previous numerical simulation of groundwater resources evaluation, the spring flow rate was not taken as a fitting factor. In fact, it is impossible for practices like MODFLOW to take the spring flow as a fitting factor, because they need to solve a coefficient that lacks physical meaning when identifying the model.
The model selects 3 typical point springs and 22 line springs (spring ditches) to simulate their water flow dynamics, and takes the measured flow as an important fitting object to obtain the hydrogeological parameters of the water-bearing system and predict the water flow dynamics of the springs. This will improve the simulation of the model.
Fourthly, the hysteresis description of atmospheric precipitation and surface water on phreatic recharge.
Groundwater infiltration recharge such as precipitation, surface water and canal system is an important source of groundwater recharge. If not handled properly, the model will be distorted. Therefore, it is of great significance to correctly describe the lag of infiltration recharge to improve the simulation of the model.
At present, in the evaluation of groundwater resources, the following two methods are mainly used to calculate the infiltration recharge of precipitation, surface water and canal system:
1) permeability coefficient method without considering lag. This method can generally be used to describe the infiltration recharge of groundwater in small buried depth areas, but in some areas, the buried depth of groundwater is large, and if the lag characteristics of infiltration recharge such as rainfall and surface water are not considered, the model will be distorted.
2) Translation lag permeability coefficient method. Some research reports regard the time difference between the peak time of groundwater level and the peak time of rainfall as a unified lag time. For example, if the time difference is 3 months, the rainfall of 65438+ 10 will permeate and replenish in April, and the rainfall of February will permeate and replenish in May. There is no rain in March, and there is no rain infiltration in June. It is also obviously unreasonable to deal with delayed replenishment in this way. Because the peak time of groundwater level is the sum of all the factors that affect the groundwater level dynamics, besides rainfall recharge, there are also the results of groundwater evaporation, irrigation infiltration, overflow, groundwater exploitation and other factors that cause the imbalance of groundwater flow. In particular, the artificial factor of groundwater exploitation has become the main control factor of groundwater dynamics in this area, which completely covers the natural dynamics of groundwater. It is not appropriate to study the lag of rainfall infiltration recharge by this method. Even if this area is not developed at all, or even the influence of other factors other than rainfall on groundwater dynamics is ignored, this method has obvious problems only from the influence of rainfall infiltration recharge on groundwater level dynamics. The rainfall process is intermittent and discontinuous, but its recharge to groundwater is continuous. This is because groundwater flow, especially unsaturated flow, has obvious regulation and storage function, which makes the recharge lag and continuous. Although there is only a few tens of days of rain in a year, there is no rain during the rest time, but in areas with deep groundwater level, recharge works every day (the situation of strong evaporation in shallow groundwater level area is more complicated), and it is impossible to describe its lag with a time difference translation.
The above-mentioned "permeability coefficient method" has certain applicability in areas with small groundwater depth. However, due to the deep groundwater depth in some areas of this area, such as the highest groundwater level at the top of the alluvial fan in Changma is 290m m, the monthly rainfall multiplied by the permeability coefficient can not completely replenish the groundwater layer below it. Therefore, ignoring the lag characteristics of rainfall recharge will lead to the distortion of the model.
In this study, the "lag weight coefficient method of rainfall infiltration recharge" (Chen Chongxi et al., 199 1, 1998) is used to describe rainfall infiltration recharge. This method can basically reflect the objective reality and is very practical.
According to the theory of groundwater unsaturated flow, the distribution curve of groundwater recharge during a rainfall is shown in Figure 5-8. This is a single peak curve. The shape of the curve depends on the buried depth of groundwater level and the lithology of vadose zone, and of course it is also related to the distribution of initial water content in vadose zone. We will identify the latter factor and discuss the role of the former.
Figure 5-8 Distribution Curve of Initial Rainfall Replenishment
Based on the unsaturated flow theory of groundwater in vadose zone, the distribution characteristics of rainfall recharge intensity in a certain period (such as month and ten days) under different groundwater depth and lithology conditions are studied. In order to cooperate with the discussion of groundwater resources evaluation numerical method, the continuous characteristic curve is changed to the form of discrete points (because the numerical method is discrete in time). Without losing generality, we take the time period as 1 month.
For the case that the buried depth of groundwater level is large and/or the permeability coefficient of vadose zone is small, the distribution of infiltration intensity is a unimodal function, and it may be a monotonic decreasing function when the buried depth of groundwater level is small (Figure 5-9). If there is a certain rainfall RA(mm) in 10, the sum of rainfall infiltration recharge in this month and subsequent months is the total recharge formed by rainfall in this period. The ratio of this value ∑REk to 10 monthly rainfall ra is the (total) infiltration recharge coefficient α. The RAtio ra of monthly infiltration to 10 rainfall is the infiltration recharge coefficient αk of monthly 10 rainfall. that is
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Figure 5-9 Duration Distribution Curve of Rainfall Infiltration Replenishment under Different Buried Depth (Lithology)
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: RA is 10 monthly rainfall, mm; REk is the rainfall infiltration recharge in k month 10 month, mm; RE is (total) supplementary quantity, mm; αk is the infiltration recharge coefficient of rainfall in k months 10; α is the (total) infiltration recharge coefficient.
If the rainfall in the last month of 10 in Figure 5-9 is 65,438+000 mm/month, the ordinate of the figure is converted into the monthly permeability coefficient α k (k = 0,65,438+0,2,3, ...) and expressed as a percentage.
The monthly infiltration recharge coefficient series can well reflect the characteristics of rainfall infiltration recharge lag. In order to make the description of the lagging feature array more universal, we define:
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
therefore
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
This formula shows that ωk satisfies the definition of weight, so we call ω k (k = 0, 1, 2, ...) as the coefficient (group) of delayed replenishment weight.
The advantage of introducing lag recharge weight coefficient (group) is that the same lag recharge weight coefficient can be used to calculate the monthly recharge quantity in different infiltration coefficient areas (different infiltration recharge quantities in the same month in different areas), which means that the lag recharge weight coefficient is more universal. So our task is to find a characteristic function that can satisfy the above conditions, namely:
1) They are unimodal functions or monotonically decreasing functions, and the peak value (or maximum value) decreases with the delay of the occurrence time of the peak value (or maximum value);
2) The sum of the weight coefficients is1;
3) It has strong adaptability and can well match the delayed infiltration recharge under different groundwater levels and different lithology conditions.
After repeated analysis and research, we decided to use the following discrete function to express the distribution of recharge weight coefficient of delayed infiltration.
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where ω(j-i) is the weight coefficient of monthly rainfall in I and monthly recharge in J; H is the depth of groundwater level, m; B is the empirical coefficient related to factors such as lithology of vadose zone, m; It can be called the interval coefficient of buried depth. This coefficient is the increasing function of permeability coefficient of vadose zone.
In this formula, H is known, and there is only one empirical coefficient B related to lithology and other factors. This coefficient can be determined by model identification, so that the recharge weight coefficient of delayed infiltration can be calculated by the above formula. The (total) infiltration coefficient of each district is also determined by model identification, from which the recharge RE(j-i) of this month and subsequent months can be calculated under the action of rainfall RA(i) in a certain month, namely
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
So far, we have analyzed the contribution of rainfall in a month to infiltration recharge in that month and subsequent months. To simplify the calculation, for the practical problem-the recharge of J month by multi-month rainfall, we approximately use the algebraic sum of the contribution of I month rainfall to the phreatic recharge of this month and previous months to calculate the total infiltration recharge obtained by J month, namely
Investigation and Evaluation on Rational Development and Utilization of Groundwater Resources in Shule River Basin of Hexi Corridor
Where: RES(j) is the total rainfall infiltration recharge in J month and previous months, mm; RE(j-i) is the rainfall infiltration recharge from J month to J (j > I) month, mm; α is the infiltration recharge coefficient; Ω (j > I) is the lag recharge right coefficient of I month rainfall to J (j > I) month; RA(i) is the rainfall in January, mm.
According to formula (5- 17), 75 weight coefficient curves with delayed rainfall infiltration recharge are obtained, corresponding to 75 buried depth intervals of groundwater level, and the buried depth difference Bm of each interval (that is, the serial number of the weight coefficient curves changes with each increase of buried depth Bm). In Figure 5- 10, only 12 typical weight coefficient curves are drawn, and between the two curves, there are still 5-7 curves not drawn.
The proposal and use of the lag weight coefficient method of rainfall recharge not only greatly improves the simulation effect, but also is very convenient to use. The lag mechanism of surface water infiltration recharge is similar to that of rainfall infiltration, and it is also described by the above method.
Fig. 5- 10 rainfall infiltration recharge delay weight coefficient curve
Parameter division of verb (abbreviation of verb)
In this model, the water-bearing medium is thick and vertically divided into six simulation layers. Even if the same lithology (such as loam) is at different burial depths, its hydrogeological parameters are not equal because of its different consolidation degree. In this case, how to give the spatial distribution of initial parameter estimation is another difficult problem. According to the conventional method, the same lithology should be divided into different parameter zones according to different horizons (different buried depths), so there are too many parameter zones, which brings a lot of workload to parameter adjustment.
In order to solve this problem, the model adopts the following methods: the Kh value of the same lithology at different burial depths is taken as the same partition, but its parameter value decays exponentially with burial depth. In this way, each parameter partition only needs to be described by 1 horizontal permeability coefficient Kh0 and 1 attenuation coefficient γ, that is, Kh=Kh0×e-γh(h is the buried depth value). The spatial distribution of other parameters such as Kz value, μd value and μs value is also given in the same way. In this way, the number of model parameters is greatly reduced, and the workload of model parameter adjustment is correspondingly reduced.
6. Drilling simulation
Artesian wells generally exist in northwest inland basins, depressions in front of alluvial fans, around large karst springs and deep confined aquifers in plain areas. Up to now, there is still a large area of groundwater artesian zone in China (although the area of artesian zone has decreased due to a large number of groundwater exploitation). Therefore, how to describe artesian well in the model is an important hydrogeological problem. However, there is no module to simulate the flow of artesian wells in the MODFLOW software launched by the US Geological Survey, which seems to be not negligence, but a certain degree of difficulty. Conventional groundwater flow model can't solve this kind of problem. So far, the artesian well in the numerical model has given the artesian flow helplessly. However, artesian wells are different from pumping wells, in which people can generally give the pumping flow artificially, and the water level in the well is its response; When the water level in the well drops to the bottom of the well, people can't pump water freely (but most of them contain numerical models of pumping wells, which can't objectively understand that "the water level in the well is the constraint of pumping quantity"). However, artesian well can't be given artificial gravity flow, because gravity flow is the response of artesian well to surrounding environment (meteorology, hydrology, groundwater exploitation and other factors), and the flow of artesian well belongs to the category of prediction, which is output information, not input information.
This project adopts "seepage-pipe flow coupling model" to simulate artesian well. As long as the known head boundary is set at the wellhead of artesian well, that is, the head is equal to the wellhead elevation, the flow of artesian well can be simulated. This method has been successfully applied to "numerical simulation of karst water resources in eastern Weibei, Shaanxi Province" (Chen Chongxi et al., 2003).
Seven, the formation of the initial head
Initial head distribution is an indispensable condition for numerical simulation of unstable groundwater flow. Usually, the initial head value of each node (grid point) is obtained by interpolating the water level of a certain number of observation wells. However, when the study area is large, there are few observation wells, the observation wells are located in different horizons, and there are mixed holes, it is difficult to get the initial head distribution of each simulation layer by interpolation.
This project adopts "parameter-initial head iteration method" (Chen Chongxi, Lin Min, etc. , 199 1; Chen Chongxi, Pei Shunping, 200 1), determine the initial head distribution of each layer. The specific treatment methods are as follows: the initial simulation time is July 2003, and the simulation starts from July 5438+0, 2006. With the initial estimated values of hydrogeological parameters, the simulated monthly head value is fitted with the measured value in July 2003, and the initial head of the simulation period (July 2003) is roughly determined; Then the initial head is used to calculate the parameters according to the conventional method, and the iterative solution is obtained. When the fitting reaches a certain level, it is necessary to comprehensively test whether the model conforms to the basic laws, otherwise it will be readjusted again until it meets the requirements. Although the above methods have a lot of work, the effect is good.