Underground water vulnerability evaluation method and system based on hydrogeological model

Abstract

The invention belongs to the technical field of hydrogeological evaluation, and particularly discloses an underground water vulnerability evaluation method and system based on a hydrogeological model.The method comprises the steps that underground water vulnerability evaluation parameters, including aquifer parameters, underground water flow parameters and pollutant characteristic parameters, of a research area are obtained; constructing a three-dimensional anisotropic hydrogeological model based on a water flow motion equation and a pollutant migration equation; calculating underground water vulnerability characteristic indexes including a pollutant arrival time index, a pollutant concentration attenuation index and a system recovery capability index; constructing a dynamic coupling groundwater vulnerability index, calculating a time-varying contribution coefficient and a nonlinear coupling index of each index, and determining a dynamic threshold value to obtain a groundwater vulnerability quantitative evaluation result of the research area; according to the method, scientificity and accuracy of underground water vulnerability evaluation are improved, and decision support is provided for underground water resource protection and management.

Description

Technical Field

[0001] The present invention relates to the technical field of hydrogeological evaluation, and in particular to a groundwater vulnerability evaluation method and system based on a hydrogeological model. Background Art

[0002] Traditional groundwater vulnerability assessment methods are mainly based on static indicator systems, such as DRASTIC, GOD and other methods. These methods usually adopt a simplified weight accumulation method, ignoring the dynamic characteristics and nonlinear coupling effects of the groundwater system, and it is difficult to accurately reflect the actual response process of the groundwater system to pollutants; in addition, traditional methods usually assume that the aquifer is a homogeneous isotropic medium, and cannot consider the heterogeneity and anisotropy characteristics under actual geological conditions, resulting in deviations between the evaluation results and the actual situation.

[0003] Therefore, there is an urgent need to develop a vulnerability assessment method based on hydrogeological models that can reflect the dynamic response characteristics of the groundwater system, so as to improve the scientificity and accuracy of groundwater vulnerability assessment and provide decision-making support for groundwater resource protection and management. Summary of the invention

[0004] The purpose of the present invention is to provide a groundwater vulnerability assessment method and system based on a hydrogeological model, which is used to solve the problems of static evaluation, ignoring the heterogeneous and anisotropic characteristics and nonlinear coupling effects between indicators in traditional groundwater vulnerability assessment methods, and to achieve the purpose of dynamic quantitative evaluation and scientific early warning of groundwater vulnerability.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions: In a first aspect, the present invention provides a groundwater vulnerability assessment method based on a hydrogeological model, and obtains a quantitative assessment result of groundwater vulnerability in a study area based on a pollutant dynamics coupling method. The assessment method comprises the following steps: Step S1, obtaining groundwater vulnerability assessment parameters of the study area, wherein the groundwater vulnerability assessment parameters include: aquifer parameters, groundwater flow parameters and pollutant characteristic parameters.

[0006] Step S2, constructing a three-dimensional anisotropic hydrogeological model based on water flow equations and pollutant transport equations according to the dynamic response characteristics of the groundwater system to pollutants.

[0007] Step S3, calculating characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, wherein the characteristic indicators include a pollutant arrival time index, a pollutant concentration attenuation index, and a system recovery capacity index.

[0008] Step S4, constructing a dynamically coupled groundwater vulnerability index, calculating the time-varying contribution coefficient and nonlinear coupling index of each indicator and determining the dynamic threshold, and obtaining the quantitative evaluation results of groundwater vulnerability in the study area.

[0009] The aquifer parameters include specific water storage rate , hydraulic conductivity tensor K(x,y,z) and volumetric water content θ; the groundwater flow parameters include hydraulic head h, groundwater velocity vector v and source-sink term q; the pollutant characteristic parameters include pollutant concentration C, hydrodynamic diffusion coefficient tensor D, pollutant attenuation coefficient λ, pollution source intensity qs and pollution source concentration Cs.

[0010] The three-dimensional anisotropic hydrogeological model is constructed using the tetrahedron finite element discretization method. The expressions of the water flow equation and pollutant transport equation in the model are: The equation of motion for water flow is: ; Pollutant transport equation: ; in, is the specific water storage rate, which indicates the amount of water that can be released or absorbed by a unit volume of aquifer under a unit head change. is the water head, which is the sum of the potential height and the pressure head, representing the energy level of groundwater. It represents the rate of change of water head over time, indicating the rate of rise and fall of groundwater level; is the spatial coordinate The heterogeneous hydraulic conductivity tensor, the larger the K value, the easier it is for water to flow through the medium; represents the variation in the spatial distribution of groundwater flow, i.e. the difference in the amount of water flowing into and out of the system; is the source-sink term, which indicates the increase or decrease of water in a unit volume per unit time; is the water content per unit volume, is the pollutant concentration, is the hydrodynamic diffusion coefficient tensor, which indicates the degree and rate of pollutant diffusion. represents the net change in pollutant mass flux caused by diffusion; is the groundwater velocity vector, is the pollutant attenuation coefficient, is the pollution source intensity, which indicates the amount of water injected into the system per unit volume per unit time. is the pollution source concentration, which indicates the concentration of pollutants injected into the water; is a vector differential operator. In the three-dimensional rectangular coordinate system (x, y, z), is defined as: .

[0011] The pollutant arrival time index (RTI) is calculated based on the time required for the pollutant concentration to reach 50% of the initial concentration. The system recovery capacity index is calculated based on the ratio of the maximum pollutant concentration at the target point to the initial pollutant concentration. Calculated based on system contamination response time and recovery time.

[0012] The pollutant arrival time index RTI calculation formula is: ,in is the time required for the pollutant concentration to reach 50% of the initial concentration, obtained by solving the pollutant transport equation. It is a scale parameter and is adjusted dynamically with the aquifer thickness and hydraulic gradient.

