Method and apparatus for identifying structure of DNAPL source zone based on bayesian monitoring design

Abstract

The application discloses a DNAPL pollution source area structure identification method and equipment based on a Bayesian monitoring design, and belongs to the technical field of pollution hydrogeological exploration. The method firstly generates an initial parameter field through a random percolation model, and runs a forward model to obtain concentration simulation data; then, relative entropy is calculated by using a Bayesian monitoring design to determine an optimal observation point and obtain measured data; subsequently, the parameter field is iteratively updated based on an ensemble smoother multi-data assimilation method, and the monitoring design and data assimilation process are cyclically executed; finally, the forward simulation and inner optimization process are repeatedly performed through an outer loop until the optimal monitoring design is obtained, and the final estimation result of the pollution source area structure is output. The application realizes more accurate pollution source area structure identification under limited cost or limited monitoring data by using the Bayesian principle, and significantly improves the identification accuracy and efficiency of the DNAPL pollution source area structure in a strong heterogeneous site.

Description

Technical Field

[0001] This invention relates to the field of pollution hydrogeological exploration technology, specifically to a DNAPL pollution source area structure identification method and equipment based on Bayesian monitoring design. Background Technology

[0002] In groundwater environments, dense non-aqueous phase liquids (DNAPLs) are characterized by high density, low viscosity, high volatility, strong penetrability, and recalcitrant degradation. Compared to light non-aqueous phase liquids (LNAPLs), DNAPLs have a wider migration range and more complex structures in the groundwater environment, significantly increasing the difficulty of remediation. Therefore, a key prerequisite for the effective remediation of DNAPL-contaminated sites is the accurate characterization of the spatial distribution structure of the pollution source area and the pollution plume, and the distribution morphology of the pollution plume is mainly determined by the internal structure of the pollution source area.

[0003] The structure of a contamination source area essentially describes the complex mass distribution characteristics of DNAPL in space, specifically involving the breadth of the contamination range, pollutant saturation, occurrence form, and spatial location distribution. This structure is significantly controlled by the heterogeneity of the aquifer medium (such as the spatial distribution of permeability coefficients), and even small changes in permeability coefficients can significantly affect the structure of the contamination source area. Therefore, accurately characterizing the structure of the contamination source area and the spatial heterogeneity of permeability coefficients has become the most critical technical aspect in the remediation of DNAPL-contaminated sites.

[0004] To achieve high-precision identification of DNAPL pollution source area structures, traditional survey methods generally rely on intensive borehole sampling. These methods are costly, time-consuming, and have low feasibility in complex real-world sites. To address the limited monitoring data, designing an effective monitoring network to obtain monitoring data that maximizes information gain is crucial. However, traditional optimization methods (such as statistical methods, spatial interpolation, optimization-based methods, and hybrid methods) are primarily designed for optimizing monitoring networks for groundwater point source pollution identification and are difficult to apply to sampling optimization for complex DNAPL pollution source area structures.

[0005] Therefore, in order to improve the accuracy of the characterization of pollution source area structure under the premise of limited cost (or limited monitoring data), it is necessary to develop an adaptive optimization sampling strategy for DNAPL pollution source area structure identification, so as to achieve optimal identification of pollution source area structure and thus provide reliable technical support for the accurate remediation of contaminated sites. Summary of the Invention

[0006] Purpose of the invention: To address the shortcomings of traditional uniform sampling schemes, such as high cost and long implementation cycle, and the difficulty of applying traditional optimization methods to the sampling optimization of DNAPL contamination source area structures with irregular shapes and significant influence from the heterogeneity of aquifer media, this invention proposes an adaptive multi-level sampling optimization method based on Bayesian monitoring design. This method aims to achieve more accurate identification of contamination source area structures with limited cost (or limited monitoring data). The adaptive multi-level sampling optimization framework based on Bayesian monitoring design of this invention can provide a reliable decision-making basis for risk assessment and remediation of DNAPL contaminated sites.

[0007] Technical solution: To achieve the above-mentioned objectives, the present invention adopts the following technical solution:

[0008] A method for identifying the structure of DNAPL pollution source regions based on Bayesian monitoring includes the following steps:

[0009] An initial saturation field and an initial permeability field are generated using a stochastic permeation model. Based on these initial saturation fields and initial permeability fields, a forward model of groundwater flow and solute transport is run to obtain simulated pollutant concentration data.

[0010] Based on the initial saturation field, initial permeability field, and the simulated pollutant concentration data, the relative entropy is calculated using the Bayesian monitoring design method. The location with the maximum relative entropy is determined as the optimal observation point, and observation data is acquired at the optimal observation point. The observation data includes saturation data, permeability data, and concentration data.

[0011] Based on the initial saturation field, the initial permeability coefficient field, and the observed data, the ESMDA (Ensemble Smoothing Multiple Data Assimilation) method is used for iterative updates until the maximum number of iterations is reached, thereby obtaining the updated saturation field and the updated permeability coefficient field.

