A simulation and prediction method for groundwater pollution in karst areas with multiple media

Abstract

The present invention belongs to the technical field of water pollution prevention and control, and specifically discloses a method for simulating and predicting groundwater pollution in karst areas with multiple media. The method conducts pumping experiments, geophysical surveys and field tracer experiments to determine the connectivity of faults and karst pipelines, water flow velocity and average water flow velocity in pipelines, pipeline connectivity and branching structure; constructs a three-dimensional non-steady-state groundwater flow equation for the matrix to determine the area and duration of the pollution source and the characteristic pollution factors and leakage concentration; calibrates hydrogeological parameters based on the groundwater level; uses the weighted variance of the monitored water level, spring flow rate and spring water solute concentration as the calibration objective function; converts the fixed recharge coefficient into the recharge intensity, and predicts the water level fluctuation of the long observation hole, the spring flow rate and the concentration of pollution factors in the spring water based on the non-steady-state flow simulation. The present invention can accurately predict the evolution of karst springs and water quality, and can provide a decision-making basis for the protection and restoration of groundwater resources in karst areas.

Description

Technical Field

[0001] The invention relates to the technical field of water pollution prevention and control, and in particular to a method for simulating and predicting groundwater pollution in a karst area with multiple media. Background Art

[0002] Karst groundwater is an important freshwater resource, and the assessment and prevention of karst groundwater pollution is an engineering problem that needs to be solved urgently. Since there are many caves, underground rivers and other advantageous flow channels inside the aquifer, the pollutants that leak into the pipeline through sinkholes, faults, etc. migrate very quickly, making karst groundwater easy to pollute and difficult to control. Therefore, using effective methods to accurately characterize the faults, fissures and karst pipeline media in the karst area and predict the pollution of karst groundwater will help protect and repair groundwater in the karst area and provide protection for human drinking water health.

[0003] At present, there are still difficulties in accurately simulating the pollution assessment and prevention and control of karst groundwater. The traditional equivalent porous media method regards the pipeline as a highly permeable medium, but it is limited by Darcy's law and cannot reflect the characteristics of rapid solute migration in the pipeline. In the existing technology, MODFLOW, an internationally used groundwater simulation software developed by the United States Geological Survey, can approximate the drainage simulation of karst pipelines with its Drain or MNW2 package, but it is only applicable to cases with low Reynolds numbers and cannot simulate spring flow and pipeline solute migration. Although MODFLOW-2005 CFP v2 is recognized as a more mature karst groundwater simulation program and is widely used in groundwater flow simulation of karst matrix-pipeline dual media, the improved CFPv2 is also only used to simulate solute migration in karst pipelines (Reimann, 2015; Xu, 2015; Chang, 2019; Kavousi, 2020). To this end, the commercial software COMSOL Multiphysics uses Stokes-Darcy coupling boundaries to directly simulate the water volume and solute exchange between pipes and matrices, but it is limited by its huge computational complexity and cannot be applied to regional scales. In addition, in the assessment of karst groundwater pollution, existing technologies rarely consider regional faults as pollutant leakage channels, and lack rapid modeling and prediction methods for karst groundwater pollution that can be directly applied on site.

[0004] Therefore, there is an urgent need for a simulation method that can accurately identify faults and pipelines in karst areas and conduct groundwater pollution, so as to achieve the gradual calibration of parameters from simple models to complex models and ensure the simulation accuracy and computational efficiency of the model, thereby accurately predicting the karst spring flow and the migration patterns of pollutants through dominant channels such as faults and pipelines, and providing technical support for the protection and restoration of water resources in karst areas. Summary of the invention

[0005] In view of the deficiencies in the prior art, the present invention provides a method for simulating and predicting groundwater pollution in karst areas with multiple media.

[0006] The present invention is implemented in the following steps: obtaining hydrogeological parameters and monitoring data of the target study area, inferring the karst pipeline structure, constructing a three-dimensional hydrogeological model, model identification and calibration, and pollution simulation prediction. The specific contents are as follows:

[0007] A. Obtain hydrogeological parameters and monitoring data of the target study area: Select appropriate boreholes from the existing geological exploration boreholes in the target area for pumping experiments, and obtain the matrix permeability coefficient based on the pumping experiments K m , and the empirical formula for steady flow single-hole pumping of submerged non-complete wells is used to calculate K m :

[0008] ,

[0009] Where: Q For a fixed pumping volume [L 3 T -1 ]; s The water level of the pumping well is lowered [L]; s 1 is the water level drawdown in the observation well [L]; L is the filter length [L]; r 1 is the distance from the observation well to the center of the pumping well [L]; r well is the radius of the pumping well [L];

[0010] Obtaining field observation data: Counting groundwater levels in boreholes in the target area, selecting representative observation holes for long-term water level observations, and counting karst spring flow and concentrations of characteristic pollutants;

[0011] B. Inferring the karst pipeline structure: Based on the distribution of karst springs and sinkholes in the target area, a geophysical survey scheme is used to infer the regional faults, joints and water-bearing karst pipeline structure in the target area; field tracer experiments are used to determine the connectivity and water flow velocity of the faults and karst pipelines in the target area, and then the average water flow velocity of the pipeline is inferred based on the time when the tracer reaches the karst spring, and the connectivity and branching structure of the pipeline are determined based on the peak characteristics of the tracer penetration curve;

[0012] C. Constructing a three-dimensional hydrogeological model: First, the karst fissures are generalized into the matrix, the regional faults are used as the dominant seepage channels according to the water conductivity, and the karst underground rivers are simplified into segmented equal-diameter seepage pipes. The matrix and the pipes are coupled through water volume and solute exchange to construct the matrix three-dimensional unsteady groundwater flow equation;

[0013] Then, based on the regional faults, joints and water-bearing karst pipeline structures inferred in step B, a three-dimensional hydrogeological model is constructed in combination with regional geological data. At the same time, the area and duration of the pollution source are determined based on the historical evolution of the pollution source. Subsequently, the leaked percolation water is sampled for full water quality analysis and testing to determine the characteristic pollution factors and leakage concentration. C 0 ;

[0014] D. Model identification and calibration: First, based on the hydrogeological parameters obtained in step A as matrix parameters, try to input the hydrogeological parameters of the fault and pipeline into the three-dimensional hydrogeological model established in step C, simulate the steady-state flow field of the flow field according to the initial conditions and boundary conditions of the three-dimensional hydrogeological model, and compare the simulated water level with the measured water level. According to the comparison results, adjust the hydrogeological parameters to reduce the simulation error;

