A method for efficiently calculating the failure probability of groundwater pollution based on thermodynamic integration

Abstract

The present invention discloses a method for efficiently calculating the failure probability of groundwater pollution based on thermodynamic integration, comprising the following steps: calculating the failure probability formula coefficients by the Monte Carlo method, selecting different parameters that determine the smoothness of the indicator function as the temperature gradient in the parallel annealing algorithm; generating the relevant initial field of parameters that obey the Gaussian distribution by sequential Gaussian; generating the recommended sample of the current permeability field according to the pre-conditioned Crank-Nicholson recommended distribution; calculating the acceptance probability within the Markov chain to determine whether the sampling is accepted within the chain; calculating the acceptance probability between Markov chains to determine whether the sampling is accepted between chains; repeating the above process to sample the Gaussian permeability field sample of the next aquifer until it is determined to terminate the sampling, and calculating the failure probability according to the thermodynamic integration formula. This method is an efficient method for calculating the failure probability of groundwater pollution, and is more efficient than the Monte Carlo and subset simulation methods in dealing with the failure probability problem of groundwater pollution.

Description

Technical Field

[0001] The invention relates to the technical field of hydrological statistics, and in particular to a method for efficiently calculating the failure probability of groundwater pollution based on thermodynamic integration. Background Art

[0002] With the development of industry and agriculture, including the use of domestic garbage, pesticides and fertilizers, industrial wastewater and other pollutants, the groundwater environment has been damaged to a certain extent. However, due to the uncertainty of hydrogeological parameters such as permeability coefficient, the groundwater system has uncertainty, which makes the calculation of groundwater pollution difficult. Although some advanced technologies can reduce this uncertainty at this stage, there are also certain probability problems, especially some toxic substances with very low content but significant harm to human health. It is necessary to strictly control their content in groundwater to avoid negative impacts on human production and life. Therefore, calculating the probability of groundwater contamination can provide certain help for groundwater quality evaluation and pollution remediation, thereby reducing the risk of groundwater pollution. At the same time, it has certain guiding significance for building a decision-making model for water environment risk assessment system.

[0003] Failure probability refers to the probability that a certain indicator in the system is not within its reasonable threshold when considering various uncertainties in the system. It is used in civil engineering, petroleum and natural gas engineering and other fields. In groundwater pollution failure probability, it refers to the probability that the concentration of pollutants in the water body exceeds a certain standard when various uncertainties in the groundwater system are considered (mainly the uncertainty of the permeability coefficient). Due to the influence of various uncertainties, there are certain difficulties in calculating the failure probability. Commonly used methods include subset simulation method and Monte Carlo method. Among them, Monte Carlo method is the most traditional and basic calculation method. The essence of its method is the law of large numbers, so its calculation results are relatively accurate. However, for low-probability events, the required sample size is large, resulting in low calculation efficiency. On this basis, relevant researchers proposed the subset simulation method to calculate the failure probability; the subset simulation method requires a series of intermediate events. Since it converts low-probability events into high-probability events and multiplies them, it can improve the calculation efficiency and is widely used in failure probability problems. Although subset simulation has advantages over Monte Carlo method in terms of calculation time efficiency, the accuracy of the calculated failure probability is relatively poor. Summary of the invention

[0004] Purpose of the invention: In order to overcome the shortcomings of the background technology, the present invention discloses a method for efficiently calculating the failure probability of groundwater pollution based on thermodynamic integration. The method couples the parallel annealing algorithm with thermodynamic integration, changes the original indicator function, sets multiple chains, and integrates the results using the thermodynamic integration formula to obtain the failure probability of groundwater pollution.

[0005] Technical solution: The method disclosed in the present invention for efficiently calculating the failure probability of groundwater pollution based on thermodynamic integration comprises the following steps:

[0006] S1. Based on the underground aquifer, the coefficients of the groundwater pollution failure probability formula are calculated by the Monte Carlo method, and the parameters that determine the smoothness of the indicator function and the weight of the Markov chain corresponding to each indicator function are selected according to the Gaussian integral criterion;

[0007] S2. Generate the initial field of groundwater aquifer permeability coefficient that obeys Gaussian distribution through sequential Gaussian, calculate the concentration of aquifer pollutants, find the response value according to the function function, and bring the response value into the indicator function to calculate the result;

[0008] S3, generate a suggested sample of the current aquifer hydraulic conductivity field according to the pre-conditioned Crank-Nicholson suggested distribution, calculate the aquifer pollutant concentration, find the response value according to the functional function, and substitute the response value into the indicator function to calculate the result;

[0009] S4, calculating the acceptance probability of the permeability field according to the result obtained by the indicator function in the Markov chain, and judging whether the sample in the chain accepts sampling;

[0010] S5, calculating the acceptance probability of the permeability field between the Markov chains according to the results obtained from the indicator function, and determining whether the samples between the chains are accepted for sampling;