[0012] The updated saturation field and permeability coefficient field are used as new initial parameter fields. Forward modeling, monitoring point selection and data assimilation are carried out again until the target number of wells is reached, the optimal monitoring design is obtained, and the corresponding estimation results of saturation field and effective permeability coefficient field are obtained.

[0013] Furthermore, the forward model for groundwater flow and solute transport is as follows:

[0014] The presence of DNAPL will be due to the effective permeability coefficient K eff = K i ·K r (S N This is considered in the groundwater flow model, where K i K represents the inherent hydraulic conductivity. r S represents the relative hydraulic conductivity. NThe saturation level is given; assuming the groundwater flow is steady, its governing equation is:

[0015]

[0016]

[0017] The boundary conditions are:

[0018]

[0019]

[0020] in, Let Γ be the gradient, q be the vector velocity, and h be the head; D For Dirichlet boundaries; h D For Γ D The given head on; Γ N This is the Neumann boundary; n represents perpendicular to Γ. N The outer unit vector;

[0021] Based on a steady groundwater flow field, assuming that DNAPL reaches local equilibrium, the DNAPL pollution source area can be considered as the Dirichlet boundary within the study area. Therefore, the DNAPL transport problem can be solved using the convection-dispersion equations for steady flow:

[0022]

[0023] The initial and boundary conditions are as follows:

[0024]

[0025]

[0026]

[0027] in, Porosity C is the diffusion coefficient, C is the concentration of DNAPL in the dissolved phase, C0 is the initial concentration of DNAPL, and t is the solute transport time. For effective porosity, q s Neumann boundary Γ N solute flux at C S For the solubility of DNAPL, For x k The preset concentration at the location, Ω represents the entire study area, and Γ represents the internal Dirichlet boundary region. D This is the DNAPL contamination source area.

[0028] Furthermore, the Bayesian monitoring design method specifically calculates the relative entropy using the following formula:

[0029]

[0030] Where m represents pollution source parameters, including the effective permeability coefficient lnK. eff Saturation S N d represents the observed values ​​of pollution source parameters, including the effective permeability coefficient lnK. eff Saturation S N Concentration C; c represents the optimal sampling scheme corresponding to the observed parameter values; Let c be the relative entropy of the sampling scheme. Let be the posterior distribution of the parameter. The prior distribution of the parameters is given; the optimal sampling position is selected from the sampling optimization scheme with the largest relative entropy.

[0031] Furthermore, the specific method for determining the optimal sampling location includes the following steps:

[0032] Initialization: Set the initial iteration state, the number of iterations j=1, and input all possible horizontal positions and uniformly selected vertical position parameters;

[0033] Selecting the optimal horizontal sampling position: Calculate the relative entropy of each horizontal position and select the position c with the maximum relative entropy. j This serves as the current optimal horizontal sampling point;

[0034] Obtaining posterior samples and inversion analysis: at the selected optimal level position c j The observations are obtained and ESMDA inversion is performed to update the posterior distribution of the parameters;

[0035] Apply spatial weight decay: Apply weight decay to the selected location and its adjacent areas;

[0036] Iterative loop: Increment the iteration count by 1, repeat the steps of selecting the optimal level sampling location, obtaining posterior samples and inversion analysis, and applying spatial weight decay until the termination condition is met, and finally output an adaptive single-level monitoring network design.

[0037] Furthermore, the specific method for determining the optimal sampling location includes the following steps:

[0038] Initialization: Set the initial iteration state, let the iteration number j=1, and input all possible horizontal and vertical sampling position parameters;

[0039] Adaptive optimization of vertical sampling position: For each horizontal position, the optimal vertical sampling point is dynamically selected from high to low based on the relative entropy index, and a weight decay mechanism is applied to the selected position and its adjacent area.

[0040] Selecting the optimal horizontal sampling position: Based on the sum of the relative entropy of all vertical points at each horizontal position, the optimal sampling position is dynamically selected from high to low, and finally the globally optimal horizontal well placement position is determined;

[0041] Iterative loop: Increment the iteration count by 1, repeat the adaptive optimization of the vertical sampling position and the selection of the optimal horizontal sampling position steps until the termination condition is met, and finally output an adaptive multi-level monitoring network design.

[0042] Furthermore, the ensemble smoother multi-data assimilation (ESMDA) method includes the following steps:

[0043] Choose the number of data assimilation times N a and the corresponding expansion coefficient for each assimilation iteration , where t=1,2,...,N a ;

[0044] Draw N samples from the prior distribution to form the initial sample parameter set;

[0045] The iteration runs starting from t=1. In each iteration, based on the model parameters updated in the previous iteration, the forward model is run to obtain the set of simulated values ​​corresponding to the samples:

[0046]

[0047] In the formula, , where is the model parameter obtained in the t-th iteration; Let f(m) be the simulated value obtained in the t-th iteration, and f(m) be the forward model.

[0048] Calculate the error covariance matrix using the following formula:

[0049]

[0050] In the formula, C MD Let C be the covariance matrix between the model parameters and the simulated values. DD The autocovariance matrix of the simulated values. The mean of the model parameters, This is the average of the simulated values; the superscript T indicates transpose.

[0051] Iterate N according to the following formula a This is used to update the sample parameter set and obtain an inverse estimate of the posterior distribution of the parameters:

[0052]

[0053] In the formula, C D Let d be the covariance matrix of the observation error. obs These are the observations after adding perturbations.

[0054] A DNAPL pollution source area structure identification system based on Bayesian monitoring design includes:

[0055] The forward modeling module is used to generate an initial saturation field and an initial permeability field through a stochastic permeation model. Based on the initial saturation field and the initial permeability field, the forward modeling model of groundwater flow and solute transport is run to obtain simulated data of pollutant concentration.

[0056] The monitoring point selection module is used to calculate the relative entropy based on the initial saturation field, the initial permeability field, and the simulated pollutant concentration data, using the Bayesian monitoring design method, to determine the location with the largest relative entropy as the optimal observation point, and to acquire observation data at the optimal observation point, the observation data including saturation data, permeability data, and concentration data;

[0057] The data assimilation and update module is used to iteratively update the initial saturation field, the initial permeability coefficient field, and the observed data using the ensemble smoother multiple data assimilation (ESMDA) method until the maximum number of iterations is reached, thereby obtaining the updated saturation field and the updated permeability coefficient field.

[0058] The iterative control and result output module is used to take the updated saturation field and permeability coefficient field as the new initial parameter fields, and re-perform forward simulation, monitoring point selection and data assimilation update until the target number of wells is reached, obtain the optimal monitoring design, and obtain the corresponding estimation results of saturation field and effective permeability coefficient field.

[0059] The present invention also provides an electronic device, comprising: one or more processors; a memory; and one or more programs, wherein the one or more programs are stored in the memory and configured to be executed by the one or more processors, wherein when the programs are executed by the processors, they implement the steps of the DNAPL contamination source region structure identification method based on Bayesian monitoring design as described above.

[0060] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described above.

[0061] The present invention also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described above.

[0062] Beneficial Effects: This invention addresses the scenario of limited monitoring data and uncertain groundwater pollution source areas by proposing a groundwater pollution source area identification method based on Bayesian monitoring design. Using relative entropy as the objective function, it quantifies the information content and uncertainty of monitoring data. By acquiring observations at the location of maximum relative entropy, it updates the understanding of the pollution source area and the heterogeneity of the aquifer medium. Combined with the ESMDA method, it inversely estimates the complex and realistic DNAPL source area structure, thus providing an efficient way to identify the pollution source area structure with limited cost (or limited monitoring data). This provides a reliable decision-making basis for the risk assessment and remediation of DNAPL contaminated sites and has significant scientific implications for groundwater pollution risk prevention and control. Attached Figure Description

[0063] Figure 1 This is a flowchart of the Bayesian monitoring design and parameter estimation process.

[0064] Figure 2 A flowchart for single-level adaptive sampling optimization;

[0065] Figure 3 A flowchart for multi-level adaptive sampling optimization;

[0066] Figure 4 Example diagrams of the effective permeability reference field and DNAPL saturation reference field in numerical experimental examples;

[0067] Figure 5 This shows the spatial distribution of monitoring wells and their corresponding sampling points in the four examples Cases 1-4 of the numerical experiment, as well as the regions of well sampling locations and saturation values ​​optimized by Bayesian monitoring design.

[0068] Figure 6 The results of permeability and saturation field inversion between the reference field and the four examples Cases 1-4;

[0069] Figure 7 The statistical distribution of the three key source region structural parameters for each embodiment;

[0070] Figure 8 Simulation results of the DNAPL contamination source region structure for each embodiment;

[0071] Figure 9 DNAPL mass decay curves between the pollution source region structure and the reference field for each embodiment. Detailed Implementation

[0072] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0073] This invention provides a method for identifying the structure of DNAPL pollution source regions based on Bayesian monitoring design, referring to... Figure 1 The method includes the following steps.

