Water quality pollution risk assessment system based on underground water

Abstract

The invention discloses a water quality pollution risk assessment system based on underground water, and particularly relates to the technical field of underground water pollution risk assessment, comprising the following steps: collecting data through a data monitoring module; the dynamic evaluation module generates a risk evaluation benchmark based on the numerical model; the entropy change calculation module and the flux calculation module carry out dynamic diagnosis from two aspects of reaction network complexity and toxicity transport intensity respectively; a game regulation and control module queries a preset game matrix according to the dynamic indexes, and generates a regulation and control strategy with micro-disturbance detection and dark response inversion functions; and the precise intervention module executes the strategy to guide pollution conversion. According to the method, dynamic early warning of hidden risks such as dark response is realized, and system abnormity is captured through entropy change and flux indexes; a self-adaptive regulation and control closed loop based on game learning is established, and risk response is converted into active guidance from passive response; finally, targeted precise repair is achieved, and secondary risks and resource waste caused by mechanism misjudgment in traditional measures are avoided.

Description

Technical Field

[0001] This invention relates to the field of groundwater pollution risk assessment technology, and more specifically, to a groundwater-based water pollution risk assessment system. Background Technology

[0002] Traditional groundwater pollution risk assessment and management systems are based on classical hydrogeology, solute transport theory, and environmental chemistry. Their core technical approach typically follows a linear model of "monitoring-simulation-assessment-management": first, pollutant concentrations and hydrogeological parameters are obtained through on-site investigation and monitoring; second, numerical models are used to simulate the spatiotemporal evolution of the pollution plume; third, risk quantification and zoning are performed based on the simulation results; and finally, monitoring or remediation plans are formulated according to the risk level. This methodology can provide effective decision support for scenarios with clear pollution sources, relatively stable hydrogeological conditions, and well-defined pollutant transformation pathways.

[0003] However, with the deepening of industrialization and the increasing complexity of pollutants, the limitations of traditional technological paradigms in addressing new and complex pollution scenarios are becoming increasingly apparent. The core challenge stems from the fact that groundwater systems are complex mega-systems characterized by open, dynamic, multiphase, and multi-component reactions, exhibiting highly nonlinear and uncertain behavior. A particularly prominent and unresolved problem is the "dark reaction and toxic evolution risk under the coupling of multiple pollutants." Specifically, when chlorinated hydrocarbon solvents from industrial sites, heavy metals from metal smelting sites, nitrogen and phosphorus nutrients from agricultural activities, and emerging organic pollutants coexist in aquifers, they do not migrate independently. In complex geochemical gradient fields (such as redox transition zones), catalyzed by indigenous microbial communities, or through interfacial reactions on mineral surfaces, these pollutants may undergo a series of undesigned and unintended chain reactions—the "dark reactions." For example, in a reducing environment rich in organic matter and sulfates, the reductive dechlorination process of trichloroethylene may stall at the highly toxic vinyl chloride stage due to competition for electron donors or changes in the microbial community structure, rather than being completely dechlorinated to harmless ethylene; and heavy metal ions may indirectly affect the degradation pathways and rates of organic pollutants by altering enzyme activity or acting as cofactors.

[0004] The existing technological system exhibits the following disconnects when addressing this fundamental challenge: At the perception level, traditional monitoring focuses on the concentration of limited target pollutants and lacks the ability to conduct high-frequency, in-situ monitoring of the full spectrum of water chemistry fingerprints and overall biotoxicity. It cannot capture early chemical signals (such as the emergence of specific intermediate products) and effect signals (such as mutations in comprehensive toxicity) that indicate the occurrence of "dark reactions".

[0005] At the cognitive and predictive level, the core reaction modules of mainstream numerical models rely on simplified, pre-defined chemical kinetic equations, which are essentially parameterizations of known reaction pathways. For unknown, "emergent" dark reaction networks, the models have neither input parameters nor corresponding mathematical expressions, so their predictions carry the risk of "unknown unknowns" and may deviate significantly from reality.

[0006] At the decision-making and intervention level, risk management strategies based on static model predictions are essentially reactive. They cannot proactively detect the system's vulnerability and reaction potential before drastic changes in toxicity occur, and they lack a "diagnostic treatment" logic that can dynamically adjust based on system feedback. When model predictions fail due to latent reactions, the entire control system loses its basis for action, often only able to respond to pollution events that have already occurred with a lag. Therefore, this invention proposes a groundwater-based water pollution risk assessment system to address the aforementioned problems. Summary of the Invention

[0007] To achieve the above objectives, the present invention provides the following technical solution: A groundwater-based water pollution risk assessment system includes the following steps: The data monitoring module is used to collect data on the concentration of multiple pollutants, hydrogeological parameters, and geochemical indicators in groundwater through a sensor network and automatic sampling devices deployed in the aquifer. The dynamic assessment module is used to simulate the release probability, transport path, spatial distribution and concentration evolution of pollutants based on the collected multidimensional data and through a pre-set numerical model coupled with groundwater flow and solute transport, and output dynamic risk assessment results including pollution plume range prediction and risk level classification. The entropy change calculation module is used to calculate and output the entropy change value of the reaction network based on the collected multidimensional chemical fingerprint data, by analyzing the diversity of chemical species, the temporal changes in concentration, and the uncertainty of reaction pathways. This is used to quantify the dynamic complexity of the underground pollutant reaction network. The flux calculation module is used to acquire the comprehensive toxicity data measured by online biotoxicity sensors deployed at key hydrological sections downstream of the pollution plume, and combine it with a pre-set groundwater seepage velocity model to calculate and output the downstream toxicity equivalent flux. This flux characterizes the comprehensive toxicity transport intensity through the monitoring section per unit time. The game-theoretic control module receives data on the entropy change of the reaction network and the downstream toxicity equivalent flux. By querying a game matrix pre-constructed based on different geochemical backgrounds and pollutant combinations, it maps and generates a dynamic control strategy. This strategy includes micro-perturbation injection commands for active detection and feedback analysis logic for inverting the dark reaction type based on the system response. The precision intervention module is used to execute the received dynamic control strategy. By injecting targeted co-metabolite matrix, catalytic material or electron donor and acceptor into the groundwater through the injection well network, it guides the pollutant transformation pathway to the direction of water safety or complete mineralization.

[0008] The sensor network includes at least a water quality sensor and a water level sensor, and an automatic sampling device is used to collect water samples according to a preset cycle.

[0009] In a preferred embodiment, in the dynamic evaluation module: The pre-set numerical model of groundwater flow and solute transport is constructed with a computational kernel based on the finite difference method or the finite element method. The coupled numerical model uses hydrogeological parameters collected by the data monitoring module as the basic input parameters for model construction, calibration and boundary condition setting, and multi-pollutant concentration data as the input of solute concentration field to simulate the release intensity and initial spatial distribution of pollution sources. The simulation execution process of the coupled numerical model is as follows: First, the partial differential equations describing the movement of groundwater are solved to obtain the head distribution and velocity field of the study area; then, the convection-diffusion-reaction equations describing the migration and transformation of pollutants are solved in a coupled manner, and the complete dataset of the dynamic changes of pollutant concentration in the spatial and temporal dimensions is iteratively calculated and output. This dataset constitutes the intermediate results of concentration distribution for risk assessment.

[0010] The dynamic risk assessment results generated by the dynamic assessment module are obtained through further processing of intermediate results. The dynamic risk assessment results consist of two parts: the first part is the prediction of the spatial extent boundary of the pollution plume within a future set time period, extracted based on intermediate results of concentration distribution; The second part involves converting the predicted exposure concentration of pollutants into a human health risk index, and classifying them according to legally mandated or pre-set risk acceptable level thresholds, ultimately generating a spatial zoning map of carcinogenic risk and a spatial zoning map of non-carcinogenic risk.

[0011] In a preferred embodiment, the partial differential equations describing the groundwater movement are solved to obtain the hydraulic head distribution and velocity field of the study area, which is achieved through the following steps: Based on hydrogeological parameters, the underground space of the study area is divided into three-dimensional grid units; On each grid cell, Darcy's law and the law of conservation of mass are applied to establish a mathematical relationship describing the flow balance of the groundwater cell; The mathematical relationships established by all grid cells in the study area are combined and substituted with the boundary conditions determined by the water level data collected by the data monitoring module to form a closed linear equation system. An iterative solution algorithm is used to numerically calculate the closed linear equation system, and the hydraulic head value at the center node of each grid cell is obtained. The set of all node hydraulic head values ​​constitutes the hydraulic head distribution. Based on the head distribution, the hydraulic conductivity coefficient in the hydrogeological parameters, and Darcy's law, the groundwater seepage velocity in the three-dimensional direction of each grid cell is calculated and output. The collection of seepage velocities of all grid cells constitutes the velocity field.