[0013] The pollutant concentration decay index The calculation formula is: ,in is the maximum concentration of pollutants at the target point, is the initial pollutant concentration, is the distance from the pollution source to the target point, is the medium adsorption attenuation coefficient; The system resilience index The calculation formula is: ,in and are the time required for the pollutant concentration to reach 90% and 10% of the initial concentration, respectively, and are obtained by solving the pollutant transport equation. is the initial value of pollution, is the final equilibrium concentration, The time required for system recovery.

[0014] The dynamic coupled groundwater vulnerability index DCG is calculated using a nonlinear coupling equation: ;in, is the time-varying contribution coefficient, It is a nonlinear coupling index, which is used to characterize the synergistic effect among various indicators.

[0015] The time-varying contribution coefficient Determined by dynamic threshold method: ,in, For the The sensitivity function of an indicator is calculated as follows: ; In the formula, is the corresponding indicator value, is the indicator variance that changes with time, representing the spatiotemporal uncertainty of the indicator.

[0016] The nonlinear coupling index Determined by the theory of self-organized criticality: ,in, For the The coefficient of variation of the indicator, is the mean of the coefficient of variation of all indicators, For indicators The fractal dimension, which characterizes the spatial self-similarity of the index, is calculated through variogram analysis.

[0017] The dynamic threshold is determined by the theory of self-organized criticality: the frequency distribution function of the DCG value in the study area is calculated ; Search The inflection point is the critical point of phase transition ;based on Determine the dynamic classification threshold of groundwater vulnerability level.

[0018] The dynamic classification thresholds of groundwater vulnerability levels include: Very low vulnerability: ; Low vulnerability: Moderate vulnerability: ; High vulnerability: ; Extremely high vulnerability: .

[0019] According to the dynamic threshold division results, the groundwater vulnerability zoning map of the study area was generated using the geographic information system.

[0020] In a second aspect, the present invention provides a groundwater vulnerability assessment system based on a hydrogeological model, which is used to implement the method of the first aspect of the right, and the system includes: a parameter acquisition module, a model construction module, an indicator calculation module and a vulnerability assessment module connected in sequence.

[0021] The parameter acquisition module is used to obtain groundwater vulnerability assessment parameters in the study area, including aquifer parameters, groundwater flow parameters and pollutant characteristic parameters.

[0022] The model building module is used to build a three-dimensional anisotropic hydrogeological model based on water flow motion equations and pollutant migration equations according to the dynamic response characteristics of the groundwater system to pollutants.

[0023] The index calculation module is used to calculate characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, including pollutant arrival time index, pollutant concentration attenuation index and system recovery capacity index.

[0024] The vulnerability assessment module is used to construct a dynamically coupled groundwater vulnerability index, calculate the time-varying contribution coefficient and nonlinear coupling index of each indicator and determine the dynamic threshold, so as to obtain a quantitative assessment result of groundwater vulnerability in the study area.

[0025] The system also includes: a result visualization module, which is used to divide the results according to dynamic thresholds, generate a groundwater vulnerability zoning map of the study area using a geographic information system, and provide a spatiotemporal evolution analysis function.

[0026] Compared with the prior art, the present invention has the following beneficial effects: The present invention overcomes the limitations of traditional evaluation methods by establishing a theoretical system for groundwater vulnerability evaluation based on a hydrogeological model, wherein the heterogeneity and anisotropy of aquifers are taken into account, a comprehensive evaluation system including pollutant arrival time index, concentration decay index and system recovery capacity index is constructed, and a nonlinear coupling equation is used to calculate the kinetically coupled groundwater vulnerability index. At the same time, the self-organized critical state theory is introduced to determine the dynamic threshold, thus realizing the automation of the evaluation process and the visualization of the results. The evaluation results can be dynamically updated according to environmental changes, thus providing a scientific basis and decision-making support tool for groundwater resource protection and pollution prevention and control. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 It is a flow chart of the groundwater vulnerability assessment method based on the hydrogeological model of the present invention; Figure 2 It is a schematic diagram of the composition of the groundwater vulnerability assessment system based on the hydrogeological model of the present invention. DETAILED DESCRIPTION

[0028] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention is described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0029] It should be noted that the present invention solves the problems existing in the existing groundwater vulnerability assessment methods, including: the traditional assessment model is too simplified, ignoring the dynamic characteristics of the groundwater system and the pollution process; failing to consider the heterogeneity and anisotropy of the geological medium; using a simple linear superposition method for evaluation, failing to reflect the nonlinear coupling relationship between indicators; the evaluation threshold setting is highly subjective and lacks scientific basis; the evaluation results lack a dynamic updating mechanism in time and space, making it difficult to adapt to environmental changes. The present invention constructs a three-dimensional anisotropic hydrogeological model, quantitatively simulates the groundwater flow and pollutant migration process, establishes a characteristic indicator system based on the dynamic response of pollutants, uses a nonlinear coupling equation to calculate the groundwater vulnerability index, and determines the dynamic threshold through the self-organized critical state theory, thereby achieving scientific evaluation and effective early warning of groundwater vulnerability.

[0030] Example 1 like Figure 1 As shown, it is a flow chart of the groundwater vulnerability assessment method based on the hydrogeological model of the present invention, and the quantitative assessment result of groundwater vulnerability in the study area is obtained based on the pollutant dynamics coupling method. The assessment method includes the following steps: Step S1, obtaining groundwater vulnerability assessment parameters of the study area, wherein the groundwater vulnerability assessment parameters include: aquifer parameters, groundwater flow parameters and pollutant characteristic parameters.

[0031] The aquifer parameters include specific water storage rate , hydraulic conductivity tensor K(x,y,z) and volumetric water content θ; the groundwater flow parameters include hydraulic head h, groundwater velocity vector v and source-sink term q; the pollutant characteristic parameters include pollutant concentration C, hydrodynamic diffusion coefficient tensor D, pollutant attenuation coefficient λ, pollution source intensity qs and pollution source concentration Cs.

[0032] Step S2, constructing a three-dimensional anisotropic hydrogeological model based on water flow equations and pollutant transport equations according to the dynamic response characteristics of the groundwater system to pollutants.

