Identification of aquifer pollution’s point sources with the reciprocity principle

The principle of reciprocity, called Maxwell–Betti theorem, initially used in mechanics in an elastic structure, establishes a relation of equality between two distinct strains under different loads. In this paper, we extend and apply this principle to flow and solute transport equations in porous media, in order to perform the pollution sources identification in aquifers. We developed general 2D expressions of the reciprocity principle for transient transport problems. This model leads to a linear equations set, with point sources coordinates, concentrations and associated water fluxes as unknowns The proposed model is then applied to the Rocky Mountain Arsenal aquifer (Konikow in Modeling Chloride Movement in the Alluvial Aquifer at the Rocky Mountain Arsenal, Colorado. Technical Report Water-Supply Paper 2044, USGS, 1979), where polluted water is injected into a well in the domain. The used inverse technique successfully recovered the position and the pollutant concentration in addition to the associated water flux. In addition, we developed and implemented the inverse method for different knowledge levels of the degrees of the aquifer contamination, i.e. more or less data available in the field. Multiple pollution point sources and noisy data situations are also developed and tested with high efficiency. The proposed method would be easy and useful to be implemented in the modeling software now widely used by researchers and groundwater managers. It can thus be applied in real case studies, to help authorities and regulators to efficiently identify the polluters and the contamination process, i.e. its location, onset, duration and the associated mass and water fluxes.

Nowadays, the industrial development and the intensification of the agricultural activities, introduce continuously in the environment new molecules, more or less carefully produced and used. This requires the utmost vigilance on behalf of the sanitary and water management authorities [1][2][3] . Indeed, a major part of the world population relies, to some extent, on groundwater for drinking, and for crops and food production. The protection of these resources should thus be of a great concern.
For many years, it was believed that the layers of soil and sediment above an aquifer act as a natural filter that retains pollutants and thus protects the groundwater. However, it has been widely recognized that the capacity of these soil layers to retain pollutants can be exceeded very quickly 4 . In addition, since an aquifer is polluted, it may become unusable for decades. The remediation of contaminated groundwater is inherently complex and expensive and can require long periods of time and sometimes centuries 5 .
Actually, decontamination of polluted groundwater is a huge challenge in Hydrogeology, in relation with the environmental and health requirements. A main challenge in any rehabilitation action is the evaluation of the degree of contamination. This includes the identification of unknown sources of contamination and the corresponding water and solute fluxes that led to the present state of pollution.
This process is important for both understanding the implementation of adequate remedial measures and for the identification of causes and responsibilities. In fact, the identification of the location and the level of pollution sources is crucial for the application of the polluter-pay principle adopted by the United Nation Conference of Rio de Janeiro Environment and Development declaration in its article 16 6 . In addition, the contaminant source identification could also be a means of dissuading potential infringements of laws for pollutant discharge and waste repository managements.
In this context, many works have considered this inverse problem in hydrogeology. Most of the existing studies concerns the recovering of the points-sources locations and/or the contaminant release histories 7,8 and/ or the number of these points-sources 9 .
Atmadja and Bagtzoglou 10 presented a review on mathematical methods that have been developed for the study of identifying sources of contamination. Authors classified these methods as deterministic or stochastic www.nature.com/scientificreports/ and associated to an optimization model 8,[11][12][13][14][15][16][17] . Other heuristic approaches based on genetic algorithm are proposed in Refs. 9,11 . In this paper, we introduce a new method based on the reciprocity principle, that allows the simultaneous identification of pollution point-sources locations and their pollutants concentrations, from the concentration's measurements in the aquifer domain. This principle, also known as the Maxwell-Betti theorem 18 , was first introduced in mechanics for linear problems 19 . It stipulates that for a linear elastic structure subject to two forces F and G, the work resulting from the application of the force F on the displacement field, yielded by the force G, is equal to the work resulting from the application of the force G on the displacement field yielded by the force F.
From a phenomenological perspective, this principle establishes strong relationships between different sets of forces and the consequent displacements applied to a given structure.
In mechanics, the reciprocity principle is usually applied to obtain displacements due to complex forces by using proxy problems with simpler forces that are more easily solved. Within the framework of groundwater flow scenario, reciprocity between two interference pumping tests was analyzed by Bruggeman 20 for Darcian flows in an unbounded, heterogeneous porous medium. Hariga et al. [21][22][23][24] have also applied the reciprocity principle in groundwater flows by using sources and boundary conditions as forcing terms and the resulting head field as a consequence.
In this research, the reciprocity principle is applied to the transport equation to recover the features of the pollutant point sources in aquifers. We show that this method allows to evaluate the position, solutes concentration and the associated injected water flux, for point solutes sources, from the knowledge of the cumulative mass flux through the boundary at any time in the considered interval and the concentrations on all the domain at the given time. The proposed method can then be considered as a "direct" one according to the Neuman classification 25 .
Obviously, in the real-world scenarios of contamination problems, since the concentrations are measured in a finite number of points, these should first be interpolated throughout the considered domain.
We show that the accuracy of the present identification method, depends on the number of available data. However, reasonable results are also obtained with few data which is a valuable insight to guide managers in the process of pollution sources identification.
The remainder of this paper is structured as follows. After a general review of the reciprocity principle ("General formulation and interpretation of the reciprocity principle" section), we derive it for the advection-diffusion equation with a constant injection during a given time ("The transport and flow equations in porous media" section). We then illustrate it in point-source pollution identification ("Illustration of the method for pollution sources identification" section) for four scenarii: with complete data, with few observations, with noisy data and with multiply point sources.