[0012] In a preferred embodiment, iteratively calculating and outputting a complete dataset of the dynamic changes in pollutant concentrations in spatial and temporal dimensions is achieved through the following steps: Using the velocity field as the input of the convection term driving pollutant transport, the hydrodynamic dispersion term of pollutants is calculated based on the dispersion and velocity field in hydrogeological parameters. Construct reaction terms that describe the biodegradation, chemical redox, and adsorption-desorption processes that occur during the migration of pollutants; At each 3D grid cell and at each simulation time step, the convection, dispersion, and reaction terms are coupled and solved simultaneously to update the pollutant concentration of the grid cell at the end of that time step. Using the updated pollutant concentration as the initial value for the next time step, the above coupled solution process is repeated to achieve iterative calculation for all simulation time steps; finally, the pollutant concentration set of all grid cells at all simulation time steps is output, which constitutes a complete dataset of pollutant concentration dynamically changing in spatial and temporal dimensions.

[0013] In a preferred embodiment, the specific process of reacting to the network entropy change is as follows: The first step is data preparation and preprocessing: receiving multidimensional chemical fingerprint data continuously collected by the data monitoring module and organizing it into a concentration matrix with time as the sequence and chemical species as the dimension; The second step is to calculate the entropy of species diversity components: For each monitoring period, the relative concentration ratio of all detected chemical species is statistically analyzed and regarded as a probability distribution. The Shannon information entropy of this distribution is calculated to obtain the first entropy value component that characterizes species richness and evenness. The third step is to calculate the dynamic entropy of the concentration process: For each chemical species, extract its concentration value over the entire monitoring time series and calculate the permutation entropy of the concentration time series; the permutation entropy quantifies the disorder and unpredictability of concentration changes by comparing the pattern complexity formed by the relative size relationship of adjacent data points in the sequence, and integrates the permutation entropy of all species to obtain the second entropy value component. The fourth step is to calculate the entropy of the reaction path inference: Based on the pre-set knowledge graph of pollutant transformation paths, the reaction network graph of all chemical species detected in the current period is reconstructed; by analyzing the degree distribution of each node in the graph (i.e., the number of possible reactions in which a species participates) and the edge weight (estimated based on the covariance of reactant and product concentrations), the uncertainty entropy of the reconstructed network graph structure is calculated to obtain the third entropy value component. Step 5: Synthesis and output of entropy change value: Input the first, second and third entropy value components into a preset nonlinear fusion function. This function simulates the synergistic and antagonistic effects of the three entropy values ​​in reflecting the degree of chaos in the system, and outputs a normalized scalar value, which is the entropy change value of the network.

[0014] In a preferred embodiment, the specific process of downstream toxicity equivalent flux is as follows: Simultaneously acquire comprehensive toxicity data measured by one or more online biotoxicity sensors deployed at key hydrological sections, as well as current groundwater seepage velocity vector field data corresponding to the location of the monitoring section provided by the dynamic assessment module; The raw comprehensive toxicity data output by the online biotoxicity sensor is converted into equivalent concentration values ​​expressed in uniform toxicity equivalent units based on a preset dose-response relationship standard curve. The key hydrological monitoring section is discretized into multiple calculation units in the vertical direction. For each calculation unit, the groundwater seepage velocity vector at its location is obtained, and the component of the velocity vector in the direction perpendicular to the monitoring section is calculated. The vertical velocity component is multiplied by the toxicity equivalent concentration value of the corresponding unit, and then multiplied by the representative area of ​​the calculation unit to obtain the toxicity equivalent flux through the unit. The toxicity equivalent flux of all computational units is summed to obtain the total downstream toxicity equivalent flux passing through the entire monitoring section, which is then output as an indicator characterizing the comprehensive risk intensity of the pollution plume transport downstream.

[0015] In a preferred embodiment, the specific process of generating a dynamic control strategy by querying the game matrix mapping in the game control module is as follows: Based on the baseline risk partition map generated by the first round of simulation in the dynamic evaluation module, initial row and column strategies are set for the game matrix. The row strategy corresponds to different potential hidden reaction type assumptions, and the column strategy corresponds to different micro-perturbation intervention actions. The value of each cell in the matrix represents the initial value of the expected comprehensive benefit score when the corresponding column action is used to intervene in the corresponding row reaction type under a specific risk level area. The entropy change value of the reaction network output by the entropy change calculation module and the downstream toxicity equivalent flux output by the flux calculation module are used as a set of joint state indices; based on the numerical range of the joint state index, a corresponding sub-matrix block in the game matrix is ​​located. Within the located sub-matrix block, find the cell with the highest current comprehensive benefit score; parse the column strategy corresponding to the cell, i.e. the specific micro-perturbation intervention action, into an executable micro-perturbation injection instruction, which at least includes the type, concentration, dosage, and injection point coordinates of the injected reagent; at the same time, use the row strategy corresponding to the cell, i.e. the assumed dark reaction type, as the core target for verification and inversion in the feedback analysis logic.

[0016] In a preferred embodiment, the baseline plume prediction map and baseline risk zoning map generated by the dynamic assessment module are used to directly guide the system initialization settings, specifically including: To guide the deployment of key hydrological monitoring sections in the flux calculation module: Based on the maximum migration distance and mainstream direction of the pollution plume indicated in the baseline pollution plume prediction map, the system plans and generates at least one virtual monitoring section perpendicular to the groundwater flow direction at a preset safe distance downstream of the pollution source; the spatial coordinates and geometric range of this section are set as the actual deployment location and monitoring range of the online biotoxicity sensor. The system is used to provide initial strategy mapping relationships for the preset game matrix in the game control module: the system analyzes the risk level of different geographical areas in the benchmark risk zoning map and uses it as the initialization condition of the game matrix; for each risk level zone, the system assigns a set of potential dark reaction type assumptions related to the typical pollutant combination in the region to the row strategies in the game matrix according to the preset rules, and assigns a set of recommended micro-perturbation basic actions that match the hydrogeochemical conditions of the region to the column strategies, thereby completing the filling of the initial values ​​of the game matrix.

[0017] The technical effects and advantages of this invention are as follows: This invention achieves quantitative perception and forward-looking early warning of the complexity of groundwater pollutant transformation processes by introducing two synergistic dynamic diagnostic indicators: the entropy change value of the reaction network and the downstream toxicity equivalent flux. Unlike traditional models that predict based solely on the concentration of limited target substances, this invention utilizes multidimensional chemical fingerprint data to calculate entropy changes, enabling sensitive identification of abnormal fluctuations and uncertainty leaps in pollutant types, concentration sequences, and reaction pathways, thereby capturing early chemical signals of "dark reactions." Simultaneously, by fusing online biotoxicity data and seepage models to calculate toxicity flux, it achieves dynamic assessment of the comprehensive hazard transport intensity of pollution plumes. This dual diagnostic mechanism allows the system to issue timely warnings during the latency period when toxic intermediates are generated in large quantities and migrate downstream, significantly overcoming the "prediction blind spot" caused by the limitations of model-preset pathways and gaining a critical time window for risk intervention.

[0018] This invention constructs an adaptive control closed loop based on game theory, encompassing "perception-diagnosis-decision-feedback," elevating risk response from static contingency plan execution to dynamic strategy optimization. The system maps entropy changes and flux data to control strategies with both "active detection" and "precise guidance" functions through a pre-set game matrix. Specifically, the system intelligently selects and implements micro-perturbations (such as injecting specific probe reagents) based on the current state, and inversely infers the types and mechanisms of latent dark reactions by analyzing the groundwater system's response to these perturbations. This "dynamic control" strategy enables the system not only to adapt to complex and ever-changing geochemical conditions but also to proactively explore unknown risks and learn to optimize, achieving a paradigm shift from passively responding to known pollution to proactively managing uncertainty, significantly improving the intelligence level and decision reliability of risk management.

[0019] This invention establishes remediation intervention on a solid foundation of dynamic diagnosis and game-theoretic decision-making, achieving a leap from "extensive remediation" to "targeted regulation." The precision intervention module executes instructions optimized by game-theoretic logic, capable of delivering highly targeted co-metabolite matrices, catalytic materials, or electron donors and acceptors to specific identified reaction bottlenecks or harmful pathways. For example, when the diagnosis indicates that the reductive dechlorination process is stalled at the vinyl chloride stage due to insufficient electron donors, the system can precisely inject slow-release lactate and specific dehalogenated bacterial activators. This "diagnosis first, treatment later" model avoids the ineffective investment, remediation rebound, and even secondary pollution risks caused by mechanistic misjudgments in traditional remediation, thereby fundamentally improving the technical effectiveness, long-term safety, and resource economy of remediation measures, providing an innovative and sustainable solution for addressing complex contaminated sites. Attached Figure Description

