Groundwater Vulnerability Assessment Method and System Based on Hydrogeological Model

Abstract

The present invention belongs to the technical field of hydrogeological evaluation, and specifically discloses a method and system for evaluating groundwater vulnerability based on a hydrogeological model. The method includes obtaining groundwater vulnerability evaluation parameters for the study area, including aquifer parameters, groundwater flow parameters, and pollutant characteristic parameters; constructing a three-dimensional anisotropic hydrogeological model based on the groundwater flow equation and the pollutant transport equation; calculating characteristic indicators of groundwater vulnerability, including pollutant arrival time index, pollutant concentration attenuation index, and system recovery ability index; constructing a kinetic coupling groundwater vulnerability index, calculating the time-varying contribution coefficients and non-linear coupling indices of each indicator and determining the dynamic threshold, so as to obtain a quantitative evaluation result of the groundwater vulnerability in the study area; this method improves the scientificity and accuracy of groundwater vulnerability evaluation and provides decision-making support for the protection and management of groundwater resources.

Description

Technical Field

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

[0002] Traditional methods for evaluating groundwater vulnerability mainly rely on static index systems, such as the DRASTIC, GOD, etc. These methods usually adopt a simplified weighted accumulation method, ignoring the dynamic characteristics and non-linear 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 and isotropic medium, unable to consider the heterogeneity and anisotropy characteristics under actual geological conditions, resulting in a deviation between the evaluation results and the actual situation.

[0003] Therefore, there is an urgent need to develop a vulnerability evaluation method based on a hydrogeological model that can reflect the dynamic response characteristics of the groundwater system, so as to improve the scientificity and accuracy of groundwater vulnerability evaluation 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 method and system for evaluating groundwater vulnerability based on a hydrogeological model, which is used to solve the problems of static evaluation, ignoring the heterogeneity and anisotropy characteristics, and non-linear coupling effects between indicators in traditional groundwater vulnerability evaluation methods, and achieve the purpose of dynamic quantitative evaluation and scientific early warning of groundwater vulnerability.

[0005] To achieve the above purpose, the present invention adopts the following technical solutions:

[0006] In the first aspect, the present invention provides a method for evaluating groundwater vulnerability based on a hydrogeological model, and obtains a quantitative evaluation result of groundwater vulnerability in the study area based on the pollutant dynamics coupling method. The evaluation method includes the following steps:

[0007] Step S1, obtain the groundwater vulnerability evaluation parameters in the study area, and the groundwater vulnerability evaluation parameters include: aquifer parameters, groundwater flow parameters, and pollutant characteristic parameters.

[0008] Step S2, based on the dynamic response characteristics of the groundwater system to pollutants, construct a three-dimensional anisotropic hydrogeological model composed of a groundwater flow equation and a pollutant transport equation.

[0009] Step S3, calculate the characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, and the characteristic indicators include pollutant arrival time index, pollutant concentration attenuation index, and system recovery ability index.

[0010] Step S4: Construct a dynamic coupling groundwater vulnerability index, calculate the time-varying contribution coefficients and non-linear coupling indices of each index, and determine the dynamic threshold to obtain the quantitative evaluation result of groundwater vulnerability in the study area.

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

[0012] The three-dimensional anisotropic hydrogeological model is constructed using the tetrahedral finite element discretization method. The expressions of the groundwater flow equation and the pollutant transport equation in the model are as follows:

[0013] Groundwater flow equation: ;

[0014] Pollutant transport equation: ;

[0015] Among them, is the specific water storage rate, which represents the amount of water that can be released or absorbed per unit volume of the aquifer under a unit change in hydraulic head. is the hydraulic head, which represents the sum of the potential energy height and the pressure head and represents the energy level of groundwater. represents the rate of change of the hydraulic head with time and represents the rising and falling speed of the groundwater level. is the spatial coordinate of the heterogeneous hydraulic conductivity tensor. The larger the K value, the easier it is for water to flow through the medium. represents the spatial distribution change of the groundwater flow, that is, the difference in the amount of water flowing into and out of the system. is the source-sink term, which represents the increase or decrease in the amount of water per unit volume per unit time. is the volumetric water content per unit volume. is the pollutant concentration. is the hydrodynamic dispersion coefficient tensor, which represents the degree and rate of pollutant diffusion. represents the net change in the pollutant mass flux caused by dispersion. is the groundwater flow velocity vector. is the pollutant decay coefficient. is the source strength, which represents the amount of water injected into the system per unit volume per unit time. is the source concentration, which represents the pollutant concentration in the injected water. is the vector differential operator. In the three-dimensional rectangular coordinate system (x, y, z), is defined as: .

[0016] The pollutant arrival time index RTI is calculated based on the time required for the pollutant concentration to reach 50% of the initial concentration, and the pollutant concentration attenuation index is calculated based on the ratio of the maximum pollutant concentration at the target point to the initial pollutant concentration. The system recovery ability index is calculated based on the system pollution response time and recovery time.

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