[0074] Step S1: The Stochastic Percolation Model (SIP) can simulate the leakage, infiltration, redistribution, and dissolution processes of DNAPL in heterogeneous aquifers. By inputting geostatistical parameters and DNAPL contamination source area parameters into the Stochastic Percolation Model (SIP), a series of saturation fields and corresponding effective permeability coefficient fields can be generated.

[0075] Step S2: Based on the parameter field generated in step S1, run the forward model of groundwater flow and solute transport to obtain pollutant concentration data.

[0076] The model for groundwater flow and solute transport is as follows:

[0077] The presence of DNAPL will be due to the effective permeability coefficient K eff = K i ·K r (S N This is considered in the groundwater flow model. Where K... i K represents the inherent hydraulic conductivity. r S represents the relative hydraulic conductivity. N Let be the saturation level. Assuming the groundwater flow is steady, its governing equation is:

[0078]

[0079]

[0080] The boundary conditions are:

[0081]

[0082]

[0083] in, Γ is the gradient, q is the vector velocity (m / s), and h is the head (m); D For Dirichlet boundaries; h D For Γ D The given head on; Γ N This is the Neumann boundary; n represents perpendicular to Γ. N The external unit vector.

[0084] Based on the aforementioned stable groundwater flow field, assuming that DNAPL reaches local equilibrium, the DNAPL pollution source area can be considered as the Dirichlet boundary within the study area. Therefore, the DNAPL transport problem can be solved using the convection-dispersion equations for steady flow:

[0085]

[0086] Solving this equation will yield the pollutant concentration.

[0087] The initial and boundary conditions are as follows:

[0088]

[0089]

[0090]

[0091] in, Porosity C is the diffusion coefficient, C is the concentration of DNAPL in the dissolved phase, C0 is the initial concentration of DNAPL, and t is the solute transport time. For effective porosity, q s Neumann boundary Γ N solute flux at C S For the solubility of DNAPL, For x k The preset concentration at the location, Ω represents the entire study area, and Γ represents the internal Dirichlet boundary region. D This is the DNAPL contamination source area.

[0092] Step S3: Based on the saturation and permeability coefficient data from Step S1 and the simulated pollutant concentration data obtained in Step S2, in the initial step of j=1, N=1000 sets of parameter samples m were collected from the prior distribution of the parameters. i (i=1,2,3,4,,,N), a Bayesian monitoring design is used to obtain the location with the maximum relative entropy, and observation data is obtained at that location.

[0093] The formula for calculating the relative entropy of sampling optimization scheme c is:

[0094]

[0095] Where m represents pollution source parameters (including: effective permeability coefficient lnK) eff Saturation S N ), d represents the observed values ​​of pollution source parameters (including: effective permeability coefficient lnK). eff Saturation S N (C), where c is the optimal sampling scheme corresponding to the observed parameter values. Let c be the relative entropy of the sampling scheme. Let be the posterior distribution of the parameter. The prior distribution of the parameter is .

