River and lake groundwater refilling coupling simulation method

Abstract

The invention relates to a river and lake groundwater refilling coupling simulation method. The method comprises the steps of through water refilling riverbank zone refilling influence research, plain area refilling simulation and deep water source replacement simulation evaluation. According to the through water replenishing riparian zone refilling influence research, through a riparian zone hydrogeological conceptual model and a riparian zone numerical model, after model identification and verification, riparian introduction seepage quantity analysis and riparian zone refilling influence analysis are carried out; plain area back-filling simulation is based on a T-PROGS three-dimensional hydrogeological structure model, a typical area groundwater flow numerical model is established based on a U-Grid river channel grid encryption technology and river channel infiltration analysis of a fuzzy mathematical method, and groundwater back-filling effect analysis and groundwater over-mining treatment effect prediction analysis are carried out; in the deep water source replacement simulation evaluation, an underground water numerical model is firstly constructed, and then deep underground water replacement effect evaluation is performed, including overall evaluation of the influence of water source replacement on the underground water level and evaluation of the influence of water source replacement in different time periods on the underground water level.

Description

Technical Field

[0001] The present invention relates to the technical field of large-scale ecological water replenishment of rivers and lakes, and in particular to a coupled simulation method for groundwater replenishment in rivers and lakes. Background Art

[0002] In recent years, before the main flood season every year, key water replenishment actions such as summer ecological water replenishment and full-line water replenishment are organized and implemented. However, there is no more reasonable simulation evaluation, which cannot provide strong support for the evaluation of water replenishment effects, the impact of groundwater infiltration and replenishment, and the study of groundwater migration laws in the riverside zone. Summary of the invention

[0003] The present invention aims to solve the deficiencies of the prior art and provide a coupled simulation method for groundwater recharge in rivers and lakes.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] A coupled simulation method for groundwater recharge in rivers and lakes, including the study of the impact of recharge on the riverbank zone, recharge simulation in plain areas, and simulation evaluation of deep water source displacement;

[0006] The study on the impact of riparian zone recharge is to conduct riparian zone infiltration analysis and riparian zone recharge impact analysis after the riparian zone hydrogeological conceptual model, riparian zone numerical model and model identification verification;

[0007] The recharge simulation in the plain area is based on the three-dimensional hydrogeological structure model of T-PROGS, the river grid encryption technology based on U-Grid, and the river infiltration analysis based on fuzzy mathematics methods. The numerical model of groundwater flow in typical areas is established to analyze the effect of groundwater recharge and the effect of groundwater over-exploitation control. The analysis of the effect of groundwater recharge includes historical flow field changes and flow field changes after water recharge. The prediction analysis of the effect of groundwater over-exploitation control includes water level changes, impact range and equilibrium analysis.

[0008] The simulation evaluation of deep water source replacement first constructs a groundwater numerical model, including the generalization of the aquifer system and hydrogeological parameters, and then evaluates the effect of deep groundwater replacement, including an overall evaluation of the impact of water source replacement on groundwater levels and an evaluation of the impact of water source replacement on groundwater levels at different times.

[0009] The hydrogeological concept model of the riparian zone is divided into two layers. The first layer of the model is the phreatic aquifer, where groundwater extraction is mainly mixed extraction; the second layer of the model is the confined aquifer, where groundwater extraction is small and there is a lack of corresponding water level monitoring data; the groundwater recharge items in the area include rainfall infiltration recharge, river seepage recharge, and agricultural irrigation infiltration recharge. Agricultural irrigation infiltration recharge includes canal leakage, canal irrigation infiltration and well irrigation regression; the groundwater discharge items in the area are phreatic evaporation discharge and groundwater extraction.

[0010] In the riparian zone numerical model:

[0011] Groundwater flow can be described by the following equation:

[0012]

[0013] Among them, K xx , K yy , K zz are the components of the permeability coefficient in the x, y, and z directions (LT -1 ); h is the water level of the aquifer (L); W represents the source and sink term (T -1 );S s For unconfined aquifers, it is the water supply rate; for confined aquifers, it is the water storage rate (L -1 );t is time (T).

[0014] The model identification and verification follows the following principles: the hydrogeological parameters in the model can reflect the hydrogeological conditions in the study area; the final groundwater flow field during the simulation period is basically consistent with the actual groundwater flow field during that period; the simulated groundwater level process line is basically similar to the measured water level process line; the simulated groundwater equilibrium is consistent with the actual water equilibrium status.

[0015] The analysis of riparian zone infiltration is based on the water balance method, using monitoring data such as rainfall, evaporation and water withdrawal to calculate the infiltration rate of each river section.

[0016] The impact analysis of riparian zone recharge uses a model to simulate the groundwater level distribution in each river section without water recharge, and compares it with the measured groundwater level under the current water recharge situation.

[0017] The numerical model of groundwater flow in a typical area uses the following differential equation to describe the well-defined problem:

[0018]

[0019] Where:

[0020] Ω-seepage area; h-water level of groundwater system (m); K-horizontal permeability coefficient of water-bearing medium (m / d); K z -vertical permeability coefficient of water-bearing medium (m / d); ε-source and sink term of aquifer (1 / d); h0-initial water level distribution of the system (m); S-water storage rate of aquifer below free surface (1 / m); Γ0-upper boundary of seepage area, i.e. free surface of groundwater; μ-gravity water supply of phreatic aquifer on phreatic surface; p-evaporation and precipitation infiltration intensity of phreatic surface, etc. (m / d); Γ1-known water level boundary; h1-known boundary water level value (m); Γ2-flow boundary of seepage area; K n-permeability coefficient in the normal direction of the boundary surface (m / d), n-normal direction of the boundary surface; q-flow rate per unit area of ​​Γ2 boundary (m / d), inflow is positive, outflow is negative, and the water-blocking boundary is 0.

[0021] The flow field changes after water replenishment are specifically evaluated in the following four aspects:

[0022] Changes in flow field morphology: Compare the dynamic changes of the flow field in the water replenishment area at different time periods;

[0023] Variation of water level in river channel: Carry out variation analysis on water level before and after water replenishment in river channel, and draw variation chart to illustrate the variation of water level before and after water replenishment;

[0024] Recharge impact range: According to the model fitting accuracy, the water level variation range of 0.5m was selected to simulate the water level changes under the conditions of "recharge" and "no recharge" in the study area, and the range of water level rise caused by river water replenishment was determined under the same other hydrogeological conditions;

[0025] Equilibrium analysis: Analyze the groundwater equilibrium in the study area under the two conditions of "water replenishment" and "no water replenishment", and judge whether water replenishment can effectively restore the groundwater level from the perspective of equilibrium changes.

[0026] The beneficial effects of the present invention are as follows: in order to clarify the impact of river and lake replenishment on groundwater, the present invention conducts a coupled simulation study on river and lake groundwater replenishment from multiple dimensions of point-line-surface, and conducts groundwater replenishment impact studies on different areas such as typical riverbank zones, river and lake water replenishment receiving areas, and deep pressure mining and replacement areas, respectively, which provides strong support for water replenishment effect evaluation, groundwater infiltration and replenishment impact, and research on groundwater migration laws in riverbank zones. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 It is a schematic diagram of the process of the present invention;

[0028] The following is a detailed description of the embodiments of the present invention with reference to the accompanying drawings. DETAILED DESCRIPTION

[0029] The principles and features of the present invention are described below in conjunction with the accompanying drawings. The embodiments are only used to explain the present invention and are not used to limit the scope of the present invention. The present invention is described in more detail by way of example with reference to the accompanying drawings in the following paragraphs. The advantages and features of the present invention will become clearer according to the following description. It should be noted that the accompanying drawings are all in a very simplified form and are not in precise proportions, and are only used to facilitate and clearly assist in explaining the purpose of the embodiments of the present invention.

[0030] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those commonly understood by those skilled in the art of the present invention. The terms used herein in the specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The term "and / or" used herein includes any and all combinations of one or more related listed items.

