(Shandong Institute of Geology and Mineral Resources, Jinan 2500 14)

About the author: Wang Shaojuan (1976-), female, engineer, mainly engaged in research on hydrology, geothermal, environment and disaster geology.

From the establishment of MODFLOW model, the generalization of aquifer and boundary conditions in calculation area, the treatment of recharge and discharge items, the treatment of model calculation items, model verification and other aspects, this paper analyzes the reliability and credibility of the model, and puts forward two different optimization schemes for Jinan spring area under the condition of ensuring spring water supply.

Keywords: MODFLOW mode; Establish; Inspection; Optimization; App application

Jinan Baoquan water supply has always been the focus of attention from all walks of life. Shandong Institute of Geology and Mineral Resources has been engaged in hydrogeology research of Jinan spring area for many years, especially the dynamic observation of groundwater level, which has laid a solid foundation for us to study Jinan spring area. In addition, on the basis of many years' research on the hydrogeological conditions of Jinan spring area, combined with the dynamic data of water level in the large-scale pumping test in the west of Jixi in 2003, the water resources exploitation potential of Jinan spring area is optimized by using MODFLOW modular three-dimensional finite difference groundwater flow model.

Brief introduction of 1 MODFLOW model

MODFLOW modular three-dimensional finite difference groundwater flow model is a set of software developed by McDonald and Harbaugh of US Geological Survey in 1980s, which is specially used for numerical simulation of three-dimensional finite difference groundwater flow in porous media. It is characterized by modular program structure, simplified discrete method and diversified solution methods.

In this work, the 5. 17 model is adopted for the numerical simulation of groundwater flow in Jinan spring area.

2 Jinan spring hydrogeological conceptual model

2. 1 simulation range

On the plane, Mashan fault is the boundary in the west; The northern part is bounded by Jinan rock mass and Carboniferous and Permian strata, including: Wu Dong fault-Liu Changshan section, bounded by -350m limestone roof elevation due to magmatic intrusion; On the first line of Weilizhuang-Mashan fault, according to the research results of coalfield geological exploration data, the contact zone between Ordovician and Carboniferous and Permian at the elevation of -350 m ~-400 m is the boundary; The east is bounded by Wu Dong fault; It is bounded by the bottom of Zhangxia Formation in the middle of Cambrian in the south, starting from Gangxinzhuang in the west, Song Cun, Huanggujing, Xiaogushan, Dagushan, Wohushan, Hanjiazhuang and Dongjiazhuang in the east, reaching Xiaofo Temple in the south and Wu Dong fault in the east. The simulated area is 847.5km2.

Raise the groundwater level from the vertical simulation range to the lowest horizontal elevation depth-160m. Because the limestone aquifer in the north of the study area is covered by igneous rocks, considering the depth of karst development, the vertical minimum calculated depth is -423 meters, while in the south of the study area, the maximum calculated thickness is 400 meters.

2.2 Generalization of aquifer structure

In this model, the aquifer in the study area is divided into two layers, which are summarized as follows:

The first layer is phreatic aquifer, and the target layer is mainly Quaternary Holocene and Pleistocene strata. The lithology of the aquifer is medium-coarse sand and gravel, and the total thickness of the aquifer is 10~40m ~ 40m, which is simplified as two-dimensional flow.

The second layer is confined water aquifer, and the target layers are Cambrian Zhangxia Formation, Fengshan Formation and Ordovician strata. The vertical seepage velocity component of fractured karst water in the model area is very small and can be ignored, which can be summarized as two-dimensional water flow. However, considering that there is overflow recharge between the first aquifer and the second aquifer, it is generalized as three-dimensional flow. K(h-e) = t is approximately regarded as a constant. Therefore, in the calculation, the south fractured karst phreatic area is also treated as confined water flow.

2.3 boundary conditions

The boundary conditions of the study area can be summarized as water blocking boundary and runoff boundary, that is, the second kind of boundary. The aquifer in the study area is divided into two layers.

The first layer is the phreatic aquifer, and its boundary is summarized as follows:

East boundary: it is bounded by the eastern edge of the alluvial fan of Yufu River. It can be summarized as a water-resisting boundary.

Southern boundary: the top of alluvial-diluvial fan is water-proof boundary.