[0015] Then, based on the results of step B, the weighted variance of the monitored water level, spring flow, and spring water solute concentration is used as the calibration objective function:

[0016] ,

[0017] Where: Φ is the optimization objective function; k represents water level, spring flow and water quality; i is the observation data sequence number; N is the total number of observation data; w i For the i The weighting coefficient of each observation data; X obs It is the experimental data set; X sim is a simulated data set; r H , r Q , r C They are the simulation deviations of water level, spring flow and water quality, respectively;

[0018] E. Pollution simulation prediction: based on the dynamic rainfall of a hydrological cycle at the meteorological station in the target area q p,t , fixed recharge coefficient according to parameter partition α Converted into supply intensity q re,t ,Right now q re,t = αq p,t, adjust the simulation time of the three-dimensional hydrogeological model to a time series, adjust the three-dimensional hydrogeological model to a non-steady-state flow, and then run the non-steady-state flow model to obtain the water level fluctuation of the Changguan hole, the spring flow rate, and the concentration of pollution factors in the spring water.

[0019] Furthermore, the specific process of step B is as follows:

[0020] B1. Geophysical survey: According to the distribution of karst springs and sinkholes in the target area, joint profiling method and / or high-density electrical sounding method are used for detection. For the detection data, individual distortion points are deleted, and then inversion and mapping are performed. The boundary with inverted pseudo-resistivity less than the specific value contour line and large gradient is taken as the low-resistance anomaly area. According to the distribution of "zero value points" and low-resistance anomaly areas, the regional faults, joints and water-bearing karst pipeline structures in the target area are inferred;

[0021] B2. Field tracer experiment: The connectivity and water flow velocity of regional faults and karst conduits are determined by field tracer experiment, in which the background value is detected at the receiving point before the tracer is released, and then the fluorescence spectrophotometer is used for automatic and continuous observation at the receiving point;

[0022] Tracer recovery:

[0023] ,

[0024] Where: M r and M tot are the recovered mass and the put-in mass [M] respectively; i is the number of recycling times; C i No. i Secondary recovery concentration [ML -3 ]; Q i No. i Secondary flow rate [L 3 T -1 ]; Δ t Monitoring interval [T];

[0025] Then the average water flow rate in the pipeline is estimated based on the time it takes for the tracer to reach the karst spring. ,in L c is the possible pipeline path length, t 0 is the arrival time of the tracer; then the connectivity and branching structure of the pipeline are determined according to the peak characteristics of the tracer penetration curve, where a single peak is generally a single karst pipeline; double peaks or multiple peaks are branching structures.

[0026] Furthermore, the measurement electrode distance of the joint profile method in step B1 is AB / 2=55~155m; the high-density electrical sounding method is used to detect the depth range of 15~180m, using the Wenner array α 1 and Schlumber arrangement α 2Two types of electrode arrangements.

[0027] Furthermore, the tracer is any one or any combination of fluorescent whitening agent, sodium fluorescein, carmine and rhodamine B.

[0028] Furthermore, the specific process of step C is as follows:

[0029] C1. Constructing mathematical model: The karst fissures are generalized as matrix, the regional faults are used as the dominant seepage channels according to the water conductivity, the karst underground river is simplified into segmented equal-diameter seepage pipes, and the matrix and the pipes are coupled through water volume and solute exchange to construct the matrix three-dimensional non-steady-state groundwater flow equation:

[0030] ,

[0031] Where: ▽ is the Nabla operator; S s is the matrix water storage rate [L -1 ]; H is the matrix head [L]; t is time [T]; K is the permeability tensor [LT -1 ]; q ex is the exchange flux per unit volume between the matrix and the pipe [T -1 ]; q s is the unit volume source and sink term [T -1 ];

[0032] The exchange flux per unit volume between the matrix and the tube is q ex Assuming a linear relationship:

[0033] ,

[0034] Where: d is the Dirac function, which is used to indicate whether there is a pipeline node in the matrix grid; c is the water exchange coefficient per unit volume between the matrix and the pipe [LT -1 ], h is the water head in the pipe [L];

[0035] The flow rate in the pipeline is calculated according to the flow pattern, using the Hagen-Poiseuille formula or the Colebrook-White empirical formula to calculate the pipeline flow rate;

[0036] According to the mass conservation at each pipeline node:

[0037] ,

[0038] Where: Q c,i The pipeline node i The inflow of the connected pipes [L 3 T -1 ]; Q ex is the flow rate from the pipeline to the substrate [L 3 T -1 ]; Q r is the direct recharge at the pipeline node through sinkholes or rainfall [L 3 T -1 ], which is usually obtained through observation or estimation;

[0039] The solute migration in the matrix mainly considers convection, diffusion and delay effects, and ignores the chemical reaction in the matrix; it is expressed as:

[0040] ,

[0041] Where: f is the matrix porosity [-]; C m is the mass concentration of the solute in the matrix [ML -3 ]; t is time [T]; oh It is used to identify the water flow direction between the matrix and the pipe [-], which is determined by the head difference. When its value is 1, it means that the groundwater is discharged from the equivalent continuous medium to the pipe, and when its value is 0, it means that the groundwater is discharged from the pipe to the matrix; ▽ is the Nabla operator; D m is the matrix diffusion [L]; v m is the matrix water flow rate [LT -1 ]; q m,s represents the volume flow of the source and sink in the matrix [T -1 ]; C m,s Represents the solute concentration in the source and sink [ML -3 ]; q ex represents the exchange rate of unit volume medium-pipe water flow [T -1 ]; dis the Dirichlet function, which is used to indicate whether there is a pipeline node in the matrix grid; C c is the solute concentration in the pipeline [ML -3 ]; R N is the reaction rate [ML -3 T -1 ];

[0042] The solute migration in the pipe ignores the hydrodynamic diffusion effect and is described by the one-dimensional convection equation:

[0043] ,

[0044] Where: V c is the volume of water in the pipe [L 3 ]; l is the length of the pipeline [L]; v c is the flow velocity in the pipeline [LT -1 ]; q c,s represents the volume flow rate of the pipeline source and sink [T -1 ]; C c,s Represents the solute concentration of the source and sink [ML -3 ].

[0045] C2. Constructing conceptual model: According to the surface watershed and river boundary, the simulation scope and aquifer type are determined, and the vertical atmospheric precipitation infiltration recharge intensity is determined according to the meteorological data of previous years; at the same time, the karst fissures and water-conducting faults are generalized as equivalent porous media, and the karst pipeline is assumed to be a cylindrical pipeline distributed in three-dimensional space;

