CN114088908B - Method for smoothly depicting groundwater simple non-point source information based on multiple data assimilation sets - Google Patents

Abstract

The invention discloses a method for smoothly depicting underground water simple surface source information based on multiple data assimilation sets, which comprises the following steps of: determining the range of a suspicious contaminated area, setting initial values of various information parameters of an initial elliptical surface source and the number of times of assimilation of observation data; assimilating observation data, and updating each surface source parameter by combining a pollutant solute transport model and utilizing a multi-data assimilation set smoothing method; and circulating the previous process until the set assimilation times. The invention can approximately depict the position, the shape, the pollutant release time, the release duration and the release concentration of a simple surface source through the rotating ellipse; only the approximate value range of each parameter of the surface source is known, and each information parameter of the pollution source can be identified by assimilating the observation data of the state variable; the feasibility of the multi-data homogenization set smoothing in the aspect of identifying the non-point source is verified.

Description

Method for smoothly depicting groundwater simple non-point source information based on multiple data assimilation sets

Technical Field

The invention belongs to the technical field of hydrology statistics, and particularly relates to a method for smoothly depicting simple non-point source information of underground water based on multiple data homogenization sets.

Background

Groundwater is an important component of water resources and is an important source of fresh water for human life, agriculture and industry, and the quality of groundwater is closely linked with human health. Once groundwater is contaminated, it is necessary to know the source of the contamination, such as where and for how long the contamination has been released. With the advance of urbanization process in China, water resource development is unreasonable, groundwater is seriously polluted, and problems of tracing and remedying pollution sources and the like need to be solved urgently. Reliable source identification is critical to understanding groundwater contamination and to improving remediation efficiency for groundwater management.

The reverse modeling method for identifying the underground water pollution source is widely concerned due to effectiveness and high efficiency, and can be divided into three categories according to the characteristics: an optimization method, a probabilistic method, and a deterministic method. Early methods for identifying sources of contamination were optimization methods, primarily characterized by minimizing the difference between the simulated concentration and the observed values; the probabilistic approach is to maximize the posterior probability of some source parameter for a given observation; the deterministic method is to solve the advection-diffusion equation over time. But these are difficult to accurately characterize a number of unknown information parameters. And the method has outstanding performance in identifying static parameters and pollution source information parameters, wherein the position and release information of a pollution source point can be identified simultaneously by the smooth multi-data assimilation set, and the method is proved in synthesis and practical cases.

From a scale point of view, groundwater sources of pollution can be simply divided into two categories: point sources and area sources. Point sources are typically caused by landfills, gas stations, factory effluents and municipal sewage, while non-point sources are typically caused by agricultural fertilizers, livestock, poultry farming manure effluents and chemical leaks from chemical plants. At present, most pollution source identification researches only relate to identifying source position and release information by a single-point source or multiple scattered point sources, or a point source is adopted to simplify a surface source, but the quality of pollutant release is underestimated by directly simplifying the surface source to the central point. Since the complexity of the contamination source identification increases as the number of surface source emissions increases, there is little work on the surface source identification. The few face source identification works are that the face source area is processed into a normalized and angular form by using a simple graph consisting of rectangles or rectangles.

Disclosure of Invention

The purpose of the invention is as follows: in order to overcome the defects of the background technology, the invention discloses a method for smoothly depicting underground water simple surface source information based on a multi-data assimilation set.

The technical scheme is as follows: the invention discloses a method for smoothly depicting simple underground water non-point source information based on multi-data assimilation set, which comprises the following steps:

s1, generating an initial logarithmic permeability coefficient field by utilizing sequential Gaussian simulation according to hydrogeological statistical information; determining a suspicious region of a pollution source according to early-stage field investigation, determining a value range of information parameters (such as the position of an ellipse central point, the length of a long semi-axis and a short semi-axis, a rotation angle, initial release time of the pollutant, the release duration and the release amount) for approximately depicting a simple area source of the pollutant, and randomly sampling in the value range to generate an initial value of each parameter; setting the times of multidata assimilation;

s2, assimilating the jth data, and predicting a global pollutant concentration field through a solute transport equation based on updated (initial) surface source information;

s3, updating each information parameter of the surface source by utilizing an updating process in a multi-data assimilation set smoothing technology based on the difference between the pollutant concentration observation data and the predicted concentration data in all assimilation time periods at the observation well;

and S4, repeating the assimilation process of S2-S3 until the set assimilation times are reached.

Wherein S1 specifically comprises:

generating an initial logarithmic permeability coefficient field lnK, an initial pollution source central point position Xs (x coordinate axis direction) and Ys (y coordinate axis direction) approximately depicted by a rotating ellipse, long and short semi-axis lengths Ra and Rb, an ellipse anticlockwise rotating angle B around a central point, initial release time Ti, release time delta T and release mass load rate M, and constructing a parameter vector S:

further, S2 specifically is:

based on the initial pollutant concentration field (state variable) C ₀ And (j-1) the surface source information parameter vectors of data assimilation iteration update

State transfer equation using solute transport

Predicting a j-th iteration global pollutant concentration distribution field

Further, S3 specifically is:

model parameter augmentation vector based on j-1 th data assimilation iteration update

Kalman gain K _j And a concentration observed value C at the observation well ^o And the predicted concentration

Difference between and observation error epsilon _j Updating the face source parameter vector after the j-th data assimilation

Kalman gain K in the update process _j Is the cross covariance D between the surface source related parameter at all times and the predicted concentration at the observation point _SC,j And autocovariance D between predicted concentrations of observation points obtained at all times _CC,j And the covariance of the observation error R _j As a function of (a) or (b),

wherein:

K _j ＝D _SC,j (D _CC,j +a _j R _j ) ^-1 。

from the above formula, the covariance of the observation error R _j And error of observationDifference epsilon _j Respectively formed by a _j And

and (4) amplifying.

Has the advantages that: compared with the prior art, the invention has the advantages that: firstly, the position, the shape, the pollutant release time, the release duration and the release concentration of a simple surface source can be approximately described through a rotating ellipse; secondly, the method can identify each information parameter of the pollution source by assimilating the observation data of the state variable only by knowing the approximate value range of each parameter of the non-point source; finally, the method verifies the feasibility of the multi-data homogenization set smoothing in the aspect of identifying the non-point source.

Drawings

FIG. 1 is a schematic diagram of a true simple surface source and an observation well position under three scenarios in an embodiment;

FIG. 2 is a flow chart of the method of the present invention;

FIG. 3 is a graph of average absolute error (AAB) and collective divergence (ESp) of information of a non-observation data assimilation and updating of area source after different times of data assimilation in three scenarios in the example;

FIG. 4 is a box plot of the initial and updated sets of simple area source information parameters computed for three scenarios in the example.

Detailed Description

The technical solution of the present invention is further described below with reference to the accompanying drawings and examples.

The invention is verified in a saturated steady-state confined aquifer, and the method can smoothly and effectively invert the simple area source information of the underground water based on multiple data assimilation sets. Setting a two-dimensional aquifer space: 800[ L ]. 800[ L ], thickness of 80[ L ], discrete as 80[ 80 ] 1 unit cells, each unit cell being 10[ L ]. 10[ L ], model east-west boundary being constant water head boundary and given water heads being 80[ L ] and 300[ L ], south-north boundary being water-proof boundary.

As shown in figure 1, the position and the shape of a real surface source, the position distribution of 30 observation wells and 2 verification wells are given, wherein a triangle represents an observation well, a square represents a verification well and an irregular areaThe fields represent different shapes of the area source. Through a sequential multi-Gaussian simulator GCOSIM3D [ G Lou mez-Hern a ndez and journal, 1993]Generating a mean value of-2ln ^-1 ]And a standard deviation of 1ln [ lt ], [ solution of [ 1 ] n ^-1 ]The maximum correlation length and the minimum correlation length are respectively 300[ L ]]And 200[ 2 ], [ L ]]A log hydraulic permeability coefficient reference field of a high-dimensional gaussian distribution of spherical variation function with an angle of 135 degrees.

The initial pollutant concentration of the aquifer is 0[ ML ] ^-3 ]And setting other groundwater flow and solute transport parameters as mean values: porosity is set to 0.3[ -]The longitudinal dispersion degree is set to be 3.0[ L ]]The ratio of the lateral and longitudinal dispersivity was set to 0.5. It is assumed that only advection and dispersion processes exist for solute transport of contaminants, and the entire solute transport process is unstable. The real area source setting is as shown in fig. 1, and is divided into three scenarios (oval, circular and simple irregular graphs), and the respective position and shape and information parameters such as mass load rate and release time released by each grid in the area of the area source are shown in table 1. For the present embodiment, the simulation total time period is set to 10950[ T ]]And equally dividing the time interval into 100 time steps, the flour source is from 985.5[ T ]]The release is started until 2299.5[ 2 ] T]The release is stopped (i.e., the release period is 1314T]) And recording the concentration change in the observation well in the first 50 time steps for inversion of the surface source parameters.

Calculating a pollutant migration equation through a groundwater solute migration model MT3DMS [ Zheng,2010 ]:

in the formula: theta represents the effective porosity [ -](ii) a C represents the concentration [ ML ] ^-3 ](ii) a T is time [ T]；