[0020] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings; Figure 1 This is a schematic diagram of the groundwater-based water pollution risk assessment system of the present invention. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0022] Reference Figure 1 The following examples were obtained: Example 1: A groundwater-based water pollution risk assessment system, comprising the following steps: The data monitoring module, through a sensor network and automatic sampling devices deployed in the aquifer, collects multi-pollutant concentration data, hydrogeological parameters, and geochemical indicators from groundwater. This module constitutes the "nerve endings" of the entire system for sensing the state of the groundwater environment. Its significance lies in acquiring multi-dimensional and heterogeneous basic data on pollutant concentrations, hydrogeological conditions, and the geochemical environment through diverse sensing and sampling methods. These real-time and near-real-time data streams are the sole starting point and factual basis for all subsequent simulations, assessments, calculations, and decisions, ensuring the system's objective perception of the dynamics of the underground environment and laying the data foundation for shifting from static assessment to dynamic control.

[0023] The dynamic assessment module, based on collected multidimensional data and a pre-set numerical model coupling groundwater flow and solute transport, simulates the release probability, transport path, spatial distribution, and concentration evolution of pollutants. It outputs dynamic risk assessment results, including predictions of the pollution plume's extent and risk level classification. The significance of this module lies in constructing and running a "digital twin" benchmark model based on first principles of physicochemical principles. Utilizing monitoring data to drive the pre-set numerical model, it simulates the entire process of pollutant transport and transformation in the underground medium, thereby generating predictions of the spatial extent of the pollution plume and classifications of potential health risks over a future period. Its output provides a "static" spatiotemporal prediction baseline for risk, serving not only as the basis for traditional risk management decisions but also as the initial reference coordinate system and strategic planning foundation for all subsequent dynamic diagnostic and control actions.

[0024] The entropy change calculation module, based on collected multidimensional chemical fingerprint data, analyzes chemical species diversity, concentration temporal changes, and reaction pathway uncertainties to calculate and output the entropy change value of the reaction network. This quantifies the dynamic complexity of the underground pollutant reaction network. The significance of this module lies in achieving a quantitative diagnosis of the inherent complexity and instability of the entire groundwater pollutant reaction network. Moving beyond a focus on the concentration of a single pollutant, it calculates a comprehensive index characterizing the system's "disorder" or "unpredictability"—the reaction network entropy change value—by analyzing the diversity, dynamic behavior, and interactions of chemical species. This index can keenly capture early signals of anomalous processes such as "dark reactions" that standard models cannot describe, representing a crucial intelligent leap from "predicting the known" to "detecting the unknown."

[0025] The flux calculation module acquires comprehensive toxicity data measured by online biotoxicity sensors deployed at key hydrological sections downstream of the pollution plume. Combined with a pre-set groundwater seepage velocity model, it calculates and outputs the downstream toxicity equivalent flux. This flux characterizes the intensity of comprehensive toxicity transport through the monitoring section per unit time. The significance of this module lies in measuring the intensity and flux of the comprehensive risk of pollution plume transport to downstream sensitive targets. By integrating online biotoxicity sensors (reflecting the comprehensive effect) and the seepage velocity model (reflecting transport capacity), the "downstream toxicity equivalent flux" is calculated. This indicator organically combines the toxicity and concentration of pollutants with groundwater flow, dynamically and intuitively characterizing the rate of risk "dose" transport, providing a direct and quantitative decision-making benchmark for assessing risk urgency and triggering different levels of emergency responses.

[0026] The game theory-based regulation module receives data on the entropy change of the reaction network and the downstream toxicity equivalent flux. By querying a pre-constructed game matrix based on different geochemical backgrounds and pollutant combinations, it maps and generates a dynamic regulation strategy. This strategy includes micro-perturbation injection commands for active detection and feedback analysis logic for inverting the type of hidden reactions based on the system response. The significance of this module lies in acting as the system's "intelligent decision-making center," realizing a paradigm shift from "passive response" to "active game theory" in regulation. It receives two dynamic indicators: entropy change (characterizing internal system uncertainty) and flux (characterizing external output risk). By querying the pre-set game matrix, it maps and generates a dynamic strategy with both "probing" and "treatment" purposes. Its core significance is that the system no longer waits for the risk to fully manifest but actively issues controllable "micro-perturbations" as probes and infers the type of hidden "dark reactions" based on the system's feedback patterns, thereby providing advanced and targeted decision-making commands for precise intervention.

[0027] The precision intervention module executes the received dynamic control strategies by injecting targeted co-metabolite matrices, catalytic materials, or electron donors and acceptors into the groundwater through the injection well network. This guides the pollutant transformation pathways towards safe water quality or complete mineralization. The significance of this module lies in its role as the system's "execution terminal," translating intelligent decisions into physicochemical actions that alter underground environmental processes. Following instructions from the game theory control module, targeted co-metabolite matrices, catalytic materials, or electron donors and acceptors are injected into specific locations through the injection well network. Its objective is not blind total reduction, but rather "guiding" or "correcting" identified specific reaction pathways. The aim is to reverse the pollutant transformation from the formation of highly toxic intermediates to safe pathways of low toxicity or complete mineralization, ultimately achieving root-cause, precise, and proactive remediation of groundwater quality risks.

[0028] In one specific implementation, the sensor network includes at least water quality sensors and water level sensors, with an automatic sampling device for collecting water samples at a preset cycle. First, a three-dimensional observation network consisting of multiple monitoring wells is deployed in the target aquifer. Each monitoring well is equipped with a water quality sensor for continuously measuring conventional physicochemical indicators, a water level sensor for recording dynamic changes in water level, and a temperature sensor for monitoring formation temperature. These sensors, at a preset acquisition frequency, transmit the acquired groundwater conductivity, dissolved oxygen, redox potential, water level depth, and temperature data to a data center via a telemetry terminal unit, forming a continuous time-series dataset. This provides fundamental hydrological dynamics and basic geochemical background information for subsequent analysis.

[0029] Automated sampling devices deployed at key locations automatically collect groundwater samples from designated depths at preset intervals, such as every 24 hours or after each precipitation event. The collected samples undergo immediate on-site filtration and stabilization pretreatment before being sent to a laboratory equipped with chromatography-mass spectrometry (GC-MS) and inductively coupled plasma mass spectrometry (ICP-MS). In the laboratory, technicians perform full-spectrum analysis of the water samples, quantitatively determining the precise concentrations of volatile organic compounds (VOCs) such as trichloroethylene and benzene compounds, heavy metal ions such as chromium, arsenic, and lead, and nutrients such as nitrates and nitrites, thereby generating a multi-pollutant concentration data report covering both organic and inorganic contaminants.

[0030] To obtain key hydrogeological parameters, the implementation team selected representative well locations within the monitoring well network to conduct multi-fall pumping tests. By recording the dynamic water level response data of the main well and observation wells during pumping, the hydraulic conductivity, storage coefficient, and effective porosity of the aquifer were calculated using existing technologies such as the Theis formula or numerical inversion methods. Simultaneously, aquifer core samples collected during drilling were sent to the laboratory for further determination of pore structure characteristics and permeability through mercury intrusion porosimetry and permeability experiments. These measured parameters collectively constitute a set of fundamental hydrogeological parameters characterizing groundwater flow and solute transport capabilities.

[0031] Geochemical indicators were acquired through a multi-pronged approach. In-situ water quality sensors deployed underground continuously monitored pH, redox potential, and dissolved oxygen concentration. Regularly collected water samples were analyzed in the laboratory using ion chromatography to determine the composition of major anions and cations, titration to determine total alkalinity, and spectroscopic methods to determine the forms and concentrations of variable-valence elements such as iron, manganese, and sulfur. These data were integrated and validated to form a comprehensive set of indicators reflecting the geochemical environment of groundwater, including redox zoning characteristics, acid-base buffering capacity, and potential mineral precipitation and dissolution tendencies, providing key environmental constraints for understanding pollutant transformation pathways.

[0032] In the dynamic assessment module: a pre-set coupled numerical model of groundwater flow and solute transport is constructed based on the finite difference method or the finite element method. The coupled numerical model uses hydrogeological parameters collected by the data monitoring module as the basic input parameters for model construction, calibration and boundary condition setting, and uses multi-pollutant concentration data as the input of solute concentration field to simulate the release intensity and initial spatial distribution of pollution sources. The simulation execution process of the coupled numerical model is as follows: first, the partial differential equations describing the groundwater movement law are solved to obtain the hydraulic head distribution and velocity field of the study area; then, the convection-dispersion-reaction equations describing the pollutant migration and transformation law are coupled and solved, and the complete dataset of the dynamic changes of pollutant concentration in the spatial and temporal dimensions is iteratively calculated and output. This dataset constitutes the intermediate result of concentration distribution for risk assessment.

[0033] In one specific implementation, a three-dimensional hydrogeological structure model of the study area is first established. Based on geological exploration borehole data, geophysical exploration results, and stratigraphic records from monitoring wells, the study area is vertically generalized into a multi-layered structure including unconfined aquifers, weakly permeable layers, and confined aquifers, and horizontally partitioned according to geomorphic units and sedimentary facies characteristics. Subsequently, the finite difference method is used as the numerical calculation kernel to discretize the three-dimensional structural model into a regular grid cell system. For example, rectangular grids are used on the horizontal plane, and the vertical plane is divided into several layers according to the stratigraphic thickness. Finally, a spatially discrete system containing hundreds of thousands to millions of active units is generated, providing a carrier for subsequent numerical calculations.