[0096] Since the difference between the posterior and prior distributions of the parameters is caused by the observed values, the greater this difference, the more information the observed values ​​carry about the unknown parameters. Therefore, the optimal sampling location is the sampling optimization scheme that maximizes the relative entropy from the prior to the posterior of the parameters.

[0097] Traditional sampling optimization schemes often employ a single-level sampling optimization process, optimizing only the horizontal position of the monitoring well. (Refer to...) Figure 2 ,include:

[0098] (1) Initialization: Set the initial iteration state (j=1) and input all possible horizontal positions and uniformly selected vertical position parameters.

[0099] (2) Select the optimal horizontal sampling position: Calculate the relative entropy of each horizontal position and select the position with the maximum relative entropy (c). j ), which is the current optimal horizontal sampling point.

[0100] (3) Obtaining posterior samples and inversion analysis: at the selected optimal level position (c j Obtain the observations and perform ESMDA inversion to update the posterior distribution of the parameters.

[0101] (4) Apply spatial weight decay: Apply a 10% weight decay to the selected location and its adjacent area (e.g., ±2 grid intervals) to avoid the subsequent sampling locations being too concentrated and to improve the rationality of spatial distribution.

[0102] (5) Iterative loop: Repeat steps (2) to (4) (j=j+1) until the termination condition (i.e. the preset number of wells) is met, and finally an adaptive single-level monitoring network design is output.

[0103] Furthermore, this invention specifically proposes an adaptive multi-level sampling strategy based on Bayesian monitoring design, which not only adaptively optimizes the location of the monitoring well but also adaptively optimizes the sampling depth. (Refer to...) Figure 3 ,include:

[0104] (1) Initialization: Set the initial iteration state (j=1) and input all possible horizontal and vertical sampling position parameters.

[0105] (2) Adaptive optimization of vertical sampling position: For each horizontal position, the system dynamically selects the optimal vertical sampling point from high to low based on the relative entropy index, and avoids excessive concentration of positions through the weight decay mechanism, thereby achieving refined depth sampling.

[0106] (3) Select the optimal horizontal sampling position: Based on the sum of the relative entropy of all vertical points at each horizontal position, the optimal sampling position is dynamically selected from high to low, and finally the globally optimal horizontal well placement position is determined.

[0107] (4) Iteration: Repeat steps (2) and (3) (j=j+1) until the termination condition is met, and finally output an adaptive multi-level monitoring network design.

[0108] The single-level optimization strategy proposed by the traditional sampling optimization scheme can significantly improve the amount of information in the acquired observations compared with the traditional uniform well layout. If the multi-level optimization proposed in this method is further adopted, the sampling depth can be simultaneously optimized based on the horizontal positioning of the monitoring wells, thereby obtaining a better monitoring network design and providing effective guidance for screening depth and post-drilling sample analysis.

[0109] Step S4: Based on the series of saturation fields and permeability coefficient fields generated in step S1, and the saturation, permeability coefficient and concentration data obtained in step S3, the saturation fields and permeability coefficient fields are iteratively updated using the ESMDA method. After each iteration, the original parameter fields are replaced with new saturation fields and permeability coefficient fields. Steps S3 to S4 are repeated until the maximum number of iterations is reached.

[0110] The ESMDA (Ensemble Smoother with Multiple Data Assimilation) method can be summarized as follows:

[0111] (1) Select the number of data assimilation times (N) a ) and the corresponding expansion coefficient for each assimilation iteration ( , where t=1,2,...,N a ).