[0011] S6. Repeat the above process to sample the Gaussian permeability coefficient field of the next aquifer until the sampling is terminated and the failure probability of groundwater contamination is calculated according to the thermodynamic integral formula.

[0012] Among them, in S1, different parameters σ that determine the smoothness of the indicator function are selected according to the Gaussian integral criterion to generate different indicator functions I σ , and then obtain the temperature gradient in the parallel annealing algorithm and the weight w of each Markov chain k , and calculate the coefficient P of the groundwater pollution failure probability formula based on the indicator function and the functional function by the Monte Carlo method a ,

[0013]

[0014]

[0015]

[0016] Where N is the number of Monte Carlo sampling, I σ is the indicator function, I a To calculate P aThe indicator function used, G(x) is the functional function, which calculates the difference between the pollutant concentration and the pollutant concentration standard through simulation, and the response value is y.

[0017] Furthermore, the failure probability of groundwater contamination in S2 is mainly determined by the uncertainty of the aquifer permeability field. The initial state time t = 0 is set, and the initial state x of the aquifer permeability field that obeys the Gaussian distribution is generated using sequential Gaussian. 0 , calculate the aquifer pollutant concentration, functional function response value and indicator function value

[0018] x 0 ~N(0,C)

[0019] Where C is the variance of Gaussian distribution.

[0020] Furthermore, the preconditioned Crank-Nicholson technique in S3 can effectively handle high-dimensional sampling problems, so the proposed sample x' of the aquifer permeability field is generated according to its proposed distribution, and the aquifer pollutant concentration, functional function response value and indicator function value are calculated.

[0021]

[0022] In the formula, β is the jump factor and ε is the interference term.

[0023] Furthermore, S4 calculates the sample acceptance probability α in the Markov chain according to the indicator function value obtained in S3, generates a random number u∈[0,1] that obeys a uniform distribution, and determines whether to accept the sample.

[0024]

[0025]

[0026] S5 For Markov chains of different temperatures, calculate the indirect acceptance probability α of the Markov chain s , use the same method as in the chain to determine whether to accept sampling

[0027]

[0028] Furthermore, S6 repeats sampling until the termination condition is reached, and the groundwater pollution failure probability value q of each Markov chain is calculated. k , and calculate the failure probability P of groundwater pollution according to the thermodynamic integral formula F ,

[0029]

[0030]

[0031]

[0032] In the formula, w k is the kth weight, q k is the parameter σ k The failure probability of the corresponding Markov chain.

[0033] Beneficial effects: Compared with the prior art, the advantages of the present invention are: first, under the premise of ensuring the accuracy of calculating the failure probability of groundwater pollution, the calculation time efficiency of the present invention is greatly improved compared with the classical Monte Carlo method; second, under the condition of roughly the same calculation time, compared with the traditional subset simulation method for calculating the failure probability, the present invention can obtain more accurate calculation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 Schematic diagram of the location of the observation wells in the embodiment, where the number of observation wells is 9, 36 and 49 respectively;

[0035] Figure 2 It is the algorithm flow chart of the present invention;

[0036] Figure 3 : is a graph showing the variation of the mean failure probability with the number of samples when the number of observation wells is 9, 36, and 49 respectively in the embodiment;

[0037] Figure 4 : is a graph showing the variation of the mean failure probability with the number of observation wells when the sample numbers are 6000, 8000 and 10000 respectively in the embodiment;

[0038] Figure 5 Graph showing how RMSE changes with the number of samples and the number of observation wells in the embodiment. DETAILED DESCRIPTION

[0039] The technical solution of the present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0040] In order to verify the present invention, it is verified in a saturated confined aquifer, the aquifer space is 5000 [L] × 5000 [L] in two dimensions, the thickness is 50 [L], and it is discretized into 50 × 50 × 1 cells, each cell is 100 [L] × 100 [L] × 50 [L] (no specific units are used, only dimensions are given, and any set of corresponding units will get corresponding results), the east and west boundaries are water-proof boundaries and the given water levels are 300 [L] and 80 [L] respectively, and the south and north boundaries are water-proof boundaries. In this embodiment, the reference field of the permeability coefficient of the groundwater aquifer obeys a logarithmic multi-Gaussian distribution. The reference field distribution satisfies the mean of 2.5ln [LT -1 ], with a standard deviation of 2ln[LT -1], the maximum and minimum correlation lengths are 300L and 200L respectively, the exponential variogram model with an angle of 135 degrees, and it is generated by the sequential simulation generator [Deutsch and Journel, 1998].