[0018] The pollutant concentration attenuation index The calculation formula is: , where is the maximum pollutant concentration 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;

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

[0020] The dynamic coupling groundwater vulnerability index DCG is calculated using a non - linear coupling equation:

[0021] ; where, is the time - varying contribution coefficient, is the non - linear coupling index, used to characterize the synergy effect between various indicators.

[0022] The time - varying contribution coefficient is determined by the dynamic threshold method:

[0023] , where, is the sensitivity function of the

[0024] ; In the formula, is the corresponding index value, is the index variance that changes with time, representing the spatio-temporal uncertainty of the index.

[0025] The non-linear coupling index is determined by the self-organized critical state theory:

[0026] , where is the coefficient of variation of the th index, is the mean value of the coefficient of variation of all indices, is the fractal dimension of the index , representing the spatial self-similarity of the index, and is calculated by variogram analysis.

[0027] The dynamic threshold is determined by the self-organized critical state theory: Calculate the frequency distribution function of the DCG value in the study area ; Find the inflection point of as the phase transition critical point ; Based on determine the dynamic division threshold of the groundwater vulnerability level.

[0028] The dynamic division threshold of the groundwater vulnerability level includes:

[0029] Very low vulnerability: ; Low vulnerability: Moderate vulnerability: ; High vulnerability: ; Extremely high vulnerability: .

[0030] According to the dynamic threshold division result, use the geographic information system to generate the groundwater vulnerability zoning map of the study area.

[0031] In the second aspect, the present invention provides a groundwater vulnerability evaluation system based on a hydrogeological model for performing the method of the first aspect of the rights. The system includes, connected in sequence: a parameter acquisition module, a model construction module, an index calculation module, and a vulnerability evaluation module.

[0032] The parameter acquisition module is used to acquire the groundwater vulnerability evaluation parameters of the study area, including aquifer parameters, groundwater flow parameters, and pollutant characteristic parameters.

[0033] The model construction module is used to construct a three-dimensional anisotropic hydrogeological model based on the groundwater flow equation and the pollutant transport equation according to the dynamic response characteristics of the groundwater system to pollutants.

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

[0035] The vulnerability evaluation module is used to construct a dynamic coupling groundwater vulnerability index, calculate the time-varying contribution coefficient and non-linear coupling index of each index and determine the dynamic threshold, so as to obtain the quantitative evaluation result of groundwater vulnerability in the study area.

[0036] The system further includes: a result visualization module, which is used to divide the results according to the dynamic threshold, generate a groundwater vulnerability zoning map of the study area by using a geographic information system, and provide a spatio-temporal evolution analysis function.

[0037] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0038] By establishing a theoretical system for groundwater vulnerability evaluation based on a hydrogeological model, the present invention overcomes the limitations of traditional evaluation methods, takes into account the heterogeneity and anisotropy characteristics of the aquifer, constructs a comprehensive evaluation system including pollutant arrival time index, concentration attenuation index and system recovery ability index, calculates the dynamic coupling groundwater vulnerability index by using a non-linear coupling equation, and introduces the self-organized critical state theory to determine the dynamic threshold, realizing the automation of the evaluation process and the visualization of the results, and being able to dynamically update the evaluation results according to environmental changes, providing a scientific basis and decision-making support tool for groundwater resource protection and pollution prevention and control. Description of the Drawings

[0039] Figure 1 It is a flowchart of the groundwater vulnerability evaluation method based on a hydrogeological model of the present invention;

[0040] Figure 2 It is a schematic diagram of the composition of the groundwater vulnerability evaluation system based on a hydrogeological model of the present invention. Detailed Embodiments

[0041] In order to make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions in the present invention will be clearly and completely described below. Apparently, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art without making creative efforts based on the embodiments of the present invention fall within the protection scope of the present invention.

[0042] It should be noted that the present invention solves the problems existing in the existing groundwater vulnerability assessment methods, including: the traditional assessment models are too simplistic, ignoring the dynamic characteristics of the groundwater system and the pollution process; the heterogeneity and anisotropy of the geological medium are not considered; a simple linear superposition method is used for assessment, failing to reflect the non-linear coupling relationship between indicators; the setting of the assessment threshold is highly subjective and lacks a scientific basis; the assessment results lack a spatio-temporal dynamic update mechanism and are difficult to adapt to environmental changes. The present invention constructs a three-dimensional anisotropic hydrogeological model, quantitatively simulates the groundwater flow and pollutant transport processes, establishes a characteristic index system based on the dynamic response of pollutants, calculates the groundwater vulnerability index using a non-linear coupling equation, and determines the dynamic threshold through the self-organized critical state theory, realizing the scientific assessment and effective early warning of groundwater vulnerability.

[0043] Example 1