[0031] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0032] A coupled simulation method for river-lake groundwater recharge, such as Figure 1 As shown, it includes the study on the impact of recharge on the riverbank of the through-flow replenishment, the recharge simulation in the plain area and the simulation evaluation of deep water source replacement; the study on the impact of recharge on the riverbank of the through-flow replenishment is to analyze the infiltration of the riverbank and the impact of recharge on the riverbank through the hydrogeological conceptual model of the riverbank, the numerical model of the riverbank and the model identification and verification; the recharge simulation in the plain area is based on the three-dimensional hydrogeological structure model of T-PROGS, the river grid encryption technology based on U-Grid and the river infiltration analysis based on fuzzy mathematics method to establish a groundwater flow numerical model in a typical area, and conduct groundwater recharge effect analysis and groundwater over-exploitation control effect prediction analysis; the groundwater recharge effect analysis includes historical flow field changes and flow field changes after water replenishment; the groundwater over-exploitation control effect prediction analysis includes water level changes, impact range and equilibrium analysis; the deep water source replacement simulation evaluation first constructs a groundwater numerical model, including the generalization of the aquifer system and hydrogeological parameters, and then evaluates the effect of deep groundwater replacement, including the overall evaluation of the impact of water source replacement on groundwater level and the evaluation of the impact of water source replacement on groundwater level at different periods. Specific embodiment 1:

[0034] Taking the Beijing-Hangzhou Grand Canal as an example, a study on the impact of replenishment on the riverbank zone was conducted.

[0035] The Beijing-Hangzhou Grand Canal is the world's longest and largest ancient canal, and it is also one of the oldest canals. Starting from Yuhang (now Hangzhou) in the south and ending at Zhuojun (now Beijing) in the north, it passes through Zhejiang, Jiangsu, Shandong, Hebei, Tianjin, and Beijing, and connects the five major water systems of Haihe River, Yellow River, Huaihe River, Yangtze River, and Qiantang River, with a total length of about 1,794 km. The section of the Grand Canal north of the Yellow River is 822 km, accounting for 45.8% of the total length. The water replenishment in this embodiment uses the section of the Grand Canal north of the Yellow River as the main through-line, starting from Dongbianmen, Beijing, in the north, through Tonghui River and Beiyunhe to Sanchahekou, Tianlu City, starting from Weishan Gate, Liaocheng City, Shandong Province in the south, through Xiaoyunhe, Weiyunhe, and Nanyunhe to Sanchahekou, Tianlu City, involving Beijing, Tianlu, Hebei, and Shandong, and flowing through 8 prefecture-level administrative regions and 31 county-level administrative regions. Since the 1970s, due to the rapid increase in urban and industrial water consumption, the river section north of the Yellow River has gradually dried up.

[0036] 1. Conceptual model of riparian zone hydrogeology

[0037] (1) Simulation scope and boundary conditions

[0038] The scope of this simulation is based on the assessment scope, covering 62 counties and districts in 5 provinces and cities, with a total area of ​​50,500 km2. The specific scope is shown in the figure.

[0039] In the horizontal direction, the simulation area is bounded by the administrative boundaries of the counties and cities, which are transformed into the universal water head boundary. In the vertical direction, the water table of the phreatic aquifer is the upper boundary, which receives the replenishment of atmospheric precipitation infiltration and river infiltration; the bottom boundary of the Quaternary loose layer is the lower boundary, and the underlying Tertiary clastic rocks are generalized as the water-blocking boundary.

[0040] (2) Hydrogeological structure model

[0041] The model is divided into two layers. The Quaternary system in the study area is divided into the Lower Pleistocene, Middle Pleistocene, Upper Pleistocene and Holocene from old to new, corresponding to four aquifer groups, and the interlayer groups are relatively independent but have certain hydraulic connections. Among them, the hydraulic connection between the Holocene and the Upper Pleistocene is relatively close, and groundwater exploitation is mainly mixed mining. This simulation transforms it into the first layer of the model, which is a phreatic aquifer; the Lower Pleistocene is buried deep, the groundwater exploitation volume is small, and the statistical data is difficult to distinguish from the Middle Pleistocene exploitation volume. At the same time, there is a lack of corresponding water level monitoring data, so these two layers are generalized as the second layer of the model, which is a confined aquifer.

[0042] (3) Groundwater sources and sinks

[0043] The main groundwater recharge items in the area include rainfall infiltration recharge, river seepage recharge, agricultural irrigation infiltration recharge (canal leakage, canal irrigation infiltration and well irrigation return). The discharge items are mainly groundwater evaporation discharge and groundwater extraction.

[0044] 2. Numerical model of riparian zone

[0045] (1) Mathematical model

[0046] The groundwater flow in this embodiment can be described by the following equation:

[0047]

[0048] Among them, K xx , K yy , K zz are the components of the permeability coefficient in the x, y, and z directions (LT -1 ); h is the water level of the aquifer (L); W represents the source and sink term (T -1 );S s For unconfined aquifers, it is the water supply rate; for confined aquifers, it is the water storage rate (L-1 );t is time (T).

[0049] (2) Simulation software

[0050] This embodiment uses the MODFLOW-USG program of the GMS software to numerically simulate the groundwater flow in the study area, and on this basis, predicts the changes in groundwater levels in the Beijing-Hangzhou Grand Canal after the water replenishment and the scope of water replenishment. Unlike conventional MODFLOW programs, MODFLOW-USG supports a variety of structured and unstructured grid types, such as nested grids and grids based on prismatic triangles, rectangles, hexagons and other unit shapes, where each grid cell can be connected to any number of adjacent cells. This flexible subdivision method can be used to focus resolution along rivers and around wells, or to subdivide individual layers to better represent hydrological stratigraphic units, etc. In this embodiment, this subdivision method facilitates the assignment of planar exploitation of the confined aquifer below the uppermost aquifer.

[0051] (3) Space-time discreteness

[0052] The simulation area is gridded as unstructured grid with a grid size of 1km×1km. The vertical model is generalized into two layers, the first layer is the phreatic aquifer and the second layer is the confined aquifer. The simulation period of the model is from January 2010 to December 2020, with one stress period per month, for a total of 144 stress periods.

[0053] (4) Initial conditions and boundary conditions

[0054] The groundwater flow field drawn from the borehole measurements in January 2010 in the North China Plain was selected, and part of the simulation range was selected as the initial flow field of the model.

[0055] (5) Hydrogeological parameters

[0056] The initial values ​​of the hydrogeological parameters in the simulation area are given according to the actual characteristics of the water-bearing rock group and the water level geological parameters of the North China groundwater flow model in combination with the drilling data. The hydrogeological parameters are subsequently identified by fitting the hydrological process lines of each observation well to determine the final parameter values.

[0057] (6) Source and sink processing

[0058] ① Rainfall infiltration recharge: The rainfall infiltration recharge in this embodiment is calculated on a county basis. The rainfall infiltration recharge is calculated by multiplying the rainfall in each district by the corresponding infiltration coefficient. The calculation formula is as follows.

[0059] Q = α × P × S;

[0060] Where Q is the monthly precipitation infiltration recharge, m 3 / M; α is the precipitation infiltration coefficient; P is the monthly precipitation, m / M; S is the area of ​​the calculation area, m 2 .

[0061] This example collects monthly rainfall data involving 62 counties and districts. The rainfall infiltration coefficient refers to the rainfall infiltration coefficient field in the North China Groundwater Sustainable Development Atlas. As shown in the figure, the rainfall infiltration recharge in 2020 is calculated to be 50.93×108m 3 / a. The rainfall infiltration is input into the model using the RCH package.

[0062] ② River seepage recharge: In this embodiment, the equilibrium method is first used to calculate the seepage amount of the river channel, and the calculation formula is as follows.

[0063] Q RL =Q1-Q2+Q in -Q Z -Q S -0.1·(EP)·S;

[0064] In the formula, Q RL The leakage of the water replenishment section is 104m 3 ; Q1 and Q2 are the water replenishment of the starting section and the water discharge of the end section respectively; Q in Q is the inflow volume of the river section; Z is the outflow volume of the river section, 104m 3 ;Q S is the channel storage variable in the river section; E and P are the evaporation and precipitation acting on the river section, mm; S is the area of ​​the corresponding river section, km 2 The leakage volume was input into the model by river section using the RIV module provided by GMS.

[0065] ③ Agricultural irrigation infiltration recharge: The agricultural irrigation infiltration recharge in the study area is divided into canal leakage, canal irrigation infiltration and well irrigation regression recharge. The canal leakage and canal irrigation infiltration are calculated based on the irrigation area, and the well irrigation regression recharge is calculated based on the county. The specific canal leakage, canal irrigation infiltration, well irrigation regression recharge and irrigation area distribution refer to the "Research Report on the Evolution Mechanism and Regulation of Groundwater in the North China Plain" and follow the corresponding part of the simulation scope. The canal leakage, canal irrigation infiltration and well irrigation regression recharge are processed into surface data with irrigation areas and counties as basic units and distributed to each month. The agricultural irrigation infiltration is input into the model using the RCH package. The well irrigation regression recharge in 2019-2010 is 2.83×108m 3 / a, the canal leakage and canal irrigation infiltration are 6.03×108m 3 / a.