Western boundary: it is bounded by the western edge of the alluvial-diluvial fan of Beishahe River, which is basically parallel to Mashan fault and perpendicular to the direction of pore diving, and is a water-resisting boundary.

Northern boundary: Because of the Quaternary gravel aquifer, it extends far to the north, and the northern calculation boundary is a permeable drainage boundary.

The second layer is confined aquifer, and the boundary is simplified as:

North boundary: the central and eastern parts are intrusive igneous rocks, the limestone roof is buried deeply, and the rock mass forms obvious water-resisting layer; Only the limestone roof in the west of the boundary is buried shallowly, and the fracture karst develops and extends outside the boundary, which is generalized as a water-resisting boundary.

Eastern boundary: Wu Dong fault is the boundary, in which Xiagelao-Xujiazhuang section is a water-blocking fault and a water-resisting boundary; The Xujiazhuang-Dashuipo section is generally blocked by water, but it is weakly permeable in the northern part of the refinery and local areas of the block plant, which is reduced to a flowing boundary.

Southern boundary: reduced to zero flow boundary. The calculation boundary is taken from the bottom interface of Zhangxia Formation of Middle Cambrian. In the south, the shale of Xuzhuang Formation, Maozhuang Formation and Shantou Formation is sandwiched with thin limestone and old strata of Taishan Group, so the southern boundary is generalized as a water-resisting boundary.

West boundary: Mashan fault is the boundary, and the section south of Xiguan in Changqing is the water-blocking fault as the water-resisting boundary; In Changqing-Qianlong section, the pumping test shows that the water levels on both sides of the fault are basically the same and have strong permeability, so it is treated as a flow boundary. Because the changxiao water source is far away from the spring group, Mashan fault can be used as a relative boundary.

Internal boundary: Qianfoshan fault runs through the whole simulation area from south to north, with waterproof boundary in the south and permeable boundary in the north of Nanjiao Hotel.

See Figure 1 and Figure 2 for the generalization of the boundary geometric conditions of the study area.

Figure 1 Calculation range of the first floor

Figure 2 Calculation range of the second layer

2.4 Treatment of groundwater recharge and discharge projects

2.4. 1 groundwater exploitation

Under the current conditions, groundwater exploitation in the study area is mainly urban tap water, self-owned mines by industrial and mining enterprises, irrigation water for rural agricultural life and water for township enterprises.

2.4.2 Infiltration and overflow recharge

Vertical infiltration and overflow recharge mainly include rainfall infiltration, surface water leakage and Quaternary pore water overflow recharge. According to the stratum structure, rock characteristics, topography and geomorphology, the parameters are partitioned, and the initial values are given respectively, which are identified and finally determined in the simulation process.

Transverse leakage

Lateral seepage refers to the lateral exchange quantity flowing into or out of the study area from the outside. In the process of simulation, the unit flow on each flow interface is given first, and then it is identified by simulation and finally determined.

After generalizing the hydrogeological conditions in the study area, the model can be summarized as a three-dimensional unsteady flow mathematical model of heterogeneous anisotropic unconfined-confined fractured karst.

3 MODFLOW model calculation

MODFLOW is mainly simulated by three-dimensional finite difference method. Its basic principle is that the flow of groundwater in three-dimensional porous media can be expressed by the following partial differential equation without considering the change of water density:

Shandong environmental geology collection

Where Kxx, Kyy and Kzz are the components of permeability coefficient in X, Y and Z directions, here we assume that the principal axis direction of permeability coefficient is consistent with the direction of coordinate axis, and the dimension is (lt-1); H is the head (l); W is the unit volume flow (T- 1), which represents the inflow or discharge of water; Sx is the water storage rate of porous media (L- 1).

Groundwater system often exchanges materials, energy and information with the environment, and the environment is always changing. In the interaction with the environment, all the elements of the aquifer (such as water level, water quantity, hydrochemical composition, water temperature, etc. ) changes with time. The change of groundwater elements is the result of the imbalance of water, salt, heat and energy budget in aquifer. For example, when the recharge of the aquifer is greater than the discharge, the water storage capacity increases and the groundwater level rises; On the contrary, when the recharge water is less than the discharge water, the water storage capacity decreases and the groundwater level drops. The factors affecting groundwater dynamics can be divided into two categories. One is the input of aquifer information by environment, such as precipitation, recharge of groundwater by surface water, artificial exploitation or recharge of groundwater, and the influence of ground stress on groundwater. The other is the factors that change the input information, which mainly involve the hydrogeological, geological and topographic conditions of groundwater.