[0034] The discrete model was assigned and calibrated using field-measured hydrogeological parameters. The hydraulic conductivity, storage coefficient, and effective porosity obtained from pumping tests for different zones, as well as those from laboratory core analysis, were assigned to the corresponding model grid cells. Using long-term observed stable water level data from the monitoring well network as the initial flow field, and river water levels or known constant head boundaries as the model's boundary conditions, the groundwater flow model was run for steady-state calibration. By adjusting the hydraulic parameters of each zone, the root mean square error between the simulated and measured water level values ​​from each monitoring well was made to meet the preset accuracy requirements, thus obtaining a reliable velocity field model that truly reflects the regional groundwater flow characteristics.

[0035] A solute transport model fully coupled with the flow model was constructed, and the pollution source was initialized using multi-pollutant concentration data. This solute transport model shares the same set of spatially discrete grids. Multi-pollutant concentration data obtained through intensive monitoring at specific historical time points were used as the initial concentration field input. For example, the concentration detection results of trichloroethylene and cis-dichloroethylene in deep soil and groundwater samples obtained from the first comprehensive survey of the site left behind after the relocation of a chemical plant were used. Kriging spatial interpolation was then used to generate the initial concentration values ​​of each pollutant in each grid cell of the model, thereby numerically reconstructing the original spatial distribution of the pollution plume as the starting point for simulating its future transport evolution.

[0036] The flow model and solute transport model were integrated, and key reaction process parameters controlling pollutant migration and transformation were defined. Based on the calibrated velocity field, longitudinal and lateral dispersion values ​​were assigned to different lithological zones according to empirical dispersion formulas. For target pollutants, such as trichloroethylene, first-order decay reactions or more complex multi-step tandem biodegradation reaction chains were defined in their transport simulations. The relevant reaction rate constants were initially set based on previous microcosm experiments or site-specific geochemical conditions (such as monitored redox potentials and pH ranges). After completing all the above steps, a calibrated coupled numerical model with well-defined initial and boundary conditions and incorporating key physicochemical processes was constructed. This model is ready for execution and prediction.

[0037] It should be noted that in the groundwater-based water pollution risk assessment system involved in this invention, the selection of the finite difference method or the finite element method as the computational kernel of the numerical model is based on a balance between mature practices in the field of groundwater simulation and the actual application requirements of this invention. These two methods differ fundamentally in their mathematical principles and spatial discretization strategies, but both effectively serve the fundamental purpose of this invention: constructing a "dynamic assessment benchmark." The core of the finite difference method lies in using difference approximation to replace differentiation, directly solving for the head or concentration values ​​at the nodes of a regular grid (such as rectangles or hexahedrons). Its advantages lie in its intuitive formula derivation, relatively mature program implementation, and high computational efficiency. In the specific implementation of this invention, if the hydrogeological structure of the target site, after reasonable generalization, can be described by a relatively regular three-dimensional grid system (e.g., a clearly layered aquifer in a sedimentary plain), and the main focus is on the overall migration trend of the pollution plume at the regional scale, then using the finite difference method to construct the model is an efficient and reliable choice. Many industry standard software programs (such as MODFLOW) are based on this method, facilitating integration and verification.

[0038] The core of the finite element method (FEM) lies in dividing the solution domain into a series of irregular but closely connected sub-elements (such as triangles and tetrahedrons), and constructing an approximate function within each element for solution. Its greatest advantage lies in its unparalleled adaptability to complex geometries and boundaries. In the implementation of this invention, if the geological conditions of the target site are extremely complex, such as including steep faults, irregular lenses, meandering river boundaries, or when a fine characterization of the local flow field near the pollution source is required, the FEM becomes a better choice due to its mesh flexibility. It can more accurately fit the actual geological structure, thus providing a more accurate local velocity field under complex conditions, which is crucial for subsequent calculations of the "downstream toxicity equivalent flux".

[0039] This invention presents both methods as optional solutions because the core innovation lies not in the specific numerical solution algorithm itself, but in constructing a closed-loop system that starts with "benchmark model prediction" and subsequently uses "reaction network entropy change" and "toxicity equivalent flux" for dynamic diagnosis and control. Regardless of the numerical method used, as long as a calibrated benchmark model capable of simulating pollutant transport can be constructed based on monitoring data, all subsequent innovative steps can be smoothly initiated and supported. Secondly, this approach provides flexibility for implementation. Implementers can select the most suitable and efficient numerical engine as the kernel of the "pre-built model" based on the geological complexity of a specific site, data foundation, and technical preferences, without altering the overall system architecture and workflow. This ensures the method of this invention possesses strong adaptability and feasibility under different geographical environments and engineering conditions, thereby broadening its protection scope and application scenarios.

[0040] The dynamic risk assessment results generated by the dynamic assessment module are obtained by further processing of the intermediate results. The dynamic risk assessment results include two parts: the first part is the prediction of the spatial range boundary of the pollution plume within a future set time period extracted based on the intermediate results of concentration distribution; the second part is to convert the predicted exposure concentration of pollutants into a human health risk index, and classify the levels according to the legal or preset acceptable risk level thresholds, and finally generate a carcinogenic risk spatial zoning map and a non-carcinogenic risk spatial zoning map.

[0041] The generation of dynamic risk assessment results begins with the analysis of intermediate concentration distribution results output by the coupled numerical model. This intermediate result is a complete dataset containing predicted concentrations of specific pollutants in each three-dimensional grid cell within the study area at various future time points. A snapshot of the three-dimensional spatial distribution of concentrations at a future predetermined time point (e.g., the end of the first year after the simulation begins) is extracted from this dataset. Subsequently, the predicted concentration value of each grid cell is compared with the concentration limits for that type of pollutant specified in the National Groundwater Quality Standard. Through spatial analysis, all adjacent grid cells with concentrations exceeding the limits are aggregated and connected. Using an isosurface extraction algorithm, a well-defined three-dimensional spatial closed envelope is automatically generated. The projected outer boundary of the spatial volume defined by this envelope on the horizontal plane is output as a predicted boundary map of the spatial extent of the pollution plume within the future predetermined time period. This map is stored in a standard geographic information vector format and can be directly used for subsequent spatial planning and management.

[0042] Within the identified pollution plume's spatial range, a conversion calculation of the human health risk index is performed. For sensitive locations such as residential areas potentially affected by the pollution plume, the pollutant concentration variation curve over time for the grid cell containing that location is extracted from intermediate results of the same concentration distribution throughout the simulation period. Based on the model and exposure parameters recommended in the national "Technical Guidelines for Soil Pollution Risk Assessment of Construction Land," and considering the actual use of groundwater in the area (e.g., as a drinking water source), the lifetime average exposure dose for residents via drinking water is calculated. Then, this exposure dose is coupled with the pollutant's toxicity parameters: for pollutants with carcinogenic effects, the exposure dose is multiplied by the substance's carcinogenic slope factor to obtain the lifetime carcinogenic risk index for that location; for pollutants with non-carcinogenic effects, such as nitrates, the exposure dose can be divided by the substance's reference dose. For better quantification, the lifetime carcinogenic risk index can be used as a reference, defined as the hazard quotient index for that location.

[0043] Human health risk indices include two categories: hazard quotient (HQ) and lifetime carcinogenic risk index. The calculated risk indices are classified according to legally acceptable risk levels. For carcinogenic risk indices, levels below 1 / 1,000,000 are typically classified as acceptable risk, between 1 / 1,000,000 and 1 / 10,000 as low risk, between 1 / 10,000 and 1 / 1,000 as medium risk, and above 1 / 1,000 as high risk. For non-carcinogenic hazard quotients, levels less than 1 are typically classified as risk-free, and levels greater than 1 are classified as risky. Using a geographic information system (GIS), the risk indices of all calculated grid cells within the study area are classified according to the above thresholds, and different color codes are assigned to different levels.

[0044] Based on the above classification results, formal risk spatial zoning maps were generated, with the carcinogenic risk level zoning map and the non-carcinogenic risk level zoning map each created as independent thematic maps. Each map includes a clear risk level legend, geographic coordinate grid, scale bar, and labels of important features. These two zoning maps, together with the aforementioned pollution plume spatial extent boundary prediction map, constitute a complete dynamic risk assessment result. They intuitively demonstrate the spatial distribution pattern and severity of potential health risks caused by groundwater pollution at a specific future time point, providing a quantitative and spatial scientific basis for risk management decisions.