[0066] ④ Phreatic evaporation discharge: Phreatic evaporation is affected by factors such as water surface evaporation intensity, phreatic depth, and lithology of the aeration zone. This model sets the phreatic evaporation limit depth to 4m, the surface elevation is used as the evaporation surface elevation, and the evaporation intensity is based on the county. The meteorological data is used to generate the evaporation intensity sequence, and the evaporation intensity data is input into the model using the EVT module of GMS. The specific calculation of phreatic evaporation is shown in the formula.

[0067]

[0068] Among them, E g is the groundwater evaporation intensity, mm; d0 is the evaporation limit buried depth, m; d is the water level buried depth, m; E0 is the water surface evaporation intensity, mm; n is the empirical coefficient, which is taken as 1 in the model.

[0069] ⑤ Groundwater extraction:

[0070] The simulation area is located in the North China Plain, which has many urban communities, dense population, developed agriculture, a large number of groundwater extraction wells and a large amount of extraction. Artificial extraction is the main way of groundwater discharge in the study area. It can be divided into agricultural, domestic and industrial extraction according to its use. Among them, agricultural extraction has obvious seasonality. Irrigation is 1-2 times during the winter wheat sowing period (October and November), and 2-3 times during the greening period-jointing and heading period-filling and maturity period (March-June). During this period, the groundwater extraction volume increased significantly, which are the two stages with the largest underground extraction volume. There is no obvious seasonal change in the extraction of industrial and domestic water, and this part of the extraction volume is evenly distributed to each month. This study collected the extraction volume data of each district and county, and allocated it to each water-using department and each month according to a fixed ratio. The extraction volume of the groundwater aquifer in the study area from 2019 to 2020 was 24.73×108m 3 / a, and the extraction volume of confined aquifer is 14.63×108m3 / a.

[0071] The WELL module and RCH module of GMS software can handle groundwater exploitation. Considering that there are many exploitation wells in the simulation area, this simulation treats groundwater exploitation as a surface exploitation input model. The RCH module can handle surface exploitation well.

[0072] 3. Model Identification Verification

[0073] The identification and verification of the model should follow the following principles: (1) the hydrogeological parameters in the model can reflect the hydrogeological conditions in the study area; (2) the final groundwater flow field during the simulation period should be basically consistent with the actual groundwater flow field during that period; (3) the simulated groundwater level process line should be basically similar to the measured water level process line; (4) the simulated groundwater equilibrium should be consistent with the actual water equilibrium status.

[0074] The measured groundwater flow field in the simulation area at the end of 2019 was collected and compared with the final groundwater flow field calculated by the groundwater model. The fitting status of the phreatic aquifer is shown in the figure. The total and bulk flow direction of the phreatic aquifer is from the north, west and south to the Bohai Bay in the east. Near the northern boundary, the simulated water level is different from the measured water level. The analysis shows that this is due to the large parameter variation at the alluvial fan in front of the mountain. Overall, the simulated flow field in the shallow layer is basically consistent with the measured flow field and the trend is consistent. The simulated flow field basically reflects the flow characteristics of groundwater.

[0075] The fitting status of the confined aquifer is shown in the figure. Since the confined aquifer is simulated in this simulation, the extraction volume is processed in a planar manner, which cannot accurately depict the concentrated extraction of groundwater. Therefore, the fitting of the funnel is not very good, and only the trend of the funnel is depicted. However, overall, the simulated flow field and the measured flow field are basically similar in form, the simulated flow field and the measured flow field trends are basically consistent, and the range of the falling funnel is roughly consistent with the actual one.

[0076] 4. Analysis of the impact of groundwater recharge in riparian zones

[0077] (1) Infiltration analysis

[0078] Based on the water balance method, the infiltration rate of each river section is calculated using monitoring data such as rainfall, evaporation and water withdrawal. From March 1 to June 1, the water inflow in the water replenishment trunk line was 155 million m 3 The water volume diverted from the area is 534 million m 3 , the new tank storage capacity is 29 million m 3 After deducting the precipitation, the evaporation loss is 119.1 million m 3 The cumulative infiltration volume is 144 million m 3 The average infiltration rate is 13%. The infiltration rate of the North Canal is 19%, the infiltration rate of the South Canal is 10%, the infiltration rate of the Wei Canal is 12%, and the infiltration rate of the Xiaoyun Canal is 7%. The specific analysis is as follows:

[0079] ①Northern Canal

[0080] The infiltration rate of the North Canal is calculated using the typical river section from Tongxian Station to Kuangergang Station. The rainfall is measured by Beiguan Weir Station and Tumenlou Station, the evaporation is measured by Tumenlou Station, and the trough storage capacity is calculated based on the water level changes at Tongxian Station and Kuangergang Station. The rock particles in the front area of ​​the North Canal are coarse and suitable for infiltration, and the middle and lower reaches have water all year round and the infiltration rate is low. Since the water replenishment, the amount of water flowing into the North Canal from Zhangjiawan Station of Liangshui River is 115 million m 3 According to statistics in recent years, the water from Beijing's sewage river flows into the river section from Kuangergang to Qujiadian hub, about 40 million m 3 The North Canal diverts about 152 million m3 of water through Niumutun, Tumenlou (Qing) and other canals. 3According to statistics in recent years, about 80 million m3 of water has been diverted through the Kangergang hub, Qujiadian hub, and agricultural water diversion along the line in Tianlu City. 3 According to the water balance calculation, the average infiltration rate of the North Canal is 19%, and the infiltration volume is 58 million m 3 .

[0081] ②South Canal

[0082] The infiltration rate of the South Grand Canal is calculated using the river section from Gengliyang Station to Jiuxuan Sluice. The rainfall is calculated using the monitoring values ​​of 11 rain gauges including Anling Station, and the evaporation is calculated using the monitoring values ​​of 4 evaporation stations including Si Nu Si Sluice and Anling Station. The trough storage capacity is calculated based on the water level changes of Gengliyang Station and Jiuxuan Sluice. The South Grand Canal has more clay soil and weaker permeability. Since the water replenishment, the water diversion volume of the South Grand Canal is 237 million m 3 According to the water balance calculation, the infiltration rate of the South Canal is 10%, and the infiltration volume is 31 million m 3 .

[0083] ③Weiyun River

[0084] The infiltration rate of Weiyun River is calculated from the section from Linqing Station to Si Nu Si Hub. The evaporation and rainfall are based on the monitoring values ​​of Linqing Station, and the trough storage capacity is calculated based on the water level changes of Linqing Station and Si Nu Si Hub. The proportion of clay in the Weiyun River channel is relatively small and has good permeability. Since the water replenishment, the water diversion volume of Weiyun River is 0.66 billion m 3 According to the water balance calculation, the infiltration rate of Weiyun Canal is 12% and the infiltration volume is 45 million m 3 .

[0085] ④Small Canal

[0086] The infiltration volume of the Xiaoyunhe section was calculated from the Dongping Lake Outlet Gate to the Qiutun Hub, and the evaporation and rainfall were respectively calculated from the monitoring values ​​of Liuqiao Gate and Liaocheng Station. The Xiaoyunhe main canal 0+000-73+860, a total of 74km, was lined and had weak permeability. According to the water balance calculation, the infiltration rate of the Xiaoyunhe was 7%, and the infiltration volume was 100 million m 3 .

[0087] (2) Analysis of the impact of backfilling

[0088] The model was used to simulate the groundwater level distribution in each river section from March to June 2023 without water replenishment, and compared with the measured groundwater level under the current water replenishment conditions.

[0089] According to preliminary analysis, the maximum impact distance of groundwater on both sides of the water replenishment trunk line is 5.11km, the average impact distance is 1.56km, and the cumulative impact area is 990.88km 2. The average impact distance of the North Canal on groundwater is 1.14km; the average impact distance of the South Canal on groundwater is 1.39km; the average impact distance of the Wei Canal on groundwater is 2.04km; and the average impact distance of the Small Canal on groundwater is 0.76km. Since March to May is the agricultural irrigation period along the coast, the groundwater level is still mainly declining. Ecological water replenishment has slowed down the decline of groundwater levels around the river to a certain extent. The specific analysis results are as follows:

[0090] Table 1 Scope and extent of groundwater recharge impact (as of June 1, 2023)

[0091]

[0092] ①Northern Canal

[0093] According to the simulation prediction of whether there is water replenishment or not, the water level in the area near the North Canal is rising. The increase of groundwater level caused by water replenishment is generally 0.3m, and the maximum increase of groundwater level caused by water replenishment is 0.7m. The maximum distance of water replenishment on groundwater is 2.40km, the average distance of water replenishment is about 1.14km, and the affected area is 144.96km 2 .

[0094] ②South Canal

[0095] According to the simulation prediction of whether there is water replenishment or not, the overall water level along the South Canal has rebounded. The groundwater level rise caused by water replenishment is generally 0.1m, and the maximum water level rise near the river bank is 0.3m. The maximum impact distance of water replenishment on groundwater is 3.52km, the average impact distance is about 1.39km, and the affected area is 364.79km 2 .

[0096] ③Weiyun River

[0097] According to the simulation prediction of whether there is water replenishment or not, the groundwater level of Weiyun River has risen significantly as a whole, with an increase of 0.9m in general, and the maximum water level in the riverside area has risen to about 1.1m. The maximum impact distance of water replenishment on groundwater is 5.11km, the average impact distance is about 2.04km, and the affected area is 394.26km 2 .

[0098] ④Small Canal

[0099] The southern section of the Xiaoyun River is lined, and the area affected by the groundwater level is only in Linqing City near the Qiutun hub. The groundwater level rise caused by water replenishment is generally 0.1m, and the maximum rise at the junction with the Weiyun River is about 0.5m. The maximum impact distance of water replenishment on groundwater is 4.72km, the average impact distance is 0.76km, and the potential impact area is 86.87km 2 . Specific embodiment 2:

[0101] A simulation study on replenishment in plain area was conducted using Daqing River as an example.

[0102] The plain area of ​​the Daqing River Basin is located in the middle of the basin, bordering the Bohai Sea to the east, the Taihang Mountains to the west, the Yongding River to the north, and the Hutuo River to the south. The administrative division includes the plain parts of Beijing, Tianlu, and Hebei provinces (cities), with an area of ​​20,000 km2. According to the genesis and morphological characteristics, it can be divided into piedmont alluvial diluvial inclined plains, central and eastern alluvial lacustrine plains, the Yellow River alluvial fan, and coastal alluvial marine plains. There are many rivers and lakes in the basin. Baiyangdian, known as the "Pearl of North China", is located in the southeastern part of the region. There are 9 rivers flowing into the lake from the south, west, and north of Baiyangdian, namely Zhulong River, Xiaoyi River, Tang River, Qingshui River, Fu River, Cao River, Pu River, Ping River, and Baigou Yinhe River.

[0103] 1. Three-dimensional hydrogeological structure model based on T-PROGS.

[0104] The solution to the problem of heterogeneity of water storage media is usually to combine and analyze the vertical lithology data and profile data of the borehole, so as to classify and merge them according to the water-bearing characteristics of the lithofacies, and finally generalize the spatial distribution structure of the aquifer, that is, vertical stratification and horizontal zoning. However, in the generalization process, due to the different distribution forms of lithofacies and the large variability in space, and the appearance of lenses will cause the lithology to mutate, the situation is more complicated at the regional scale. It is generally difficult to truly and accurately reflect the spatial distribution of aquifers and impermeable layers. There are also human factors in the description of hydrogeological parameters on the plane, and the setting of parameter zoning is a difficult point in the study. Therefore, when establishing a three-dimensional groundwater numerical simulation model, all these factors will directly or indirectly affect the accuracy of the simulation model, and then affect the prediction and analysis of the groundwater flow field.

[0105] The geostatistics method with transition probability as the core algorithm identifies the lithology distribution in three-dimensional space based on discrete drilling data, forming a three-dimensional discrete block structure of the hydrogeological structure. Through Markov-based simulation, a mathematical model of the groundwater system is established with hydrogeological units, the simulation results of the water storage medium are characterized by lithology parameters, and the aquifer structure is constructed in a random model. This solves certain medium heterogeneity and provides basic support for the establishment of a numerical simulation model for groundwater flow. In addition, this project combines three-dimensional geological structures with fuzzy mathematical methods to calculate the segmented infiltration coefficient of the river section, providing quantitative data for the recharge model and river recharge capacity analysis.

[0106] (1) Drilling information

[0107] The establishment of a 3D geological structure model of the plain area of ​​the Daqing River Basin requires a large amount of geological drilling data in the area as raw data. According to previous work results and the collected drilling holes in the area, there are 83 holes in total. From the perspective of revealing the strata, they can be roughly divided into drilling holes that expose the Quaternary System, drilling holes that expose the Quaternary System, and shallow drilling holes. The description of lithology in each category is different. The drilling holes that expose the Quaternary System can reach more than 400m, and the shallow drilling is only ten to tens of meters. In terms of quantity, the drilling data is relatively rich, meeting the requirement of 3 holes / 100km. 2 The standard is met, but it is still necessary to virtual expand according to the profile data and geological data so that the distribution of the boreholes in the whole area can be more uniform, which is conducive to better establishing a three-dimensional geological model after interpolation. After expansion, the number of boreholes reached 95.

[0108] (2) Drilling generalization and arrangement

[0109] According to the needs of structural simulation, the contents of the borehole data are organized, including borehole number, geographic coordinates (latitude and longitude coordinates or geodetic coordinates), borehole depth, borehole elevation, and borehole lithology description. For the 83 real boreholes established in the Daqing River Basin and 12 virtual boreholes, a total of 95 boreholes were organized and cataloged. According to the borehole lithology description, the borehole lithology is uniformly planned as: medium-coarse sand, fine sand, silt and clay. Among them, medium-coarse sand and fine sand are the main water-bearing media in the study area, while silt and clay are the main relatively water-repellent media. According to the borehole data, it is discretized into 2m in the vertical direction, and 1-4 is used as the number of the four lithologies for transfer probability calculation, and the Kriging interpolation equation group is established.

[0110] After the original data enters the program, it will be used to statistically calculate the average lithology content, vertical average thickness and transfer probability matrix. The simulation process is as follows: ① The original drilling data is tabulated and catalogued, including unified coordinates, stratigraphic age, lithology code, lithology bottom boundary depth, top boundary depth, ground elevation and lithology description; ② The lithology is generalized and classified, mainly based on the particle size, generally not more than 5 categories; ③ The borehole spatial coordinates, lithology segment elevation and lithology code are extracted; ④ The lithology is processed as an indicator variable, and the lithology is expressed in vector form to obtain the data file used for calculation.

[0111] (3) Regional division

[0112] After repeated experiments, with simulation accuracy as the main factor, it was finally established that the study area was discreted into 2km×2km in plane and 2m in vertical direction. The depth involved was greater than 550m, which was able to reveal the Quaternary system. Finally, 885,780 unit blocks were obtained.

[0113] (4) Implementation of random simulation of hydrogeological structure

[0114] After gridding the study area, the transfer probability curves and Markov transfer probability curves of the four generalized lithologies of the borehole were obtained through T-PROGS calculation. The vertical simulation effect of the constructed hydrogeological structure model was determined by fitting the two. After debugging and identification, it was extended to the entire area, and finally a three-dimensional geological structure reflecting the hydrogeological structure of the study area was formed.

[0115] According to the simulation results, the sedimentary structure of the entire Daqing River Basin plain area has a certain cyclicity. From top to bottom, the Quaternary grain size changes from fine to coarse and then to fine, forming a relatively complete sedimentary cycle. From the front of the mountain to the coastal area, the lithology changes from coarse to fine. Most of the piedmont area is zcs with large grain size, which transitions to fine sand in the central plain, and fine grains are increasingly developed, such as ft and fn. In the near coastal area, the content of fine grains reaches the maximum.

[0116] 2. U-Grid river grid encryption technology

[0117] This work will use the quadtree method to grid densify the four rivers in the study. Considering the characteristics of river sediment changes along the river and the calculation efficiency, the piedmont area is dominated by coarse-grained sediments and the mainstream field changes slowly, so the densification degree is relatively low. In the central and eastern parts of the relevant plains, the groundwater level changes more dramatically, and the densification degree is also increased accordingly.