There are many factors affecting groundwater dynamics, and the influencing mechanism is complex. In the numerical simulation of groundwater flow mechanics, all influencing factors must be fully considered, and the specific treatment is as follows:

(1) Division of calculation area. In this model, the total area of Jinan is divided into 80 rows, 1 15 columns and 9200 units, each unit is a 500m×500m square. There are 3390 units in the study area, with a total area of 847.5km2.

(2) Treatment of aquifer. When the MODFLOW model simulates groundwater flow, the total aquifer can be divided into the following two types: 0 indicates that the aquifer belongs to confined aquifer; 1 indicates that the aquifer is strictly pressureless.

(3) Treatment of boundary conditions. In the MODFLOW model, the boundary is treated as follows: for the active concentration battery, the assignment is1; For the inactive concentration battery, the assignment is 0; For units with constant concentration, the assignment is-1. In this simulation work, all units outside the study area are regarded as inactive units, and the assigned value is 0; For the units in the study area, the value is 1 according to the active unit, among which, because the south of T village of Qianfoshan fault in the area is a water-blocking fault, the value of this unit is also treated as an inactive unit, with the value of 0. For the above runoff boundary and recharge boundary, the model adopts pumping and water injection methods to generalize respectively.

(4) Division of time periods. MODFLOW provides users with the input of simulation time period, including time period, time period length, time step and other options. In the process of this model, different divisions are made according to the specific simulation situation, and the time unit is days.

(5) Initial water level. The initial water level of this simulation is input according to the actual water level elevation. According to the water level data of each long measuring point and observation point, the initial isobath is drawn first, and then the data is input into the model, and the unit is m.

(6) Precipitation replenishment treatment. The main source of groundwater recharge in the study area is atmospheric precipitation recharge. According to previous work and experience, the precipitation infiltration coefficients of each layer in the study area are assigned as follows: the infiltration coefficients of Cambrian Zhangxia Formation and Fengshan Formation are 0.33, the infiltration coefficients of Cambrian Gushan Formation and Changshan Formation are 0.00 1, the infiltration coefficient of Ordovician concealed area is 0.39, and the infiltration coefficient of Ordovician limestone exposed area is 0. The infiltration coefficient of Quaternary in the western suburbs is 0.3 ~ 0.32, in which the infiltration coefficient in this area is adjusted to 0. 15 due to the urbanization development of Changqing County, and the precipitation unit is converted into head height m.

(7) Treatment of wells. In MODFLOW, wells are divided into injection wells and production wells, in which injection wells are represented by positive numbers and production wells are represented by negative numbers, and the unit is m3/d. ..

(8) Agricultural development, agricultural supply and domestic water treatment. In this simulation process, the agricultural water consumption in the study area is treated as follows: (1) the formal expression of the Quaternary agricultural and domestic water consumption intensity in the western suburbs; Karst water for agriculture and domestic use is expressed in the form of mining wells in this simulation process, and agricultural irrigation regression is directly deducted from agricultural exploitation according to the coefficient of 0. 15 according to previous research results.

4 Numerical simulation and model updating

The whole significance of numerical simulation of groundwater flow system is that mathematical models can "truly" reflect objective entities to explain the past, explain the present and predict the future. Only the mathematical model and its solution process have been completed, and the mathematical simulation is not over yet. Identifying and testing model parameters is an essential and important step in mathematical simulation.

4. 1 model identification

The identification of model parameters is to solve the inverse problem, that is, to calculate hydrogeological parameters and verify the mathematical model according to pumping test data or actual production and water level data. Because of the ill-posed nature of inverse problem, it has long been a problem that modelers cannot solve well. Parameter partition is a key step in the process of parameter identification. According to the geological structure, rock characteristics and karst development degree of the simulation study area, combined with the long-term research results of the evolution law of fractured karst water in northern China, the simulation study area is divided into 13 horizontal permeability coefficient, 32 hydraulic conductivity coefficients, 16 water supply, 32 water storage coefficients and 18 cross-flow coefficient.