[0046] Then, based on the regional faults, joints and water-bearing karst pipeline structures inferred in step B, a three-dimensional hydrogeological model is constructed in combination with regional geological data; at the same time, the pollution source area and duration are determined based on the historical evolution of the pollution source; then the pollution source is generalized: the leaked percolation water is sampled for full water quality analysis and testing to determine the characteristic pollution factors and leakage concentration C 0 .

[0047] Furthermore, in the step C1, in order to determine the pipeline flow state, the value range of the Reynolds number Re is defined [Re min ,Re max ]:①When Re>Re max ② When the flow rate decreases, the flow state changes from laminar flow to turbulent flow; <Re min The flow state changes from turbulent flow to laminar flow; the above method is used to achieve a smooth transition of the pipeline flow state to ensure the stability of the numerical calculation;

[0048] The flow rate in the pipeline is calculated based on the flow pattern:

[0049] (i) Under laminar flow conditions, the Hagen-Poiseuille formula is used to calculate the pipeline flow rate:

[0050] ,

[0051] Where: h is the pipe head [L]; Δ l is the local pipeline length [L]; d is the pipe diameter [L]; r is the density of water [ML -3 ], g is the gravitational acceleration [LT -2 ]; m is the dynamic viscosity of water [ML -1 T -1 ]; t is the pipeline tortuosity [-]; Δ h is the head difference in the local pipe, Δ h / t Δ l is the hydraulic slope of the pipeline.

[0052] (ii) Under turbulent flow conditions, the Colebrook-White empirical formula is used to calculate the pipeline flow rate:

[0053] ,

[0054] Where: ξ is the roughness height of the pipe wall [L]; A is the cross-sectional area of ​​the pipe [L 2 ]; v is the average flow velocity in the pipe [LT -1 ].

[0055] Furthermore, after the preliminary calibration of the hydrogeological parameters in step D, the time for the pollutants to reach the karst spring is determined based on the average water flow rate in the pipeline calculated based on the tracer experiment in step B on the basis of the steady-state flow field. t 0 Then, the three-dimensional hydrogeological model established in step C is used for preliminary simulation, that is, to simulate the rapid migration time of sewage from the pollution source along the pipeline. t 0 ; The above t 0 Compare the simulated results of the karst spring water quality with the first data of the measured concentration of the karst spring in step A, and adjust the pipe diameter d , pipeline wall roughness height ξ, pipeline tortuosity t Matrix diffusivity Dm Matrix solute delay factor R m , so that the above t 0 The simulated concentration of pollutants at the spring point at time t matches the measured value, which is used as the initial state of the unsteady flow model.

[0056] Furthermore, the E step also includes prediction evaluation: using simulation error and parameter sensitivity to evaluate the reliability of the prediction results of the unsteady flow model simulation; wherein the parameter sensitivity is evaluated by mature software PEST or UCODE, and the simulation error is evaluated by using the Nash coefficient NSE and the Kling-Gupta efficiency coefficient KGE, and the expression is:

[0057] ,

[0058] ,

[0059] ,

[0060] Where: is the test average value; is the simulated average value; r is the Pearson linear correlation coefficient, α = s sim / s obs is the ratio of the simulated and measured standard deviations, β = is the ratio of the simulated and measured average values.

[0061] Furthermore, based on the reliability evaluation of the simulation prediction results of the model, the changes in the pollution of karst springs over time are further predicted based on the non-steady-state flow model. Then, according to the hydrogeological parameters obtained in step A and the karst pipeline structure in step B, the regional water diversion tunnel and curtain grouting are designed to prevent and control the rapid migration of pollution in the karst area.

[0062] The beneficial effects of the LSTM invention are:

[0063] 1. The present invention can accurately identify the regional karst hydrogeological structure, especially the dominant karst groundwater channel, through pumping tests, geophysical surveys and field tracer experimental techniques; and simulate and predict pollution migration based on the accurate karst hydrogeological structure, and can provide effective groundwater pollution risk assessment results, which has good guiding significance for groundwater pollution control in karst areas.

[0064] 2. The present invention generalizes the multiple media of faults, fissures and pipelines in karst areas, establishes a mathematical model, and couples the matrix and pipeline through water volume and solute exchange theory, considers the pollution diffusion of karst fissure media and the rapid migration effect of pipeline pollutants, and uses long-term series monitoring data to calibrate the model parameters. The model parameters are evaluated based on the Nash coefficient and the Kling-Gupta efficiency coefficient, which ensures the reliability of the model and improves the practicability of the mechanism model in actual engineering.

[0065] 3. The present invention uses the calibrated model to predict the pollution prevention and control effect of karst groundwater, which helps to reduce engineering costs and propose an optimized karst groundwater pollution control plan. It not only provides a decision-making basis for the protection and restoration of groundwater resources in karst areas, but also avoids the problem of low efficiency of groundwater control plans in karst areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 is a flow chart of the present invention;

[0067] Figure 2 A three-dimensional hydrogeological model in an embodiment of the present invention;

[0068] Figure 3 is the calibration result of the steady-state flow model in the embodiment of the present invention;

[0069] Figure 4 The figure is a comparison between the simulation result of the karst spring flow in the embodiment of the present invention and the measured value;

[0070] Figure 5 The figure is a comparison between the simulation results and the measured values ​​of the characteristic pollution factors of the karst spring in the embodiment of the present invention;

[0071] Figure 6 The figure is a comparison of the pollution control effects of curtain grouting in the embodiments of the present invention. DETAILED DESCRIPTION

[0072] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0073] like Figure 1 to 6 As shown, the present invention includes the steps of obtaining hydrogeological parameters and monitoring data of the target study area, inferring the karst pipeline structure, constructing a three-dimensional hydrogeological model, model identification and calibration, and pollution simulation prediction, and the specific contents are as follows:

[0074] A. Obtain hydrogeological parameters and monitoring data of the target study area: Select appropriate boreholes from the existing geological exploration boreholes in the target area for pumping experiments, and obtain the matrix permeability coefficient based on the pumping experimentsK m , and the empirical formula for steady flow single-hole pumping of submerged non-complete wells is used to calculate K m :

[0075] ,

[0076] Where: Q For a fixed pumping volume [L 3 T -1 ]; s The water level of the pumping well is lowered [L]; s 1 is the water level drawdown in the observation well [L]; L is the filter length [L]; r 1 is the distance from the observation well to the center of the pumping well [L]; r well is the radius of the pumping well [L];

