According to the above heterogeneous, horizontally isotropic and vertically variable three-dimensional space structure and unstable thermal fluid system, it can be described by the following differential equation definite solution problem:
Figure 5- 16 Ng group 3D solid model modeling demonstration diagram
Exploration, development and utilization of geothermal fields in sedimentary basins
Where: ω is the seepage area; H is the water level elevation of hot fluid (m); Kx, KY and KZ are the permeability coefficients (m/d) in X, Y, Z, Y and Z directions in thermal fluid reservoirs, respectively; Kn is the permeability coefficient (m/d) of thermal fluid reservoir in the normal direction of boundary surface; Kz→ vertical permeability coefficient of thermal fluid reservoir within the boundary (m/d); S is the water storage rate (1/m); σ is the resistance coefficient of the common boundary, σ=L/K 1, l is the horizontal distance (m) from the model boundary to the common boundary, and K 1 is the average permeability coefficient (m/d) of the thermal fluid reservoir between the model boundary and the common boundary; ε is the source-sink term (1/d); H0 is the initial water level distribution (m), and h0 = h0 (x, y, z); Hb is the water level distribution on the common head boundary (m); Is the normal direction of the boundary surface; Is the normal direction of the interface between the top and bottom of the thermal storage; γ 1 is the general head boundary; γ 2 is the lateral water-resisting boundary of seepage zone; γ 3 is the upper and lower boundary of seepage zone;
(2) Numerical model
Most mathematical models cannot be solved by analytical methods. Digitization is to transform a mathematical model into a solvable numerical model. Commonly used numerical finite element method and finite difference method.
According to the characteristics of geothermal geological conceptual model of thermal reservoir, the richness and accuracy of actual data, the thermal reservoir of Guantao Formation is evaluated by finite difference method.
Modflow (Modular Fine Difference Groundwater Flow Model), a computer program that simulates the characteristics of groundwater flow and pollutant migration, is selected as the model software. The system is a three-dimensional simulation and evaluation software for groundwater flow and solute transport.
(3) Dividing calculation area units
Before the dynamic simulation of hot fluid, the simulation area must be divided into units. In addition to following the general division principle, the following actual conditions should be fully considered: ① calculation area boundary, lithologic division boundary, administrative division boundary, etc. (2) Observation wells should be placed on the nodes of subdivision units as far as possible; ③ In concentrated mining area, due to the large hydraulic gradient and the changing trend of flow field, the cutting should be properly densified. Therefore, the non-equidistant rectangular partition method is adopted in the plane, and the partition spacing is appropriately increased in the concentrated area of production wells. The outer and missing areas are defined as inactive units along the boundary. It is divided into 60 rows and 55 columns, including 1580 invalid cells. The maximum subdivision unit is about 1.8× 1.8km2 square unit, and the minimum subdivision unit is 0.28×0.38km2 rectangular unit. Vertically, thermal reservoirs are divided into NG ⅰ, middle weak permeable layers NG ⅱ and NG ⅲ. Then Ng has (60×55- 1580)×3 three-layer subdivision units, one of which is 5 160, as shown in Figure 5- 17 and Figure 5- 18.
(4) Boundary conditions
Natural gas thermal storage can be summarized into three boundary conditions: ① natural flow boundary; ② Zero-flow water-resisting boundary; ③ Internal recharge boundary formed along Cangdong fault (Figure 5- 19). Different types of boundary conditions are treated as follows:
A. Natural head boundary
The natural flow boundary is treated by the general head boundary in variable head boundary package (GHB) in MODFLOW software. The purpose of using this boundary condition is to avoid the unnecessary outward expansion of the model area to meet the water head affected by factors in the model. Therefore, the variable head boundary conditions are usually specified outside the model domain. The amount of hot fluid entering or leaving the computing unit from an external source is directly proportional to the difference between the head of the computing unit and the head of the external source.
On the universal head boundary element, the amount of hot fluid entering or leaving the computing unit from the hot fluid outside the element is directly proportional to the difference between the head H of the computing unit and the head hb of the universal head boundary hot fluid. Therefore, a linear relationship can be established between the hot fluid quantity of the calculation unit and the pressure head of the calculation unit, namely:
Fig. 5- 17 Ng plane element profile
Fig. 5- 18 Ng3D unit section
Figure 5- 19 Ng Calculation Boundary Condition Diagram
Exploration, development and utilization of geothermal fields in sedimentary basins
Where: Qb is the flow entering the calculation unit from the common head boundary (m3/d); Cb is the hydraulic conductivity coefficient between the general head boundary and the calculation unit (m2/d); (L×W) is the cross-sectional area (m2) of the boundary and grid element; K is the average permeability coefficient (m/d) of the aquifer material that separates the source/sink from the model grid; D is the distance (m) between the boundary source/sink and the model grid; Hb is the hot fluid pressure head (m) on the common pressure head boundary; H is the heat flow head (m) of the calculation unit.
When Cb is very large, the function of GHB subroutine package is the same as that of constant head boundary. Therefore, the water level and hydraulic conductivity coefficient Cb on the common head boundary determine the inflow and outflow of the boundary Qb.
B. Zero-flow water blocking boundary
In addition to the western boundary of Neogene Guantao Formation and some missing areas, the zero-flow water-resisting boundary is defined, and the top and bottom plates of thermal reservoirs are also defined as water-resisting boundaries. Undoubtedly, for the water-resisting boundary, the inflow and outflow Qf values of the calculation unit and the external domain are zero.
C. Internal recharge boundary formed along Cangdong fault
According to the characteristics of Cangdong fault, a channel connecting the upper and lower thermal reservoirs with geothermal fluid is formed along Cangdong fault. Overexploitation of natural gas thermal fluid in Tanggu District and Dagang District of Binhai will lead to excessive reduction of thermal fluid head in this layer, thus stimulating the replenishment of its lower thermal reservoir. Therefore, the vertical recharge of other aquifers should be considered along the Cangdong fault. In this paper, the river subroutine package provided by the numerical simulation software MODFLOW is used to simulate the hot fluid supply of Cangdong fault to the calculation area. The calculation formula is as follows:
Exploration, development and utilization of geothermal fields in sedimentary basins
Where: Qfau is the flow (m3/d) from the crack into the calculation unit, with a positive value indicating entry and a negative value indicating outflow; Cfau is the hydraulic conductivity coefficient (m2/d) between the external source of the fault zone and the fault zone calculation unit; L and w are the length and width (m) of the fault passing through the grid of the calculation unit; K is the average permeability coefficient (m/d) of the calculation unit between the fracture zone and the lower thermal reservoir; M is the thickness of calculation unit between fracture zone and lower thermal reservoir (m); Hfau is the thermal fluid head (M) of the lower thermal reservoir between fracture zones in the calculation unit; H is the head of calculation unit between fracture zones (m).
The head of bedrock thermal fluid reservoir in the lower part of Cangdong fault zone is generally higher than that of ng thermal fluid, and the seepage recharge of fault zone is generally positive. That is, the thermal fluid of the lower bedrock is replenished to the upper Ng thermal reservoir through the uplift of the fault zone. However, due to the limitation of original data, the accuracy of inflow and outflow of Qfau is rough and cannot be given accurately. Only in the flow field simulation can the calibration be made according to the simulation results of the flow field.
(5) Selection of simulation period and time dispersion
The numerical simulation starts from 2002 1 month and ends in 20061month, and ***5 years1month is taken as the whole simulation time length of this layer model identification. According to the actual characteristics of geothermal fluid development and utilization, a mining cycle is divided into two stress periods (time periods): heating period (10 to March 25th10, *** 135 days) and non-heating period (March 26th to/kloc).
(6) recharge, runoff and discharge conditions
For geothermal fluid, due to the deep burial, the corresponding recharge items are relatively simple: the recharge sources mainly include lateral runoff recharge, water release from thermal reservoir rock skeleton caused by the pressure change of thermal fluid in thermal reservoir, recharge from the same group of internal stratification levels and artificial recharge. Overexploitation of Neogene Guantao Formation may stimulate some bedrock to be replenished through Cangdong fault. The drainage methods are mainly lateral runoff drainage of downstream units, manual mining and overflow replenishment of other over-exploited layers in the same group.
The treatment of lateral runoff and internal recharge boundary has been introduced before, so I won't go into details here. The following mainly analyzes the processing of other supplementary items.
A. elastic water release
The main reason of water release from rock skeleton is the pressure change of hot fluid in thermal reservoir. The water release from rock skeleton is not only affected by the pressure change, but also the elastic water release coefficient of thermal reservoir is a crucial parameter. Its value is calculated by the following formula
qi = l×w×S *×δh 5-25
Where: Qi is the unit elastic water release (m3); L is the unit length (m); W is the unit width (m); S* is the elastic water release coefficient (dimensionless) of the hot fluid enrichment section where the unit is located; Δ h is the change in head per unit hot fluid (m)
B. Replenishment and discharge of interlayer overflow in the same group
Because of the three-dimensional flow model, the interlayer runoff is automatically processed by MODFLOW model software according to the given relevant parameters.
C. Mudstone water release
Mudstone water release refers to the bound water commonly existing in soil and the water released by capillary water due to soil compression. Mudstone drainage will cause land subsidence. Obviously, the land subsidence caused by geothermal mining is limited, and the amount of mudstone released will not be too large, but it can not be ignored.
According to the research on geothermal fluid replenishment in Tanggu District, Tianjin (Li Minglang, 1998) in June of 5438+0998, the water released by cohesive soil accounts for about 1/5 of the total released water (elastic water release+clay water release). When the total water discharge is high, it accounts for about 40% of the hot fluid output. However, with the development of resources, the pressure density of cohesive soil increases after water release, and its water release should gradually decrease. With the decrease of thermal storage pressure, the total water release will gradually decrease.
D. Artificial mining and reinjection treatment
For manual mining and replenishment, both are replenished in the form of mining wells in the model, with positive mining and negative replenishment.
(7) Hydrogeological parameters
Hydrogeological parameters mainly refer to permeability coefficient K and elastic water release coefficient S*, the former K reflects the permeability of thermal reservoir fluid; The latter S* reflects the water storage (water discharge) capacity of thermal reservoirs. The selection of hydrogeological parameters is based on the calculation needs of three-dimensional model, and the method of partition assignment is mainly adopted. According to the buried distribution and thermal control structure distribution characteristics of ng thermal reservoir, the collected depressurization test data are used for comprehensive statistical analysis, and preliminary zoning is carried out according to the required hydrogeological parameters (see Figure 5-20). Conditional regional projects supplement the preliminary depressurization test of geothermal wells, and blank areas give estimated values according to the data of adjacent areas and the characteristics of stratigraphic structural changes.