The fitting time of the model began on June 6th, 2003 and ended on June 6th, 2004. June 2003 ~ June 2004 is a typical hydrological year, and the annual precipitation in the study area is 870.8mm, which is a wet year. In addition, since September 6th, 2003, the 4th Da Chun delegation has made a full comeback. Choosing this hydrological year as the mathematical simulation stage can better reflect the dynamic changes of groundwater system in Jinan Spring Group Research Area. In model identification, * * * divided 72 time periods, and the length of each time period was divided according to the observation time of long measuring points, each time period ranged from 4 days to 6 days, and the time step was 1 day. In the whole fitting process, the actual observed water level in the observation well in the study area is used as the basis for model identification.

The parameters identified by the model are basically consistent with the results of hydrogeological investigation.

According to the water level observed on June 6, 2003, the initial flow field of karst water in Jinan spring area described by computer simulation (see Figure 3) basically reflects the actual situation of groundwater flow field in the research area at that time.

Fig. 3 Plan view of groundwater flow field at the initial stage of the second layer simulation.

Figure 4 Plan of Groundwater Flow Field on June 26th, 2003

Model identification results: During the numerical simulation period, there are 49 observation wells, of which 1 1 belongs to Quaternary diving observation wells in the western suburbs; Among the 38 observation wells, there are 23 karst water observation well in the western suburb, 9 observation well in the urban area, 3 observation well in the eastern suburb, and the remaining 3 observation wells are distributed in the southern region. According to the model calculation, the simulated flow field on June 26th, 2003 is shown in Figure 4, which shows that the mathematical model has high reliability.

See Figure 5, Figure 6 and Figure 7 for the calculation results of model identification in the study area.

Fig. 5 Fitting curve of J65 observation well

Fig. 6 fitting curve of baotuquan 1 observation well.

Figure 7 A3-4 Fitting Curve of Observation Well

During the simulation period, the total recharge of karst groundwater is greater than the total discharge, and accordingly, the water level of the observed well generally shows an upward trend. The above results are reflected in the simulation period, and the model is in good agreement with the actual situation.

4.2 Model test

After model parameters are identified, these parameters need to be tested. In this model test, the water level data of the pumping observation well in Jixi Pumping Station from June 6, 2003 to July 28, 2003 were used for model correction, which lasted for 52 days and was divided into 20 periods, the length of each period was divided according to the observation time of the pumping observation well in Jixi Pumping Station, and the time step was 1 day. In the model inspection stage, the observation wells are consistent with those in the model fitting stage. See Figure 8 ~ Figure 10 for the calculation results of model identification in the study area.

Fig. 8 Fitting curve of J65 observation well

Fig. 9 Fitting curve of Baotuquan 1

Figure 10 A3-4 Fitting Curve of Observation Well

Due to the model modification, during the pumping period in western Jilin Province, large-scale pumping tests were carried out in Qiaozili, Gucheng and Lengzhuang, and the water level dropped. In the model test, the total recharge of karst groundwater is less than the discharge, and accordingly, the water level of the observed well generally shows a downward trend. At the same time, due to the continuous precipitation in Jinan after mid-July, the groundwater in the spring area has been replenished and the water level has risen obviously, which is in good agreement with the model inspection process, further explaining the accuracy of the model.

4.3 model reliability analysis

As can be seen from the figure, the simulation fitting effect is good, and the points with large errors generally appear between 1 1 and 1. The reason is mainly caused by the model itself and the structure of the observation well. At this stage, the rainy season of one year has just passed, the actual groundwater has been replenished to a certain extent, and the water level is slowly decreasing. However, during the simulation period, due to insufficient precipitation, the precipitation replenishment of the model.

Generally speaking, the fitting error of karst water level in the western suburbs from June 6 to July 27 is slightly larger, but the error effect is still satisfactory. The error analysis is as follows: the simulation time of this model revision is from June 6, 2003 to July 28, 2003. During this period, due to the increase of water production in the ancient city in the western suburbs, Lengzhuang and Qiaozili 19.2× 104 m3/d, it was generalized as a mining mine in the model, and the simulated water level of the observation well near the water plant suddenly dropped.