The reciprocity principle applied to the advection-diffusion equation
General formulation and interpretation of the reciprocity principle. Built on the Maxwell-Betti principle 18 , for the sake of simplicity we express the reciprocity principle in a general mathematical framework for linear elliptic problems. Let V be a Hilbert space associated to a domain Ω , a a bilinear form on V , assumed to be symmetric, continuous and coercive and l i a linear form defined on V for i = 1, 2 assumed to be continuous. Then we define the following variational problem: With (u i , ϕ) successively equal to (u 1 , u 2 ) and (u 2 , u 1 ) and using the symmetry of operator a , the reciprocity principle can be expressed by the identity: The reciprocity principle is fundamentally derived from the symmetry and the bi-linearity of the form a . From a physical perspective, it relates the responses to different external and internal forcing terms (source/sink terms, boundary conditions) of a given phenomenon on a fixed structure represented by form a . Finally, the reciprocity principle is similar to Green's second identity giving way to the boundary element methods.
The crucial point of the method is the relevant choice of the test functions. The test functions should ideally be closely related to the initial problem but should also lead to much simpler and, if possible, analytical solutions.

The transport and flow equations in porous media. The advection-diffusion equation governing
solute transport in a domain over a time interval [t0 , tf] is: with C(x, y, t), the solute concentration [ML −3 ] ; ω the aquifer porosity; D the hydrodynamic diffusion-dispersion tensor [L 2 T −1 ] ; V the Darcy velocity [LT −1 ] ; C D the prescribed concentration at the Dirichlet boundary Γ D ; Φ N the prescribed flux at the remaining Neuman boundary Γ N ; C 0 the initial concentration distribution and Q is the source term.
In this work, we consider the case of point-sources pollutant injected during a finite time, so that Q is expressed by the following equation: (2) l 1 (u 2 ) = l 2 (u 1 ).  −3 ] which starts at instant t i and stops at l i + t i ; so l i represents the time during which the pollutant is injected in the aquifer. δ(x − S i ) is the Dirac function which is non zero only in the point-source located at S i = (x i , y i ) and Π is a rectangular function defined as: with H the Heaviside step function. The transport problem (3) is coupled with the groundwater flow model via the Darcy velocity: Equation (6) is an expression of Darcy's law, integrated over the thickness of the aquifer to lead to the classical horizontal 2D representation of aquifers 26 . The term T represents the transmissivity and is equivalent to the integral of the hydraulic conductivity over the thickness of the aquifer under the Dupuit assumption.
In a stationary 2D case, the hydraulic head h is solution of the following problem: with N f the number of flow point-source (N f > N p ) ; T(x, y) the transmissivity field, Q sj the fluid volume flux rate at the j th point-source [T −1 ] , h D the prescribed heads at the flow Dirichlet boundary Γ ′ D , Q N the prescribed flux at the remaining flow Neuman boundary Γ ′ N . Note that systems (3) and (7) are coupled via the Darcy velocity (6) and that the pollutant point-sources form a sub-set of the pumping point-sources.
Note that the assumption of a steady state for a certain period of time is widely adopted in hydrogeology and hydrogeological modeling. It concerns the period during which we can consider that the aquifer is in a little disturbed regime, with inflows and outflows that balance each other. This often corresponds to the aquifer before its exploitation or with a still low exploitation level 26 . When building a hydrogeological model of a given aquifer, we search in the available database, the moment that corresponds to a significant increase or decrease in the exploitation and/or piezometer, apart from seasonal variations if any. This moment is then considered as the beginning of the transient regime. The inertia of the hydrogeological systems being such that this moment is expressed in year and the various fluxes entering/leaving the system are considered constant for the period of time from the infinite to this date. The transient regime then begins and the various fluxes, and more rarely the boundary conditions, are calculated for each period of time chosen according to the reactivity/inertia of the aquifer: often monthly or seasonal or even annual.
Let us recall that the considered inverse problem's unknowns are: C si the concentration of the pollutant released during the time interval [t i , t i + l i ] and the corresponding point source's position S i = (x i , y i ) . In the other hand we have hydraulic head's measurements and pollution concentration's measurements in some points of the domain Ω.
We hereafter establish the reciprocity expression for the transient transport equation in a generic 2D domain.
Reciprocity principle with the advection-diffusion equation. The reciprocity principle can be applied to the advection-diffusion equation using test functions ϕ that verify: Multiplying the first equations of systems (3) by φ and integrating it over Ω then applying two times Green's first identity, lead to the following equations: Then, using the fact that the test function verify equation (8), leads to: in Ω. where C f = C(x, y, t f ).
Note that for a given test function ϕ , the left hand side of Eq. (10) is known. Then we can determine the point-source concentrations ( C si ) and locations ( x i , y i ) by solving a linear system constitutes with Eq. (10) with different functions ϕ , as we will show it in the next example.
So in conclusion, the reciprocity principle relates the concentration and flux values at the boundary, the concentration values in the domain at any time with the pollution point-sources parameters.
However, since the wells water flow are also the pollutant point-source injection and as the flow equations are steady ones, it will be easier to apply the reciprocity principle to the system (7) as a first step of the identification procedure to recover the wells position (see [21][22][23] ). Then the Eq. (10) is exploited to recover the pollutant concentration. This identification process will be illustrated in the following examples.