[0112] (2) Draw N samples from the prior distribution to form the initial sample parameter set.

[0113] (3) Start iteratively from t=1, and in each iteration, update the model parameters m obtained in the previous iteration. i Run the forward model to obtain the set of simulated values ​​corresponding to the sample:

[0114]

[0115] In the formula, , where is the model parameter obtained in the t-th iteration; f(m) is the forward model.

[0116] (4) Calculate the error covariance matrix according to the following formula:

[0117]

[0118] In the formula, C MD Let C be the covariance matrix between the model parameters and the simulated values. DD The autocovariance matrix of the simulated values. The mean of the model parameters, This is the average of the simulated values, and the superscript T indicates transpose.

[0119] (5) Iterate N according to the following formula a This is used to update the sample parameter set and obtain an inverse estimate of the posterior distribution of the parameters:

[0120]

[0121] In the formula, C D Let d be the covariance matrix of the observation error. obs These are the observations after adding perturbations.

[0122] Step S5: Based on the saturation field and permeability coefficient field finally updated in step S4, replace the parameter field generated in step S1, repeat steps S2 to S4 until the target number of wells is reached, obtain the optimal monitoring design, and obtain the corresponding estimation results of the saturation field and effective permeability coefficient field.

[0123] The following numerical experiments further illustrate the feasibility of Bayesian monitoring design for identifying groundwater pollution source areas. The example is a two-dimensional heterogeneous confined aquifer with an area of ​​40m × 20m, divided into 80 × 40 = 3200 grids, each grid being 0.5m in length. In this example, the DNAPL pollution source area consists of a single component, DNAPL. Figure 4 (a) in the figure shows the effective permeability reference field for this example. Figure 4 (b) shows the DNAPL saturation reference field for this example. The parameter values ​​used in the model are shown in Table 1. The symbols in parentheses of the parameters indicate units.

[0124] Table 1. Reference model parameter settings in numerical experiments

[0125]

[0126] I x ,I z Let X represent the correlation lengths in the x and z directions, respectively. 2Let x represent the variance, K be the permeability coefficient, x and x' be the coordinates of two locations in space, |x - x'| represent the Euclidean distance between these two points, and I be the correlation length.

[0127] To better verify the optimization effect of the adaptive multi-level sampling strategy based on Bayesian monitoring design, four embodiments were designed to highlight its advantages over traditional uniform sampling strategies and single-level adaptive sampling strategies. Embodiment 1 (Case 1) and Embodiment 2 (Case 2) employ uniform sampling designs with 10 and 3 monitoring wells, respectively; Embodiment 3 (Case 3) adopts a single-level adaptive sampling strategy, optimizing the horizontal well placement of the 3 monitoring wells. In the above schemes, all monitoring wells are arranged with 10 sampling points at 1-meter intervals along the vertical direction (Z direction). In contrast, Embodiment 4 (Case 4) applies a multi-level adaptive sampling strategy, simultaneously optimizing the horizontal spatial distribution and vertical sampling depth of the 3 monitoring wells. The adaptive sampling strategy of Embodiment 3 is as follows: Figure 2 As shown, the adaptive sampling strategy used in Example 4 is as follows: Figure 3 As shown, the number of monitoring wells and the layout parameters of sampling points in each embodiment are summarized in Table 2.

[0128] Table 2. Number of monitoring wells and sampling point layout parameters in each embodiment.

[0129]