[0118] The grid density of the upper reaches of Juma River is 500m×500m, the grid density of the middle reaches is 250m×250m, and the grid density of the lower reaches is 100m×100m. The grid density of the upper reaches of Tang River is 500m×500m, the grid density of the middle reaches is 250m×250m, and the grid density of the lower reaches is 100m×100m. The grid density of the upper reaches of Dasha River is 500m×500m, the grid density of the middle reaches is 250m×250m, and the grid density of the lower reaches is 100m×100m. The grid density of the Yongding River is 250m×250m within about 1 km of the river channel, 100m×100m within about 500 meters from the river channel, 50m×50m within about 200 meters from the river channel, and 25m×25m in the river channel.

[0119] 3. River infiltration analysis based on fuzzy mathematics method

[0120] (1) Basic information of the river

[0121] This work mainly calculates the recharge infiltration coefficients of the Yongding River, Juma River, Tang River and Dashaha River, and then analyzes the infiltration capacity of their different sections.

[0122] The Juma River in the study area is 122.47 km long. According to the profile characteristics and the phreatic level, the thickness of the infiltration section is determined to be 3-22 m. The Tang River is 117.75 km long. According to the profile characteristics and the phreatic level, the thickness of the infiltration section is determined to be 2-47 m. The Dasha River is 191.68 km long. According to the profile characteristics and the phreatic level, the thickness of the infiltration section is determined to be 5-39 m.

[0123] (2) Fuzzy mathematics method

[0124] This study integrates the geological exploration data of lithology borehole stratification in the study area, and uses the generalized borehole lithology stratification data (BoreholesData) to construct a lithology joint profile based on the GMS software platform. At the same time, the above data is used to build a three-dimensional geological structure model of the study area based on Tins universal surface grid division and Solid solid filling generation function, and the T-PROGS transfer probability model of the study area is calculated based on MCMOD and TSIM tools. Comprehensively analyze the joint profile of a single river channel and the information of lithology stratification on the cross section obtained by cutting along the river channel of the two structural models, draw the river channel lithology profile in CorelDraw high-precision drawing software, organize the intuitive information of the lithology profile and combine it with fuzzy mathematics method, and establish fuzzy mathematical models of each river channel.

[0125] Fuzzy mathematical model: For a given domain: Y {rock layers of different lithology}, there is a fuzzy subset: A {influence on river infiltration rate}. There is a certain number μ A (δ)∈[a,b], then there is a mapping:

[0126] μ A (δ): Y→[a, b], δ→μ A (δ);

[0127] In the formula, μ A (δ) is called the membership function of A, μ is called the element in the domain Y, and the membership degree of δ to the fuzzy subset A indicates the degree to which δ belongs to A.

[0128] According to the lithology profile data of the study area, the lithology of the study area is divided into four categories: clay, silt, fine medium sand, medium coarse sand and above. The vertical permeability coefficients of the four lithologies are given according to the average value of the general experience of the region. After standardization, the permeability evaluation values ​​Q of the four lithologies are obtained. i (i=1,2,3,4:Q i ∈A); the river channel is divided into n sections according to the generalized shallow aquifer permeability, and the area S of the four types of lithology in the profile ni Calculate the proportion P of the four types of rock formations in the nth segment by the total area Sn of the segment ni :

[0129]

[0130] Then calculate the comprehensive permeability evaluation value Q of the nth rock layer n :

[0131] Q n =∑P ni Q i ;

[0132] According to relevant data, the maximum and minimum infiltration rates of the river channel are max , min , establish the river infiltration rate and permeability evaluation value Q n The linear equation between:

[0133]

[0134] Solve the equation to get the segmental infiltration rate λ of each river section n .

[0135] (3) Analysis of calculation results of infiltration rate of river channel recharge segment in the study area

[0136] The Yongding River has the strongest infiltration capacity in the Yuanbo Lake-Lugouqiao section, with an infiltration rate of 0.8, which is relatively weak at 0.2 in the Tianlu Wuqing District-Beichen District section, and basically distributed between 0.4-0.7 in other areas. The Juma River has the strongest infiltration capacity in segment 1, with an infiltration rate of 0.6, which is relatively weak at 0.4 in segment 2, and basically at 0.47 in other areas. The Tang River has the strongest infiltration capacity in segment 2, with an infiltration rate of 0.6, which is relatively weak at 0.41 in segment 4, and basically at 0.48 in other areas. The Dasha River has the strongest infiltration capacity in segment 1, with an infiltration rate of 0.53, which is relatively weak at 0.43 in segment 4, and basically at 0.48 in other areas. It can also be seen from the calculation results that the infiltration capacity in the front of the mountain area is relatively strong, and begins to weaken in the plain area. The eastern part of the plain is basically the area with the lowest infiltration capacity. This feature is also consistent with the sedimentary characteristics of alluvial fan and alluvial lake plain strata.

[0137] Table 2 Infiltration rate of river sections

[0138]

[0139] 4. Numerical model of groundwater flow in typical areas

[0140] (1) Mathematical model

[0141] According to the hydrogeological conceptual model of the study area, the well-defined problem can be described by the following differential equation:

[0142]

[0143] Where: Ω-seepage area; h-water level of groundwater system (m); K-horizontal permeability coefficient of water-bearing medium (m / d); K z -vertical permeability coefficient of water-bearing medium (m / d); ε-source and sink term of aquifer (1 / d); h0-initial water level distribution of the system (m); s-water storage rate of aquifer below free surface (1 / m); Γ 0- The upper boundary of the seepage area is the free surface of the groundwater; μ-gravity water supply of the phreatic aquifer on the phreatic surface; p-evaporation and precipitation infiltration intensity of the phreatic surface (m / d); Γ1-known water level boundary; h1-known boundary water level value (m); Γ2-flow boundary of the seepage area; K n -permeability coefficient in the normal direction of the boundary surface (m / d), n-normal direction of the boundary surface; q-flow rate per unit area of ​​Γ2 boundary (m / d), inflow is positive, outflow is negative, and the water-blocking boundary is 0.

[0144] (2) Model structure

[0145] ① Spatial Discretization

[0146] According to the aquifer structure, boundary conditions and groundwater flow characteristics of the study area, the Daqing River Basin was divided into 2km×2km grids, of which there are 5031 effective cells in the first layer and 4015 effective cells in the second and third layers.

[0147] ② Determination of simulation period

[0148] According to the project data collection and the actual situation of groundwater measurement, the simulation period of the Daqing River numerical simulation model is from June 2010 to December 2019, with one stress period per month, a total of 114 stress periods, and the time step is set to 1 day. The intensity of external source and sink terms remains unchanged in each stress period. The flow field in June 2010 is used as the initial flow field, and the flow field in December 2019 is used as the fitting verification flow field.

[0149] ③Determined solution condition processing

[0150] Initial conditions: The Daqing River model uses the groundwater level data in June 2010, and uses interpolation and extrapolation methods to obtain the initial water levels of shallow and deep groundwater. There is less water level data that can distinguish the third layer, so its initial flow field is set the same as the second layer. Boundary conditions: The inflow and outflow of each boundary calculated in Chapter 4 are brought into the model, and appropriate adjustments are made through flow field fitting.

[0151] ④ Source and sink item processing

[0152] The composition of source and sink items in the study area is relatively complex, the amount of data processing is large, and the data sources are widely provided in the form of GIS and database. The surface source data include calculation parameters such as precipitation infiltration coefficient and well irrigation regression coefficient. At the same time, the mining intensity and rainfall can also be given in the form of surface data; the linear data include boundary inflow and river infiltration; the mining is given in the form of point expansion to surface.

[0153] Area data: Convert the area source and sink data into intensity discretization to the grid, then superimpose the calculation, convert it to the well flow on the grid, generate a GMS external well file, and bring it into the model. Linear data: First topology to the center of the grid, convert it to a single grid well flow, add it to the GMS external well file and bring it into the model. For example, if a river is set as a flow boundary, the relationship between the river and groundwater can be characterized by injection or pumping wells. Point data: Similar to linear data, first topology to the center of the grid, and then bring it into the model calculation in the form of a pumping well or injection well.

[0154] ⑤Hydrogeological parameter processing

[0155] The parameter partition values ​​are brought into the model interface as graphic data and input into the model through the internal algorithm. By fitting the dynamic curves of the groundwater flow field and typical pores, various hydrogeological parameters of the aquifer are identified.

[0156] 5. Model Identification and Verification