[0077] Obtaining field observation data: Counting groundwater levels in boreholes in the target area, selecting representative observation holes for long-term water level observations, and counting karst spring flow and concentrations of characteristic pollutants;

[0078] B. Inferring the karst pipeline structure: Based on the distribution of karst springs and sinkholes in the target area, a geophysical survey scheme is used to infer the regional faults, joints and water-bearing karst pipeline structure in the target area; field tracer experiments are used to determine the connectivity and water flow velocity of the faults and karst pipelines in the target area, and then the average water flow velocity of the pipeline is inferred based on the time when the tracer reaches the karst spring, and the connectivity and branching structure of the pipeline are determined based on the peak characteristics of the tracer penetration curve;

[0079] C. Constructing a three-dimensional hydrogeological model: First, the karst fissures are generalized into the matrix, the regional faults are used as the dominant seepage channels according to the water conductivity, and the karst underground rivers are simplified into segmented equal-diameter seepage pipes. The matrix and the pipes are coupled through water volume and solute exchange to construct the matrix three-dimensional unsteady groundwater flow equation;

[0080] Then, based on the regional faults, joints and water-bearing karst pipeline structures inferred in step B, a three-dimensional hydrogeological model is constructed in combination with regional geological data. At the same time, the area and duration of the pollution source are determined based on the historical evolution of the pollution source. Subsequently, the leaked percolation water is sampled for full water quality analysis and testing to determine the characteristic pollution factors and leakage concentration. C 0 ;

[0081] D. Model identification and calibration: First, based on the hydrogeological parameters obtained in step A as matrix parameters, try to input the hydrogeological parameters of the fault and pipeline into the three-dimensional hydrogeological model established in step C, simulate the steady-state flow field of the flow field according to the initial conditions and boundary conditions of the three-dimensional hydrogeological model, and compare the simulated water level with the measured water level. According to the comparison results, adjust the hydrogeological parameters to reduce the simulation error;

[0082] Then, based on the results of step B, the weighted variance of the monitored water level, spring flow, and spring water solute concentration is used as the calibration objective function:

[0083] ,

[0084] Where: Φ is the optimization objective function; k represents water level, spring flow and water quality; i is the observation data sequence number; N is the total number of observation data; w i For the i The weighting coefficient of each observation data; X obs It is the experimental data set; X sim is a simulated data set; r H , r Q , r C They are the simulation deviations of water level, spring flow and water quality, respectively;

[0085] E. Pollution simulation prediction: based on the dynamic rainfall of a hydrological cycle at the meteorological station in the target area q p,t , fixed recharge coefficient according to parameter partition α Converted into supply intensity q re,t ,Right now q re,t = αq p,t , adjust the simulation time of the three-dimensional hydrogeological model to a time series, adjust the three-dimensional hydrogeological model to a non-steady-state flow, and then run the non-steady-state flow model to obtain the water level fluctuation of the Changguan hole, the spring flow rate, and the concentration of pollution factors in the spring water.

[0086] The specific process of step B is as follows:

[0087] B1. Geophysical survey: According to the distribution of karst springs and sinkholes in the target area, joint profile method and / or high-density electrical sounding method are used for detection. For the detection data, individual distortion points are deleted, and then inversion and mapping are performed. The boundary with the inverted pseudo-resistivity less than the contour line of a specific value (determined according to the actual drilling data combined with geophysical data) and a large gradient is taken as the low-resistance anomaly area. According to the distribution of "zero value points" and low-resistance anomaly areas, the regional faults, joints and water-bearing karst pipeline structures in the target area are inferred;

[0088] B2. Field tracer experiment: The connectivity and water flow velocity of regional faults and karst conduits are determined by field tracer experiment, in which the background value is detected at the receiving point before the tracer is released, and then the fluorescence spectrophotometer is used for automatic and continuous observation at the receiving point;

[0089] Tracer recovery:

[0090] ,

[0091] Where: M r and M tot are the recovered mass and the put-in mass [M] respectively; i is the number of recycling times; C i No. i Secondary recovery concentration [ML -3 ]; Q i No. i Secondary flow rate [L 3 T -1 ]; Δ t Monitoring interval [T];

[0092] Then the average water flow rate in the pipeline is estimated based on the time it takes for the tracer to reach the karst spring. ,in L c is the possible pipeline path length, t 0 is the arrival time of the tracer; then the connectivity and branching structure of the pipeline are determined according to the peak characteristics of the tracer penetration curve, where a single peak is generally a single karst pipeline; double peaks or multiple peaks are branching structures.

[0093] The measurement electrode distance of the combined profiling method in step B1 is AB / 2=55~155m; the high-density electrical sounding method is used to detect the depth range of 15~180m, using the Wenner array α 1 and Schlumber arrangement α 2Two types of electrode arrangements.

[0094] The tracer is any one or any combination of fluorescent whitening agent, sodium fluorescein, carmine and rhodamine B.

[0095] The specific process of step C is as follows:

[0096] C1. Constructing mathematical model: The karst fissures are generalized as matrix, the regional faults are used as the dominant seepage channels according to the water conductivity, the karst underground river is simplified into segmented equal-diameter seepage pipes, and the matrix and the pipes are coupled through water volume and solute exchange to construct the matrix three-dimensional non-steady-state groundwater flow equation:

[0097] ,

[0098] Where: ▽ is the Nabla operator; S s is the matrix water storage rate [L -1 ]; H is the matrix head [L]; t is time [T]; K is the permeability tensor [LT -1 ]; q ex is the exchange flux per unit volume between the matrix and the pipe [T -1 ]; q s is the unit volume source and sink term [T -1 ];

[0099] The exchange flux per unit volume between the matrix and the tube is q ex Assuming a linear relationship:

[0100] ,

[0101] Where: d is the Dirac function, which is used to indicate whether there is a pipeline node in the matrix grid; c is the water exchange coefficient per unit volume between the matrix and the pipe [LT -1 ], h is the water head in the pipe [L];

[0102] The flow rate in the pipeline is calculated according to the flow pattern, using the Hagen-Poiseuille formula or the Colebrook-White empirical formula to calculate the pipeline flow rate;

[0103] According to the mass conservation at each pipeline node:

[0104] ,

[0105] Where: Q c,iThe pipeline node i The inflow of the connected pipes [L 3 T -1 ]; Q ex is the flow rate from the pipeline to the substrate [L 3 T -1 ]; Q r is the direct recharge at the pipeline node through sinkholes or rainfall [L 3 T -1 ], which is usually obtained through observation or estimation;