[0130] Figure 5 Figures (a)-(d) show the spatial distribution of monitoring wells and their corresponding sampling points in Cases 1-4, respectively. Based on the Bayesian monitoring design, the optimal sampling locations for Cases 3 and 4 are shown below. Figure 5 (c) and Figure 5 As shown in (d) in the figure. The well sampling locations optimized by Bayesian monitoring design are those with saturation values ​​greater than 0 (S). N Regions with values ​​greater than 0 show significant spatial correlation. Figure 5 (e)-(f)). Among them, Case 4 (see Figure 5 (e) in this paper outperforms Case 3 in accurately identifying high-saturation regions (see [reference]). Figure 5 (f) in the figure captures the structural distribution of pollution source areas at different aquifer depths through vertical optimization sampling. This multi-level sampling optimization strategy obtains more information than sampling strategies that only optimize monitoring well locations and traditional uniform sampling strategies, providing key observational data for efficiently estimating the DNAPL saturation field (SZA) and permeability field.

[0131] Figure 6The permeability and saturation field inversion results were compared between the reference field and Cases 1-4. The results show that both Case 1 and Case 4 can effectively characterize the DNAPL contamination source region structure similar to the reference field. Case 1 shows a higher similarity in characterizing the saturation field because it selects more well locations, maximizing the use of observational data from the reference field for inversion estimation; however, more well locations also mean a significant increase in cost. Compared to Case 1, the cost of Cases 2-4 is significantly reduced. However, with the decrease in cost, the estimation accuracy of Case 2, with its uniform sampling, also decreases significantly. Figure 6 As shown in the third row, Case 2 failed to accurately reconstruct the spatial distribution of the saturation field, and the estimated saturation value was significantly lower than the reference data. Furthermore, in the permeability field, Case 2's ability to identify high and low permeability regions was also very limited. Compared to Case 2, Case 3 roughly captured the structure of the saturation field, with a higher saturation value than Case 2 and closer to the reference field, but there is still room for improvement. Case 4 performed better than Case 2 and Case 3 in reconstructing the pollution source area structure, and its estimated saturation value was closest to Case 1 and the reference saturation field. In addition, Case 4 successfully identified low-permeability regions in the permeability field. These results demonstrate that the multi-level adaptive sampling strategy based on Bayesian monitoring design has significant advantages in both efficiency and accuracy.

[0132] To further quantify and evaluate the accuracy of the saturation field estimation in the four examples, the source region projected area A and the average relative permeability at the local scale were introduced as structural indicators of the pollution source region. and the average local scale DNAPL saturation This is used to evaluate the accuracy of the estimated pollution source region structure. A represents the area obtained by projecting the source region onto a control plane perpendicular to the flow direction, normalized to the total area of ​​the control plane. This index can reasonably approximate the contribution of the saturation domain to the total mass flux. It represents the average relative permeability within the source region, which is positively correlated with flow velocity and mass transfer rate, thereby promoting the passage of more fluid through the source region. This represents the average DNAPL saturation level within the source region. Both local relative permeability and DNAPL saturation are only affected by S... N >5% of the numerical units were used for calculation. All three variables were normalized based on their respective initial values ​​to quantify the change in the source region relative to its initial state. Furthermore, this invention uses the ganglia-to-pool ratio (GTP) to quantify the ratio of the NAPL volume at or below residual saturation to the NAPL volume above residual saturation within the pollution source region, expressed as:

[0133]

[0134] In the formula, The maximum residual saturation is represented by this index, which distinguishes the distribution of DNAPL in the ganglia and pool phases. In the early stages of source region evolution, mass flux dynamics are controlled by the dissolution process of the ganglia; subsequently, the source region gradually transitions from being dominated by ganglia to being dominated by pools. Therefore, the initial source region structure may have a significant impact on mass flux. This ratio, by predicting mass conversion efficiency and long-term source region behavior, can provide crucial information for the formulation of remediation strategies.

[0135] Figure 7 The statistical distribution of three key source region structural parameters for each embodiment is presented in box plot form. Each blue horizontal dashed line in the figure represents an index value calculated from the reference SZA; each box plot corresponds to one embodiment and shows the distribution of the corresponding index in all posterior implementations of that embodiment. Figure 7 In the diagram, (a) represents the projected area A of the source region, and (b) represents the average relative permeability at the local scale. (c) represents the average DNAPL saturation at the local scale. Furthermore, the results of GTP from different embodiments were compared, see [link / reference]. Figure 7 (d) Case 2 uses only 3 uniformly distributed monitoring wells to estimate the unknown S. N The field is not fully represented, thus failing to reproduce the true GTP. In Case 3, compared to Case 2, the source region structure parameters are closer to the reference values, but the results still differ somewhat from Case 4. Unlike the passive uniform sampling of Case 2, Case 4 employs an adaptive sampling framework driven by Bayesian experimental design, strategically prioritizing well locations with the highest information value. This adaptive decision-making process maximizes relative entropy and iteratively optimizes the monitoring network configuration, thereby systematically assimilating observational data from the most informative sampling points and achieving a high-fidelity reconstruction of the pollution source region structure. In summary, from the perspective of characterizing the pollution source region structure, Case 4 can effectively characterize the pollution source region structure within a certain cost, showing a significant advantage over other schemes.

[0136] To assess the impact of inversion accuracy on pollution source area lifetime prediction, this invention uses a DNAPL migration and dissolution model to predict the natural decay process of pollution source areas based on reference values ​​and the permeability field and DNAPL saturation field obtained from inversion estimation. Figure 8The simulation results of the DNAPL contamination source region structure are presented. The results of the reference source region in the first row indicate that all residual DNAPL will be essentially depleted within 30 years, with only a few DNAPL pooling areas remaining; after approximately 50 years, the entire contamination source region will completely dissipate. A comparison between the reference results and the simulation results shows that Case 1 reconstructed the decay process of the contamination source region quite well, and Case 4 also roughly reproduced this process; in contrast, Case 2 and Case 3 underestimated the total mass of DNAPL and failed to accurately reproduce the actual evolution of the contamination source region over 50 years.

[0137] To quantitatively compare the ability of different sampling strategies to assess the lifetime of DNAPL contamination source areas, Figure 9 The DNAPL mass decay curves estimated based on various sampling strategies and compared to the reference field were analyzed. The results showed that Case 1 established a performance baseline, with its predictions closest to the reference field. Case 2 incorrectly predicted complete decay of the source region in approximately 45 years, while Case 3, although successfully predicting a source region lifetime of 50 years, showed a significant discrepancy between its estimated mass and the reference value. In contrast, Case 4 achieved 99.7% accuracy in initial mass estimation compared to the Case 1 baseline and successfully predicted the pollutant residue at 50 years, demonstrating a 44.5% improvement in prediction accuracy compared to Case 2. This provides a more reliable initial state for DNAPL repair prediction.

[0138] This invention also provides a DNAPL pollution source region structure identification system based on Bayesian monitoring design, comprising:

[0139] The forward modeling module is used to generate an initial saturation field and an initial permeability field through a stochastic permeation model. Based on the initial saturation field and the initial permeability field, the forward modeling model of groundwater flow and solute transport is run to obtain simulated data of pollutant concentration.

[0140] The monitoring point selection module is used to calculate the relative entropy based on the initial saturation field, the initial permeability field, and the simulated pollutant concentration data, using the Bayesian monitoring design method, to determine the location with the largest relative entropy as the optimal observation point, and to acquire observation data at the optimal observation point, the observation data including saturation data, permeability data, and concentration data;

[0141] The data assimilation and update module is used to iteratively update the initial saturation field, the initial permeability coefficient field, and the observed data using the ensemble smoother multiple data assimilation (ESMDA) method until the maximum number of iterations is reached, thereby obtaining the updated saturation field and the updated permeability coefficient field.

[0142] The iterative control and result output module is used to take the updated saturation field and permeability coefficient field as the new initial parameter fields, and re-perform forward simulation, monitoring point selection and data assimilation update until the target number of wells is reached, obtain the optimal monitoring design, and obtain the corresponding estimation results of saturation field and effective permeability coefficient field.

[0143] The DNAPL pollution source area structure identification system based on Bayesian monitoring design can realize all the technical solutions in the above method embodiments. The functions of each functional module can be specifically implemented according to the methods in the above method embodiments. The specific implementation process can be referred to the relevant descriptions in the above method embodiments, which will not be repeated here.

[0144] The present invention also provides an electronic device, comprising: one or more processors; a memory; and one or more programs, wherein the one or more programs are stored in the memory and configured to be executed by the one or more processors, wherein when the programs are executed by the processors, they implement the steps of the DNAPL contamination source region structure identification method based on Bayesian monitoring design as described above.

[0145] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described above.

[0146] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, computer devices, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0147] This invention is described with reference to a flowchart of a method according to embodiments of the invention. It should be understood that each step in the flowchart and combinations thereof can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing device to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing device, generate instructions for implementing the process. Figure 1 A system that specifies functions in one or more processes.