[0041] For this example, the initial contaminant concentration in the aquifer is 0 [ML -3 ], except for the water head at the western boundary, the initial water head of the aquifer is set to 200 [L], and other groundwater flow and pollutant transport parameters are set to the mean: porosity is set to 0.3 [-], longitudinal dispersion is set to 3 [L], and the ratio of transverse dispersion to longitudinal dispersion is set to 0.1. The total simulation time is set to 9125 [T] and refined into 25 time periods. The real pollutants are set at time 1095 [T] to 1000 [MT -1 ] is released at a mass loading rate of 1825[T] and the release stops at time 1825[T].

[0042] The number and distribution of pollutant concentration observation wells are as follows: Figure 1 shown.

[0043] The following groundwater unsteady flow equation is calculated using the groundwater flow simulator MODFLOW [McDonald and Harbaugh, 1988]:

[0044]

[0045] Where: S s is the water storage coefficient [L -1 ]; t is time [T]; represents the divergence operator; K is the permeability coefficient [LT -1 ]; is the Labla operator; H is the water level value [L]; q is the volume flow rate per unit volume of the aquifer [T -1 ]; W is the source and sink per unit volume [T -1 ].

[0046] The pollutant transport equation is calculated using the groundwater solute transport simulator MT3DMS [Zheng, 2010].

[0047]

[0048] Where: θ represents the effective porosity in the aquifer [-]; C represents the concentration [ML -3 ]; D m represents the molecular diffusion coefficient [L 2 T -1 ]; α is the diffusion tensor [L]; v is the velocity vector [LT -1 ];q sis the volume flow rate per unit volume of the aquifer [T -1 ]; C s is the source-sink flux concentration.

[0049] TIPT: Thermodynamic Integral Coupled Parallel Annealing Algorithm

[0050] SS: Subset Simulation

[0051] MC: Monte Carlo method

[0052] like Figure 2 The method of the present invention shown comprises the following steps:

[0053] S1. Based on the underground aquifer, the coefficients of the groundwater pollution failure probability formula are calculated by the Monte Carlo method, and the parameters that determine the smoothness of the indicator function and the weight of the Markov chain corresponding to each indicator function are selected according to the Gaussian integral criterion;

[0054] S2. Generate the initial field of groundwater aquifer permeability coefficient that obeys Gaussian distribution through sequential Gaussian, calculate the concentration of aquifer pollutants, find the response value according to the function function, and bring the response value into the indicator function to calculate the result;

[0055] S3, generate a suggested sample of the current aquifer hydraulic conductivity field according to the pre-conditioned Crank-Nicholson suggested distribution, calculate the aquifer pollutant concentration, find the response value according to the functional function, and substitute the response value into the indicator function to calculate the result;

[0056] S4, calculating the acceptance probability of the permeability field according to the result obtained by the indicator function in the Markov chain, and judging whether the sample in the chain accepts sampling;

[0057] S5, calculating the acceptance probability of the permeability field between the Markov chains according to the results obtained from the indicator function, and determining whether the samples between the chains are accepted for sampling;

[0058] S6. Repeat the above process to sample the Gaussian permeability coefficient field of the next aquifer until the sampling is terminated and the failure probability of groundwater contamination is calculated according to the thermodynamic integral formula.

[0059] When implementing:

[0060] According to the Gaussian integral criterion, different parameters σ that determine the smoothness of the indicator function are selected to generate different indicator functions I σ , and then obtain the temperature gradient in the parallel annealing algorithm and the weight w of each Markov chain k , and calculate the coefficient P of the groundwater pollution failure probability formula based on the indicator function and the functional function by the Monte Carlo method a ,

[0061]

[0062]

[0063]

[0064] Where N is the number of Monte Carlo sampling, I σ is the indicator function, I a To calculate P a The indicator function used, G(x) is the functional function, which calculates the difference between the pollutant concentration and the pollutant concentration standard through simulation, and the response value is y.

[0065] The failure probability of groundwater contamination is mainly determined by the uncertainty of the aquifer permeability field. The initial state time t = 0 is set, and the initial state x of the aquifer permeability field that obeys the Gaussian distribution is generated using sequential Gaussian. 0 , calculate the aquifer pollutant concentration, functional function response value and indicator function value

[0066] x 0 ~N(0,C)

[0067] Where C is the variance of Gaussian distribution.

[0068] The preconditioned Crank-Nicholson technique can effectively handle high-dimensional sampling problems. Therefore, according to its proposed distribution, the proposed sample x' of the aquifer permeability field is generated, and the aquifer pollutant concentration, the functional function response value and the indicator function value are calculated.

[0069]

[0070] In the formula, β is the jump factor and ε is the interference term.

[0071] According to the indicator function value obtained above, the sample acceptance probability α in the Markov chain is calculated, and a random number u∈[0,1] that obeys a uniform distribution is generated to determine whether to accept the sample.

[0072]

[0073]

[0074] For Markov chains of different temperatures, calculate the indirect acceptance probability α of the Markov chain s , use the same method as in the chain to determine whether to accept sampling