[0045] To obtain the hydraulic head distribution and velocity field of the study area, a system of partial differential equations describing groundwater movement is solved through the following steps: Based on hydrogeological parameters, the underground space of the study area is divided into three-dimensional grid cells; for each grid cell, Darcy's law and the law of conservation of mass are applied to establish a mathematical relationship describing the flow balance of the groundwater unit; the mathematical relationships established for all grid cells in the study area are combined and substituted with boundary conditions determined by water level data collected by the data monitoring module to form a closed linear equation system; an iterative solution algorithm is used to numerically calculate the closed linear equation system to obtain the hydraulic head value at the central node of each grid cell, and the set of all node hydraulic head values ​​constitutes the hydraulic head distribution; based on the hydraulic head distribution, the hydraulic gradient and the hydraulic conductivity coefficient in the hydrogeological parameters, and according to the calculation formula of Darcy's law, the groundwater seepage velocity in the three-dimensional direction of each grid cell is calculated and output, and the set of seepage velocities of all grid cells constitutes the velocity field.

[0046] Specifically, the process begins with three-dimensional spatial discretization of the target area based on previously acquired hydrogeological parameters to construct a numerical computation grid. For example, for a site affected by chemical waste, the geological survey report generalizes the site from top to bottom into fill layers, silty unconfined aquifers, clay impermeable layers, and gravel confined aquifers. Based on the thickness and distribution of each layer, a 10m x 10m rectangular grid is used for horizontal subdivision, and the unconfined aquifers are subdivided into three layers and the confined aquifers into two layers vertically. This results in a three-dimensional structured grid system containing tens of thousands of active cells, each assigned a typical hydraulic conductivity coefficient and effective porosity value for that zone obtained from field tests.

[0047] The fundamental physical laws of groundwater movement are applied to each three-dimensional grid cell to establish governing equations. For each cell, according to Darcy's law, the groundwater flow through each face is proportional to the hydraulic conductivity of that face, the flow area, and the head difference between adjacent cells. Simultaneously, according to the law of conservation of mass, under steady-state flow conditions, the total flow into the cell must be equal to the total flow out of the cell. Combining these two laws, a linear algebraic equation can be written for each grid cell, with the head of adjacent cells as unknowns. For example, for an internal cubic cell, the equation describes a balance relationship where the sum of water flows from the six adjacent cells (top, bottom, left, right, front, and back) is zero.

[0048] The linear algebraic equations established by all grid cells within the study area are integrated and combined with boundary conditions to form a closed large-scale system of linear equations. The boundary conditions are determined by monitored well water level data; for example, the monitored water level of a river in the northern part of the site is set as a constant head boundary, and the groundwater watershed in the southern part of the site is set as a zero flux boundary. Using these boundary values, the corresponding terms in the system of equations are assigned values, thus ensuring that the system of equations has a mathematically unique solution. Formally, this system of equations can be expressed as the product of the coefficient matrix and the unknown head vector equals the constant vector on the right-hand side.

[0049] An iterative algorithm is used to solve the aforementioned large-scale linear equation system to obtain the overall head distribution, and then the three-dimensional velocity field is calculated. Specifically, efficient iterative algorithms such as the preconditional conjugate gradient method can be used to solve the problem on a computer until the change in the calculated head value of all grid nodes is less than one millionth of a meter between two consecutive iterations. After solving, each grid node obtains an accurate head value, and the head values ​​of all nodes constitute the spatial head distribution. Subsequently, according to Darcy's law, the hydraulic gradient is calculated based on the head difference between adjacent nodes, and then multiplied by the hydraulic conductivity coefficient of the medium contained in that grid cell to obtain the seepage velocity components of that cell in the X, Y, and Z directions. The set of velocity vectors of all cells in the entire study area ultimately forms the three-dimensional groundwater velocity field used to simulate the transport of pollutants.

[0050] The iterative calculation and output of a complete dataset showing the dynamic changes of pollutant concentrations in both spatial and temporal dimensions is achieved through the following steps: The velocity field is used as the input for the convection term driving pollutant transport; based on the dispersion and velocity field in hydrogeological parameters, the hydrodynamic dispersion term of the pollutants is calculated; a reaction term describing the biodegradation, chemical oxidation-reduction, and adsorption-desorption processes occurring during pollutant migration is constructed; at each 3D grid cell and each simulation time step, the convection, dispersion, and reaction terms are coupled and solved simultaneously to update the pollutant concentration of that grid cell at the end of that time step; the updated pollutant concentration is used as the initial value for the next time step, and the above coupled solution process is repeated to achieve iterative calculation for all simulation time steps; finally, the set of pollutant concentrations for all grid cells at all simulation time steps is output, which constitutes a complete dataset showing the dynamic changes of pollutant concentrations in both spatial and temporal dimensions.

[0051] In one specific implementation, the iterative calculation process for solute transport begins with acquiring and processing three-dimensional steady-state velocity field data generated from groundwater flow simulations. This velocity field provides the seepage velocity components in the X, Y, and Z directions for each cell in a gridded form. Taking trichloroethylene as an example, when simulating its transport from the pollution source to the surrounding area, the convection term needs to be calculated first. At each grid cell, the convection term quantifies the net flux of pollutants carried by groundwater flow, and its value is obtained by calculating the difference in pollutant mass flowing into and out of each face of that cell. Specifically, the pollutant mass flowing into a face is equal to the pollutant concentration at the center point of the upstream adjacent cell multiplied by the groundwater volumetric flow rate through that face, and the volumetric flow rate is determined by the surface area, normal direction, and seepage velocity vector at that point. Numerically, an upstream weighted method is typically used to determine the interface concentration to ensure the stability of the calculation.

[0052] Based on the hydrogeological dispersion parameters of the aquifer and the local velocity field, the hydrodynamic dispersion term is calculated. Dispersion is contributed by both mechanical dispersion and molecular diffusion, and its magnitude is characterized by the dispersion tensor. The dispersion tensor depends on the pore-average velocity and the longitudinal and transverse dispersion. In each grid cell, the dispersion term describes the diffusion and dispersion flux of pollutants due to the concentration gradient. Its numerical calculation involves a discrete approximation of the spatial second derivative of the concentration at the cell center. For example, in a regular finite difference grid, the dispersion flux component along the X direction can be obtained by calculating the concentration difference between adjacent cells (east and west cells) and multiplying it by the dispersion coefficient derived from the local velocity and dispersion.

[0053] A reaction term describing the transformation of pollutants during migration is constructed. For the biodegradation of trichloroethylene under anaerobic conditions, the reaction term can be expressed as a first-order kinetic decay. Its numerical calculation is the trichloroethylene concentration within the grid cell multiplied by a first-order decay rate constant pre-determined through site microcosm experiments or literature data. For example, if the trichloroethylene concentration in a cell is 2 mg / L and the set first-order decay rate constant is 0.1 mg / year, then the mass concentration of trichloroethylene removed by the reaction in that cell within one year can be obtained by multiplying the concentration by the rate constant, resulting in 0.2 mg / L per year. For more complex reaction networks, such as those involving multi-step degradation to generate intermediate products, the reaction terms will form a simultaneous system of equations describing the ebb and flow of concentrations of various species.

[0054] The convection, dispersion, and reaction terms are explicitly and implicitly coupled in time and solved simultaneously in space to update the concentration field at the end of each time step. Within each simulation time step (e.g., one day), the mass conservation relationship of pollutants within each grid cell is expressed as a difference equation: the increase in pollutant mass within the cell equals the net flux input by convection and dispersion, minus the amount consumed by reaction. The difference equations of all grid cells are combined to form a large system of linear equations with the new time-instance concentrations of all cells as unknowns. An efficient iterative solver (such as the generalized minimum residual method) is used to solve this system of equations, thereby obtaining the updated concentration distribution of the entire computational domain at the end of the new time step in one go. The updated concentration field is used as the initial condition for the next time step, and the above process of calculating and coupling the convection, dispersion, and reaction terms is repeated until the simulation of the entire set duration is completed. Finally, the concentration values ​​of all grid cells at all time points are output, forming a complete spatiotemporal dataset for risk assessment.

[0055] In one specific implementation, the calculation of the entropy change value of the reaction network first involves the systematic organization and preprocessing of the monitoring data. Groundwater samples continuously collected by the monitoring network are analyzed in the laboratory to generate reports containing the concentrations of various organic and inorganic pollutants. These data are arranged chronologically to construct a concentration matrix with the monitoring period as rows and each specific detected chemical substance as columns. For example, in twelve consecutive months of monitoring of a chemically contaminated site, this matrix may contain twelve time points (rows) and the average concentrations (columns) of dozens of chemical substances such as trichloroethylene, cis-dichloroethylene, vinyl chloride, benzene, and hexavalent chromium in different monitoring wells. Missing values ​​are appropriately filled using interpolation methods to ensure data continuity.