[0157] After the model is run and calculated, the spatiotemporal distribution of groundwater levels under given hydrogeological parameters and equilibrium conditions based on the conceptual model can be formed. By fitting the actual flow field and the duration curves of typical groundwater monitoring wells, the specific quantities of hydrogeological parameters and source and sink items can be identified, so that the model can be combined with the hydrogeological conditions of the study area to more accurately and quantitatively study the recharge and discharge relationship in the simulation area and predict the changes in groundwater flow field and water equilibrium under given recharge and pressure extraction schemes.

[0158] (1) Flow field fitting analysis

[0159] Judging from the fitting of the shallow flow field, the simulated flow field can basically reflect the flow law and trend of groundwater, and the overall fitting is good. The fitting effect of the piedmont zone in the northwest is better, reflecting the flow direction and change characteristics of shallow groundwater. The fitting error in the southwest is slightly larger, but it can basically reflect the flow characteristics of groundwater. In the Baoding area, a groundwater drop funnel was formed due to the high mining intensity. To the northeast and southeast, in the Xiongxian-Bazhou-Dai and Gaoyang-Lixian-Qingyuan areas, due to the problem of groundwater over-exploitation, a large groundwater funnel was formed, and the simulation results can reflect the morphological characteristics of the flow field in the funnel area. Therefore, it can be said that the simulation effect of the shallow flow field in the study area is good and can reflect the actual groundwater flow field.

[0160] The overall deep flow field is also consistent with the measured flow field shape. The deep groundwater flows from west to east as a whole, and the flow characteristics are obvious. In short, the calculated flow field of the simulated study area can reflect the actual groundwater flow characteristics, and the fitting effect is good.

[0161] (2) Water level process line fitting analysis

[0162] According to the distribution of water level observation points in the model area, 158 observation points were selected, and the absolute and relative errors between the simulated water level and the observed water level at each fitting point were calculated for analysis.

[0163] There are 5600 data points in total at 158 ​​observation points, with an average error of 1.16m. Among them, there are 70 points with an error less than 1m, 61 points with an error between 1 and 2m, and 27 points with an error greater than 2m. The analysis results are shown in Table 3. It can be seen from the fitting that some of the fitting results are good, and the difference between the simulated water level and the actual observed water level is small, which can reflect the dynamic change trend of the water level at this point; the difference between the simulated water level and the measured water level is large, but the change trend is basically the same.

[0164] Table 3 Fitting error distribution table

[0165] Mean absolute error (m) ＜1 1-2 ＞2 Number of observation holes 70 61 27 Percentage of the total (%) 44.3 38.6 17.1 Number of shallow observation holes 55 49 21 Number of deep observation holes 15 12 6

[0166] The model grid area is 4km 2 The water level at the center of the grid represents the entire 4km 2 The changes in regional upper water levels can better characterize the movement patterns and trends of groundwater in the entire study area and meet the requirements of calculation accuracy.

[0167] (3) Identification and analysis of hydrogeological parameter fields

[0168] According to the analysis of the hydrogeological conditions in the study area, the hydrogeological parameter division is basically the same as the groundwater system boundary. The main hydrogeological parameters of the phreatic aquifer are: the permeability coefficient of the phreatic aquifer and the water supply degree of the vadose zone. The main hydrogeological parameters of the confined aquifer are: permeability coefficient and water release coefficient.

[0169] The hydraulic conductivity can be calculated using the following formula (for each calculation unit):

[0170]

[0171] Where: K 1，g , the permeability coefficient of the lithofacies corresponding to the hydrogeological unit; m 1，g , the horizontal length of the lithofacies corresponding to the hydrogeological unit in the calculation unit; 1, the lithofacies corresponding to each hydrogeological unit in the calculation unit; g, each hydrogeological unit in the calculation unit.

[0172] The numerical simulation model of groundwater in the North China Plain based on the block structure assigns different hydrogeological parameters to different lithofacies.

[0173] 6. Analysis of groundwater recharge effect

[0174] (1) Historical flow field changes

[0175] The shallow groundwater flow field in the plain area of ​​the Daqing River Basin in 1959, 1984, 2005 and 2019 and the deep groundwater flow field in 1980, 2005, 2010 and 2019 were mapped, forming 8 periods of high-precision stratified groundwater flow fields.

[0176] Comprehensive changes in shallow groundwater flow fields show that the shallow groundwater flow field in the Daqing River Plain area has changed from its natural state in the late 1950s. With the continuous increase in the development and utilization of groundwater, local groundwater funnels have gradually formed. Around the 1980s, some small groundwater funnels had appeared in the central and eastern parts of the plain. The mining areas were mainly distributed in Gaoyang, Lixian, Qingyuan and Xiongxian. The groundwater level line dropped to about 0m, and the front part of the mountain was basically able to maintain its natural form. By the end of the 20th century, the amount of groundwater extraction increased sharply. In the eastern region, local small funnels had developed into composite funnels, and the groundwater level in the funnel area generally dropped to about -20m. By 2019, with the implementation of various over-exploitation management measures, the groundwater flow field has changed further, and some areas have rebounded, but the water level in most areas is still low.

[0177] (2) Flow field changes after water replenishment

[0178] The effect of river replenishment was evaluated and analyzed by comparing the groundwater flow field state and water balance composition in the plain area of ​​the Daqing River Basin before and after river replenishment. The following four aspects were evaluated specifically:

[0179] ① Changes in flow field morphology

[0180] The change in flow field morphology is the most direct manifestation of the impact of river recharge on groundwater. This work will illustrate the impact of water recharge on the groundwater flow field morphology by comparing the dynamic changes of the flow field in the recharge area in different time periods.

[0181] ② Water level fluctuation in river channel

[0182] Carry out variation analysis on the water level before and after water replenishment in the river, and draw a variation diagram to illustrate the change of water level before and after water replenishment.

[0183] ③ Impact scope of backfill

[0184] According to the model fitting accuracy, a water level variation range of 0.5m was selected to simulate the water level changes in the study area under "water replenishment" and "no water replenishment" conditions, and to determine the range of water level rise caused by river water replenishment under the same other hydrogeological conditions.

[0185] ④ Balance analysis

[0186] The groundwater balance in the study area under the two conditions of "water replenishment" and "no water replenishment" was analyzed, and whether water replenishment can effectively restore the groundwater level was judged from the perspective of equilibrium changes.

[0187] The overall flow line of the Juma River has tightened, with greater changes in the north of Zhuozhou City and more obvious changes in the water level in Rongcheng County. As time goes on, the water flow gradually expands to both sides of the river. By the end of the water replenishment period, the flow field is basically stable, with an overall recovery range of 1-2m.

[0188] The water level in the Tanghe River channel from the eastern part of Dingzhou City to Qingyuan County has risen significantly, and the water level has risen to varying degrees during the recharge process. One year after the recharge, the water level has stabilized year by year, and the water level has risen by about 1m.

[0189] The areas where Dashaha River is most significantly affected by river recharge are in the middle and lower reaches of the river. From Shenze County to Gaoyang County in the east, as water recharge continues, the groundwater level in this section gradually rises, and the overall water level rises by about 1m.

[0190] VII. Prediction and analysis of the effect of comprehensive management of groundwater overexploitation

[0191] (1) Backfill and pressure mining plans

[0192] According to the "Action Plan for Comprehensive Management of Groundwater Overexploitation in North China", the long-term goal is to be achieved by 2035. The management measures are mainly to reduce overexploitation and replenish water from rivers and lakes. The replenishment volume in 2022 is 50-80 million m3 of Juma River. 3 , Tanghe River replenished 70-90 million m 3 , waterfall river 0.6-0.8 billion m 3 , Dasha River 0.8-100 million m 3 The Hebei region has achieved an overall reduction in groundwater extraction of 1.66 billion m 3 .

[0193] (2) Scenario Forecast

[0194] ①Water level changes

[0195] The above data were brought into the model to extend the simulation period to 2035. The rainfall was set according to the 10-year recurrence pattern, and other conditions remained unchanged. After simulation, the groundwater flow field in the study area in 2035 was obtained. By 2035, the groundwater will flow from northwest to southeast as a whole. There will be no large funnels in the study area, but the overall water level in Suning County is relatively low.

[0196] Except for Renqiu, Dacheng and Wen'an in the east, Xingtang, Xinle, Zhengding and Wuji in the south, and Suning and Hejian in the southeast, the water level in the study area has rebounded as a whole compared with 2017. The rebound in most areas can reach more than 5m. In Laishui, Xiong'an New Area, Fangshan, Mancheng, Gaoyang and Lixian, the water level has rebounded by more than 20m. The regional water level rebound effect is obvious.