Illustration of the method for pollution sources identification
First, we check the methodology with a single point-source pollution in the case of a hydrogeologic configuration inspired from the work by Konikow 27 who predicted long term pollutants dispersion in groundwater flow due to a leaky chemical pond under the Rocky Mountain arsenal in Colorado, US. The model setup is directly inspired from a test case in SUTRA code developed by Voss 28,29 .
In 1984 Voss considers the Rocky Mountain Arsenal to demonstrate some of the capabilities of SUTRA software modelling 28 . This example serves him to demonstrate the applicability of SUTRA to an areal constant density solute transport problem.
Here, in our study, we think that this case is perfectly suited to our objectives insofar as it includes production and injection wells, all the characteristics of which must be found by the inverse method developed in this work.
We considered the same types of boundary conditions as the original problem and the geometric and hydrodynamic data used are those used in the works of Hariga et al. 21,23 . The results of the inverse calculation are therefore compared to those of the direct calculation for the chosen dataset. The porous media is supposed with isotropic heterogeneous properties. The geometry, with boundary conditions data, is sketched on Fig. 1 The objective is to find the injected concentration of pollution C s and the position ( x s , y s ) of the point source,i.e. the pond, from the piezometric heads and the concentrations in the domain and the concentration flux over the boundaries at different times. As the problem has three unknowns, the reciprocity method requires only three virtual fields ϕ 1 , ϕ 2 and ϕ 3 . We choose three simple polynomial functions: Their derivatives in the direction of the normal to the boundary of the domain are given by: where n = (n x , n y ) the outward normal to the boundary.
Applied to these test functions, Eq. (10) and using the fact that: leads to the following equations: Ω ∂φ 1 ∂n = 0, ∂φ 2 ∂n = n x and ∂φ 3 ∂n = n y .  (11), we note that we directly obtain the injected concentration C s from the knowledge of the cumulative mass flux through the boundary at any time in the interval [ t 0 , t f ] and the concentration on all the domain Ω at the final time t f . Then we replace this value in Eq. (12) to obtain x p and in Eq. (13) in order to obtain y p . In Eqs. (12) and (13), we need to have the concentrations field at each time to identify the point-source position. However, as in this example, the pollution point-source is also a source point for the flow equation, we can apply the reciprocity principle to the stationnary problem (7) as done by Hariga et al in [21][22][23] to identify the point-source position. It's more easy and it consists in multiplying the first equation of system (7) by simple test   www.nature.com/scientificreports/ functions ( φ 1 = x and φ 2 = y ) and then using the Green formula twice times. So that at the end, we find the following expression for x p and y p : For the numerical study we define the relative errors as: for position's identification, where S the vector position and L 2 the euclidean norm. And: for fluxes' identification, where f is the injected pollution's concentration and L 1 the absolute value.
Identification with complete data over the domain. For the injected concentration recovering we use Eq. (11), whereas for the point-source position identification we compare the two methodologies: the transport one by the use of Eqs. (12) and (13) and the stationary flow one by applying Eq. (14). We consider t 1 = 25 days and l 1 = 5 days. Figures 3, 4 and 5 show, respectively, the hydraulic head, the concentration distribution and the total flux at the pond over the domain.
On Table 1 and on Fig. 6, we show the injected concentration recovered values for different period ( t f = 45, t f = 60, t f = 90 and t f = 180 days) and note that when the ratio t f t 1 increases the parameter recovering is hard.
On Table 2 we give the different point-source location identification. We note that the recovering with stationary flow reciprocity (Eq. (14)) is better than the one with reciprocity principle applied to the transport (14)  www.nature.com/scientificreports/ problem (Eqs. (12), (13)). This is due to the presence of time integration in the two last equations which lead to numerical errors.
Pollution source identification with few observation points. In this section, we perform the identification methodology with the sole knowledge of some 'measurements' for the hydraulic heads at the initial   Kriging is an interpolation process that produces an optimum, linear, and unbiased estimate of the property under examination based on the available data, with the least amount of error. The advantages of kriging over more traditional interpolation methods are that kriging integrates the spatial structure of the data in the form of a variogram model in its estimate process. Moreover, it is an exact interpolator because the surface created goes across the experimental points (unless a nugget effect is incorporated). This is why Kriging interpolation has been used in hydrogeology for many years, since the work of Delhomme 32 , to estimate hydraulic parameters throughout a complete domain, optimize recognitions, simulate interfaces, and so on.
With interpolated heads, the resulting point-source position remains close to their reference at around 10% (identified location is (2489.7 m, 4470.7 m)). For the injected concentration identification with kriged concentrations, the error values are shown in Table 3 and on Fig. 7, for different number of retained observations. We note that as expected the error increases when the number of observations decreases until it reaches 27% for only 20 observations. Noise sensibility pollution source identification. In actual field cases, errors may affect the collected data. There are many reasons for these errors, including the limited accuracy of equipment and the influence of sampling conditions. The quality and credibility of a given measurement depends on these errors, which must be properly identified and evaluated. So sensitivity to noise is performed on the last problem.  We have tested the effect of different levels of noise on the identification process, by adding an uniform white noise, with zero mean, to the pollution concentration measurement . The noise levels from 2 to 8% that we tested, are of the same orders of magnitude as those found in similar studies (between 5 and 10% 33 ). Table 4 shows the relative errors for different noise levels in the case with complete data. We note that the error remains acceptable (max.23%) until 8% of noise.
Multiply point-source pollution identification. In this case, we consider that in additional to the pond, pumping wells are also point-source pollution. We change the volumetric flux's sign for the wells and affect them the following pollution concentration: C 1 s = 0.5 kg/m 3 , C 2 s = 0.1 kg/m 3 and C 3 s = 0.5 kg/m 3 (exact positions and concentrations ares summarized on Table 5). So the studied inverse problem is to identify the four concentrations from the knowledge of the cumulative mass flux through the boundary at any time in the considered interval and the concentration on all the domain at the final time considered. We identify the positions concentrations by using the reciprocity principle respectively with the darcean equations and the advection-diffusion. As shown on Table 6, positions as well as concentrations are identified with a satisfying errors which don't exceed 8%.