[0106] in, Q ex It is automatically calculated from the head difference and permeability coefficient between the pipe and the matrix, which is different for each pipe location; Qr sinkhole recharge is generally estimated through observation and then manually input into the model;

[0107] The solute migration in the matrix mainly considers convection, diffusion and delay effects, and ignores the chemical reaction in the matrix; it is expressed as:

[0108] ,

[0109] Where: f is the matrix porosity [-]; C m is the mass concentration of the solute in the matrix [ML -3 ]; t is time [T]; oh It is used to identify the water flow direction between the matrix and the pipe [-], which is determined by the head difference. When its value is 1, it means that the groundwater is discharged from the equivalent continuous medium to the pipe, and when its value is 0, it means that the groundwater is discharged from the pipe to the matrix; ▽ is the Nabla operator; D m is the matrix diffusion [L]; v m is the matrix water flow rate [LT -1 ]; q m,s represents the volume flow of the source and sink in the matrix [T -1 ]; C m,s Represents the solute concentration in the source and sink [ML -3 ]; q ex represents the exchange rate of unit volume medium-pipe water flow [T -1 ]; d is the Dirichlet function, which is used to indicate whether there is a pipeline node in the matrix grid; C c is the solute concentration in the pipeline [ML -3 ];R N is the reaction rate [ML -3 T -1 ];

[0110] The solute migration in the pipe ignores the hydrodynamic diffusion effect and is described by the one-dimensional convection equation:

[0111] ,

[0112] Where: V c is the volume of water in the pipe [L 3 ]; l is the length of the pipeline [L]; v c is the flow velocity in the pipeline [LT -1 ]; q c,s represents the volume flow rate of the pipeline source and sink [T -1 ]; C c,s Represents the solute concentration of the source and sink [ML -3 ].

[0113] C2. Constructing conceptual model: According to the surface watershed and river boundary, the simulation scope and aquifer type are determined, and the vertical atmospheric precipitation infiltration recharge intensity is determined according to the meteorological data of previous years; at the same time, the karst fissures and water-conducting faults are generalized as equivalent porous media, and the karst pipeline is assumed to be a cylindrical pipeline distributed in three-dimensional space;

[0114] Then, based on the regional faults, joints and water-bearing karst pipeline structures inferred in step B, a three-dimensional hydrogeological model is constructed in combination with regional geological data; at the same time, the pollution source area and duration are determined based on the historical evolution of the pollution source; then the pollution source is generalized: the leaked percolation water is sampled for full water quality analysis and testing to determine the characteristic pollution factors and leakage concentration C 0 .

[0115] In the step C1, in order to determine the pipeline flow state, the value range of the Reynolds number Re is defined [Re min ,Re max ]:①When Re>Re max ② When the flow rate decreases, the flow state changes from laminar flow to turbulent flow; <Re min The flow state changes from turbulent flow to laminar flow; the above method is used to achieve a smooth transition of the pipeline flow state to ensure the stability of the numerical calculation;

[0116] The flow rate in the pipeline is calculated based on the flow pattern:

[0117] (i) Under laminar flow conditions, the Hagen-Poiseuille formula is used to calculate the pipeline flow rate:

[0118] ,

[0119] Where: h is the pipe head [L]; Δ l is the local pipeline length [L]; d is the pipe diameter [L]; r is the density of water [ML -3 ], g is the gravitational acceleration [LT -2 ]; m is the dynamic viscosity of water [ML -1 T -1 ]; t is the pipeline tortuosity [-]; Δ h is the head difference in the local pipe, Δ h / t Δ l is the hydraulic slope of the pipeline.

[0120] (ii) Under turbulent flow conditions, the Colebrook-White empirical formula is used to calculate the pipeline flow rate:

[0121] ,

[0122] Where: ξ is the roughness height of the pipe wall [L]; A is the cross-sectional area of ​​the pipe [L 2 ]; v is the average flow velocity in the pipe [LT -1 ].

[0123] After the preliminary calibration of the hydrogeological parameters in step D, the time for the pollutants to reach the karst spring is determined based on the average water flow rate in the pipeline calculated based on the tracer experiment in step B on the basis of the steady-state flow field. t 0 Then, the three-dimensional hydrogeological model established in step C is used for preliminary simulation, that is, to simulate the rapid migration time of sewage from the pollution source along the pipeline. t 0 ; The above t 0 Compare the simulated results of the karst spring water quality with the first data of the measured concentration of the karst spring in step A, and adjust the pipe diameter d , pipeline wall roughness height ξ, pipeline tortuosity t Matrix diffusivity D m Matrix solute delay factor R m , so that the above t 0The simulated concentration of pollutants at the spring point at time t matches the measured value, which is used as the initial state of the unsteady flow model.

[0124] The E step also includes prediction evaluation: using simulation error and parameter sensitivity to evaluate the reliability of the prediction results of the unsteady flow model simulation; wherein the parameter sensitivity is performed using mature software PEST or UCODE, and the simulation error is evaluated using the Nash coefficient NSE and the Kling-Gupta efficiency coefficient KGE, the expression is:

[0125] ,

[0126] ,

[0127] ,

[0128] Where: is the test average value; is the simulated average value; r is the Pearson linear correlation coefficient, α = s sim / s obs is the ratio of the simulated and measured standard deviations, β = is the ratio of the simulated and measured average values.

[0129] It should be noted that the general NSE value range is (-∞,1]. When NSE ≥ 0.5, it means that the model is acceptable; when NSE ≥ 0.65, it means that the model is good; when NSE ≥ 0.75, it means that the model is very good. However, the Nash coefficient is not easy to be used alone as an indicator to judge the quality of the model. The KGE value range is the same as that of NSE, but compared with NSE, KGE can provide a more comprehensive model performance evaluation.

[0130] Based on the reliability evaluation of the simulation prediction results of the model, the present invention further predicts the change of karst spring pollution over time based on the non-steady-state flow model, and then, according to the hydrogeological parameters obtained in step A and the karst pipeline structure in step B, designs the regional water diversion tunnel and curtain grouting to prevent and control the rapid migration of pollution in the karst area.

[0131] Example 1

[0132] S100: Select appropriate boreholes from existing geological exploration boreholes in the target area to conduct pumping experiments, and obtain the matrix permeability coefficient based on the pumping experiments K m , and the empirical formula for steady flow single-hole pumping of submerged non-complete wells is used to calculate K m :

[0133] ,