[0033] The dynamic response characteristics of the groundwater system to pollutants are mainly manifested as nonlinear laws in the process of temporal and spatial evolution. This dynamic response includes: (1) There is a hysteresis effect in the migration of pollutants in groundwater, that is, there is a certain time delay from the release of the pollution source to the observation point; (2) The concentration of pollutants is non-uniformly distributed in space, which is affected by the heterogeneity of the aquifer and the dominant flow path; (3) During the migration of pollutants, pollutants are affected by various physical and chemical processes such as adsorption, desorption, and degradation, showing concentration attenuation characteristics; (4) The groundwater system has a certain self-purification ability. After the pollution source is removed, the system will gradually return to its initial state, but the recovery process is time-dependent. By considering the directional differences in the hydraulic conductivity tensor, the three-dimensional anisotropic hydrogeological model can accurately reflect the different characteristics of groundwater flow and pollutant migration in different directions, thereby more realistically simulating the diffusion law of pollutants under complex geological conditions.

[0034] The three-dimensional anisotropic hydrogeological model is constructed using the tetrahedral finite element discretization method, which is a numerical method suitable for modeling complex three-dimensional geological structures. Compared with the traditional finite difference method, it has better spatial adaptability and can accurately express the geometric characteristics of irregular boundaries and internal geological structures. In this method, the study area is divided into a large number of tetrahedral units, and linear interpolation functions are used in each unit to describe the spatial distribution of unknown variables (head and concentration). After discretization, a large sparse matrix equation system is formed, which is solved by efficient solvers such as the preconditioned conjugate gradient method. An implicit backward difference format is used for the time term to ensure numerical stability. When meshing, local encryption is performed around geological interfaces, faults and pollution sources to improve calculation accuracy.

[0035] The expressions of the water flow equation and pollutant transport equation in the model are: The equation of motion for water flow is: ; Pollutant transport equation: ; in, The relative water storage rate indicates the amount of water that can be released or absorbed by a unit volume of aquifer under unit head change. The unit is usually 1 / m. It reflects the elastic water storage characteristics of the aquifer. The larger the value, the stronger the water storage capacity of the aquifer. The hydraulic head is the sum of the potential height and the pressure head, representing the energy level of groundwater. It is measured in meters (m). The hydraulic head determines the flow direction of groundwater (from high head to low head). It represents the rate of change of water head over time, indicating the rate of rise and fall of groundwater level; is the spatial coordinate The heterogeneous hydraulic conductivity tensor describes the three-dimensional spatial distribution characteristics of the water conductivity of the geological medium. As a tensor, it means that the water conductivity may be different in different directions (non-isotropic), with units of m / day or m / s. A higher K value means that water flows more easily through the medium. represents the variation in the spatial distribution of groundwater flow, i.e. the difference in the amount of water flowing into and out of the system; It is a source-sink term, which indicates the increase or decrease of water in a unit volume per unit time; such as rainfall infiltration (source) or pumping (sink), with units of 1 / s or 1 / day. It is the water content per unit volume, which represents the volume ratio of water in the unit volume of the medium. It is dimensionless and ranges from 0 to 1. In the saturated zone, it is approximately equal to the porosity.

[0036] It is the pollutant concentration, which indicates the mass of pollutants per unit volume of water, usually expressed in mg / L or g / m³. It indicates the rate of change of pollutant mass in unit volume of medium per unit time. It is the hydrodynamic diffusion coefficient tensor, which indicates the degree and rate of pollutant diffusion and describes the combined effect of pollutant diffusion and mechanical dispersion in groundwater. It is expressed as a tensor, indicating that the diffusion characteristics are different in different directions. represents the net change in pollutant mass flux caused by diffusion; It is the groundwater velocity vector, which describes the speed and direction of groundwater flow. Its unit is m / day or m / s. It is the actual groundwater flow velocity (not Darcy velocity) and is related to porosity.

[0037] Indicates the rate of concentration reduction of pollutants due to degradation, decay, etc., in units of 1 / day or 1 / s. The larger the value, the faster the pollutant degradation; It indicates the total amount of pollutant degradation or decay in unit volume per unit time; is the pollutant attenuation coefficient, The pollution source intensity indicates the amount of water injected into the system per unit volume per unit time, in units of 1 / day or 1 / s. It is the pollution source concentration, which indicates the concentration of pollutants injected into the water, and its unit is the same as C. is a vector differential operator. In the three-dimensional rectangular coordinate system (x, y, z), is defined as: Together, these two equations form a coupled system that describes groundwater flow and contaminant transport, taking into account the interaction of hydraulic processes (flow) and material transport processes (diffusion, convection, degradation). The anisotropic representation enables the model to more accurately simulate the behavior of groundwater systems under complex geological conditions.

[0038] Step S3, calculating characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, wherein the characteristic indicators include a pollutant arrival time index, a pollutant concentration attenuation index, and a system recovery capacity index.

[0039] The pollutant arrival time index (RTI) is calculated based on the time required for the pollutant concentration to reach 50% of the initial concentration. The system recovery capacity index is calculated based on the ratio of the maximum pollutant concentration at the target point to the initial pollutant concentration. Calculated based on system contamination response time and recovery time.

[0040] The pollutant arrival time index RTI calculation formula is: ,in is the time required for the pollutant concentration to reach 50% of the initial concentration, obtained by solving the pollutant transport equation. It is a scale parameter and is adjusted dynamically with the aquifer thickness and hydraulic gradient. is a scale parameter, and its physical meaning is the characteristic time scale, which is closely related to the physical properties and hydraulic conditions of the aquifer. The dynamic adjustment equation of the parameters is: is the reference value under reference conditions (usually 100 days), M is the thickness of the aquifer, is the reference thickness (taken as 10 meters), J is the hydraulic gradient, is the reference hydraulic gradient (taken as 0.01). The exponents 0.5 and 0.8 respectively reflect the influence of aquifer thickness and hydraulic gradient on pollutant migration time. These index values ​​are obtained through a large number of numerical simulations and experimental data fitting. In thicker aquifers, the vertical migration distance increases, resulting in longer pollutant arrival time; while the increase in hydraulic gradient will accelerate the groundwater flow rate, thereby reducing the arrival time of pollutants.