Conclusion
In this study, we developed and applied the idea of reciprocity, which is taken from the field of mechanics and is generally applicable to linear problems, to the identification of pollution point sources in the advection-diffusion equation for the transport of solutes in aquifers.
This theory generates a strong link between the forcing terms and the resulting fields for two different forcing sets in the example of pollution transport in a Darcian flow, as illustrated in this study.
We have limited the scope of our work to the conservative transport of contaminants in the context of this research and for the sake of relevance and simplification. However, situations involving adsorption and/or degradation can be established. Furthermore, we believe that the method described here may be easily integrated into hydrogeological modeling codes such as SUTRA, MODFLOW, FEFLOW, and others, which are now widely used by water resource and environmental managers all over the world. The outcomes of the identification technique are satisfactory even with little information and available data on the state of aquifer contamination.
The simplicity of the subsequent identification process is indeed the main attractiveness of the reciprocity principle. Since the computations are limited to solving small linear equations where coefficients are given by some numerical evaluation of integrals, the computational cost of this direct method is quite minimal. The main drawback is the large amount of data required, which includes full information of concentration and fluxes at the boundaries and interfaces (which may not be available). However, as shown in Ref. 34 , using incomplete data sets is also still possible.
Other advanced parameter identification methods can be utilized in conjunction with these basic direct identification methods. They can serve as a preliminary assessment before embarking on more expensive field and laboratory investigations. These can also be effective persuasion techniques for decision-makers and competent authorities.