[0075]

[0076] Repeat the sampling until the termination condition is reached, and calculate the groundwater pollution failure probability value q of each Markov chain k , and calculate the failure probability P of groundwater pollution according to the thermodynamic integral formulaF ,

[0077]

[0078]

[0079]

[0080] In the formula, w k is the kth weight, q k is the parameter σ k The failure probability of the corresponding Markov chain.

[0081] like Figure 3 As shown in the figure, the mean failure probability of groundwater pollution calculated by thermodynamic integration coupled with parallel annealing method is closer to Monte Carlo method and converges faster than subset simulation method. Figure 4 As shown in the figure, the mean failure probability of 6000, 8000, and 10000 samples are plotted against the number of observation wells. It can be seen that as the number of observation wells increases, the mean failure probability is getting closer and closer to the result obtained by the Monte Carlo method, and compared with the subset simulation method, the result of this method is closer to the Monte Carlo method. Figure 5 As shown in the figure, the RMSE changes with the number of samples and the number of observation wells. Overall, the RMSE of the results obtained by this method is smaller than that of the subset simulation method. It can be seen that the thermodynamic integral coupled parallel annealing algorithm has obvious advantages in computing time efficiency and accuracy compared with traditional methods in calculating the failure probability of groundwater pollution.

[0082] Through the above analysis, the thermodynamic integral coupled parallel annealing algorithm can efficiently calculate the failure probability of groundwater pollution.

Claims

1. A method for efficiently calculating the failure probability of groundwater pollution based on thermodynamic integration, characterized in that: The following steps are involved: S1. Based on the underground aquifer, the coefficients of the groundwater pollution failure probability formula are calculated by the Monte Carlo method, and the parameters that determine the smoothness of the indicator function and the weight of the Markov chain corresponding to each indicator function are selected according to the Gaussian integral criterion; S2. Generate the initial field of groundwater aquifer permeability coefficient that obeys Gaussian distribution through sequential Gaussian, calculate the concentration of aquifer pollutants, find the response value according to the function function, and bring the response value into the indicator function to calculate the result; S3, generate a suggested sample of the current aquifer hydraulic conductivity field according to the pre-conditioned Crank-Nicholson suggested distribution, calculate the aquifer pollutant concentration, find the response value according to the functional function, and substitute the response value into the indicator function to calculate the result; S4, calculating the acceptance probability of the permeability field according to the result obtained by the indicator function in the Markov chain, and judging whether the sample in the chain accepts sampling; S5, calculating the acceptance probability of the permeability field between the Markov chains according to the results obtained from the indicator function, and determining whether the samples between the chains are accepted for sampling; S6, repeat the above process to sample the next aquifer Gaussian permeability coefficient field sample until the sampling is determined to be terminated, and calculate the groundwater contamination failure probability according to the thermodynamic integral formula; Among them, in S1, different parameters σ that determine the smoothness of the indicator function are selected according to the Gaussian integral criterion to generate different indicator functions I σ , and then obtain the temperature gradient in the parallel annealing algorithm and the weight w of each Markov chain k , and calculate the coefficient P of the groundwater pollution failure probability formula based on the indicator function and the functional function by the Monte Carlo method a , Where N is the number of Monte Carlo sampling, I σ is the indicator function, I a To calculate P a The indicator function used, G(x) is a functional function, which calculates the difference between the pollutant concentration and the pollutant concentration standard through simulation, and the response value is y; The failure probability of groundwater contamination in S2 is mainly determined by the uncertainty of the aquifer permeability field. The initial state time t = 0 is set, and the initial state x0 of the aquifer permeability field that obeys the Gaussian distribution is generated using sequential Gaussian to calculate the aquifer pollutant concentration, functional function response value, and indicator function value. x0~N(0,C) Where C is the variance of Gaussian distribution; The preconditioned Crank-Nicholson technique in S3 can effectively handle high-dimensional sampling problems. Therefore, according to its recommended distribution, the recommended sample x' of the aquifer permeability field is generated, and the aquifer pollutant concentration, functional function response value and indicator function value are calculated. In the formula, β is the jump factor and ε is the interference term; S4 calculates the sample acceptance probability α in the Markov chain based on the indicator function value obtained in S3, generates a random number u∈[0,1] that obeys a uniform distribution, and determines whether to accept the sample S5 For Markov chains of different temperatures, calculate the indirect acceptance probability α of the Markov chain s , use the same method as in the chain to determine whether to accept sampling S6 repeats sampling until the termination condition is reached, and calculates the groundwater pollution failure probability value q of each Markov chain k , and calculate the failure probability P of groundwater pollution according to the thermodynamic integral formula F , In the formula, w k is the kth weight, q k is the parameter σ k The failure probability of the corresponding Markov chain.

Images