[0148] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including an instruction set implemented in a process. Figure 1 The function specified in one or more processes.

[0149] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 Steps of a specified function in one or more processes.

Claims

1. A method for identifying the structure of DNAPL pollution source regions based on Bayesian monitoring design, characterized in that, Includes the following steps: An initial saturation field and an initial permeability field are generated using a stochastic permeation model. Based on these initial saturation fields and initial permeability fields, a forward model of groundwater flow and solute transport is run to obtain simulated pollutant concentration data. Based on the initial saturation field, initial permeability field, and the simulated pollutant concentration data, the relative entropy is calculated using the Bayesian monitoring design method. The location with the maximum relative entropy is determined as the optimal observation point, and observation data is acquired at the optimal observation point. The observation data includes saturation data, permeability data, and concentration data. Based on the initial saturation field, the initial permeability coefficient field, and the observed data, the ESMDA (Ensemble Smoothing Multiple Data Assimilation) method is used for iterative updates until the maximum number of iterations is reached, thereby obtaining the updated saturation field and the updated permeability coefficient field. The updated saturation field and permeability coefficient field are used as new initial parameter fields to perform forward modeling, monitoring point selection and data assimilation updates again until the target number of wells is reached, the optimal monitoring design is obtained, and the corresponding estimation results of saturation field and effective permeability coefficient field are obtained. Specifically, the Bayesian monitoring design method calculates the relative entropy using the following formula: ， Where m represents pollution source parameters, including the effective permeability coefficient K. eff Saturation S N d represents the observed values ​​of pollution source parameters, including the effective permeability coefficient K. eff Saturation S N Concentration C; c represents the optimal sampling scheme corresponding to the observed parameter values; Let c be the relative entropy of the sampling scheme. Let be the posterior distribution of the parameter. The prior distribution of the parameters is used; the optimal sampling position is selected as the sampling optimization scheme with the largest relative entropy. The ESMDA (Ensemble Smoothing Multiple Data Assimilation) method includes the following steps: Choose the number of data assimilation times N a and the corresponding expansion coefficient for each assimilation iteration , where t=1,2,...,N a ; Draw N samples from the prior distribution to form the initial sample parameter set; The iteration runs starting from t=1. In each iteration, based on the model parameters updated in the previous iteration, the forward model is run to obtain the set of simulated values ​​corresponding to the samples: ， In the formula, , where is the model parameter obtained in the t-th iteration; Let f(m) be the simulated value obtained in the t-th iteration, and f(m) be the forward model. Calculate the error covariance matrix using the following formula: ， In the formula, C MD Let C be the covariance matrix between the model parameters and the simulated values. DD The autocovariance matrix of the simulated values. The mean of the model parameters, This is the average of the simulated values; the superscript T indicates transpose. Iterate N according to the following formula a This is used to update the sample parameter set and obtain an inverse estimate of the posterior distribution of the parameters: ， In the formula, C D Let d be the covariance matrix of the observation error. obsi These are the observations after adding perturbations.