Represents the divergence operator; d _m Represents the molecular diffusion coefficient [ L ² T ^-1 ](ii) a α is the dispersion tensor [ L](ii) a v is the flow velocity vector [ LT ^-1 ]Passing through the head H

Associating;

is Laplace operator; q. q.s _s Volumetric flow rate per unit volume of aquifer ^-1 ]；C _s Is the source sink flux concentration.

The head H at steady state conditions is obtained by calculating the groundwater flow steady state equation using a groundwater flow model MODFLOW [ McDonald and Harbaugh,1988 ]:

in the formula: k is the hydraulic permeability coefficient [ LT ^-1 ](ii) a W is the source and sink quantity LT of the aquifer per unit volume ^-1 ]。

For the present example, as shown in table 1, the number of iterations for smoothly assimilating state variable observation data using multiple data assimilation sets is assumed to be 1 time, 2 times, 4 times, and 6 times.

Table 1 defines parameters of real simple area sources under three scenarios

Context	S1	S2	S3
				Number of assimilation iterations	1，2，4，6	1，2，4，6	1，2，4，6
Shape of pollution source	Oval shape	Circular shape	Irregular figure
				Xs
	150	150	/
				Ys	540	540	/
Ra	80	60	/
				Rb	40	60	/
B	30	/	/
				Ti	985.5	985.5	985.5
ΔT	2299.5	2299.5	2299.5
				M	1000	1000	1000

For the present embodiment, as shown in table 2, 500 sample sets are randomly selected from the uniform distribution of the surface source central point position, the long and short semi-axis length, the rotation angle, the initial release time, the release duration, and the release concentration in the suspicious range approximately depicted by the rotation ellipse.

Table 2 shows the suspicious ranges of the parameters of the area source under the three scenarios of the embodiment

Parameter(s)	Suspicious region
		Xs	110-210
Ys	460-560
		Ra	40-140
Rb	10-80
		B	0-90
Ti	0-3175.5
		ΔT	1204.5-6679.5
M	950-1200

As shown in fig. 2, the method for smoothly depicting groundwater simple non-point source information based on multiple data assimilation sets of the invention comprises the following steps:

s1, generating an initial logarithmic permeability coefficient field by utilizing sequential Gaussian simulation according to hydrogeological statistical information; determining a suspicious area of a pollution source according to early-stage field investigation, determining a value range of information parameters (such as the position of an ellipse central point, the length of a long semi-axis and a short semi-axis, a rotation angle, initial release time of the pollutant, release duration and release amount) for approximately depicting a simple area source of the pollutant, and randomly sampling in the value range to generate an initial value of each parameter; the number of times of multidata assimilation is set.

S2, assimilating the jth data, and predicting a global pollutant concentration field through a solute transport equation based on updated (initial) surface source information.

And S3, updating each information parameter of the surface source by utilizing an updating process in a multi-data assimilation set smoothing technology based on the difference between the pollutant concentration observation data and the predicted concentration data in all assimilation time periods at the observation well.

And S4, repeating the assimilation process of S2-S3 until the set assimilation times are reached.

In the specific implementation:

(1) Generating an initial logarithmic permeability coefficient field lnK, an initial pollution source central point position Xs (x coordinate axis direction) and Ys (y coordinate axis direction) approximately depicted by a rotating ellipse, long and short semi-axis lengths Ra and Rb, an ellipse anticlockwise rotating angle B around a central point, initial release time Ti, release time delta T and release mass load rate M, and constructing a parameter vector S:

(2) Based on the initial pollutant concentration field (state variable) C ₀ And (j-1) each information parameter vector of face source updated by data assimilation iteration

State transition equation using solute transport

Predicting a j-th iteration global pollutant concentration distribution field

(3) Model parameter augmentation vector based on j-1 th data assimilation iteration update

Kalman gain K _j And a concentration observed value C at the observation well ^o And the predicted concentration

Difference between and observation error e _j Updating the face source parameter vector after the j-th data assimilation

Kalman gain K in the update process _j Is the cross covariance D between the surface source related parameter at all times and the predicted concentration at the observation point _SC,j And autocovariance D between predicted concentrations of observation points obtained at all times _CC,j And the covariance of the observation error R _j Is a function of (a) a function of (b),

wherein

K _j ＝D _SC,j (D _CC,j +a _j R _j ) ^-1

From the above formula, the covariance of the observation error R _j And the observation error ε _j Respectively formed by a _j And

and (4) amplifying.

(4) The above processes are repeated until the set data assimilation times are reached.