[0041] For example: a homogeneous aquifer with a thickness of 10 meters, a pollution source at coordinates (0,0,0), and an initial pollutant concentration = 100 mg / L; we need to calculate the time required for the pollutant to reach the coordinate (100,0,0) and reach 50 mg / L; for this aquifer: volumetric water content θ = 0.3, longitudinal diffusion coefficient DL = 10 m² / day; transverse diffusion coefficient DT = 1m² / day; vertical diffusion coefficient DV = 0.1 m² / day; groundwater velocity v = 0.5 m / day (x direction); pollutant attenuation coefficient λ = ; In this simplified one-dimensional case, an analytical solution can be used.

[0042] For an instantaneous point source pollution, the analytical solution in one dimension is: ; It is necessary to find the time t that satisfies C(100,t) = 50 mg / L; Here, a numerical method is used to solve the complete three-dimensional pollutant transport equation.

[0043] Build a 3D finite element model of the aquifer and randomly sample parameters (e.g. DL may vary in the range of 8-12 m² / day); for each set of parameters, calculate the change in pollutant concentration over time; determine the concentration of pollutants at the target point / 2 time; spatial step Δx = 1m; time step Δt = 0.1 day; in each time step, we calculate the pollutant concentration at each grid point; the following results are obtained through this numerical simulation: At t = 180 days, the concentration at the target point (100,0,0) is 48 mg / L; At t = 185 days, the concentration at the target point was 52 mg / L; By linear interpolation, we can determine that t50 ≈ 182.5 days. Once we have t50 = 182.5 days, we can calculate the RTI: Assuming β = 100 (a scaling parameter related to aquifer thickness and hydraulic gradient): ; This RTI value reflects the resistance of the point to pollutants, with lower values ​​indicating a longer arrival time for pollutants and therefore lower vulnerability.

[0044] The pollutant concentration decay index The calculation formula is: ,in is the maximum concentration of pollutants at the target point, is the initial pollutant concentration, is the distance from the pollution source to the target point, is the medium adsorption attenuation coefficient; the medium adsorption attenuation coefficient α is a comprehensive parameter that characterizes the adsorption and natural attenuation capacity of geological media for pollutants, and its calculation formula is: in is the distribution coefficient (L / kg), which indicates the distribution ratio of the pollutant between the solid and liquid phases; is the volume density of the medium (kg / L); is the effective porosity; is the pollutant degradation rate constant (1 / day); is the actual groundwater flow rate (m / day); The larger the value, the stronger the medium's ability to adsorb and degrade pollutants, and the more significant the attenuation of pollutants with distance. The values ​​vary greatly, for example, The values ​​are mainly affected by adsorption, while organic pollutants are affected by both adsorption and degradation.

[0045] The system resilience index The calculation formula is: ,in and are the time required for the pollutant concentration to reach 90% and 10% of the initial concentration, respectively, and are obtained by solving the pollutant transport equation. is the initial value of pollution, is the final equilibrium concentration, The time required for system recovery.

[0046] System recovery time The determination is based on the time course of the pollutant concentration returning to the background value or environmental standard; in actual calculations, an exponential decay model is used to fit the concentration change curve over time after the pollution source is removed: ; where k is the decay rate constant, determined by nonlinear regression. System recovery time Defined as the concentration from the peak Decrease to the final equilibrium concentration The time required for 95% of ; in the molecule The term is used to standardize the comparison under different initial pollution conditions, making the index dimensionless.

[0047] Step S4, constructing a dynamically coupled groundwater vulnerability index, calculating the time-varying contribution coefficient and nonlinear coupling index of each indicator and determining the dynamic threshold, and obtaining the quantitative evaluation results of groundwater vulnerability in the study area.

[0048] The dynamic coupled groundwater vulnerability index DCG is calculated using a nonlinear coupling equation: ;in, is the time-varying contribution coefficient, It is a nonlinear coupling index, which is used to characterize the synergistic effect among various indicators.

[0049] The time-varying contribution coefficient Determined by dynamic threshold method: ,in, For the The sensitivity function of an indicator is calculated as follows: ; In the formula, is the corresponding indicator value, is the indicator variance that changes with time, representing the spatiotemporal uncertainty of the indicator.

[0050] The nonlinear coupling index Determined by the theory of self-organized criticality: ,in, For the The coefficient of variation of the indicator, is the mean of the coefficient of variation of all indicators, For indicators The fractal dimension, which characterizes the spatial self-similarity of the index, is calculated through variogram analysis.

[0051] Nonlinear coupling index The determination of the self-organized critical state theory is based on the theory of self-organized critical state, which holds that complex systems will spontaneously tend to critical states during the evolution process. It is obtained by calculating the ratio of the standard deviation to the mean of the indicator in the study area, reflecting the spatial variation of the indicator.

[0052] Fractal dimension Characterize the complexity and self-similarity of the spatial distribution of indicators, and calculate through variogram analysis: First, calculate the semi-variogram at different distance intervals h , and then fit in the logarithmic coordinate system and The relationship between the slope of the curve and the fractal dimension D is: slope = 3-D (in two-dimensional space) or slope = 4-D (in three-dimensional space); The larger the value, the more complex the spatial distribution of the indicator and the more significant the corresponding nonlinear coupling effect.

[0053] The dynamic threshold is determined by the theory of self-organized criticality: the frequency distribution function of the DCG value in the study area is calculated ; Search The inflection point is the critical point of phase transition (calculate The first-order derivative and the second-order derivative of are used to determine the inflection point of the frequency distribution function by solving the point where the first-order derivative is zero and the second-order derivative is negative. The DCG value corresponding to the inflection point is the critical point of phase transition); based on Determine the dynamic classification threshold of groundwater vulnerability level.