[0044] As Figure 1 shown, it is a flow chart of the groundwater vulnerability assessment method based on a hydrogeological model of the present invention. Based on the pollutant dynamics coupling method, a quantitative assessment result of the groundwater vulnerability of the study area is obtained. The assessment method includes the following steps:

[0045] Step S1, obtain the groundwater vulnerability assessment parameters of the study area. The groundwater vulnerability assessment parameters include: aquifer parameters, groundwater flow parameters, and pollutant characteristic parameters.

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

[0047] Step S2, according to the dynamic response characteristics of the groundwater system to pollutants, construct a three-dimensional anisotropic hydrogeological model composed of a groundwater flow equation and a pollutant transport equation.

[0048] The dynamic response characteristics of the groundwater system to pollutants are mainly manifested as non-linear laws in the spatio-temporal evolution process. This dynamic response includes: (1) There is a lag 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 pollutant concentration shows non-uniform distribution in space, affected by the heterogeneity of the aquifer and the preferential flow path; (3) The pollutants are affected by various physico-chemical processes such as adsorption, desorption, and degradation during the migration process, showing the attenuation characteristics of concentration; (4) The groundwater system has a certain self-purification ability. After the pollution source is removed, the system will gradually return to the initial state, but the recovery process is time-dependent. The three-dimensional anisotropic hydrogeological model can accurately reflect the different characteristics of groundwater flow and pollutant migration in different directions by considering the directional differences of the hydraulic conductivity tensor, thus more realistically simulating the pollutant diffusion law under complex geological conditions.

[0049] The three-dimensional anisotropic hydrogeological model is constructed by 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 elements, and the spatial distribution of unknown variables (head and concentration) is described by linear interpolation functions within each element. After discretization, a large sparse matrix equation system is formed and solved by efficient solvers such as the preconditioned conjugate gradient method. The implicit backward difference format is used for the time term to ensure numerical stability. During grid division, local refinement is carried out around geological interfaces, faults, and pollution sources to improve the calculation accuracy.

[0050] The expressions of the flow motion equation and the pollutant transport equation in the model are as follows:

[0051] Flow motion equation: ;

[0052] Pollutant transport equation: ;

[0053] Where, is the specific storage rate, which represents the amount of water that can be released or absorbed by a unit volume of aquifer under a 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. is the head, which represents the sum of the potential energy height and the pressure head and represents the energy level of groundwater. It is in meters (m). The head determines the flow direction of groundwater (from high head to low head). represents the rate of change of the head with time, indicating the rising and falling speed of the groundwater level; is the spatial coordinate The heterogeneous hydraulic conductivity tensor describes the three-dimensional spatial distribution characteristics of the water-conducting ability of geological media. As a tensor representation, it means that the water-conducting ability may be different in different directions (anisotropy), with the unit of m / day or m / s. A higher K value indicates that water flow through the medium is easier. Represents the spatial distribution change of groundwater flow, that is, the difference in the amount of water flowing into and out of the system; Is the source-sink term, representing the increase or decrease in the amount of water per unit volume per unit time; such as rainfall infiltration (source) or pumping (sink), with the unit of 1 / s or 1 / day. Is the water content per unit volume, representing the volume ratio of water in the medium per unit volume, dimensionless, ranging from 0 to 1; in the saturated zone, it is approximately equal to the porosity.

[0054] Is the pollutant concentration, representing the mass of pollutants in unit volume of water, usually expressed in mg / L or g / m³. Represents the change rate of the mass of pollutants in unit volume of the medium per unit time. Is the hydrodynamic dispersion coefficient tensor, representing the degree and rate of pollutant diffusion, describing the combined effect of pollutant diffusion and mechanical dispersion in groundwater. As a tensor representation, it shows that the dispersion characteristics are different in different directions. Represents the net change in the pollutant mass flux caused by dispersion; Is the groundwater flow velocity vector, describing the velocity and direction of groundwater flow, with the unit of m / day or m / s. It is the actual groundwater flow velocity (rather than the Darcy velocity) and is related to the porosity.

[0055] Represents the concentration reduction rate of pollutants due to processes such as degradation and decay, with the unit of 1 / day or 1 / s. A larger value indicates faster pollutant degradation; Represents the total amount of pollutants degraded or decayed per unit volume per unit time; Is the pollutant attenuation coefficient, Is the source strength of pollutants, representing the amount of water injected into the system per unit volume per unit time, with the unit of 1 / day or 1 / s, Is the source concentration of pollutants, representing the pollutant concentration in the injected water, with the same unit as C. Is the vector differential operator. In a three-dimensional rectangular coordinate system (x, y, z), Is defined as: . These two equations together constitute a coupled system for describing groundwater flow and pollutant transport, considering the interaction between the hydraulics process (flow) and the mass transfer process (diffusion, convection, degradation). The anisotropic representation enables the model to more accurately simulate the behavior of groundwater systems under complex geological conditions.

