Evaluating pollution sources from limited concentration data of their contaminant plumes in groundwater systems is an inverse problem that presents computational challenges because of its ill-posedness. In this work, the Green element method (GEM) is used to solve inverse contaminant transport problems in 2-D orthotropic homogeneous aquifers. The GEM discretization of the differential equation produces an over-determined, ill-conditioned global matrix that is decomposed by the singular value decomposition method and solved by the least square method with Tikhonov regularization. Five test cases of pollution sources of unknown strengths are used to evaluate the performance of GEM. The results indicate that the current methodology is capable of correctly predicting the strength of pollution sources and the historical concentration plumes that they produce. However, the computation accuracy is influenced by the location of the observation points in relation to the source, data errors of observed concentrations, the transport mode, and the value of the Courant number used in the simulations. Data errors at observation points influence more significantly the GEM prediction of the pollution source strength than the concentration plume, while numerical artefacts of dispersion and oscillations observed in direct solutions of advection-dominant transport cases are also evident in inverse modelling.
INTRODUCTION
The quality of groundwater, which accounts for the water source and productive activities of a significant population of the world, can degrade due to the release of contaminants into aquifers from pollution sources. These sources include landfill sites, mine tailings, septic tanks, underground repositories of toxic and radioactive materials, unintentional spills of oils, hydrocarbons, and toxic substances on land surfaces, pesticides and fertilizers applied on agricultural lands, and seepages from polluted streams. Contamination may result from a single source or a combination of sources. In most instances, when contamination is detected, the locations and characteristics of the pollution sources are unknown. Therefore, a paramount step in addressing groundwater contamination is precise identification of pollution sources, their strength, spatial and temporal distributions, and this in essence constitutes an inverse groundwater contaminant transport problem.
There are various classes of inverse groundwater contaminant transport problems that have been presented in the literature. These include identifying the contaminant source location (Gorelick et al. 1983; Wagner 1992; Mahar & Datta 2000), estimating the number of possible pollution sources (Ayvaz 2010), estimating the release history or source strength (Skaggs & Kabala 1994), and recovering the historical spatial distribution of contaminant plumes (Michalak & Kitanidis 2004). In this work, we seek to simultaneously recover the signal of boundary pollution sources and the resulting spatial and temporal distributions of the plume.
Unlike in direct groundwater contaminant transport where the problem is well-posed with known boundary and initial conditions, inverse problems may seek to determine these input variables and their solutions may exhibit non-uniqueness, non-existence and instability (Sun 1996). Generally, for an inverse groundwater contaminant transport problem to be solvable additional information is provided by concentration measurements at some observation points.
Various numerical methods and their variants have previously been applied to inverse groundwater contaminant transport problems. Examples include the finite difference method (Mahar & Datta 2000; Neupauer & Wilson 2004; Ayvaz 2010; Prakash & Datta 2014), the finite element method (FEM) (Gorelick et al. 1983; Wagner 1992), the boundary element method (BEM) (Lesnic et al. 1998; Katsifarakis et al. 1999) and the adjoint method (Neupauer & Wilson 1999; Michalak & Kitanidis 2004). The methodology proposed in this study is based on the Green element method (GEM) which is founded on the singular integral theory of BEM, but implements the theory in an element-by-element manner as in FEM. By this approach, the resultant coefficient matrix is banded, the singular integral theory of BEM is amenable to solving nonlinear problems in heterogeneous domains, and point and distributed singularities are easier to handle because of their being aligned with the free-space Green's function which is in itself singular. GEM has successfully been previously applied to inverse modelling of heat conduction problems (Taigbenu 2015), and this is the first time it is being applied to inverse groundwater contaminant transport problems.
The GEM formulation followed in this work was recently developed (Taigbenu 2012), in which a second-order polynomial expression is used to approximate the internal normal fluxes so that only the solution for the concentration or the flux is calculated at external nodes and the concentration at internal nodes. Although this formulation is slightly less accurate than a flux-based one, its ability to generate fewer degrees of freedom per node than the latter makes it more attractive for inverse modelling. The resultant coefficient matrix from applying this formulation to the inverse contaminant transport is ill-conditioned, and is decomposed by the singular value decomposition (SVD) method and solved by the least square method with the Tikhonov regularization (Tikhonov & Arsenin 1977).
The objectives of this study therefore are: (1) to demonstrate the use of the GEM and a regularization approach to recover a contaminant's release history and historical distribution from limited observation data; (2) to investigate the influence of the Peclet and Courant numbers on source identification; and (3) to evaluate the influence of the location of observation points and measurement errors on source identification.
INVERSE DISPERSION–ADVECTION–REACTION EQUATION
RESULTS AND DISCUSSION
Test case 1
The influence of the location of the observation point in relation to the pollution source at x = 0 is assessed with this test example using the parameters: u = 1.0 km/yr, D = 0.025 km2/yr, L = 1.0 km, Δx = 0.025 km and Δt = 0.0125 yr. This corresponds to a marginal advection transport mode with local Peclet number Pe = uΔx/D = 1 and a Courant number Cr = uΔt/Δx = 0.5. Forty uniform rectangular elements are used in the GEM simulation, and each simulation is run until t = 0.5 yr with two observation points (γ1 = γ2) that are chosen at nodes with the same x-coordinate and located at the top and bottom boundaries of the computational domain.
Mean error of estimated source strength in relation to observation point location for Test case 1.
Mean error of estimated source strength in relation to observation point location for Test case 1.
GEM prediction of the pollution source strength for various locations of observation points.
GEM prediction of the pollution source strength for various locations of observation points.
The influence of two transport modes on the source identification is assessed with this test case. The first transport mode (i) represents an advection dominant transport with Pe = 50, and the second example (ii) is one of marginal advection transport with Pe = 1. The simulation parameters for the two transport cases are given in Table 1, while the fully implicit scheme with θ = 1 is used in all the simulations. The influence of the time step on the numerical solutions is assessed for three values of the Courant number, Cr=uΔt/Δx, namely Cr= 0.1, 0.5 and 1.0. Additional concentration data are provided at two observation points γ1 and γ2 that are located along x= 0.025 km for both transport modes.
Simulation parameters for Test case 1
Transport case . | Pe . | D (km2/yr) . | u (km/yr) . | Δx (km) . | L (km) . |
---|---|---|---|---|---|
Advection dominant (i) | 50.0 | 0.0005 | 1.0 | 0.025 | 1.0 |
Marginal advection (ii) | 1.0 | 0.025 | 1.0 | 0.025 | 1.0 |
Transport case . | Pe . | D (km2/yr) . | u (km/yr) . | Δx (km) . | L (km) . |
---|---|---|---|---|---|
Advection dominant (i) | 50.0 | 0.0005 | 1.0 | 0.025 | 1.0 |
Marginal advection (ii) | 1.0 | 0.025 | 1.0 | 0.025 | 1.0 |
Exact and inverse GEM solutions for the pollution source strength Cs for Pe= 50.
Exact and inverse GEM solutions for the pollution source strength Cs for Pe= 50.
Exact and inverse GEM solutions for the pollution source strength Cs for Pe= 1.
Source strength for Pe= 1 and Cr= 1.0 of Test case 1 for various noise levels in the observation data.
Source strength for Pe= 1 and Cr= 1.0 of Test case 1 for various noise levels in the observation data.
Test case 2
Test case 3
Mean error of estimated source strength in relation to location of observation points for Test case 3.
Mean error of estimated source strength in relation to location of observation points for Test case 3.
Exact and GEM calculated source strength along x= 0 for Test case 3.
Exact and GEM concentration (g/L) contours of Test case 3 – solid line is exact and dotted is GEM with no noise in the data.
Exact and GEM concentration (g/L) contours of Test case 3 – solid line is exact and dotted is GEM with no noise in the data.
Error plots of solutions for Test case 3 against regularization parameter.
Test case 4
It should be noted that the release history is unrestricted in time in contrast to that used by Huang et al. (2008) and Li & Mao (2011) who both used global space-time optimization techniques in their inverse modelling of this example which places a limitation on the time dimension to enable the size of the matrix that is solved to be manageable. GEM inverse solution to this example does not suffer from this limitation in time, as it allows for the numerical scheme to march forward to any time level. The direct solution procedure uses the expression for q(x,t) that is given by Equation (33) to compute the concentrations at the observation points (P = 17) that are located at a depth of 2.5 m from the top of the aquifer over a simulation period of 1 day. The inverse model then considers these concentration values at the observation points as the actual ones to calculate the release history of the pollution source.
Error plots of the inverse GEM solutions for the pollution source of Test case 4 for various noise levels of observation data.
Error plots of the inverse GEM solutions for the pollution source of Test case 4 for various noise levels of observation data.
Error plots of the inverse GEM solutions for the spatial concentration distribution of Test case 4 for various noise levels of observation data.
Error plots of the inverse GEM solutions for the spatial concentration distribution of Test case 4 for various noise levels of observation data.
Test case 5
Release history of pollution source q(x,t) g/L−1 day−1 m−1 of Test case 5 by GEM with σ= 0%.
Release history of pollution source q(x,t) g/L−1 day−1 m−1 of Test case 5 by GEM with σ= 0%.
Error plots of the inverse GEM solutions for the pollution flux q(x,t) of Test case 5 for various noise levels of observation data.
Error plots of the inverse GEM solutions for the pollution flux q(x,t) of Test case 5 for various noise levels of observation data.
Error plots of the inverse GEM solutions for the spatial concentration distribution of Test case 5 for various noise levels of observation data.
Error plots of the inverse GEM solutions for the spatial concentration distribution of Test case 5 for various noise levels of observation data.
CONCLUSION
This paper has presented the Green element solutions of the inverse dispersion–advection contaminant transport equation in 2-D orthotropic homogeneous aquifers. The GEM discretization of the differential equation gives rise to an over-determined, ill-conditioned global matrix that is solved by the least square method with Tikhonov regularization, and aided by the SVD of the matrix. Using five examples, the influences on the GEM simulation results are evaluated for various transport modes, values of the Courant number, location of observation points in relation to the pollution source, and errors in observation data. We found the following:
(1) The prediction capability of the pollution source strength by GEM is enhanced by having observation points located in close proximity to the source. This confirms the results of previous numerical investigations and what should be expected intuitively, that the integrity of pollution source recovery depends on that of measured concentrations at observation points which in part is dependent on how close those points are to the source.
(2) The GEM prediction of the source strength is generally better with more dispersion dominant transport than with the advection dominant one; for transient transport, oscillations are observed in the numerical solution at early simulation times but they die out at later times. The challenge in the GEM modelling of the advection-dominant transport problem lies with the inability of the steep concentration profiles associated with such problems to be correctly approximated by the linear interpolation that is used in the current formulation. The use of higher interpolation functions like quadratic and Hermitian functions provides better approximation of the steep concentration fronts and, hence, improved solutions (Taigbenu 1998).
(3) The two-fold influence of the Courant number Cr on the GEM solutions, namely numerical stability and frequency of available sampling data, makes it difficult to determine which influence is dominant in any particular transport case. However, small values of Cr in the GEM simulations generally produce oscillatory solutions at early simulation times when concentration information is not available at observation points, but give accurate solutions at later times when they are available frequently at those points. The oscillations are more pronounced in the GEM solutions for advection-dominant transport than for dispersion-dominant transport.
(4) For the simulated noise levels in the data at observation points, their influence on the recovery of the source strength is significant and, as expected, the influence on the inverse solutions increases with increased noise level in the data.
(5) Errors in observation measurements have greater influence on the recovery of the pollution source strength than on the historical concentration distribution.
(6) GEM is capable of identifying the active period of pollution sources as well as their strengths when they are active.
The simulations of all the test cases were carried out on a single CPU desktop 64-bit Intel® Core™ i7-2600, 3.40 GHz, 8 GB memory computer, and the simulation times with respect to spatial discretization element and time step ranged from 1.63 × 10−4 seconds for the simple 1-D example to 4.3 × 10−4 seconds for the more challenging 2-D example. These small simulation times are indicative that the proposed methodology is computationally efficient.