[0197] (2) Scope of impact

[0198] Through continuous water replenishment and groundwater pressure extraction, the total area of ​​water level rise of 0.5m will reach 14403.84km2 by 2035, accounting for 72% of the total area of ​​the study area. The rise in most areas in the central and northern regions basically reaches more than 4m.

[0199] The water levels in Dacheng, Hejian and Suning in the east, which only adopted the policy of reducing water supply but did not implement recharge, have not been effectively restored. The water levels in Xingtang, Xinle, Gaocheng and Zhengding in the middle and upper reaches of Dasha River, where the recharge conditions are poor, have also failed to recover.

[0200] (3) Equilibrium analysis

[0201] According to the forecast analysis of groundwater balance in the study area from 2022 to 2035, the study area is in positive balance as a whole, with a recharge of 28.61×108m 3 / a, river recharge 4×108m 3 / a, accounting for about 14% of the total recharge. The difference in recharge is nearly 6×108m 3 / a, the overall water level rose sharply. However, some areas are still in a negative balance state. These areas include Anping County, Anguo City, Xushui County, Suning County, Gaobeidian City, Hejian City, Dingxing County, Gu'an County, Xian County, Xiong County, Gaocheng City, Raoyang County, Boye County, Baoding City and Zhengding County. Most counties and cities are in the range where rivers cannot be recharged and the recharge conditions in the middle and upper reaches of Dasha River are poor.

[0202] Table 4 Groundwater balance in the study area from 2022 to 2035

[0203] Specific embodiment 3:

[0205] Simulation and evaluation of deep water source replacement

[0206] The implementation of water replenishment through the Beijing-Hangzhou Grand Canal during the spring irrigation period provides important guarantees for agricultural irrigation in deep over-exploited areas along the route, such as Cangzhou and Hengshui. At the same time, part of the water diversion can replace groundwater as a source of agricultural irrigation, reducing groundwater exploitation. By constructing a water source replacement model for typical water receiving areas such as Cangzhou and Hengshui, a technical basis is provided for the evaluation of deep groundwater replacement effects.

[0207] 1. Construction of groundwater numerical model

[0208] (1) Generalization of aquifer systems

[0209] ①Vertical generalization

[0210] Cangzhou: This simulation covers the entire Cangzhou area, with a total of about 14,300 km 2 Based on actual drilling data, geological information and past experience, this study generalized the aquifer into four aquifer groups, from top to bottom, namely, the I (bottom plate elevation -10m to -60m), II (bottom plate elevation -140m to 210m), III (bottom plate elevation -310m to -450m), and IV (bottom plate elevation -400m to 650m) aquifers. Among them, the I and II aquifers are shallow aquifers, and the III and IV aquifers are deep aquifers.

[0211] Hengshui City: The simulation range of the groundwater numerical model of Hengshui City is the entire territory of Hengshui City, a total of about 8836km 2 Based on actual drilling data, geological information and past experience, this study generalized the aquifers into four aquifer groups, namely, aquifers I, II, III and IV from top to bottom. Among them, aquifers I and II are shallow aquifers, and aquifers III and IV are deep aquifers.

[0212] ② Generalization of boundary conditions

[0213] According to the initial distribution characteristics of isowater levels in the study area, the boundary is generalized as a universal head boundary to depict the recharge or discharge boundary; when the isowater level is perpendicular to the boundary, it is set as a water-isolating boundary.

[0214] ③ Initial conditions

[0215] January 1, 2019 was taken as the initial moment of the model, and the initial groundwater level contour map of the study area was drawn using Kriging interpolation.

[0216] (2) Space-time Discreteness

[0217] ① Spatial Discretization

[0218] Cangzhou City: 1km×1km grid is used in the horizontal direction, and the vertical direction is divided into 4 layers according to the aquifer structure. The model has a total of 56,432 activated grids, with 14,108 grids in each layer.

[0219] Hengshui City: 1km×1km grid is used in the horizontal direction, and it is divided into 4 layers vertically according to the aquifer structure. The model has a total of 35,632 activated grids, with 8,908 grids in each layer.

[0220] ②Time Discrete

[0221] The groundwater level observation data at the Dian scale from January 2019 to June 2023 were collected, of which the daily groundwater level observation data were collected from March to June 2023. Taking into account the intra-annual and inter-annual variations of groundwater levels in the study area, the simulation period was determined to be from January 2019 to June 2023, mainly including the calibration period and verification period of the model.

[0222] A total of 160 stress periods are divided from January 2019 to the end of June 2023, and each scale is a stress period. Since the second Beijing-Hangzhou Grand Canal water replenishment time is from March to June 2023, the time step is set to 1 day from March to June each year. Among them, January 2019 to February 2023 is the model's rate period, and March to June 2023 is the model's validation period.

[0223] (3) Hydrogeological parameters

[0224] Cangzhou City: The main hydrogeological parameters involved in the model are the permeability coefficient, water supply degree and storage coefficient of the aquifer. According to the hydrogeological conditions of the simulation area and the lithology of the aquifer, the hydrogeological parameters of the aquifer are divided into zones. Finally, the hydrogeological parameters that meet the conditions of the study area are determined through the calibration and verification of the model. The model is divided into four aquifer groups, with a total of 24 zones.

[0225] Hengshui City: The groundwater level observation wells and hydrogeological parameters in Hengshui City are divided into 10 zones.

[0226] (4) Evaluation of deep groundwater replacement effect

[0227] ① Overall assessment of the impact of water source replacement on groundwater level

[0228] The Beijing-Hangzhou Grand Canal will start full-line water replenishment on March 1, 2023. As of June 1, the water diversion volume of Cangzhou City (Cangxian, Botou City, Qingxian, Wuqiao County, Huanghua City, Dongguang County and Nanpi County) and Hengshui City (Jingxian and Gucheng County) reached 262.65 million m 3 , of which Cangzhou City diverted 181 million m 3 , Hengshui City diverts 81.65 million m 3 .

[0229] Early March to the end of May is an important period for wheat irrigation in Cangzhou and Hengshui. March to mid-April is the wheat greening and seedling pulling period, and mid-April to the end of May is the wheat filling period. Most of the agricultural irrigation areas in Cangzhou and Hengshui have carried out two agricultural irrigations in succession.

[0230] As of June 1, the irrigation area of ​​Cangzhou and Hengshui that have implemented water source replacement is 751,600 mu, of which Cangzhou has replaced 552,600 mu and Hengshui has replaced 199,000 mu. According to relevant data, the water consumption for spring wheat irrigation is about 55-90m3 / mu, most irrigation areas have completed two irrigations, and the calculated groundwater replacement volume is about 113.877 million m 3 The replacement water volume accounts for 43% of the diverted water volume, of which Cangzhou’s replacement water volume is about 79.515 million m 3 The water displacement in Hengshui is about 34.362 million m 3 The statistical table of deep groundwater replacement water volume in each county is shown in the table below.

[0231] Table 5 Statistics of irrigation area and groundwater replacement volume of the Beijing-Hangzhou Grand Canal water replenishment replacement

[0232]

[0233] ②Evaluation of the impact of water source replacement on groundwater level at different time periods

[0234] To further evaluate the effect of implementing water source replacement, this study compared the groundwater levels under the two schemes of implementing water source replacement and not implementing water source replacement at different times on April 1, May 1 and June 1, and refined the impact of the water replenishment process on deep groundwater replacement since March. According to the difference in groundwater levels in the third aquifer of Cangzhou City and Hengshui City when water source replacement was implemented and not implemented at different times, the water level difference between the two schemes was continuously increased, indicating that the implementation of water source replacement can effectively alleviate the rate of groundwater level decline under the existing water use conditions. Continuous water source replacement reduces the amount of deep groundwater extraction, which plays a positive role in alleviating regional water pressure and restoring groundwater levels.

[0235] The water level changes in each receiving county under the two scenarios of Beijing-Hangzhou Grand Canal through-filling and without water source replacement from March 1 to June 1, 2023 were statistically analyzed, as shown in Table 6. The Beijing-Hangzhou Grand Canal through-filling has a significant impact on the deep groundwater level during the water replenishment period. Among them, if water source replacement is not implemented, the water level drop in most of the water source replacement counties in Nanpi County and Botou County of Cangzhou City will increase by more than 2m; the water level drop in most of the water source replacement counties in Jing County and Gucheng County of Hengshui City will increase by more than 6m.