[0056] The first entropy component characterizing the complexity of species diversity is calculated. For each monitoring period (e.g., each month), the concentration data of all detected chemical substances within that period are statistically analyzed. The concentration of each species is divided by the total concentration of all species in that period to obtain its relative abundance proportion, treating the relative abundance of all species as a complete probability distribution. Subsequently, the Shannon information entropy formula is applied, which is the sum of probabilities multiplied by their logarithms and then negative. For example, if only three pollutants are detected in a certain month with relative abundances of 0.5, 0.3, and 0.2, the calculated first entropy component is approximately 1.03. The higher this value, the richer and more evenly distributed the chemical composition of groundwater during that period.

[0057] Parallel computation is used to compute the second entropy component, representing temporal disorder, and the third entropy component, representing network structural uncertainty. For the second component, for each chemical species (e.g., vinyl chloride), its concentration sequence over the entire monitoring period (e.g., twelve months) is extracted. The permutation entropy of this sequence is calculated: defining the embedding dimension (e.g., 3) and time delay (e.g., 1), the sequence is divided into multiple overlapping subsequences, each containing three consecutive concentration values; the three values ​​in each subsequence are sorted according to size to obtain an ordinal pattern (e.g., "small-medium-large" corresponds to pattern 1); the frequency of all possible ordinal patterns in the entire sequence is counted to form a pattern probability distribution, and the Shannon entropy of this distribution is calculated, which is the permutation entropy of the concentration sequence of that species. The arithmetic mean of the permutation entropies of all species is calculated to obtain the second entropy component. For the third component, based on a pre-defined knowledge graph, such as the known that trichloroethylene can be progressively degraded into cis-dichloroethylene and vinyl chloride under reducing conditions, a potential reaction network graph among the detected species in the current monitoring period is constructed. In the graph, nodes represent species, and edges represent possible transformation relationships. The edge weights are estimated by calculating the concentration covariance of the two connected species (e.g., trichloroethylene and cis-dichloroethylene) over the entire monitoring time series. Specifically, the calculation formula can be the absolute value of the Pearson correlation coefficient between their concentrations. Simultaneously, the degree (number of connected edges) of each node is calculated. Then, based on the weight distribution of all edges and the degree distribution of all nodes, their Shannon entropy is calculated. These two entropy values ​​are then added together with preset weights to obtain a third entropy component characterizing the uncertainty of the network structure.

[0058] The three entropy components are nonlinearly fused to generate the final entropy change value of the reaction network. The calculated first, second, and third entropy components are input into a predefined synthesis function, which can be expressed as a well-defined nonlinear polynomial. This function is designed to capture the interactions between the three components; for example, when the species diversity entropy (first component) is low but the process dynamics entropy (second component) is high, it may indicate that a few species are undergoing drastic and unpredictable transformations. In this case, the function will output a higher overall entropy value to warn of system instability. The output of the synthesis function is an unnormalized scalar, which is normalized by dividing it by the maximum reference entropy value calculated on historical "quiet period" data, so that the final result falls between 0 and 1. This normalized scalar is the entropy change value of the reaction network, and its tendency to 1 indicates that the subsurface pollutant reaction network is in a highly complex, dynamic, and unpredictable state.

[0059] In one specific implementation, the pre-defined nonlinear fusion function used is essentially a carefully designed multivariate polynomial. The independent variables of this polynomial are the calculated first entropy component (denoted as H1, representing species diversity), the second entropy component (denoted as H2, representing process dynamics), and the third entropy component (denoted as H3, representing network structure). Its core design lies in quantitatively characterizing the nonlinear superposition, cancellation, or amplification effects that may exist when the three entropy values ​​reflect the overall chaos level of the system by introducing higher-order terms and cross-product terms of the independent variables.

[0060] The first optional polynomial form focuses on expressing the individual significance of entropy components and their fundamental synergistic effects. Its general expression can be described as follows: First, assign different basic weight coefficients to the three entropy components and perform a linear weighted sum; then, add the squared term of each entropy component to capture the accelerating or decelerating trend of its own influence; finally, add pairwise product terms for every two entropy components to quantify the synergistic effect between any two. For example, when both species diversity entropy and process dynamics entropy are high, their product term will produce a large positive value, significantly boosting the overall output, which describes the extremely high system uncertainty represented by "numerous species with drastic individual changes." The coefficients of all terms need to be calibrated through historical data regression analysis or based on expert experience with physical meaning.

[0061] The second polynomial form aims to characterize specific antagonistic effects and higher-order interactions. Building upon linear, quadratic, and pairwise product terms, it further introduces cubic terms with three entropy components to describe the "S"-shaped inflection point characteristic of its influence. Simultaneously, a ternary product term with three entropy components is added. This ternary product term is crucial, simulating the "emergent" or saturation inhibition effects that may arise when three sources of uncertainty act together, exceeding the pairwise superposition. For example, when network structure entropy is extremely high, it may mean that the response path is extremely ambiguous. Even if species diversity entropy and process dynamics entropy are high, the overall risk may enter a state of "chaotic saturation" due to the system's excessive complexity. The ternary product term can be designed with negative coefficients to moderately lower the final output and avoid over-warning.

[0062] The third form is a piecewise polynomial whose coefficients can be dynamically adjusted according to the different intervals in which the entropy components are located. For example, multiple thresholds can be pre-set based on historical data or theory to divide the value range of each entropy component into three intervals: "low," "medium," and "high." Different sets of polynomial coefficients are used for different "state domains" formed by different combinations of the three entropy components. When monitoring data indicates that the system is in the state domain of "H1 high, H2 medium, H3 low," the system calls the coefficient set specifically calibrated for this state domain for calculation. This coefficient set may significantly amplify the weights of the H1 square term and the H1-H2 product term because this state indicates an unstable situation with complex species composition but relatively clear transformation pathways. This form, through domain-specific fitting, can more flexibly and precisely approximate complex nonlinear relationships.

[0063] Regardless of the specific form used, the final output of this nonlinear fusion function is an unnormalized real value. To obtain the final reaction network entropy change value in the [0, 1] interval, this real value can be input into a subsequent normalization function. This normalization function is typically a sigmoid function, which maps the input real number to the (0, 1) interval. By adjusting the center point and slope parameters of the sigmoid function, it can be ensured that the output value is close to 0 in historical "stable baseline" states and close to 1 in historical "extreme chaotic" states, thus ensuring consistent interpretability and comparability of the final reaction network entropy change value.

[0064] In one specific implementation, the calculation of downstream toxicity equivalent flux begins with the synchronous acquisition of multi-source data from key hydrological sections. Online biotoxicity sensors, such as those based on the principle of luminescent bacteria inhibition, are deployed in a cluster of monitoring wells downstream of the pollution source and perpendicular to the groundwater flow direction to measure the raw data of the overall inhibition rate of the water body on organisms. Simultaneously, from the current-moment three-dimensional groundwater velocity field model generated during the dynamic assessment process, the seepage velocity vector at each point in the spatial location of the monitoring section is extracted and interpolated. This vector contains both magnitude and direction information, collectively forming the section velocity field dataset. For example, in a section downstream of a site contaminated with petroleum hydrocarbons, the sensors report a luminescence inhibition rate of 35%, while the numerical model provides an average pore water velocity of 0.5 meters per day in the central region of the section, pointing due east.

[0065] Raw biotoxicity data are converted into toxicity equivalent concentrations with uniform dimensions. Based on a pre-defined dose-response standard curve, the overall inhibition rate output by the sensor is converted into an equivalent concentration of a specific reference toxicant. This standard curve is pre-established by testing a standard solution of the reference toxicant (such as benzo[a]pyrene) in the laboratory using the same type of sensor. For example, according to the calibration curve, a 35% luminescence inhibition rate may correspond to an equivalent concentration of 0.8 micrograms per liter of benzo[a]pyrene. This step normalizes all overall toxicity effects to equivalent concentrations in micrograms per liter of benzo[a]pyrene toxicity equivalents, thus enabling quantitative comparisons of overall toxicity at different times and with different pollutant compositions.

[0066] The monitoring section was spatially discretized, and unit flux calculations were performed. Vertically, the entire monitoring section was uniformly discretized into ten horizontal calculation strips, each one meter thick, based on the aquifer thickness and the arrangement of the monitoring well filters. For each calculation strip, the seepage velocity vector at the center point of the strip was first obtained from the corresponding velocity field data. The component of this velocity vector perpendicular to the normal direction of the monitoring section was calculated, representing the effective transport velocity perpendicular to the section. Next, this effective transport velocity was multiplied by the benzo[a]pyrene equivalent concentration value representing the overall toxicity of the strip obtained in the previous step, and then multiplied by the representative area of ​​the calculation strip (the product of the strip thickness and the section width) to obtain the toxicity equivalent flux per unit time through the calculation strip, with dimensions in micrograms per day (TCDs).