[0056] Step S3: Calculate the characteristic indexes of groundwater vulnerability based on the anisotropic hydrogeological model. The characteristic indexes include the pollutant arrival time index, the pollutant concentration attenuation index, and the system recovery ability index.

[0057] 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 pollutant concentration attenuation index is calculated based on the ratio of the maximum pollutant concentration at the target point to the initial pollutant concentration. The system recovery ability index is calculated based on the system pollution response time and the recovery time.

[0058] The calculation formula for the pollutant arrival time index RTI is: where is the time required for the pollutant concentration to reach 50% of the initial concentration, which is obtained by solving the pollutant transport equation, is the scale parameter, which is dynamically adjusted according to the aquifer thickness and the hydraulic gradient. is the scale parameter, whose physical meaning is the characteristic time scale and is closely related to the physical properties and hydraulic conditions of the aquifer. The dynamic adjustment equation for the parameter is: is the reference value under the reference condition (usually taken as 100 days), M is the aquifer thickness, is the reference thickness (taken as 10 m), 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 degrees of the aquifer thickness and the hydraulic gradient on the pollutant transport time. These exponent values are obtained by fitting a large amount of numerical simulation and experimental data. In an aquifer with a larger thickness, the vertical migration distance increases, resulting in an extended pollutant arrival time; while an increase in the hydraulic gradient will accelerate the groundwater flow velocity, thus reducing the pollutant arrival time.

[0059] For example: a homogeneous aquifer with a thickness of 10 m, the pollution source is at the coordinate (0, 0, 0), and the 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: the volumetric water content θ = 0.3, the longitudinal dispersion coefficient DL = 10 m² / day; the transverse dispersion coefficient DT = 1 m² / day; the vertical dispersion coefficient DV = 0.1 m² / day; the groundwater flow velocity v = 0.5 m / day (in the x direction); the pollutant decay coefficient λ = ; in this simplified one-dimensional case, an analytical solution can be used.

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

[0061] A three-dimensional finite element model of the aquifer is established, and random sampling of parameters is carried out (for example, DL may vary in the range of 8 - 12 m² / day); For each set of parameters, calculate the change of pollutant concentration with time; Determine the time when the pollutant reaches / 2 at the target point; The spatial step Δx = 1m; The time step Δt = 0.1 day; Within each time step, we calculate the pollutant concentration at each grid point; The following results are obtained through this numerical simulation:

[0062] At t = 180 days, the concentration at the target point (100, 0, 0) is 48 mg / L;

[0063] At t = 185 days, the concentration at the target point is 52 mg / L;

[0064] Through linear interpolation, we can determine that t50 ≈ 182.5 days. Once t50 = 182.5 days is obtained, the RTI can be calculated:

[0065] Assume β = 100 (a scale parameter related to the aquifer thickness and hydraulic gradient): ; This RTI value reflects the resistance of this point to pollutants. A lower value indicates a longer pollutant arrival time and thus lower vulnerability.

[0066] The pollutant concentration attenuation index The calculation formula is: where is the maximum pollutant concentration 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 characterizing the adsorption and natural attenuation ability of the geological medium by pollutants, and its calculation formula is: where is the distribution coefficient (L / kg), indicating the distribution ratio of pollutants between the solid and liquid phases; is the medium bulk density (kg / L); is the effective porosity; is the pollutant degradation rate constant (1 / day); is the actual groundwater flow velocity (m / day); The larger the value, the stronger the adsorption and degradation ability of the medium to pollutants, and the more significant the attenuation of pollutants with distance. For different types of pollutants, the values vary greatly. For example, the value of heavy metals is mainly affected by adsorption, while that of organic pollutants is affected by both adsorption and degradation.

[0067] The system recovery ability index is calculated by the following formula: , where and are the times required for the pollutant concentration to reach 90% and 10% of the initial concentration respectively, obtained by solving the pollutant transport equation. is the initial pollution value, is the final equilibrium concentration, is the time required for system recovery.

[0068] The system recovery time is determined based on the time process of the pollutant concentration recovering to the background value or environmental standard; in actual calculation, an exponential decay model is used to fit the curve of the concentration change with time after the pollution source is removed: ; where k is the decay rate constant, determined by the nonlinear regression method. The system recovery time is defined as the time required for the concentration to drop from the peak to 95% of the final equilibrium concentration , that is ; the term in the numerator is used to standardize the comparison under different initial pollution conditions, making the index dimensionless.

[0069] Step S4: Construct the kinetic coupling groundwater vulnerability index, calculate the time-varying contribution coefficient and nonlinear coupling index of each index and determine the dynamic threshold, and obtain the quantitative evaluation result of groundwater vulnerability in the study area.

[0070] The kinetic coupling groundwater vulnerability index DCG is calculated by the nonlinear coupling equation:

[0071] ; where is the time-varying contribution coefficient, is the nonlinear coupling index, used to characterize the synergy effect between each index.