[0134] Where: Q For a fixed pumping volume [L 3 T -1 ]; s The water level of the pumping well is lowered [L]; s 1 is the water level drawdown in the observation well [L]; L is the filter length [L]; r 1 is the distance from the observation well to the center of the pumping well [L]; r well is the radius of the pumping well [L].

[0135] A one-year field observation was conducted to obtain groundwater observation data from 35 areas in the target area, and 210 days of long-term observation of borehole water level, spring flow and water quality monitoring data.

[0136] S200: The specific process is as follows:

[0137] S210: According to the distribution of karst springs and sinkholes in the target area, high-density electrical sounding method is used for detection (high-density electrical sounding method is used to detect the depth range of 15~180m, using Wenner array α 1 and Schlumber arrangement α 2 two types of electrode arrangement). For the detection data, delete individual distortion points, then perform inversion and mapping, and use the inverted pseudo-resistivity ≤1000Ω·m contour and the boundary with a large gradient as the low-resistance anomaly area. According to the distribution of "zero value points" and low-resistance anomaly areas, infer the regional faults, joints and water-bearing karst pipeline structures in the target area. Finally, the geophysical survey divided a total of 34 low-resistance anomaly areas.

[0138] S220: Field tracer experiments are used to determine the connectivity and water flow velocity of regional faults and karst pipelines. The background value is detected at the receiving point before sodium fluorescein is released, and then a fluorescence spectrophotometer is used to automatically and continuously observe at the receiving point.

[0139] Tracer recovery:

[0140] ,

[0141] Where: M r and M tot are the recovered mass and the put-in mass [M] respectively; i is the number of recycling times; C i No. i Secondary recovery concentration [ML-3 ]; Q i No. i Secondary flow rate [L 3 T -1 ]; Δ t Monitoring interval [T].

[0142] Field tracer experiments showed that the sodium fluorescein in the karst spring reached a peak at around 131 hours, and the recovery rate of sodium fluorescein was about 10%. Based on the arrival time of sodium fluorescein, the average water flow velocity was estimated to be 18.8 m / h.

[0143] like Figure 2 As shown in Figure 1, a multi-media hydrogeological model based on faults, fissures and conduits in the karst area was finally constructed. The parameters used in the model are shown in Table 1.

[0144] Table 1 Model hydrogeological parameters

[0145]

[0146] S300: The specific process is as follows:

[0147] S310: The karst fissures are generalized as the matrix, the regional faults and joints are used as the dominant seepage channels according to the water conductivity, the karst underground rivers and caves are simplified as segmented equal-diameter seepage pipes, and the matrix and the pipes are coupled through water volume and solute exchange to construct the matrix three-dimensional unsteady groundwater flow equation:

[0148] ,

[0149] Where: ▽ is the Nabla operator; S s is the matrix water storage rate [L -1 ]; H is the matrix head [L]; t is time [T]; K is the permeability tensor [LT -1 ]; q ex is the exchange flux per unit volume between the matrix and the pipe [T -1 ]; q s is the unit volume source and sink term [T -1 ].

[0150] The exchange flux per unit volume between the matrix and the tube is q ex Assuming a linear relationship:

[0151] ,

[0152] Where: d is the Dirac function, which is used to indicate whether there is a pipeline node in the matrix grid; c is the water exchange coefficient per unit volume between the matrix and the pipe [LT -1 ]; h is the water head in the pipe [L].

[0153] In order to determine the flow state in the pipeline, the value range of the Reynolds number Re is defined [Re min ,Re max ]:①When Re>Re max ② When the flow rate decreases, the flow state changes from laminar flow to turbulent flow; <Re min The flow state changes from turbulent flow to laminar flow; the above method is used to achieve a smooth transition of the pipeline flow state to ensure the stability of the numerical calculation.

[0154] The flow rate in the pipeline is calculated based on the flow pattern:

[0155] (i) Under laminar flow conditions, the Hagen-Poiseuille formula is used to calculate the pipeline flow rate:

[0156] ,

[0157] Where: h is the pipe head [L]; Δ l is the local pipeline length [L]; d is the pipe diameter [L]; r is the density of water [ML -3 ], g is the gravitational acceleration [LT -2 ]; m is the dynamic viscosity of water [ML -1 T -1 ]; t is the pipeline tortuosity [-]; Δ h is the head difference in the local pipe, Δ h / t Δ l is the hydraulic slope of the pipeline.

[0158] (ii) Under turbulent flow conditions, the Colebrook-White empirical formula is used to calculate the pipeline flow rate.

[0159] ,

[0160] Where: ξ is the roughness height of the pipe wall [L]; A is the cross-sectional area of ​​the pipe [L 2 ]; v is the average flow velocity in the pipe [LT -1 ].

[0161] According to the mass conservation at each pipeline node:

[0162] ,

[0163] Where: Q c,i The pipeline node i The inflow of the connected pipes [L 3 T -1 ]; Q ex is the flow rate from the pipeline to the substrate [L 3 T -1 ]; Q r is the direct recharge at the pipeline node through sinkholes or rainfall [L 3 T -1 ], usually obtained through observation or estimation.

[0164] in, Q ex It is automatically calculated through the head difference and permeability coefficient between the pipe and the matrix, which is different for each pipe location; Qr sinkhole recharge is generally estimated through observation and then manually input into the model.

[0165] The solute migration in the matrix mainly considers convection, diffusion and delay effects, and ignores the chemical reaction in the matrix; it is expressed as:

[0166] ,

[0167] Where: f is the matrix porosity [-]; C m is the mass concentration of the solute in the matrix [ML -3 ]; t is time [T]; oh It is used to identify the water flow direction between the matrix and the pipe [-], which is determined by the head difference. When its value is 1, it means that the groundwater is discharged from the equivalent continuous medium to the pipe, and when its value is 0, it means that the groundwater is discharged from the pipe to the matrix; ▽ is the Nabla operator; D m is the matrix diffusion [L]; v m is the matrix water flow rate [LT -1 ]; q m,s represents the volume flow of the source and sink in the matrix [T -1 ]; C m,s Represents the solute concentration in the source and sink [ML -3 ]; q exrepresents the exchange rate of unit volume medium-pipe water flow [T -1 ]; d is the Dirichlet function, which is used to indicate whether there is a pipeline node in the matrix grid; C c is the solute concentration in the pipeline [ML -3 ]; R N is the reaction rate [ML -3 T -1 ];

[0168] The solute migration in the pipe ignores the hydrodynamic diffusion effect and is described by the one-dimensional convection equation:

[0169] ,

[0170] Where: V c is the volume of water in the pipe [L 3 ]; l is the length of the pipeline [L]; v c is the flow velocity in the pipeline [LT -1 ]; q c,s represents the volume flow rate of the pipeline source and sink [T -1 ]; C c,s Represents the solute concentration of the source and sink [ML -3 ].

[0171] S320: According to the surface watershed and river boundaries, the simulation scope, area and aquifer type are determined, and the vertical atmospheric precipitation infiltration recharge intensity is determined based on historical meteorological data; at the same time, karst fissures and water-conducting faults are generalized as equivalent porous media, and karst pipelines are assumed to be cylindrical pipelines distributed in three-dimensional space.

[0172] Then, based on the regional faults, joints and water-bearing karst pipeline structures inferred by S200, a three-dimensional hydrogeological model was constructed in combination with regional geological data. At the same time, the pollution source area and duration were determined based on the historical evolution of the pollution source. Subsequently, the leaked percolation water was sampled for full water quality analysis and testing to determine the characteristic pollution factors and leakage concentration. C 0 .

[0173] In this example, the pollution sources are the slag field and the ammonia nitrogen waste liquid pool, and the main characteristic pollution factors are PO 4 3- NH 4 + 、F - According to the results of the full water quality analysis, the slag field leachate PO 4 3- and F- The leakage concentrations of the ammonia nitrogen waste liquid pool F were 2550mg / L and 48.7mg / L respectively; - and NH 4 + The leakage concentrations were 825mg / L and 116mg / L respectively.

[0174] S400: Figure 3 As shown, firstly, based on the hydrogeological parameters obtained in S100 as matrix parameters, try to input the hydrogeological parameters of the fault and the pipeline into the three-dimensional hydrogeological model established in S300, simulate the steady-state flow field of the flow field according to the initial conditions and boundary conditions of the three-dimensional hydrogeological model, and compare the simulated water level with the measured water level. According to the comparison results, the hydrogeological parameters are adjusted to reduce the simulation error.

[0175] On the basis of the steady-state flow field, the steady-state flow model is adjusted to a non-steady state according to the dynamic rainfall recharge and the adjustment of the simulation time: according to the average water flow rate in the pipeline calculated based on the tracer experiment in S200, the time when the pollutants reach the karst spring is determined t 0 Then, the three-dimensional hydrogeological model established by S300 was used for preliminary simulation to simulate the rapid migration time of sewage from the pollution source along the pipeline. t 0 The solute concentration in the spring water after t 0 Compare the simulated results of the karst spring water quality at this moment with the first data of the measured concentration of the karst spring in step S100, and adjust the pipe diameter d , pipeline wall roughness height ξ, pipeline tortuosity t Matrix diffusivity D m and matrix solute delay factor R m , so that the above t 0 The simulated concentration of pollutants at the spring point at time t is matched with the measured value, which is used as the initial state of the non-steady-state model.

[0176] Then, based on the results of S200, the weighted variance of the monitored water level, spring flow and spring water solute concentration was used as the calibration objective function:

[0177] ,

[0178] Where: Φ is the optimization objective function; k represents water level, spring flow and water quality; i is the observation data sequence number; N is the total number of observation data; w i For the iThe weighting coefficient of each observation data; X obs It is the experimental data set; X sim is a simulated data set; r H , r Q , r C They are the simulation deviations of water level, spring flow and water quality, respectively.

[0179] S500: The specific process is as follows:

[0180] S510: The concentration field 131 hours after the simulated pollutant leakage is used as the initial concentration field. According to the dynamic rainfall of the meteorological station in the target area in a hydrological cycle, the fixed recharge coefficient of the parameter partition is converted into the recharge intensity, and then the water level fluctuation of the long observation hole, the spring flow and the concentration of the pollution factors in the spring water (such as Figure 4 and Figure 5 as shown).

[0181] S520: The reliability of the model simulation prediction results is evaluated by simulation error and parameter sensitivity; the parameter sensitivity is evaluated by mature software PEST or UCODE, and the simulation error is evaluated by Nash coefficient NSE and Kling-Gupta efficiency coefficient KGE, the expression is:

[0182] ,

[0183] ,

[0184] ,

[0185] Where: is the test average value; is the simulated average value; r is the Pearson linear correlation coefficient; α = s sim / s obs is the ratio of the simulated and measured standard deviations; β = is the ratio of the simulated and measured average values.

[0186] In this embodiment, in the simulation of karst spring flow, NSE is 0.72 and KGE is 0.83; in the simulation of characteristic pollutants, NSE is greater than 0.2 and KGE is greater than 0.5, which proves that the model is relatively reliable.

[0187] S530: Based on the reliability evaluation of the simulation prediction results of the model, the pollution changes of karst springs over time are further predicted based on the unsteady flow model. Then, according to the hydrogeological parameters obtained by S100 and the karst pipeline structure inferred by S200, a curtain grouting scheme is designed at the fault and pipeline connection to block the migration channel of pollutants in the karst area (such as Figure 6 as shown).