For the method of the present invention, as shown in fig. 3, an average absolute error (AAB) between each parameter value of the area source and the true value updated in three scenarios without observed data assimilation and after different data assimilation times and an aggregate divergence (ESp) diagram for measuring aggregate accuracy achieved by updating through aggregate variance square root in the embodiment are shown, and it can be seen that, when data assimilation is not performed, AAB and ESp of each parameter of the area source in three scenarios are both at a relatively high value, and as the number of iterations increases, AAB and ESp gradually decrease, and all scenarios except AAB of scenario S3 after 6 iterations tend to converge to a smaller value. This result demonstrates that as the number of iterations increases, the parameters of the area source get closer to the true values and the uncertainty also decreases significantly.

As shown in fig. 4, the box line plot calculated by the initial and updated sets of source information parameters under three scenarios in the embodiment. In the figure, the solid line in the rectangular box represents the median number, the upper and lower sides of the rectangular box represent the upper and lower quartiles respectively, and the black horizontal dotted line represents the true value of each parameter of the area source. As can be seen from the boxplot, compared with the huge uncertainty of each parameter value of the initial surface source, when the number of assimilation iterations increases, the uncertainty of all updated source parameter values decreases sharply, the median value after 4 iterations almost coincides with the true value, and only the value of the mass load rate M of the scenario S3 is smaller than the true value. The reason why M of scenario S3 is underestimated is that there are more nodes that are updated to approximate the elliptical coverage of the irregular source region than the true irregular source region, resulting in M being diluted. For the scenario S2, the rotation angle B of the final iteration of the circular surface source has no influence on the graphic depiction; for scenario S3, when the true contamination source is irregular in shape, there are no true source shape parameter values (Xs, ys, ra, rb, and B), but the updated elliptical shape parameter values may converge to a certain value quickly, and an ellipse described according to these parameters may approximately describe the irregular surface source. Therefore, the ellipse shape and the related pollution source information (initial release time, release duration and pollutant release amount) which are described by the ellipse parameters updated after the ES-MDA assimilation pollutant concentration observation data are utilized can not only effectively describe regular real area source (such as ellipse or circle) information, but also approximately describe real irregular simple area source information by using the ellipse.

Through the analysis, the method can effectively depict the information of each parameter of the simple surface source in the aquifer after assimilating the observation data of the iteration for enough times by using the smoothing technology based on the multi-data assimilation set.

Claims (1)

1. A method for smoothly depicting underground water simple surface source information based on multiple data assimilation sets is characterized by comprising the following steps:

s1, generating an initial logarithmic permeability coefficient field by utilizing sequential Gaussian simulation according to hydrogeological statistical information;

according to the field investigation of earlier stage confirm the suspicious area of pollution source and confirm that is used for approximately depicting the information parameter value range of simple area source of pollutant, include: the positions of the central points of the ellipses, the lengths of the long and short semi-axes, the rotation angle, the initial release time of pollutants, the release duration and the release amount are randomly sampled in the range to generate the initial values of all the parameters;

setting the times of multidata assimilation;

s2, assimilating the jth data, and predicting a global pollutant concentration field through a solute transport equation based on updated surface source information;

s3, updating each information parameter of the surface source by utilizing an updating process in a multi-data assimilation set smoothing technology based on the difference between pollutant concentration observation data and predicted concentration data in all assimilation time periods at an observation well;

s4, repeating the assimilation process of S2-S3 until the set assimilation times are reached;

s1 specifically comprises the following steps:

generating an initial logarithmic permeability coefficient field lnK, initial pollution source central point coordinate positions Xs and Ys approximately depicted by a rotation ellipse, long and short semi-axis lengths Ra and Rb, an anticlockwise rotation angle B of the ellipse around a central point, initial release time Ti, release time delta T and release mass load rate M, and constructing a parameter vector S:

s2 specifically comprises the following steps:

based on the initial contaminant concentration field C ₀ And j-1 times data assimilation iteration updated surface source parameter vectors

State transfer equation using solute transport

Predicting a j-th iteration global pollutant concentration distribution field

S3 specifically comprises the following steps:

model parameter augmentation vector based on j-1 th data assimilation iteration update

Kalman gain K _j Concentration observed value C at observation well ^o And predicted concentration

Difference between and observation error e _j Updating the face source parameter vector after the j-th data assimilation

Kalman gain K in the update process _j Is the cross covariance D between the surface source related parameter at all times and the predicted concentration at the observation point _SC,j And autocovariance D between predicted concentrations of observation points obtained at all times _CC,j And the covariance of the observation error R _j As a function of (a) or (b),

wherein:

K _j ＝D _SC,j (D _CC,j +a _j R _j ) ^-1

observation error covariance R _j And the observation error ε _j Respectively formed by a _j And

and (4) amplifying.

Images