[0072] The time-varying contribution coefficient is determined by the dynamic threshold method:

[0073] , where is the The sensitivity function of an index, with the calculation formula being:

[0074] ; where, is the corresponding index value, is the variance of the index changing with time, characterizing the spatio-temporal uncertainty of the index.

[0075] The non-linear coupling index is determined through the self-organized critical state theory:

[0076] , where, is the coefficient of variation of the -th index, is the mean value of the coefficient of variation of all indices, is the index 's fractal dimension, characterizing the spatial self-similarity of the index, and is calculated through variogram analysis.

[0077] The non-linear coupling index is determined using the self-organized critical state theory, which holds that a complex system will spontaneously tend towards a critical state during the evolution process. The coefficient of variation of the index is obtained by calculating the ratio of the standard deviation to the mean value of the index within the study area, reflecting the spatial variation degree of the index.

[0078] The fractal dimension characterizes the complexity and self-similarity of the spatial distribution of the index, and is calculated through variogram analysis: First, calculate the semi-variogram at different distance intervals h, and then fit the relationship curve between and in a logarithmic coordinate system. 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

[0079] value, the more complex the spatial distribution of the index, and the more significant the corresponding non-linear coupling effect. ; Search for 's inflection point as the phase transition critical point (calculate the first derivative and second derivative of , and determine the inflection point of the frequency distribution function by solving the point where the first derivative is zero and the second derivative is negative. The DCG value corresponding to this inflection point is the phase transition critical point); Based on determine the dynamic division threshold for the groundwater vulnerability level.

[0080] The dynamic division threshold for the groundwater vulnerability level includes:

[0081] Very low vulnerability: ; Low vulnerability: Medium vulnerability: ; High vulnerability: ; Extremely high vulnerability: .

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

[0083] To further illustrate the technical effects of the present invention, this method is applied taking the area around a certain industrial park as an example. The area of this region is about 25 square kilometers, the groundwater is the main drinking water source, there are multiple potential pollution sources, and the geological conditions are complex. First, the groundwater vulnerability evaluation parameters of the study area are obtained through on-site investigation: the aquifer parameters include the specific storage rate range , 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 volumetric 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; for the pollutant characteristic parameters, nitrate is selected as the indicator pollutant, and the longitudinal , transverse , vertical components of the hydrodynamic dispersion coefficient tensor, the attenuation coefficient λ is , and the pollution source intensities of the industrial area and the landfill are 0.08 m^3 / d and 0.05 m^3 / d respectively, and the concentrations are 120 mg / L and 85 mg / L respectively.

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

[0085] By calculating the non - linear coupling index ( = 2.44, = 2.88, = 2.72) and the time - varying contribution coefficient ( = 0.28, = 0.42, = 0.30), the dynamic coupling groundwater vulnerability index DCG of monitoring well R8 is obtained as 0.482. Based on the distribution of DCG values in the study area, the phase - change critical point DCGcr = 0.525 is identified, and the vulnerability levels are 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 in the study area is 12.6% for the extremely low vulnerability area, 24.8% for the low vulnerability area, 31.5% for the medium vulnerability area, 21.2% for the high vulnerability area, and 9.9% for the extremely high vulnerability area; the areas around industrial zones and landfills are mainly high and extremely high vulnerability, and the vulnerability shows a gradual change characteristic in the downstream area of the groundwater flow.

[0086] Comparing this method with the traditional DRASTIC method, the prediction accuracy for 15 pollution verification points is increased 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 dry 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 risks in the mutation areas. The actual application results show that the method of the present invention has significant advantages in evaluation accuracy, dynamic response ability, adaptability to heterogeneity and anisotropy, decision - making support ability, and early warning ability, providing a scientific basis for the protection and management of groundwater resources.

[0087] Example 2

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

[0089] 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 collection and processing functions, supports the import and parsing of original data such as borehole data, pumping tests, and laboratory permeability tests, and constructs a parameter spatial distribution field through geostatistical methods. Taking the application in an industrial park as an example, this module uses the multi-point geostatistical method 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, this 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.

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

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

[0092] The said vulnerability evaluation module is used to construct the dynamic coupling groundwater vulnerability index, calculate the time-varying contribution coefficient and non-linear coupling index of each index, and determine the dynamic threshold to obtain the quantitative evaluation result of the groundwater vulnerability in the study area.

[0093] The vulnerability evaluation module innovatively introduces the complex system theory and non-linear coupling mechanism, breaking through the limitations of traditional evaluation methods. This module uses the global sensitivity analysis method to calculate the time-varying contribution coefficient, which can identify the key control factors in different periods and different regions and realize the dynamic adjustment of the evaluation. For the non-linear relationship between indexes, the module has developed a coupling index calculation method based on the self-organized critical state theory. Through the analysis of the coefficient of variation and fractal dimension, it quantifies the spatial variation characteristics and self-similarity of the indexes. The dynamic threshold determination function uses a data-driven method. Through the inflection point detection algorithm of the frequency distribution function, it automatically identifies the critical point of the system state transition and generates a scientific and reasonable grading standard accordingly. Compared with the traditional fixed threshold, this method better adapts to the geological and hydrogeological 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 inter-annual changes, seasonal fluctuations, and extreme events on groundwater vulnerability, realizing the methodological leap from static evaluation to dynamic evaluation.