[0236] Table 6 Statistics of water level changes in each receiving area and county under the two scenarios of implementing and not implementing water source replacement (unit: m)

[0237]

[0238] Research conclusions and application effects of the present invention:

[0239] (1) The study constructed an impact assessment model for the water replenishment of the Beijing-Hangzhou Grand Canal, and compared the water level changes of each river section without water replenishment and with actual water replenishment. A multi-period high-precision stratified groundwater flow field was drawn and the flow field evolution was analyzed. Based on the water balance method, the infiltration rate of each river section was calculated using monitoring data such as rainfall, evaporation, and water withdrawal. The cumulative infiltration volume was 144 million m 3 The average infiltration rate is 13%. According to preliminary analysis, the rock particles in the front area of ​​the North Canal are coarse and suitable for infiltration, the South Canal has more clay soil and weaker permeability, the Wei Canal has better permeability, and the main canals of the Small Canal are mostly lined and have weak permeability. The maximum impact distance of groundwater on both sides of the water replenishment main line is 5.11km, the average impact distance is 1.56km, and the cumulative impact area is 990.88km 2 The average impact distances of the North Canal, South Canal, Wei Canal and Small Canal on groundwater are 1.14km, 1.39km, 2.04km and 0.76km respectively. Since March to May is the agricultural irrigation period along the coast, the groundwater level is still mainly declining. Ecological water replenishment has slowed down the decline of groundwater levels around the river to a certain extent.

[0240] (2) Based on T-PROGS, U-Grid, and fuzzy mathematics methods, a recharge model for the plain area of ​​the Daqing River Basin was constructed, and the groundwater level recovery and water balance changes under different recharge and pressure extraction scenarios were simulated and predicted. By comparing the flow field and balance differences before and after the water replenishment of the Juma River, Tang River, and Dasha River, the effect of river recharge was analyzed and predicted. According to the Action Plan, through the continuous implementation of water replenishment and pressure extraction, by 2035, the water level will rise by 0.5m in areas accounting for 72% of the total area of ​​the study area, and the recovery range in most areas in the central and northern parts will basically reach more than 4m.

[0241] (3) Taking Cangzhou and Hengshui as typical areas, a numerical model of deep water source replacement in the water receiving area of ​​the Beijing-Hangzhou Grand Canal was constructed, and the replacement effect was evaluated. The simulation results show that by reducing the amount of deep groundwater extraction through water source replacement, the rate of groundwater level decline during the spring irrigation period can be effectively alleviated and reduced under the existing water use conditions. The impact of water source replacement on groundwater levels in different periods from March to June was evaluated. According to the simulation results, if water source replacement is not implemented, the water level drop in most water source replacement counties in Nanpi County and Botou County of Cangzhou City will increase by more than 2 meters; the water level drop in most water source replacement counties in Jing County and Gucheng County of Hengshui City will increase by more than 6 meters.

[0242] The present invention is described above by way of example in conjunction with the accompanying drawings. It is obvious that the specific implementation of the present invention is not limited to the above-mentioned method. As long as various improvements are made using the method concept and technical solution of the present invention, or are directly applied to other occasions without improvement, they are all within the protection scope of the present invention.

Claims

1. A coupled simulation method for groundwater recharge in rivers and lakes, characterized in that: It includes the study of the impact of recharge on the riverbank zone of the through-flow recharge, the simulation of recharge in the plain area and the simulation evaluation of deep water source displacement; The study on the impact of recharge on the riparian zone of the through-flow replenishment project is to conduct the riparian zone infiltration analysis and riparian zone recharge impact analysis after the riparian zone hydrogeological conceptual model, riparian zone numerical model and model identification verification; The recharge simulation in the plain area is based on the three-dimensional hydrogeological structure model of T-PROGS, the river grid encryption technology based on U-Grid, and the river infiltration analysis based on fuzzy mathematics methods. The numerical model of groundwater flow in typical areas is established to analyze the effect of groundwater recharge and the effect of groundwater over-exploitation control. The analysis of the effect of groundwater recharge includes historical flow field changes and flow field changes after water recharge. The prediction analysis of the effect of groundwater over-exploitation control includes water level changes, impact range and equilibrium analysis. The simulation evaluation of deep water source replacement first constructs a groundwater numerical model, including the generalization of the aquifer system and hydrogeological parameters, and then evaluates the effect of deep groundwater replacement, including an overall evaluation of the impact of water source replacement on groundwater levels and an evaluation of the impact of water source replacement on groundwater levels at different times.

2. A river and lake groundwater recharge coupling simulation method according to claim 1, characterized in that: The hydrogeological concept model of the riparian zone is divided into two layers. The first layer of the model is the phreatic aquifer, where groundwater extraction is mainly mixed extraction; the second layer of the model is the confined aquifer, where groundwater extraction is small and there is a lack of corresponding water level monitoring data; the groundwater recharge items in the area include rainfall infiltration recharge, river seepage recharge, and agricultural irrigation infiltration recharge. Agricultural irrigation infiltration recharge includes canal leakage, canal irrigation infiltration and well irrigation regression; the groundwater discharge items in the area are phreatic evaporation discharge and groundwater extraction.

3. A river and lake groundwater recharge coupling simulation method according to claim 2, characterized in that: In the riparian zone numerical model: Groundwater flow can be described by the following equation: Among them, K xx , K yy , K zz are the components of the permeability coefficient in the x, y, and z directions (LT -1 ); h is the water level of the aquifer (L); W represents the source and sink term (T -1 );S s For unconfined aquifers, it is the water supply rate; for confined aquifers, it is the water storage rate (L -1 );t is time (T).

4. A coupled simulation method for groundwater recharge in rivers and lakes according to claim 3, characterized in that: The model identification and verification follows the following principles: the hydrogeological parameters in the model can reflect the hydrogeological conditions in the study area; the final groundwater flow field during the simulation period is basically consistent with the actual groundwater flow field during that period; the simulated groundwater level process line is basically similar to the measured water level process line; The simulated groundwater balance is consistent with the actual water balance status.

5. A river and lake groundwater recharge coupling simulation method according to claim 4, characterized in that: The analysis of riparian zone infiltration is based on the water balance method, using monitoring data such as rainfall, evaporation and water withdrawal to calculate the infiltration rate of each river section.

6. A coupled simulation method for groundwater recharge in rivers and lakes according to claim 5, characterized in that: The impact analysis of riparian zone recharge uses a model to simulate the groundwater level distribution in each river section without water recharge, and compares it with the measured groundwater level under the current water recharge situation.

7. A river and lake groundwater recharge coupling simulation method according to claim 6, characterized in that: The numerical model of groundwater flow in a typical area uses the following differential equation to describe the well-defined problem: Where: Ω - percolation area; h-water level elevation of the groundwater system (m); K-horizontal permeability coefficient of water-containing medium (m / d); K z -Vertical permeability coefficient of water-containing medium (m / d); ε - source and sink term of aquifer (1d); h0-initial water level distribution of the system (m); S-water storage rate of aquifer below free surface (1 / m); Γ0 - upper boundary of the seepage zone, i.e., the free surface of groundwater; μ - gravity water supply of the unconfined aquifer at the water table; p-evaporation and precipitation infiltration intensity of the water table (m / d); Γ1-known water level boundary; h1-known boundary water level value (m); r2 - flow boundary of the seepage region; K n - Permeability coefficient in the normal direction of the boundary surface (m / d), n- normal direction of the boundary surface; q-Γ2 is the flow rate per unit area (m / d), which is positive for inflow, negative for outflow, and 0 for water-blocking boundaries.

8. A method for coupling simulation of groundwater recharge in rivers and lakes according to claim 7, characterized in that: The flow field changes after water replenishment are specifically evaluated in the following four aspects: Changes in flow field morphology: Compare the dynamic changes of the flow field in the water replenishment area at different time periods; Variation of water level in river channel: Carry out variation analysis on water level before and after water replenishment in river channel, and draw variation chart to illustrate the variation of water level before and after water replenishment; Impact range of backfill: According to the model fitting accuracy, the water level variation range of 0.5m was selected to simulate the water level changes under the conditions of "replenishment" and "no refill" in the study area, and the range of water level rise caused by river water replenishment was determined under the same other hydrogeological conditions; Equilibrium analysis: Analyze the groundwater equilibrium in the study area under the two conditions of "water replenishment" and "no water replenishment", and judge whether water replenishment can effectively restore the groundwater level from the perspective of equilibrium changes.

Images