[0188] The above is only a preferred specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by any technician familiar with the technical field within the technical scope disclosed by the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A method for simulating and predicting groundwater pollution in a karst area with multiple media, characterized in that: It includes obtaining hydrogeological parameters and monitoring data of the target study area, inferring the karst pipeline structure, building a three-dimensional hydrogeological model, model identification and calibration, and pollution simulation and prediction steps. The specific contents are as follows: A. Obtain hydrogeological parameters and monitoring data of the target study area: Select appropriate boreholes from the existing geological exploration boreholes in the target area for pumping experiments, and obtain the matrix permeability coefficient based on the pumping experiments K m , and the empirical formula for steady flow single-hole pumping of submerged non-complete wells is used to calculate K m : ， Where: Q For a fixed pumping volume [L 3 T -1 ]; s The water level of the pumping well is lowered [L]; s 1 is the water level drop in the observation well [L]; L is the filter length [L]; r 1 is the distance from the observation well to the center of the pumping well [L]; r well is the radius of the pumping well [L]; Obtaining field observation data: Counting groundwater levels in boreholes in the target area, selecting representative observation holes for long-term water level observations, and counting karst spring flow and concentrations of characteristic pollutants; B. Inferring the karst pipeline structure: Based on the distribution of karst springs and sinkholes in the target area, a geophysical survey scheme is used to infer the regional faults, joints and water-bearing karst pipeline structure in the target area; field tracer experiments are used to determine the connectivity and water flow velocity of the faults and karst pipelines in the target area, and then the average water flow velocity of the pipeline is inferred based on the time when the tracer reaches the karst spring, and the connectivity and branching structure of the pipeline are determined based on the peak characteristics of the tracer penetration curve; C. Constructing a three-dimensional hydrogeological model: First, the karst fissures are generalized into the matrix, the regional faults are used as the dominant seepage channels according to the water conductivity, and the karst underground rivers are simplified into segmented equal-diameter seepage pipes. The matrix and the pipes are coupled through water volume and solute exchange to construct the matrix three-dimensional unsteady groundwater flow equation; Then, based on the regional faults, joints and water-bearing karst pipeline structures inferred in step B, a three-dimensional hydrogeological model is constructed in combination with regional geological data. At the same time, the area and duration of the pollution source are determined based on the historical evolution of the pollution source. Subsequently, the leaked percolation water is sampled for full water quality analysis and testing to determine the characteristic pollution factors and leakage concentration. C 0; D. Model identification and calibration: First, based on the hydrogeological parameters obtained in step A as matrix parameters, try to input the hydrogeological parameters of the fault and pipeline into the three-dimensional hydrogeological model established in step C, simulate the steady-state flow field of the flow field according to the initial conditions and boundary conditions of the three-dimensional hydrogeological model, and compare the simulated water level with the measured water level. According to the comparison results, adjust the hydrogeological parameters to reduce the simulation error; Then, based on the results of step B, the weighted variance of the monitored water level, spring flow, and spring water solute concentration is used as the calibration objective function: ， Where: Φ is the optimization objective function; κ represents water level, spring flow and water quality; i is the observation data sequence number; N is the total number of observation data; w i For the i The weighting coefficient of each observation data; X obs It is the experimental data set; X sim is a simulated data set; r H , r Q , r C They are the simulation deviations of water level, spring flow and water quality, respectively; E. Pollution simulation prediction: based on the dynamic rainfall of a hydrological cycle at the meteorological station in the target area q p,t , fixed recharge coefficient according to parameter partition α Converted into supply intensity q re,t ,Right now q re,t = αq p,t , adjust the simulation time of the three-dimensional hydrogeological model to a time series, adjust the three-dimensional hydrogeological model to a non-steady-state flow, and then run the non-steady-state flow model to obtain the water level fluctuation of the Changguan hole, the spring flow rate, and the concentration of pollution factors in the spring water.

2. The method for simulating and predicting groundwater pollution in karst areas with multiple media according to claim 1, characterized in that: The specific process of step B is as follows: B1. Geophysical survey: According to the distribution of karst springs and sinkholes in the target area, joint profiling method and / or high-density electrical sounding method are used for detection. For the detection data, individual distortion points are deleted, and then inversion and mapping are performed. The boundary with inverted pseudo-resistivity less than the specific value contour line and large gradient is taken as the low-resistance anomaly area. According to the distribution of "zero value points" and low-resistance anomaly areas, the regional faults, joints and water-bearing karst pipeline structures in the target area are inferred; B2. Field tracer experiment: The connectivity and water flow velocity of regional faults and karst conduits are determined by field tracer experiment, in which the background value is detected at the receiving point before the tracer is released, and then the fluorescence spectrophotometer is used for automatic and continuous observation at the receiving point; Tracer recovery: ， Where: M r and M tot are the recovered mass and the put-in mass [M] respectively; i is the number of recycling times; C i No. i Secondary recovery concentration [ML -3 ]; Q i No. i Secondary flow rate [L 3 T -1 ]; Δ t Monitoring interval [T]; Then the average water flow rate in the pipeline is estimated based on the time when the tracer reaches the karst spring. ,in L c is the possible pipeline path length, t 0 is the tracer arrival time; then the connectivity and branching structure of the pipeline are determined according to the peak characteristics of the tracer penetration curve, where a single peak generally represents a single karst pipeline; double or multiple peaks represent a branching structure.

3. The method for simulating and predicting groundwater pollution in karst areas with multiple media according to claim 2, characterized in that: The measurement electrode distance of the combined profiling method in step B1 is AB / 2 = 55~155 m; the high-density electrical sounding method is used to detect the depth range of 15~180m, using the Wenner array α 1 and Schlumber arrangement α 2Two types of electrode arrangements.

4. The method for simulating and predicting groundwater pollution in karst areas with multiple media according to claim 2, characterized in that: The tracer is any one or any combination of fluorescent whitening agent, sodium fluorescein, carmine and rhodamine B.

5. The method for simulating and predicting groundwater pollution in karst areas with multiple media according to claim 1, characterized in that: After the preliminary calibration of the hydrogeological parameters in step D, the time for the pollutants to reach the karst spring is determined based on the average water flow rate in the pipeline calculated based on the tracer experiment in step B on the basis of the steady-state flow field. t 0, and then use the three-dimensional hydrogeological model established in step C to perform preliminary simulation, that is, simulate the rapid migration time of sewage from the pollution source along the pipeline t 0; t Compare the simulated results of the karst spring water quality at time 0 with the first data of the measured concentration of the karst spring in step A, and adjust the pipe diameter d , pipeline wall roughness height ξ, pipeline tortuosity τ Matrix diffusivity D m Matrix solute delay factor R m , so that the above t The simulated concentration of pollutants at the spring point at time 0 matches the measured value, which is used as the initial state of the unsteady flow model.

6. The method for simulating and predicting groundwater pollution in karst areas with multiple media according to any one of claims 1 to 5, characterized in that: The E step also includes prediction evaluation: using simulation error and parameter sensitivity to evaluate the reliability of the prediction results of the unsteady flow model simulation; wherein the parameter sensitivity is evaluated by mature software PEST or UCODE, and the simulation error is evaluated by using the Nash coefficient NSE and the Kling-Gupta efficiency coefficient KGE, the expression is: ， ， ， Where: is the test average value; is the simulated average value; r is the Pearson linear correlation coefficient, α = σ sim / σ obs is the ratio of the simulated and measured standard deviations, β = is the ratio of the simulated and measured average values.

7. The method for simulating and predicting groundwater pollution in karst areas with multiple media according to claim 6, characterized in that: Based on the reliability evaluation of the simulation prediction results of the model, the change of karst spring pollution over time is further predicted based on the non-steady-state flow model. Then, according to the hydrogeological parameters obtained in step A and the karst pipeline structure in step B, the regional water diversion tunnel and curtain grouting are designed to prevent and control the rapid migration of pollution in the karst area.

Images