[0094] The said system also includes: a result visualization module, which is used to generate the groundwater vulnerability zoning map of the study area using the geographic information system according to the dynamic threshold division result and provide the spatio-temporal evolution analysis function; a decision support module, which is used to provide suggestions on groundwater protection area division, pollution risk warning, and management countermeasures based on the groundwater vulnerability evaluation result, combined with land use planning and environmental management requirements.

[0095] The vulnerability assessment module innovatively introduces the complex system theory and the non-linear coupling mechanism, breaking through the limitations of traditional assessment methods. This module uses the global sensitivity analysis method to calculate the time-varying contribution coefficients, which can identify the key control factors in different periods and regions, and achieve dynamic adjustment of the assessment. For the non-linear relationship between indicators, the module develops a coupling index calculation method based on the self-organized critical state theory. Through the analysis of the coefficient of variation and fractal dimension, it quantifies the spatial variability characteristics and self-similarity of the indicators. The dynamic threshold determination function adopts a data-driven method. Through the inflection point detection algorithm of the frequency distribution function, it automatically identifies the critical points of the system state transition, and generates a scientific and reasonable grading standard accordingly. Compared with the traditional fixed threshold, this method better adapts to the geological and hydrogeological conditions of different regions, and the assessment results are more objective and accurate. In addition, the module also supports multi-time scale assessment, which can analyze the impacts of inter-annual changes, seasonal fluctuations and extreme events on groundwater vulnerability, achieving a methodological leap from static assessment to dynamic assessment.

[0096] The decision-making support module transforms the scientific assessment results into practical management measures, which reflects the application value of the system. Based on the multi-criteria decision-making theory, this module integrates the groundwater vulnerability assessment results and socio-economic factors to generate an optimized groundwater protection area division plan. The protection area grading function considers various factors such as vulnerability levels, pollution source distributions, current land use status and planned development needs. Through spatial overlay analysis, it automatically generates the boundaries of different protection levels and suggestions for control measures. The pollution risk early warning subsystem realizes 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 persons. The management countermeasure suggestion function, based on the case library and knowledge graph, provides targeted prevention and control measure suggestions for different vulnerability regions and pollution types, such as optimizing the land use structure, adjusting the industrial layout, and improving the agricultural fertilization method. In addition, this module also supports the scenario simulation function, which can evaluate the effects and cost-benefits of different management measures, assist decision-makers in selecting the optimal management strategy, and achieve the sustainable utilization and effective protection of groundwater resources.

[0097] Combined with the application case of Embodiment 1, each functional module of the groundwater vulnerability assessment system based on the hydrogeological model is introduced in detail.

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

[0099] The parameter acquisition module is responsible for collecting and processing various parameter data of the research area. Taking the industrial park in Embodiment 1 as an example, this module integrates multi-source data acquisition interfaces and can simultaneously access borehole data, hydrogeological parameters, groundwater monitoring data, and pollutant characteristic parameters. For parameters such as specific 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 the parameters through spatial statistical methods. For example, aiming at the heterogeneity characteristics of the hydraulic conductivity in this area, the system uses 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.

[0100] The model construction module is the core component of the system and is responsible for constructing a three-dimensional anisotropic hydrogeological model. Aiming at the complex geological conditions in Embodiment 1, this module adopts adaptive grid meshing technology to locally refine the grids at key geological interfaces and around pollution sources. The total number of tetrahedral elements reaches 128,000, effectively ensuring the calculation accuracy under complex geological structures. The model solving engine adopts parallel computing technology, and the solving efficiency of the water flow motion equation and the pollutant transport equation is increased by about 3.5 times. The model calibration function realizes the automatic optimization of hydraulic parameters by integrating the PEST parameter estimation program, reducing the mean absolute error between the simulated water head and the measured water head to 0.22m, and controlling the relative error of pollutant concentration simulation within 8.5%. The model prediction ability is significantly enhanced.

[0101] The index calculation module calculates the characteristic indexes of groundwater vulnerability based on the calibrated hydrogeological model. Aiming at the monitoring well network in Embodiment 1, this module calculates the pollutant arrival time index RTI, pollutant concentration decay index CDI, and system recovery ability 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 index values of RTI = 0.487, CDI = 0.732, and SRI = 0.308. For the entire research area, the system calculates the spatial distribution of each index at a grid density of 25m×25m and generates a high-precision index distribution map, intuitively showing the spatial heterogeneity of groundwater vulnerability.