[0054] The dynamic classification thresholds of groundwater vulnerability levels include: Very low vulnerability: ; Low vulnerability: Moderate vulnerability: ; High vulnerability: ; Extremely high vulnerability: .

[0055] According to the dynamic threshold division results, the groundwater vulnerability zoning map of the study area was generated using the geographic information system.

[0056] To further illustrate the technical effect of the present invention, the method is applied to the surrounding area of ​​a certain industrial park as an example. The area of ​​this area is about 25 square kilometers. Groundwater is the main source of drinking water. There are multiple potential pollution sources and the geological conditions are complex. First, the groundwater vulnerability assessment parameters of the study area are obtained through field investigation: the aquifer parameters include specific water storage rate scope The main direction values ​​of the hydraulic conductivity tensor K(x,y,z) are 5.2, 3.8, and 0.7 m / d, respectively, and the volume water content θ ranges from 0.18 to 0.32; the groundwater flow parameters include the head distribution of 138.5-152.3 m and the average flow velocity of 0.25 m / d; the pollutant characteristic parameters select nitrate as the indicator pollutant, and the hydrodynamic diffusion coefficient tensor longitudinal , horizontal , vertical , the attenuation coefficient λ is The pollution source intensities of industrial areas and landfills were 0.08 m^3 / d and 0.05 m^3 / d, respectively, and the concentrations were 120 mg / L and 85 mg / L, respectively.

[0057] A three-dimensional anisotropic hydrogeological model containing about 128,000 units was constructed based on the tetrahedron finite element method, and the model was calibrated using the head and concentration observation data of 16 monitoring wells. The groundwater vulnerability characteristic index was calculated for monitoring well R8: the pollutant arrival time index RTI = 0.487, where = 245 days, β = 163.8; pollutant concentration decay index CDI = 0.732, where =120 mg / L, =58.2 mg / L, d=430 m, α=0.0035 ; System Recovery Index SRI = 0.308, where =182 days, =412 days, =58.2 mg / L, =3.5 mg / L, trec=485 days.

[0058] By calculating the nonlinear coupling index ( =2.44, =2.88, =2.72) and the time-varying contribution coefficient ( =0.28, =0.42, =0.30), and the dynamic coupled groundwater vulnerability index DCG of monitoring well R8 was obtained to be 0.482. Based on the distribution of DCG values ​​in the study area, the phase change critical point DCGcr=0.525 was identified, and the vulnerability levels were divided accordingly: extremely low (<0.105), low (0.105-0.263), medium (0.263-0.525), high (0.525-0.788) and extremely high (≥0.788). The research results show that monitoring well R8 belongs to the medium vulnerability level; the overall vulnerability distribution of the study area is 12.6% in the extremely low vulnerability area, 24.8% in the low vulnerability area, 31.5% in the medium vulnerability area, 21.2% in the high vulnerability area, and 9.9% in the extremely high vulnerability area; the areas around industrial areas and landfills are mainly high and extremely high in vulnerability, and the vulnerability of groundwater flowing to the downstream area shows a gradual change.

[0059] Compared with the traditional DRASTIC method, the prediction accuracy of 15 pollution verification points was improved by 24.3% (86.7% vs 62.4%). During seasonal changes, this method can reflect the dynamic changes of vulnerability in the rainy season and drought period, while the DRASTIC method only gives static evaluation results. In heterogeneous mutation areas, this method can accurately identify the vulnerability mutation areas caused by heterogeneous mutations, while the DRASTIC method tends to smooth these changes and underestimate the vulnerability risk of mutation areas. The actual application results show that the method of the present invention has significant advantages in evaluation accuracy, dynamic response capability, heterogeneity and anisotropy adaptability, decision support capability and early warning capability, providing a scientific basis for the protection and management of groundwater resources.

[0060] Example 2 like Figure 2 As shown, it is a schematic diagram of the composition of the groundwater vulnerability assessment system based on the hydrogeological model of the present invention, which is used to execute the method of Example 1. The system includes: a parameter acquisition module, a model construction module, an index calculation module and a vulnerability assessment module connected in sequence.

[0061] The parameter acquisition module is used to obtain the groundwater vulnerability assessment parameters of the study area, including aquifer parameters, groundwater flow parameters and pollutant characteristic parameters. The parameter acquisition module integrates multi-source data acquisition and processing functions, supports the import and analysis of raw data such as drilling data, pumping tests, indoor permeability tests, and constructs parameter spatial distribution fields through geostatistical methods. Taking the application of an industrial park as an example, the module uses multi-point geostatistical methods to process the lithology and hydrogeological parameters of 42 boreholes, generates a three-dimensional distribution field of hydraulic conductivity, and ensures that the spatial interpolation accuracy is within the 95% confidence interval through cross-validation, realizing the parameter scale conversion from point scale to regional scale, and providing a high-quality parameter basis for model construction. In addition, the module also has a parameter sensitivity analysis function to identify the key parameters that have the most significant impact on the model results, thereby optimizing the on-site monitoring and data collection strategies and improving the pertinence and efficiency of data acquisition.

[0062] The model building module is used to build a three-dimensional anisotropic hydrogeological model based on the water flow motion equation and the pollutant transport equation according to the dynamic response characteristics of the groundwater system to pollutants. The model building module adopts advanced computational fluid dynamics theory and finite element numerical method, which can accurately characterize the groundwater flow and pollutant transport process under complex geological conditions. The module supports a variety of grid generation strategies, including structured grids, unstructured grids and hybrid grids, and can adaptively adjust the grid density according to the complexity of geological structures. For practical engineering applications, the module has developed a dedicated boundary condition processing tool to support the setting and verification of common first-class, second-class, third-class boundary conditions and well boundary conditions. The model solution engine adopts matrix sparse storage technology and parallel computing framework, which greatly improves the computing efficiency, so that the solution time of large-scale models (more than one million grid nodes) is controlled within an acceptable range. The model calibration function integrates a variety of automatic parameter estimation algorithms, such as PEST and SCE-UA, which can optimize and adjust the model parameters based on monitoring data to ensure the accuracy and reliability of the model.

