The calculation area is located at the southern foot of Tianshan Mountain and the northern edge of Tarim Basin, starting from Longkou, Gan Wei in the north, flooding plain of Tarim River in the south, Santamu Farm in Xinhe County in the west and Halahatang Township in Kuqa County in the east. Geographical coordinates: east longitude 8215' ~ 8315', north latitude 40 45' ~ 4145', covering an area of about 5530km2. The administrative divisions belong to Kuqa County, Xinhe County and Shaya County in Aksu Prefecture, Xinjiang (Figure 5- 1).
Second, the conceptual model of hydrogeology
(1) model boundary and internal structure
Quaternary pore water-bearing system has multi-layer structure, mainly distributed in Gan Wei River basin. In this study, a three-dimensional flow model with six-layer structure is established to describe it.
On the plane, the northwest boundary of the model is missing the boundary between Neogene bedrock and Quaternary sediments before Letage Mountain, which may be a zero flux boundary. Due to the lack of data in the south of Tarim River, the southern boundary is taken at the middle line of Tarim River as the first boundary. Other lateral boundaries (northeast boundary, east boundary and west boundary) are all extensions of Quaternary water-bearing system, and can only be used as artificial boundaries: the streamline is determined according to the initial groundwater level contour map, the boundary delineated along the streamline is the second boundary, and the boundary delineated along the isobath in some sections is the first boundary (Figure 5-2).
The specific boundary conditions are determined by model identification. The total area of all models is about 5530km2.
Vertically, the top boundary of the model is the phreatic surface, except the Tarim River as the first boundary, the rest are the infiltration recharge boundary or phreatic evaporation discharge boundary. The bottom boundary is the bottom of the third confined aquifer, which is a zero flow boundary.
In order to connect with the stratification of multi-layer structure plain area (phreatic aquifer, three confined aquifers and two weak permeable layers between them), the single phreatic aquifer in front of the mountain is also divided into six corresponding layers vertically. In this way, the water-bearing medium in the whole region is generalized into a six-layer structure, and a unified mathematical description is carried out.
(2) Treatment of source and sink projects
1. Treatment of mixed well
The "seepage-pipe flow coupling model" is used to describe the mixed well: "seepage" describes the movement of groundwater and "pipe flow" describes the water flow in the wellbore. This method solves the simulation problem of mixed wells and greatly improves the simulation of the model.
2. Immersion evaporation treatment
According to the law of phreatic water evaporation and four empirical formulas summarized in Xinhe balance test site, the problem of phreatic water evaporation is dealt with:
(1) Arbilla's Nove formula
ε=ε0( 1-D/D0)b
Figure 5- 1 Calculation of regional traffic location map
Figure 5-2 Calculation Area Range Diagram
② Power function formula
ε=ε0 aD-b
③ Exponential formula
ε=ε0 ae-bD
④ Tsinghua University empirical formula.
ε=εmax( 1-e-bε0/εmax)
Where ε is the phreatic evaporation intensity; ε0 is the evaporation intensity of water surface; εmax is the ultimate evaporation intensity related to the diving depth (that is, the maximum diving evaporation intensity that may occur at this depth).
ε max = f (d) d is the diving depth; D0 is the limit depth of diving (related to soil quality); A and b are empirical coefficients (related to soil quality).
According to different conditions of soil quality and groundwater depth, appropriate empirical formulas are selected for simulation improvement: Formulas ② and ④ are more suitable for fine-grained soil, but only when the groundwater depth is more than 0.2m; For coarse-grained soil, formula ① or ③ should be selected.
3. Treatment of infiltration recharge diving such as rainfall and surface water.
We use the "lag weight coefficient method of rainfall infiltration recharge" to describe rainfall infiltration recharge. This method can not only reflect the actual situation of delayed rainfall infiltration recharge, but also has certain practicability.
The lag mechanism of surface water infiltration recharge is similar to that of rainfall infiltration, and it is also described by the above method.
Third, the mathematical model.
According to the hydrogeological conceptual model of Gan Wei River basin, a three-dimensional unsteady groundwater flow mathematical model of mixed wells is established, which is described as follows:
Groundwater exploration in Tarim basin
Where: h is the head function (m) of aquifer or weakly permeable layer; H0 is the initial head function (m) of the calculation area; H 1 is the known head function (m) of the first boundary of the calculation area; Kh and Kz are horizontal and vertical permeability coefficients (m/d) of aquifer or weak permeable layer; Qw and Vw are the production of the production well and the volume of the borehole working section; W is the algebraic sum (m/d) of the infiltration recharge intensity of atmospheric rainfall, rivers, reservoirs, canal systems and field irrigation and the evaporation intensity of groundwater level; μs is the unit water storage coefficient of aquifer or weakly permeable layer (1/m); μd is the gravity specific yield of phreatic aquifer; V is the known seepage velocity function (m/d) of the second boundary of aquifer; B 1 is the first boundary of the calculation area; B2 is the second boundary of the calculation area; D is the distribution range of the calculation area.
Fourthly, modify the mathematical model.
The arbitrary polygon mesh finite difference method is used to solve the problem. Each simulation layer in the calculation area is divided into 662 nodes and 1236 triangular elements on the plane (Figure 5-3). The summary points of the six simulation layers are 662× 6 = 3972, and the total number of units is 1236× 6 = 74 16.
The starting and ending time of long-term observation data of groundwater level in the calculation area is 1 1 in 2000 to1/in 2006. Take this period as the model identification period, and take the time step Δ t = 3D (Δ t starts to take 30 days, and is changed to 3d after running the model due to the strong nonlinear problem of well mixing).
The simulation layer is divided into six layers, and each layer is divided into 14 parameter partitions.
The correction results are as follows.
(1) Hydrogeological parameters
The corrected model has six layers, each layer is divided into 14 parameter partitions (see Figure 5-4), and the parameter values of each parameter partition are shown in Table 5- 1.
(2) Groundwater resources
Considering the periodicity of groundwater dynamics, the first anniversary (1 February: 20001day to 20061day, 65438+30 days in 2006) is selected as the balance period of water balance calculation. See Table 5-2, Table 5-3 and Table 5-3 for the specific calculation results.
From the water balance (Table 5-2), it can be seen that the groundwater recharge in the calculation area mainly comes from the infiltration of surface water, and the sum of the recharge from rivers, reservoirs, canal systems and field irrigation accounts for 96.3% of the total recharge, while the recharge from rainfall infiltration to groundwater in this area is very limited (less than1.5%). The exploitation of groundwater in the calculation area only accounts for about 1% of the total emissions, while phreatic evaporation consumes a lot of water resources (97.9%), of which ineffective evaporation accounts for 80.25% of the total emissions.
Figure 5-3 Plane sectional view
Table 5- 1 Statistical Table of Parameter Values of Hydrogeological Parameters Zoning
Figure 5-4 Parameter Zoning Diagram of Calculation Area 1 ~ 6 Layer
Table 5-2 Calculation results of water balance in the calculation area from 2000 1 October1to 20061October 301.
Table 5-3 Summary of Calculation Data of Groundwater Model Replenishment Project
Table 5-4 Summary of Calculation Data of Groundwater Model Discharge Project
Note: Field irrigation infiltration recharge, rainfall infiltration recharge and evaporation are related to groundwater depth and lithology, which are constantly changing in the model; The values in this table are the average values of the vehicle identification period from June 2006 1 65438+1October1to June 30, 2006.
(3) Error analysis
Analysis of fitting error of 1. deceleration field
Finally, the fitting curve between the simulated water level and the measured water level in the observation well is shown in Figure 5-5. Among them, the solid big dot represents the measured value of the observation well water level (one value per month), and the hollow small dot represents the simulated head value (one value every three days).
In order to illustrate the fitting of water head, table 5-5 counts the fitting error of observation well head, and tables 5-6 and 5-7 count the fitting error of observation well head in each aquifer. The statistical data is the absolute value of the difference between the simulated head and the measured head of each fitting observation well.
Error statistics show that δ H ≤ 0.5m accounts for 62. 16% of the total contrast, and δ H ≤ 1.0m accounts for 86.63%. Generally speaking, the fitting effect is good under the current data. The water level of some observation wells is affected by random factors, and the fitting error is large.
2. Error analysis of gradient field fitting
Figure 5-6 shows the contour fitting diagrams of Layer 2 and Layer 4 in June 200 1+0 1 year (both drawn by Kriging interpolation method). In the figure, the solid line is the measured constant head (m), and the dotted line is the simulated constant head (m).
Verb (abbreviation of verb) verifies mathematical model
Except for the time period determined by the model (165438+2000 10 to165438+2000 10), there are no statistical data of groundwater exploitation and groundwater level dynamics in other time periods in the calculation area, which cannot be used. However, in the model identification stage, the data used are all measured data, and the time series of the data is long (1a), so the observation well is well controlled in the plane, and the fitting results are good whether it is deceleration field or gradient field, and the identified model is reliable and can be used for prediction.
Prediction of intransitive verbs
(A) the determination of the forecast scheme
The first scheme: maintain the existing layout of mining wells (as shown in Figure 5-7) and groundwater exploitation. There are 32 wells in total, and the total production is 0.1537×108m3/a. ..
The second scheme: based on the 32 production wells in the first scheme, 70 grid production wells (surface wells) are added according to the planned water source layout (see Figure 5-8) and production capacity, with a total output of 3.57×108m3/a.
Figure 5-5 Fitting curve of water level of observation well
Table 5-5 Absolute Error between Simulated Head and Measured Head of Observation Well
Table 5-6 Absolute fitting error of observation wellhead in the second layer
Table 5-7 Absolute fitting error of observation wellhead in Layer 4
Scheme 3: design a large-scale exploitation of the second floor groundwater in the phreatic evaporation area, and intercept the ineffective evaporation of groundwater by exploiting groundwater (see Figure 5-9). The mining intensity distribution ε [unit:104m3/(km2 a)] is obtained by practical trial and error method. The premise of calculating mining intensity by trial-and-error method is that the groundwater depth should meet the demand of ecological water use, that is, the groundwater depth in vegetation area should not be greater than 4.5m under mining conditions (the empirical value of groundwater depth in Gan Wei basin is 4.0 ~ 4.5 m due to lack of data), so as to ensure that plant roots can absorb necessary soil moisture.
The quality of groundwater in the desert area in the southwest of Xinhe County and some areas in the southeast corner of Shaya County is poor, the TDS of confined water is generally greater than 3g/L, and the TDS of phreatic water is higher, which is not suitable for mining. The northwest boundary is a zero-flow boundary, so drilling wells is not allowed in the bed of Tarim River, and groundwater in other areas can be developed and utilized.
(2) The boundary conditions and the treatment of source and sink items in the forecast period.
In the forecast scheme, the annual dynamic changes of rainfall, evaporation, river and reservoir infiltration, canal system water diversion and various infiltration coefficients are consistent with the model identification stage.
As a basin-level numerical model, the planned water source in the second scheme can only be treated as a surface well, and the production of each water source is distributed to all grid points within the water source range, and the production of each grid point is distributed in proportion to the area.
Fig. 5-6 Measured and simulated isohead diagrams of layers 2 (a) and 4 (b) of 2001165438.
Figure 5-7 Location Map of Production Well in Scheme I
Figure 5-8 Water Source Distribution Map of the Second Plan
Figure 5-9 Isogram of Mining Intensity Distribution in Scheme 3
Fig. 5- 10 Scheme I 2011+0/isohead diagram of each aquifer.
(3) Forecast results
The forecast time of the first and second schemes is 10 year (forecast to 201year1month), and the forecast time of the third scheme is 20 years (forecast to 20021year).
Scheme 1: it is predicted to be 20 1 1, and the isobar of each aquifer is shown in Figure 5- 10.
Under the current mining conditions, the water level of observation wells will remain basically unchanged in the next 65,438+00 years, and the water level of some wells will increase slightly, such as B2 and B45 wells, which may be caused by unreasonable recharge factors (such as canal flow distribution) and discharge factors (such as water surface evaporation intensity) in the forecast period.
See Table 5-8 for the multi-year average groundwater balance of 10 predicted by this plan.
Table 5-8 Annual Average Water Balance of 10 Predicted by Scheme I
Scheme 2: It is predicted that by June of 20 1 year, the isobar of each aquifer is shown in Figure 5- 1 1.
Fig. 5- 1 1 Scheme II 2011+/isohead diagram of aquifer.
This scheme increases the amount of groundwater exploitation. For observation wells far away from water sources, increasing groundwater exploitation has little effect on the head dynamics of observation wells (such as C67, C59, C8, C9, B30, B3 1, etc.). This phenomenon seems to be different from other models in the past, which can be explained by the following two factors: First, this model is a large-scale model at the basin level, and the observation well is far from the centralized water source, so the hydraulic conductivity of the aquifer is not very large. Secondly, groundwater discharge is mainly phreatic evaporation. After the centralized exploitation of groundwater, the reduced phreatic evaporation due to the drop of phreatic water level in the water source area and its nearby small area is enough to balance the exploitation amount, so that the head drop in the distance is not obvious.
For observation wells located in or near the water source area, the water head has an obvious downward trend after increasing mining. Generally speaking, the drop on the second floor is greater than that on the fourth and sixth floors.
See Table 5-9 for the multi-year average groundwater balance of 10 predicted by this plan.
Table 5-9 Annual average water balance of 10 predicted by the second scheme
The third scheme: See Figure 5- 12 and Figure 5- 13 for the prediction of groundwater level and groundwater depth in each aquifer.
Fig. 5- 12 Isohead Diagram of Scheme III 20211aquifer.
Fig. 5- 13 Isogram of Groundwater Depth in the Fourth Scheme 202 1 10
When the total production increases to 10.29× 108 m3/a, the groundwater depth in the calculation area is still 4.0~4.5m, which can ensure the ecological water demand. Compared with scheme 2, the groundwater exploitation in scheme 3 is increased by 6.7 1× 108 m3/a, while the water storage consumption is only increased by about 0.195×108m3/a. The increased exploitation is mainly due to evaporation conversion, and the phreatic evaporation in scheme 3 is reduced by 6.65×.
The prediction results of the third scheme show that laying large-scale mining wells in the shallow groundwater area in the south of the calculation area and increasing the exploitation of groundwater can not only intercept a large number of ineffective evaporation, but also alleviate environmental geological problems such as soil salinization after the groundwater level drops.
As long as the layout is reasonable, groundwater resources can be most effectively developed and utilized without affecting ecological development.
See table 5- 10 (from 12, 1 202 1 10) for the groundwater balance in the 20th year of the third scheme.
Table 5- 10 Water balance tables of the third scheme in 2020 12 to 10, 2021/kloc-0, 30.