[0067] The fluxes of all discrete computational units are integrated and summed to generate a total cross-sectional flux index. The unit-time toxicity equivalent fluxes calculated from each of the ten computational strips are then summed to obtain the total downstream toxicity equivalent flux passing through the entire monitoring section. Using example data, if the sum of the fluxes from all strips is 12,000 toxicity equivalent micrograms per day, this value is output as a core risk indicator. This indicator dynamically and comprehensively reflects the "toxicity mass" transported downstream by the pollution plume per unit time. A sharp increase or sustained high value directly indicates an increasing exposure risk to sensitive downstream targets, providing the most direct and quantitative decision-making basis for subsequently determining whether and when to initiate proactive intervention measures.

[0068] In one specific implementation, the generation of dynamic control strategies begins with the initial construction of the game matrix. This process, based on a baseline risk zoning map generated through prior simulation, associates different risk level regions defined in the map with the row and column strategies of the matrix. Row strategies are defined as a series of potential dark reaction type hypotheses based on geochemical principles, such as "anaerobic reduction dehalogenation," "sulfate reduction coupled with oxidation," or "iron-manganese oxide-mediated oxidation." Column strategies correspond to a series of available micro-perturbation interventions, such as "injecting sodium lactate as a co-metabolite substrate," "injecting slow-release calcium peroxide as an oxidant," or "injecting zero-valent iron as a reducing agent and electron donor." Each cell in the matrix is ​​assigned an initial comprehensive benefit score, which, based on expert knowledge, estimates the combined expected positive benefits (such as promoting complete mineralization) and risk control (such as inhibiting the formation of highly toxic intermediates) that can be obtained when using the corresponding column action to intervene in the corresponding row reaction type in a specific risk region (such as a high-risk region).

[0069] Two core dynamic indicators obtained from monitoring and analysis—the reaction network entropy change and the downstream toxicity equivalent flux—are jointly processed to pinpoint the decision point. These two indicators are considered together as a two-dimensional state index. For example, when the reaction network entropy change is 0.7 (indicating high system complexity) and the downstream toxicity equivalent flux is 15,000 micrograms of benzo[a]pyrene equivalent per day (indicating high toxicity transport intensity), based on a preset threshold range, this state index is located in a specific sub-matrix block within the game matrix marked as "high uncertainty - high output risk." This sub-matrix block contains pre-set row and column strategy combinations and their scores specifically designed to address this severe situation.

[0070] Within the identified sub-matrix block, strategy optimization is performed. The system traverses all cells within the block and selects the one with the highest overall benefit score. Assuming the row strategy for this cell is "The anaerobic reduction dehalogenation process has stalled, leading to the accumulation of the intermediate vinyl chloride," and the column strategy is "Inject a mixed solution of sodium lactate and vitamin B12 in a specific ratio to stimulate the activity of dehalogenation microorganisms and provide necessary cofactors," the system parses this column strategy into a set of precise, executable micro-perturbation injection instructions. These instructions specify the type of injected reagent (sodium lactate and vitamin B12 solution), concentration (e.g., 10 grams of sodium lactate per liter, 1 milligram of vitamin B12 per liter), total dosage (e.g., 10 cubic meters), and the three-dimensional coordinates of the injection point determined by a hydrogeological model. Simultaneously, the corresponding row strategy—"Verify whether a stalled reduction dehalogenation process exists"—is designated as the core scientific question requiring focused attention and inversion verification in subsequent monitoring and analysis.

[0071] The generated micro-perturbation injection instructions are encapsulated into a standard operational instruction set, covering reagent preparation, pumping parameters, injection duration, and subsequent monitoring protocols. This instruction set is directly sent to the field execution unit. Meanwhile, the hypotheses regarding the type of dark reaction requiring verification are simultaneously pushed into the data analysis workflow to guide subsequent targeted analysis of monitoring data (especially post-perturbation species concentration change patterns), thus completing a full decision-making cycle from "state-based diagnostic strategy generation" to "guided intervention and hypothesis testing."

[0072] In one specific implementation, the process of guiding system initialization based on the baseline prediction map begins with the analysis and parameter extraction of the digitized map. The baseline plume prediction map generated by the first round of numerical simulation is stored in a geographic information system vector format, which includes the spatial boundary polygons of the plume at the set prediction period (e.g., one year, five years, etc.) and the groundwater flow direction arrows calculated through concentration gradients. Key parameters are extracted from this map, including the maximum possible migration distance of the plume along the main groundwater flow direction and the azimuth of the plume's central axis. For example, when dealing with a chlorinated hydrocarbon contaminated site, the maximum migration distance of the plume is read from the map as 320 meters downstream, with the main flow direction being 15 degrees east of south.

[0073] Based on the extracted parameters, the location of key monitoring sections is automatically planned on the digitized site map. Using the maximum migration distance obtained in the previous step as a base, a preset safety factor (e.g., 0.6) is multiplied to obtain the recommended distance for section placement. The calculation formula is: section distance equals the distance between the pollution source and the endpoint of the maximum migration distance multiplied by the safety factor. Continuing the previous example, the calculated section should be placed approximately 192 meters downstream of the pollution source. At this location, the software automatically generates a virtual line segment perpendicular to the flow direction of 15 degrees east of south as the monitoring section. The length of this line segment is determined based on the maximum width of the pollution plume in the vertical flow direction plus a certain margin (e.g., 20 meters). The spatial coordinates and geometric parameters of this virtual section (including the coordinates of the start and end points and the depth range) are then output as a technical document, directly guiding the subsequent installation location of online biotoxicity sensors in monitoring wells and the screening range of the monitoring wells.

[0074] The software analyzes the risk level zones in the baseline risk zoning map, identifying red areas as high-risk zones and yellow areas as medium-risk zones. For each risk level zone, it invokes pre-defined knowledge base rules. For example, for a zone identified as medium-risk with trichloroethylene as the primary historical pollutant, the pre-defined rules retrieve a set of relevant potential dark reaction type hypotheses from the knowledge base, such as "incomplete reduction dechlorination leads to vinyl chloride accumulation" or "covalent binding with natural organic matter under oxidative conditions," and assign these hypotheses as row strategies in the game matrix for that region. Simultaneously, based on the region's geochemical background data (e.g., groundwater under neutral or reducing conditions), it matches recommended basic actions from the intervention measures library, such as "injecting slow-release sodium lactate to promote reduction" or "injecting zero-valent iron to enhance reduction and adsorption," and assigns these actions as column strategies in the game matrix for that region.

[0075] Finally, after completing the initial value filling and overall initialization settings of the matrix, the system assigns an initial comprehensive benefit score to each pair of row and column strategies (i.e., a cell in the matrix). This score is derived from historical data or expert experience stored in the knowledge base regarding the expected effects of intervention measures on this type of pollutant under similar geochemical environments. Once all cells are filled with initial values, an initial game matrix that matches the specific risk pattern and pollution characteristics of the site and is ready for immediate use is formed. Simultaneously, the deployment parameters for monitoring sections have also been issued. The completion of these two tasks signifies that the entire dynamic control system possesses clear physical monitoring targets and a logical decision-making framework, laying a precise spatial and rule-based foundation for subsequent data-driven "perception-decision-intervention" closed-loop operation.

[0076] The above algorithms or formulas are all dimensionless and numerical calculations, and the results are obtained by software simulation based on a large amount of collected data to obtain the most recent real-world results. The preset parameters are set by those skilled in the art according to the actual situation.

[0077] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0078] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0079] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the devices and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0080] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A groundwater-based water pollution risk assessment system, characterized in that, Includes the following steps: The data monitoring module is used to collect data on the concentration of multiple pollutants, hydrogeological parameters, and geochemical indicators in groundwater through a sensor network and automatic sampling devices deployed in the aquifer. The dynamic assessment module is used to simulate the release probability, transport path, spatial distribution and concentration evolution of pollutants based on the collected multidimensional data and through a pre-set numerical model coupled with groundwater flow and solute transport, and output dynamic risk assessment results including pollution plume range prediction and risk level classification. The entropy change calculation module is used to calculate and output the entropy change value of the reaction network based on the collected multidimensional chemical fingerprint data, by analyzing the diversity of chemical species, the temporal changes in concentration, and the uncertainty of reaction pathways. This is used to quantify the dynamic complexity of the underground pollutant reaction network. The flux calculation module is used to acquire the comprehensive toxicity data measured by online biotoxicity sensors deployed at key hydrological sections downstream of the pollution plume, and combine it with a pre-set groundwater seepage velocity model to calculate and output the downstream toxicity equivalent flux. This flux characterizes the comprehensive toxicity transport intensity through the monitoring section per unit time. The game-theoretic control module receives data on the entropy change of the reaction network and the downstream toxicity equivalent flux. By querying a game matrix pre-constructed based on different geochemical backgrounds and pollutant combinations, it maps and generates a dynamic control strategy. This strategy includes micro-perturbation injection commands for active detection and feedback analysis logic for inverting the dark reaction type based on the system response. The precision intervention module is used to execute the received dynamic control strategy. By injecting targeted co-metabolite matrix, catalytic material or electron donor and acceptor into the groundwater through the injection well network, it guides the pollutant transformation pathway to the direction of water safety or complete mineralization.