[0063] The indicator calculation module is used to calculate the characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, including the pollutant arrival time index, the pollutant concentration attenuation index and the system recovery capacity index. The indicator calculation module has developed a complete set of post-processing analysis tool chains, which can extract key features from the spatiotemporal simulation results of the hydrogeological model and calculate vulnerability indicators. The module uses a particle tracking algorithm to calculate the migration path and arrival time of pollutants, and supports complex scenario simulations that consider processes such as adsorption and degradation. In order to improve the calculation efficiency, the module implements the batch processing function of indicator calculation, which can simultaneously process the indicator values ​​of multiple evaluation points or the entire regional grid. For different types of pollutants, the module has a rich built-in parameter library, including parameters such as the attenuation coefficient and distribution coefficient of common organic pollutants and heavy metals, which simplifies user input. In addition, the module also has a Monte Carlo simulation function. Through multiple random sampling simulations, it evaluates the impact of parameter uncertainty on the indicator calculation results, generates the probability distribution and confidence interval of the indicator, and provides a probabilistic basis for risk assessment.

[0064] The vulnerability assessment module is used to construct a dynamically coupled groundwater vulnerability index, calculate the time-varying contribution coefficient and nonlinear coupling index of each indicator and determine the dynamic threshold, so as to obtain a quantitative assessment result of groundwater vulnerability in the study area.

[0065] The vulnerability assessment module innovatively introduces complex system theory and nonlinear coupling mechanism, breaking through the limitations of traditional evaluation methods. The module uses a global sensitivity analysis method to calculate the time-varying contribution coefficient, which can identify key control factors in different periods and regions and realize dynamic adjustment of the evaluation. In view of the nonlinear relationship between indicators, the module has developed a coupling index calculation method based on the self-organized critical state theory, and quantified the spatial variation characteristics and self-similarity of indicators through coefficient of variation and fractal dimension analysis. The dynamic threshold determination function adopts a data-driven method to automatically identify the critical point of system state transition through the inflection point detection algorithm of the frequency distribution function, and generates a scientific and reasonable classification standard based on this. Compared with the traditional fixed threshold, this method better adapts to the geological and hydrological conditions of different regions, and the evaluation results are more objective and accurate. In addition, the module also supports multi-time scale evaluation, which can analyze the impact of interannual changes, seasonal fluctuations and extreme events on groundwater vulnerability, and realizes a methodological leap from static evaluation to dynamic evaluation.

[0066] The system also includes: a result visualization module, which is used to generate a groundwater vulnerability zoning map of the study area based on the dynamic threshold division results using a geographic information system, and provide a spatiotemporal evolution analysis function; a decision support module, which is used to provide groundwater protection zone division, pollution risk warning and management countermeasures suggestions based on the groundwater vulnerability assessment results, combined with land use planning and environmental management needs.

[0067] The vulnerability assessment module innovatively introduces complex system theory and nonlinear coupling mechanism, breaking through the limitations of traditional evaluation methods. The module uses a global sensitivity analysis method to calculate the time-varying contribution coefficient, which can identify key control factors in different periods and regions and realize dynamic adjustment of the evaluation. In view of the nonlinear relationship between indicators, the module has developed a coupling index calculation method based on the self-organized critical state theory, and quantified the spatial variation characteristics and self-similarity of indicators through coefficient of variation and fractal dimension analysis. The dynamic threshold determination function adopts a data-driven method to automatically identify the critical point of system state transition through the inflection point detection algorithm of the frequency distribution function, and generates a scientific and reasonable classification standard based on this. Compared with the traditional fixed threshold, this method better adapts to the geological and hydrological conditions of different regions, and the evaluation results are more objective and accurate. In addition, the module also supports multi-time scale evaluation, which can analyze the impact of interannual changes, seasonal fluctuations and extreme events on groundwater vulnerability, and realizes a methodological leap from static evaluation to dynamic evaluation.

[0068] The decision support module transforms scientific evaluation results into practical management measures, which reflects the application value of the system. Based on multi-criteria decision-making theory, this module integrates groundwater vulnerability evaluation results with socio-economic factors to generate an optimized groundwater protection zone division plan. The protection zone classification function takes into account multiple factors such as vulnerability level, pollution source distribution, land use status and planning and development needs. Through spatial overlay analysis, it automatically generates boundaries and control measures for different protection levels. The pollution risk early warning subsystem implements an automatic early warning mechanism based on threshold triggering. When the monitoring data or model prediction results exceed the preset threshold, the system automatically generates early warning information and pushes it to the relevant responsible person. The management countermeasures suggestion function is based on the case library and knowledge map, and provides targeted prevention and control measures for different vulnerable areas and pollution types, such as optimizing land use structure, adjusting industrial layout, and improving agricultural fertilization methods. In addition, this module also supports scenario simulation functions, which can evaluate the effects and cost-effectiveness of different management measures, assist decision makers in selecting the optimal management strategy, and achieve sustainable utilization and effective protection of groundwater resources.

[0069] In conjunction with the application case of Example 1, each functional module of the groundwater vulnerability assessment system based on the hydrogeological model is introduced in detail.

[0070] This system adopts a modular design architecture, and realizes scientific evaluation and dynamic early warning of groundwater vulnerability through the collaborative work of components such as parameter acquisition module, model construction module, indicator calculation module, vulnerability assessment module and result visualization module.

[0071] The parameter acquisition module is responsible for collecting and processing various parameter data of the study area; taking the industrial park in Example 1 as an example, the module integrates a multi-source data acquisition interface, which can simultaneously access drilling data, hydrogeological parameters, groundwater monitoring data, and pollutant characteristic parameters. For parameters such as the specific water storage rate, hydraulic conductivity, and volumetric water content of the aquifer in the industrial park, the system completes outlier detection and processing through data preprocessing functions, interpolates and estimates missing values, and generates a three-dimensional distribution field of parameters through spatial statistical methods. For example, in response to the heterogeneous characteristics of the hydraulic conductivity in the area, the system uses the Kriging interpolation method to generate a high-precision three-dimensional hydraulic conductivity field with a resolution of 5m×5m×1m, providing an accurate parameter basis for subsequent model construction.