2. The method according to claim 1, characterized in that, The forward model for groundwater flow and solute transport is as follows: The presence of DNAPL will be due to the effective permeability coefficient K eff = K i ·K r (S N This is considered in the groundwater flow model, where K i K represents the inherent hydraulic conductivity. r S represents the relative hydraulic conductivity. N The saturation level is given; assuming the groundwater flow is steady, its governing equation is: ， ， The boundary conditions are: ， ， in, Let Γ be the gradient, q be the vector velocity, and h be the head; D For Dirichlet boundaries; h D For Γ D The given head on; Γ N This is the Neumann boundary; n represents perpendicular to Γ. N The outer unit vector; Based on a steady groundwater flow field, assuming that DNAPL reaches local equilibrium, the DNAPL pollution source area can be considered as the Dirichlet boundary within the study area. Therefore, the DNAPL transport problem can be solved using the convection-dispersion equations for steady flow: ， The initial and boundary conditions are as follows: ， ， ， in, Porosity C is the diffusion coefficient, C is the concentration of DNAPL in the dissolved phase, C0 is the initial concentration of DNAPL, and t is the solute transport time. For effective porosity, q s Neumann boundary Γ N Solute flux at C S For the solubility of DNAPL, For x k The preset concentration at the location, Ω represents the entire study area, and Γ represents the internal Dirichlet boundary region. D This is the DNAPL contamination source area.

3. The method according to claim 1, characterized in that, The specific method for determining the optimal sampling location includes the following steps: Initialization: Set the initial iteration state, the number of iterations j=1, and input all possible horizontal positions and uniformly selected vertical position parameters; Selecting the optimal horizontal sampling position: Calculate the relative entropy of each horizontal position and select the position c with the maximum relative entropy. j This serves as the current optimal horizontal sampling point; Obtaining posterior samples and inversion analysis: at the selected optimal level position c j The observations are obtained and ESMDA inversion is performed to update the posterior distribution of the parameters; Apply spatial weight decay: Apply weight decay to the selected location and its adjacent areas; Iterative loop: Increment the iteration count by 1, repeat the steps of selecting the optimal level sampling location, obtaining posterior samples and inversion analysis, and applying spatial weight decay until the termination condition is met, and finally output an adaptive single-level monitoring network design.

4. The method according to claim 1, characterized in that, The specific method for determining the optimal sampling location includes the following steps: Initialization: Set the initial iteration state, let the iteration number j=1, and input all possible horizontal and vertical sampling position parameters; Adaptive optimization of vertical sampling position: For each horizontal position, the optimal vertical sampling point is dynamically selected from high to low based on the relative entropy index, and a weight decay mechanism is applied to the selected position and its adjacent area. Selecting the optimal horizontal sampling position: Based on the sum of the relative entropy of all vertical points at each horizontal position, the optimal sampling position is dynamically selected from high to low, and finally the globally optimal horizontal well placement position is determined; Iterative loop: Increment the iteration count by 1, repeat the adaptive optimization of the vertical sampling position and the selection of the optimal horizontal sampling position steps until the termination condition is met, and finally output an adaptive multi-level monitoring network design.

5. A DNAPL pollution source region structure identification system based on Bayesian monitoring design, used to implement the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described in any one of claims 1-4, characterized in that, include: The forward modeling module is used to generate an initial saturation field and an initial permeability field through a stochastic permeation model. Based on the initial saturation field and the initial permeability field, the forward modeling model of groundwater flow and solute transport is run to obtain simulated data of pollutant concentration. The monitoring point selection module is used to calculate the relative entropy based on the initial saturation field, the initial permeability field, and the simulated pollutant concentration data, using the Bayesian monitoring design method, to determine the location with the largest relative entropy as the optimal observation point, and to acquire observation data at the optimal observation point, the observation data including saturation data, permeability data, and concentration data; The data assimilation and update module is used to iteratively update the initial saturation field, the initial permeability coefficient field, and the observed data using the ensemble smoother multiple data assimilation (ESMDA) method until the maximum number of iterations is reached, thereby obtaining the updated saturation field and the updated permeability coefficient field. The iterative control and result output module is used to take the updated saturation field and permeability coefficient field as the new initial parameter fields, and re-perform forward simulation, monitoring point selection and data assimilation update until the target number of wells is reached, obtain the optimal monitoring design, and obtain the corresponding estimation results of saturation field and effective permeability coefficient field.

6. An electronic device, characterized in that, include: One or more processors; Memory; And one or more programs, wherein the one or more programs are stored in the memory and configured to be executed by the one or more processors, wherein when the programs are executed by the processors, they implement the steps of the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described in any one of claims 1-4.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described in any one of claims 1-4.

8. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the DNAPL pollution source region structure identification method based on Bayesian monitoring design as described in any one of claims 1-4.

Images