[0102] The vulnerability assessment module is the decision-making core of the system, responsible for constructing the dynamic coupling groundwater vulnerability index and determining the dynamic threshold. For the industrial park in Embodiment 1, this module first calculated the coefficient of variation and fractal dimension of each index, obtaining the non-linear coupling indices γ1 = 2.44, γ2 = 2.88, γ3 = 2.72; then determined the time-varying contribution coefficients w1 = 0.28, w2 = 0.42, w3 = 0.30 based on sensitivity analysis; and finally calculated the dynamic coupling groundwater vulnerability index DCG. The system innovatively introduced the self-organized critical state theory to determine the dynamic threshold. Through the inflection point analysis of the frequency distribution function, the phase transition critical point DCGcr = 0.525 was identified, and based on this, a five-level vulnerability classification standard was automatically generated. Compared with the traditional fixed threshold method, the dynamic threshold method improved the accuracy of vulnerability identification in the study area by 18.7%, especially significantly enhancing the identification ability in areas with sudden geological condition changes.

[0103] The result visualization module provides rich visualization and analysis functions. For the evaluation results in Embodiment 1, this module generated a groundwater vulnerability zoning map of the study area, using five colors to identify different vulnerability levels and supporting multi-scale spatial browsing. The spatio-temporal evolution analysis function can simulate the vulnerability changes under different seasonal conditions. For example, during the rainy season, the system identified a moderately vulnerable area on the east side of the industrial zone that would temporarily transform into a highly vulnerable area, providing a scientific basis for seasonal management measures. In addition, the system also supports scenario simulation functions. For example, simulating the impact of newly added pollution sources or changes in land use patterns on groundwater vulnerability. The evaluation results show that if a chemical plant is newly built in the northwest of the study area, it will cause the vulnerability level in a 3.2-square-kilometer area downstream to increase by 1 - 2 levels, providing a quantitative reference for regional planning and environmental impact assessment.

[0104] The decision support module combines the evaluation results with management requirements to provide scientific decision-making suggestions. For the highly vulnerable and extremely highly vulnerable areas (accounting for 31.1% of the total area) identified in Embodiment 1, the system automatically generated a hierarchical protection plan, suggesting that the extremely highly vulnerable area (9.9%) be designated as a first-level protection area, prohibiting any activities that may pollute groundwater; and the highly vulnerable 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 proposed a differential monitoring plan based on the spatio-temporal evolution characteristics of vulnerability, deploying 12 new monitoring wells in the highly vulnerable area with a monitoring frequency of once a week; and 8 monitoring wells in the moderately vulnerable area with a monitoring frequency of once a month, significantly improving the pertinence and economy of the monitoring network.

[0105] Through the collaborative work of the above modules, the system has been successfully applied to the groundwater vulnerability assessment of the industrial park in Embodiment 1. Compared with the traditional assessment method, the assessment accuracy has been improved by 24.3%, and it can reflect the dynamic changes in vulnerability caused by seasonal variations and human activities. The automation level and user-friendliness of the system have increased the assessment work efficiency by about 5 times, especially showing significant advantages in the assessment ability and dynamic response ability under complex geological conditions. The system has been successfully applied to the groundwater resource management of multiple industrial parks, cities and basins, providing a scientific basis and technical support for the delineation of groundwater protection areas, pollution risk early warning and the formulation of management countermeasures.