2. The groundwater-based water pollution risk assessment system according to claim 1, characterized in that, The sensor network includes at least a water quality sensor and a water level sensor, and an automatic sampling device is used to collect water samples according to a preset cycle.

3. The groundwater-based water pollution risk assessment system according to claim 2, characterized in that, In the dynamic evaluation module: The pre-set numerical model of groundwater flow and solute transport is constructed with a computational kernel based on the finite difference method or the finite element method. The coupled numerical model uses hydrogeological parameters collected by the data monitoring module as the basic input parameters for model construction, calibration and boundary condition setting, and multi-pollutant concentration data as the input of solute concentration field to simulate the release intensity and initial spatial distribution of pollution sources. The simulation execution process of the coupled numerical model is as follows: First, the partial differential equations describing the movement of groundwater are solved to obtain the head distribution and velocity field of the study area; then, the convection-diffusion-reaction equations describing the migration and transformation of pollutants are solved in a coupled manner, and the complete dataset of the dynamic changes of pollutant concentration in the spatial and temporal dimensions is iteratively calculated and output. This dataset constitutes the intermediate results of concentration distribution for risk assessment.

4. The groundwater-based water pollution risk assessment system according to claim 3, characterized in that, The dynamic risk assessment results generated by the dynamic assessment module are obtained through further processing of intermediate results. The dynamic risk assessment results consist of two parts: the first part is the prediction of the spatial extent boundary of the pollution plume within a future set time period, extracted based on intermediate results of concentration distribution; The second part involves converting the predicted exposure concentration of pollutants into a human health risk index, and classifying them according to legally mandated or pre-set risk acceptable level thresholds, ultimately generating a spatial zoning map of carcinogenic risk and a spatial zoning map of non-carcinogenic risk.

5. The groundwater-based water pollution risk assessment system according to claim 4, characterized in that, Solving the partial differential equations describing groundwater movement to obtain the hydraulic head distribution and velocity field of the study area is achieved through the following steps: Based on hydrogeological parameters, the underground space of the study area is divided into three-dimensional grid units; On each grid cell, Darcy's law and the law of conservation of mass are applied to establish a mathematical relationship describing the flow balance of the groundwater cell; The mathematical relationships established by all grid cells in the study area are combined and substituted with the boundary conditions determined by the water level data collected by the data monitoring module to form a closed linear equation system. An iterative solution algorithm is used to numerically calculate the closed linear equation system, and the hydraulic head value at the center node of each grid cell is obtained. The set of all node hydraulic head values ​​constitutes the hydraulic head distribution. Based on the head distribution, the hydraulic conductivity coefficient in the hydrogeological parameters, and Darcy's law, the groundwater seepage velocity in the three-dimensional direction of each grid cell is calculated and output. The collection of seepage velocities of all grid cells constitutes the velocity field.

6. The groundwater-based water pollution risk assessment system according to claim 5, characterized in that, The iterative calculation and output of a complete dataset showing the dynamic changes of pollutant concentrations in both spatial and temporal dimensions is achieved through the following steps: Using the velocity field as the input of the convection term driving pollutant transport, the hydrodynamic dispersion term of pollutants is calculated based on the dispersion and velocity field in hydrogeological parameters. Construct reaction terms that describe the biodegradation, chemical redox, and adsorption-desorption processes that occur during the migration of pollutants; At each 3D grid cell and at each simulation time step, the convection, dispersion, and reaction terms are coupled and solved simultaneously to update the pollutant concentration of the grid cell at the end of that time step. Using the updated pollutant concentration as the initial value for the next time step, the above coupled solution process is repeated to achieve iterative calculation for all simulation time steps; finally, the pollutant concentration set of all grid cells at all simulation time steps is output, which constitutes a complete dataset of pollutant concentration dynamically changing in spatial and temporal dimensions.

7. The groundwater-based water pollution risk assessment system according to claim 6, characterized in that, The specific process of the entropy change in the reaction network is as follows: The first step is data preparation and preprocessing: receiving multidimensional chemical fingerprint data continuously collected by the data monitoring module and organizing it into a concentration matrix with time as the sequence and chemical species as the dimension; The second step is to calculate the entropy of species diversity components: For each monitoring period, the relative concentration ratio of all detected chemical species is statistically analyzed and regarded as a probability distribution. The Shannon information entropy of this distribution is calculated to obtain the first entropy value component that characterizes species richness and evenness. The third step is to calculate the dynamic entropy of the concentration process: For each chemical species, extract its concentration value over the entire monitoring time series and calculate the permutation entropy of the concentration time series; the permutation entropy quantifies the disorder and unpredictability of concentration changes by comparing the pattern complexity formed by the relative size relationship of adjacent data points in the sequence, and integrates the permutation entropy of all species to obtain the second entropy value component. The fourth step is to calculate the entropy of the reaction path inference: Based on the pre-set knowledge graph of pollutant transformation paths, the reaction network graph of all chemical species detected in the current period is reconstructed; by analyzing the degree distribution and edge weight of each node in the graph, the uncertainty entropy of the reconstructed network graph structure is calculated to obtain the third entropy value component. Step 5: Synthesis and output of entropy change value: Input the first, second and third entropy value components into a preset nonlinear fusion function. This function simulates the synergistic and antagonistic effects of the three entropy values ​​in reflecting the degree of chaos in the system, and outputs a normalized scalar value, which is the entropy change value of the network.

8. The groundwater-based water pollution risk assessment system according to claim 6, characterized in that, The specific process of downstream toxicity equivalent flux is as follows: Simultaneously acquire comprehensive toxicity data measured by one or more online biotoxicity sensors deployed at key hydrological sections, as well as current groundwater seepage velocity vector field data corresponding to the location of the monitoring section provided by the dynamic assessment module; The raw comprehensive toxicity data output by the online biotoxicity sensor is converted into equivalent concentration values ​​expressed in uniform toxicity equivalent units based on a preset dose-response relationship standard curve. The key hydrological monitoring section is discretized into multiple calculation units in the vertical direction. For each calculation unit, the groundwater seepage velocity vector at its location is obtained, and the component of the velocity vector in the direction perpendicular to the monitoring section is calculated. The vertical velocity component is multiplied by the toxicity equivalent concentration value of the corresponding unit, and then multiplied by the representative area of ​​the calculation unit to obtain the toxicity equivalent flux through the unit. The toxicity equivalent flux of all computational units is summed to obtain the total downstream toxicity equivalent flux passing through the entire monitoring section, which is then output as an indicator characterizing the comprehensive risk intensity of the pollution plume transport downstream.

9. The groundwater-based water pollution risk assessment system according to claim 8, characterized in that, The specific process of generating dynamic control strategies by querying the game matrix mapping in the game control module is as follows: Based on the baseline risk partition map generated by the first round of simulation in the dynamic evaluation module, initial row and column strategies are set for the game matrix. The row strategy corresponds to different potential hidden reaction type assumptions, and the column strategy corresponds to different micro-perturbation intervention actions. The value of each cell in the matrix represents the initial value of the expected comprehensive benefit score when the corresponding column action is used to intervene in the corresponding row reaction type under a specific risk level area. The entropy change value of the reaction network output by the entropy change calculation module and the downstream toxicity equivalent flux output by the flux calculation module are used as a set of joint state indices; based on the numerical range of the joint state index, a corresponding sub-matrix block in the game matrix is ​​located. Within the located sub-matrix block, find the cell with the highest current comprehensive benefit score; parse the column strategy corresponding to the cell, i.e. the specific micro-perturbation intervention action, into an executable micro-perturbation injection instruction, which at least includes the type, concentration, dosage, and injection point coordinates of the injected reagent; at the same time, use the row strategy corresponding to the cell, i.e. the assumed dark reaction type, as the core target for verification and inversion in the feedback analysis logic.

10. The groundwater-based water pollution risk assessment system according to claim 9, characterized in that, The baseline plume prediction map and baseline risk zoning map generated by the dynamic assessment module are used to directly guide the system initialization settings, specifically including: To guide the deployment of key hydrological monitoring sections in the flux calculation module: Based on the maximum migration distance and mainstream direction of the pollution plume indicated in the baseline pollution plume prediction map, the system plans and generates at least one virtual monitoring section perpendicular to the groundwater flow direction at a preset safe distance downstream of the pollution source; the spatial coordinates and geometric range of this section are set as the actual deployment location and monitoring range of the online biotoxicity sensor. The system is used to provide initial strategy mapping relationships for the preset game matrix in the game control module: the system analyzes the risk level of different geographical areas in the benchmark risk zoning map and uses it as the initialization condition of the game matrix; for each risk level zone, the system assigns a set of potential dark reaction type assumptions related to the typical pollutant combination in the region to the row strategies in the game matrix according to the preset rules, and assigns a set of recommended micro-perturbation basic actions that match the hydrogeochemical conditions of the region to the column strategies, thereby completing the filling of the initial values ​​of the game matrix.

Images