Generally speaking, the model reflects the changing trend of groundwater flow field in the study area, which is in good agreement with the actual situation and the simulation results are reasonable, indicating that the model can be used for groundwater dynamic prediction.

Prediction and analysis of numerical model of groundwater flow

The start time of this forecast is June 17, 2004, and the target time is May 3, 20091day, * * * 18+00.

5. 1 Determination of calculation scheme

The model calculation needs input data, including precipitation, groundwater exploitation and hydrogeological parameters in the future.

5. 1. 1 Determination of precipitation

Precipitation is given by the predicted value. In the forecast precipitation series from 2004 to 2009, according to the trend of "four droughts and one abundant" and the precipitation from1999 to 2003, the overall trend is given. See table 1 for the division of specific time periods and precipitation.

Table 1 forecast period and forecast precipitation

5. 1.2 Determination of agricultural exploitation, agricultural supply and domestic water

Quaternary agricultural exploitation and domestic water are represented by evaporation intensity, and agricultural exploitation and domestic water in karst areas are represented by mining wells, and the agricultural replenishment coefficient is taken as 0. 15. In the forecast period, the agricultural exploitation amount is determined by comparing with the precipitation in previous years and combining with the agricultural exploitation amount in that year.

5. 1.3 Determination of groundwater resources exploitation scheme in the study area

At present, there are two methods to evaluate the allowable exploitation of groundwater by numerical method. One is to calculate the exploitation amount according to the depth of water level decline; The second is to give the mining plan of each well in advance, calculate the depth of water level decline, and then adjust it constantly to find a suitable plan. The first method is adopted in this simulation, that is, according to the water level constraint, the groundwater exploitation amount in Jinan city is controlled, and the exploitation amount is formulated with full consideration of the future water supply potential and water demand requirements.

5.2 Analysis of forecast results

The current layout of the first scheme is optimized without considering the supplementary source.

The first plan is based on increasing water plants in the western suburbs and reducing the number of industrial self-provided wells in the eastern suburbs. The calculation result is as shown in figure 1 1, and the total exploitation amount is 17.806× 104m3/d, of which Emei Mountain is 1.0× 104m3/d and dayangzhuang is 0.8×/d.

It can be seen from the calculation results that the groundwater level of Baotu Spring is 27.5m higher than the critical water level in the whole prediction period, and the four spring groups will be able to spew all the year round.

Figure 1 1 baotuquan 1 observation well prediction curve

Figure 12 Prediction Curve of Baotuquan 1

The second scheme is to replenish the source and optimize the mining.

Because human activities have changed the hydrogeological conditions, the groundwater recharge has decreased. In the second scheme, it is considered to optimize mining after recharging, and make full use of Jinan passenger water resources and water storage in southern mountainous areas. In addition to the precipitation infiltration recharge, this recharge is designed to recharge yufu river1/kloc-0 /×104m3/d, and Beishahe river is 9.0× 104m3/d, increasing the groundwater recharge by * * * 25×104m3.

According to the calculation results, 19× 104m3/d, 6.0× 104m3/d, 2.0× 104m3/d and 0.5 were mined in the western suburbs, Emei, dayangzhuang, Lashan, urban areas and eastern suburbs respectively.

From the calculation results, it can be seen that the water level of Baotu Spring is controlled above 27.5m, and the four spring groups can ensure perennial gushing. At the same time, although the exploitation in the western suburbs has increased, the area of underground funnel has not increased obviously due to the daily supply of 2.0× 105m3 from Yufu River and Beisha River.

In short, through the comparison of the above two mining schemes, it can be seen that recharge is a more effective engineering measure to keep the spring water gushing all the year round in the year of dry precipitation.

6 conclusion

Compared with previous numerical models, MODFLOW model is more modular in program structure, simpler in discrete method and more diversified in solving methods, which is more suitable for numerical simulation calculation of Jinan spring area. Through this application, no matter from the establishment, inspection or prediction of the model, it is relatively accurate and basically conforms to the hydrogeological conditions of Jinan spring area, and should be popularized in the field of numerical calculation of Jinan spring area.