[0106] The specific embodiments described above further elaborate on the purpose, technical solutions and beneficial effects of the present invention. It should be understood that the above are only specific embodiments of the present invention and are not used to limit the protection scope of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for evaluating groundwater vulnerability based on a hydrogeological model, which obtains a quantitative evaluation result of the groundwater vulnerability in the study area based on a pollutant dynamics coupling method, is characterized in that The evaluation method includes the following steps: Step S1: Obtain the groundwater vulnerability evaluation parameters of the study area. The groundwater vulnerability evaluation parameters include aquifer parameters, groundwater flow parameters, and pollutant characteristic parameters; Step S2: According to the dynamic response characteristics of the groundwater system to pollutants, construct a three-dimensional anisotropic hydrogeological model composed of a groundwater flow equation and a pollutant transport equation; The three-dimensional anisotropic hydrogeological model is constructed using the tetrahedral finite element discretization method. The expressions of the groundwater flow equation and the pollutant transport equation in the model are respectively: Equation of water flow motion: ; Pollutant transport equation: ; Among them, is the specific yield, which represents the amount of water that can be released or absorbed by a unit volume of aquifer under a unit change in hydraulic head, is the hydraulic head, which represents the sum of the potential energy height and the pressure head and represents the energy level of groundwater, represents the rate of change of hydraulic head with time and represents the rising and falling speed of the groundwater level; is the spatial coordinate is the heterogeneous hydraulic conductivity tensor. The larger the value of K, the easier it is for water to flow through the medium; represents the spatial distribution change of groundwater flow, that is, the difference in the amount of water flowing into and out of the system; is the source-sink term, which represents the increase or decrease in the amount of water in a unit volume per unit time; is the moisture content per unit volume, is the pollutant concentration, is the hydrodynamic dispersion coefficient tensor, which represents the degree and rate of pollutant diffusion, represents the net change in the pollutant mass flux caused by dispersion; is the groundwater flow velocity vector, is the pollutant decay coefficient, is the source strength of pollution, which represents the amount of water injected into the system per unit volume per unit time, is the source concentration of pollution, which represents the pollutant concentration in the injected water; is the vector differential operator. In a three-dimensional rectangular coordinate system (x, y, z), is defined as: ; Step S3: Calculate the characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model. The characteristic indicators include pollutant arrival time index, pollutant concentration attenuation index, and system recovery ability index; The pollutant arrival time index RTI is calculated based on the time required for the pollutant concentration to reach 50% of the initial concentration, and the pollutant concentration attenuation index is calculated based on the ratio of the maximum pollutant concentration at the target point to the initial pollutant concentration, and the system recovery ability index is calculated based on the system pollution response time and recovery time; The calculation formula for the pollutant arrival time index RTI is as follows: , where is the time required for the pollutant concentration to reach 50% of the initial concentration, obtained by solving the pollutant transport equation, is the scale parameter, dynamically adjusted according to the aquifer thickness and hydraulic gradient; The pollutant concentration decay index The calculation formula is as follows: , where is the maximum pollutant concentration at the target point, is the initial pollutant concentration, is the distance from the pollution source to the target point, is the medium adsorption decay coefficient; The system recovery ability index The calculation formula is as follows: , where and are the times required for the pollutant concentration to reach 90% and 10% of the initial concentration respectively, obtained by solving the pollutant transport equation, is the initial pollution value, is the final equilibrium concentration, is the time required for system recovery; Step S4: Construct a kinetic coupling groundwater vulnerability index, calculate the time-varying contribution coefficients and non-linear coupling index of each indicator, and determine the dynamic threshold to obtain the quantitative evaluation result of groundwater vulnerability in the study area; The kinetic coupling groundwater vulnerability index DCG is calculated using a non-linear coupling equation: ; wherein, is a time-varying contribution coefficient, is a non-linear coupling index used to characterize the synergy between various indicators; The time-varying contribution coefficient is determined by a dynamic threshold method: , where is the sensitivity function of the th index, and the calculation formula is: ; where, is the corresponding index value, is the index variance varying with time, representing the spatio-temporal uncertainty of this index; The non-linear coupling index is determined by the self-organized critical state theory: , where is the coefficient of variation of the th index, is the mean of the coefficients of variation of all indices, is the fractal dimension of the index , which characterizes the spatial self-similarity of the index and is calculated through variogram analysis.

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

3. The groundwater vulnerability assessment method based on a hydrogeological model according to claim 2, wherein, The dynamic threshold is determined by the self-organized critical state theory: calculate the frequency distribution function of the DCG values in the study area ; Find the inflection point of as the phase transition critical point ; Based on determine the dynamic division threshold for the groundwater vulnerability level; The dynamic division thresholds of the groundwater vulnerability levels include: Extremely low vulnerability: ; Low vulnerability: Medium vulnerability: ; High vulnerability: ; Extremely high vulnerability: ; According to the dynamic threshold division result, use a geographic information system to generate a groundwater vulnerability zoning map of the study area.

4. A groundwater vulnerability assessment system based on a hydrogeological model for performing the method according to any one of claims 1-3, characterized in that, The system includes, connected in sequence: a parameter acquisition module, a model construction module, an index calculation module, and a vulnerability evaluation module; The parameter acquisition module is used to obtain the groundwater vulnerability evaluation parameters of the study area, including aquifer parameters, groundwater flow parameters, and pollutant characteristic parameters; The model construction module is used to construct a three-dimensional anisotropic hydrogeological model composed of a groundwater flow equation and a pollutant transport equation according to the dynamic response characteristics of the groundwater system to pollutants; The index calculation module is used to calculate the characteristic indicators of groundwater vulnerability based on the anisotropic hydrogeological model, including pollutant arrival time index, pollutant concentration attenuation index, and system recovery ability index; The vulnerability evaluation module is used to construct a kinetic coupling groundwater vulnerability index, calculate the time-varying contribution coefficients and non-linear coupling index of each indicator, and determine the dynamic threshold to obtain the quantitative evaluation result of groundwater vulnerability in the study area.

5. The groundwater vulnerability assessment system based on a hydrogeological model according to claim 4, characterized in that, The system further includes: a result visualization module, which is used to generate a groundwater vulnerability zoning map of the study area using a geographic information system according to the dynamic threshold division result, and provide a spatio-temporal evolution analysis function.

Images