[0072] The model building module is the core component of the system, responsible for building a three-dimensional anisotropic hydrogeological model. In view of the complex geological conditions in Example 1, the module uses adaptive meshing technology to achieve local encryption of the grid around key geological interfaces and pollution sources, and the total number of tetrahedral units reaches 128,000, which effectively ensures the calculation accuracy under complex geological structures. The model solving engine uses parallel computing technology, and the efficiency of solving the water flow motion equation and the pollutant transport equation is improved by about 3.5 times. The model calibration function realizes the automatic optimization of hydraulic parameters by integrating the PEST parameter estimation program, so that the average absolute error between the simulated head and the measured head is reduced to 0.22m, the relative error of the pollutant concentration simulation is controlled within 8.5%, and the model prediction ability is significantly enhanced.

[0073] The indicator calculation module calculates the characteristic indicators of groundwater vulnerability based on the calibrated hydrogeological model. For the monitoring well network in Example 1, the module calculates the pollutant arrival time index RTI, pollutant concentration decay index CDI and system resilience index SRI of each monitoring point through an automated workflow. Taking monitoring well R8 as an example, the system obtains the pollutant concentration time series curve through numerical simulation, automatically identifies key time points such as t50, t90 and t10, and calculates the indicator values ​​of RTI=0.487, CDI=0.732 and SRI=0.308. For the entire study area, the system calculates the spatial distribution of each indicator with a grid density of 25m×25m, generates a high-precision indicator distribution map, and intuitively shows the spatial heterogeneity of groundwater vulnerability.

[0074] The vulnerability assessment module is the decision-making core of the system, responsible for constructing the dynamically coupled groundwater vulnerability index and determining the dynamic threshold. For the industrial park in Example 1, the module first calculated the coefficient of variation and fractal dimension of each indicator, and obtained the nonlinear coupling index γ1=2.44, γ2=2.88, γ3=2.72; then based on the sensitivity analysis, the time-varying contribution coefficients w1=0.28, w2=0.42, w3=0.30 were determined; finally, the dynamically coupled groundwater vulnerability index DCG was calculated. The system innovatively introduced the self-organized critical state theory to determine the dynamic threshold, and through the inflection point analysis of the frequency distribution function, the phase change critical point DCGcr=0.525 was identified, and the five-level vulnerability classification standard was automatically generated accordingly. Compared with the traditional fixed threshold method, the dynamic threshold method has increased the accuracy of vulnerability identification in the study area by 18.7%, especially in the area where the geological conditions change suddenly, the identification ability is significantly enhanced.

[0075] The result visualization module provides rich visualization and analysis functions. For the evaluation results in Example 1, the module generates a groundwater vulnerability zoning map for the study area, uses five colors to identify different vulnerability levels, and supports multi-scale spatial browsing. The spatiotemporal evolution analysis function can simulate the changes in vulnerability under different seasonal conditions. For example, in the rainy season, the system identifies that a medium-vulnerability area on the east side of the industrial zone will temporarily turn into a high-vulnerability area, providing a scientific basis for seasonal management measures. In addition, the system also supports scenario simulation functions, such as simulating the impact of new pollution sources or changes in land use on groundwater vulnerability. The evaluation results show that if a new chemical plant is built in the northwest of the study area, the vulnerability level of the downstream 3.2 square kilometers will increase by 1-2 levels, providing a quantitative reference for regional planning and environmental impact assessment.

[0076] The decision support module combines the evaluation results with management needs to provide scientific decision-making recommendations. For the high-vulnerability and extremely high-vulnerability areas identified in Example 1 (accounting for 31.1% of the total area), the system automatically generates a graded protection plan, suggesting that the extremely high-vulnerability area (9.9%) be designated as a first-level protection area, prohibiting any activities that may pollute groundwater; and that the high-vulnerability area (21.2%) be designated as a second-level protection area, restricting high-risk industrial activities and the use of pesticides and fertilizers. The system also proposes a differentiated monitoring plan based on the spatiotemporal evolution characteristics of vulnerability, with 12 new monitoring wells deployed in high-vulnerability areas with a monitoring frequency of once a week; and 8 monitoring wells deployed in medium-vulnerability areas with a monitoring frequency of once a month, which significantly improves the pertinence and economy of the monitoring network.

[0077] Through the collaborative work of the above modules, this system was successfully applied to the groundwater vulnerability assessment of the industrial park in Example 1. Compared with the traditional evaluation method, the evaluation accuracy was improved by 24.3%, and it can reflect the dynamic changes in vulnerability caused by seasonal changes and human activities. The automation level and user-friendliness of the system increase the efficiency of the evaluation work by about 5 times, especially in terms of evaluation capabilities and dynamic response capabilities under complex geological conditions. This system has been successfully applied to groundwater resource management in multiple industrial parks, cities and river basins, providing a scientific basis and technical support for the demarcation of groundwater protection areas, pollution risk warnings and the formulation of management countermeasures.

[0078] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A groundwater vulnerability assessment method based on a hydrogeological model is used to obtain a quantitative assessment result of groundwater vulnerability in the study area based on a pollutant dynamics coupling method, which is characterized by: The evaluation method comprises the following steps: Step S1, obtaining groundwater vulnerability assessment parameters of the study area, wherein the groundwater vulnerability assessment parameters include: aquifer parameters, groundwater flow parameters and pollutant characteristic parameters; Step S2, constructing a three-dimensional anisotropic hydrogeological model based on water flow equations and pollutant transport equations according to the dynamic response characteristics of the groundwater system to pollutants; Step S3, calculating characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, wherein the characteristic indicators include a pollutant arrival time index, a pollutant concentration attenuation index, and a system recovery capacity index; Step S4, constructing a dynamically coupled groundwater vulnerability index, calculating the time-varying contribution coefficient and nonlinear coupling index of each indicator and determining the dynamic threshold, and obtaining the quantitative evaluation results of groundwater vulnerability in the study area.

2. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 1 is characterized in that: The aquifer parameters include specific water storage rate , hydraulic conductivity tensor K(x,y,z) and volumetric water content θ; the groundwater flow parameters include hydraulic head h, groundwater velocity vector v and source-sink term q; the pollutant characteristic parameters include pollutant concentration C, hydrodynamic diffusion coefficient tensor D, pollutant attenuation coefficient , pollution source intensity qs and pollution source concentration Cs.

3. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 2 is characterized in that: The three-dimensional anisotropic hydrogeological model is constructed using the tetrahedron finite element discretization method. The expressions of the water flow equation and pollutant transport equation in the model are: The equation of motion for water flow is: ; Pollutant transport equation: ; in, is the specific water storage rate, which indicates the amount of water that can be released or absorbed by a unit volume of aquifer under a unit head change. is the water head, which is the sum of the potential height and the pressure head, representing the energy level of groundwater. It represents the rate of change of water head over time, indicating the rate of rise and fall of groundwater level; is the spatial coordinate The heterogeneous hydraulic conductivity tensor, the larger the K value, the easier it is for water to flow through the medium; represents the variation in the spatial distribution of groundwater flow, i.e. the difference in the amount of water flowing into and out of the system; is the source-sink term, which indicates the increase or decrease of water in a unit volume per unit time; is the water content per unit volume, is the pollutant concentration, is the hydrodynamic diffusion coefficient tensor, which indicates the degree and rate of pollutant diffusion. represents the net change in pollutant mass flux caused by diffusion; is the groundwater velocity vector, is the pollutant attenuation coefficient, is the pollution source intensity, which indicates the amount of water injected into the system per unit volume per unit time. is the pollution source concentration, which indicates the concentration of pollutants injected into the water; is a vector differential operator. In the three-dimensional rectangular coordinate system (x, y, z), is defined as: .

4. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 3 is characterized in that: The pollutant arrival time index (RTI) is calculated based on the time required for the pollutant concentration to reach 50% of the initial concentration. The system recovery capacity index is calculated based on the ratio of the maximum pollutant concentration at the target point to the initial pollutant concentration. Calculated based on system contamination response time and recovery time; The pollutant arrival time index RTI calculation formula is: ,in is the time required for the pollutant concentration to reach 50% of the initial concentration, obtained by solving the pollutant transport equation. is a scale parameter that is dynamically adjusted with aquifer thickness and hydraulic gradient; The pollutant concentration decay index The calculation formula is: ,in is the maximum concentration of pollutants at the target point, is the initial pollutant concentration, is the distance from the pollution source to the target point, is the medium adsorption attenuation coefficient; The system resilience index The calculation formula is: ,in and are the time required for the pollutant concentration to reach 90% and 10% of the initial concentration, respectively, and are obtained by solving the pollutant transport equation. is the initial value of pollution, is the final equilibrium concentration, The time required for system recovery.

5. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 4 is characterized in that: The dynamic coupled groundwater vulnerability index DCG is calculated using a nonlinear coupling equation: ;in, is the time-varying contribution coefficient, It is a nonlinear coupling index, which is used to characterize the synergistic effect among various indicators.

6. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 5 is characterized in that: The time-varying contribution coefficient Determined by dynamic threshold method: ,in, For the The sensitivity function of an indicator is calculated as follows: ; In the formula, is the corresponding indicator value, is the indicator variance that changes with time, representing the spatiotemporal uncertainty of the indicator.

7. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 5 is characterized in that: The nonlinear coupling index Determined by the theory of self-organized criticality: ,in, For the The coefficient of variation of the indicator, is the mean of the coefficient of variation of all indicators, For indicators The fractal dimension, which characterizes the spatial self-similarity of the index, is calculated through variogram analysis.

8. The groundwater vulnerability assessment method based on the hydrogeological model according to claim 1 is characterized in that: The dynamic threshold is determined by the theory of self-organized criticality: the frequency distribution function of the DCG value in the study area is calculated ; Search The inflection point is the critical point of phase transition ;based on Determine the threshold for dynamic classification of groundwater vulnerability levels; The dynamic classification thresholds of groundwater vulnerability levels include: Very low vulnerability: ; Low vulnerability: Moderate vulnerability: ; High vulnerability: ; Extremely high vulnerability: ; According to the dynamic threshold division results, the groundwater vulnerability zoning map of the study area was generated using the geographic information system.

9. A groundwater vulnerability assessment system based on a hydrogeological model, used to implement the method according to any one of claims 1 to 8, characterized in that: The system comprises: a parameter acquisition module, a model construction module, an index calculation module and a vulnerability assessment module connected in sequence; The parameter acquisition module is used to obtain groundwater vulnerability assessment parameters of the study area, including aquifer parameters, groundwater flow parameters and pollutant characteristic parameters; The model building module is used to build a three-dimensional anisotropic hydrogeological model based on water flow motion equations and pollutant transport equations according to the dynamic response characteristics of the groundwater system to pollutants; The index calculation module is used to calculate characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, including pollutant arrival time index, pollutant concentration attenuation index and system recovery capacity index; The vulnerability assessment module is used to construct a dynamically coupled groundwater vulnerability index, calculate the time-varying contribution coefficient and nonlinear coupling index of each indicator and determine the dynamic threshold, so as to obtain a quantitative assessment result of groundwater vulnerability in the study area.

10. The groundwater vulnerability assessment system based on the hydrogeological model according to claim 9 is characterized in that: The system also includes: a result visualization module, which is used to divide the results according to dynamic thresholds, generate a groundwater vulnerability zoning map of the study area using a geographic information system, and provide a spatiotemporal evolution analysis function.

Images