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

*C*is the contaminant concentration [ML

^{−3}],

**V**=

**i**

*u*+

**j**

*v*is the velocity vector of the transporting medium [LT

^{−1}],

*D*and

_{x}*D*are the hydrodynamic dispersion coefficients in the spatial dimensions [L

_{y}^{2}T

^{−1}],

*R*is the retardation factor,

*μ*is the first-order decay rate coefficient that captures that reaction process [T

^{−1}],

*t*is the time dimension [T] and

*x*and y are the Cartesian coordinates [L]. In forward modelling where the concentration distribution in the domain is sought, Equation (1) applies to a domain Ω over which the initial and boundary conditions are given as: where

*n*and

_{x}*n*are the direction cosines of the outward pointing normal

_{y}**n**to the surface Γ

_{2}in the

*x*and

*y*directions. When the conditions in Equations (2) and (3) are specified with known medium and flow variables, the solution is that of the direct problem which is tractable and can be solved by conventional numerical methods. However, the inverse problem addressed in this paper has an additional boundary segment Γ

_{3}where neither the concentration nor the flux is specified, and it is expected that the release history of the pollutant from this boundary be evaluated, making use of the concentration plume that it produces in the aquifer. In essence, the domain Ω has a boundary Γ that comprises three parts, Γ

_{1}, Γ

_{2}and Γ

_{3}. To determine the release history of the unknown boundary information, additional data of the concentration are required at observation points

*γ*= (

_{b}*x*,

_{b}*y*) within the aquifer, and these are represented as: Essentially the inverse solution of Equation (1) is required in the domain Ω with the boundary Γ = Γ

_{b}_{1}Γ

_{2}Γ

_{3}(see Figure 1).

*P*. In many real life circumstances, these observed concentration measurements are tainted with measurement errors, and it is expected that the influence of these errors on the numerical solutions is assessed. The actual measured concentrations can be expressed as: where

*σ*is the error or noise level on the observed concentration C(

*γ*

_{b}) and

*RAN*are random numbers that are generated from a standard normal distribution.

*et al.*1984): where

*r*= (

_{i}*x*,

_{i}*y*) is the source or collocation point,

_{i}*λ*= nodal angle at

*r*,

_{i}*q*is the normal contaminant flux which has an expression similar to Equation (4). Equation (7) consists of integrals along the boundary Γ with the integration variable

*s*and integrals over the domain Ω with integration variable

*A*. The first three terms in Equation (7) represent the dispersion transport or the terms on the left-hand side of Equation (1), while the terms in the domain integral represent the temporal variation, advection and reaction transport or the right-hand side of Equation (1). The variable G in Equation (7) is the fundamental solution of: in the infinite space that has the expression: and

*G** is the normal derivative of

*G*that is expressed as:

*C*,

*q*, and

**V**by basis functions of the Lagrange family, that is

*C*≈

*N*(

_{j}*r*)

*C*(

_{j}*t*) (where

*N*(

_{j}*r*) are linear interpolating functions in space and

*j*denotes the nodes of the element). The discrete element equation for each element Ω

*with boundary Γ*

^{e}*is: in which the indices*

^{e}*i*,

*j*and

*k*in Equation (11) denote nodes of the element Ω

*, while the elemental matrices are given as:*

^{e}*S*representing the first two terms,

_{ij}*L*the third term,

_{ij}*W*the first two terms in the domain integrals, and

_{ij}*U*and

_{ijk}*Y*account, respectively, for the advection in the

_{ijk}*x*and

*y*directions. The elemental integrations over rectangular elements in Equation (12) are done analytically and their expressions are presented for the first time in the Appendix (available with the online version of this paper). The discrete element equations are then aggregated over all elements used to discretize the computational domain, and this gives a matrix equation of the form: in which . At internal boundaries between the elements, the normal flux

*q*is approximated by a second order polynomial function that is expressed in terms of the concentration as described in Taigbenu (2012). This results in two variables at the boundary nodes (

*C*and

*q*) and only the concentration

*C*at the interior nodes. The temporal derivative term is approximated by a finite difference expression:

*dC/dt*≈ [

*C*

^{(2)}–

*C*

^{(1)}]/Δ

*t*that is evaluated at

*t*=

*t*

_{1}+

*θ*Δ

*t*, where 0 ≤

*θ*≤ 1, is the difference weighting factor, and Δ

*t*is the time step between the current time

*t*

_{2}and the previous one

*t*

_{1}. Introducing the approximation for the temporal derivative into Equation (13) and weighting the other terms in the equation by

*θ*, gives: where

*β*= 1−

*θ*and the bracketed superscripts 1 and 2 represent the previous and current times, respectively. Incorporating the known initial, boundary data and concentration data at observation points into Equation (14) yields the matrix equation: where in which

**A**is an

*M*×

*N*matrix, with

*M*being the number of nodes in the computational domain which is the same as the number of discrete equations generated and

*N*is the number of unknowns which are

*C*and/or

*q*at external nodes, and

*C*at internal nodes. The vector

**w**is

*N*× 1 in dimension and represents the unknown quantities that should be calculated, and

**b**is an

*M*× 1 vector of known quantities which consists of the terms on the right-hand side of Equation (14), as well as the contributions from the available concentration measurements that are given by Equation (5). For inverse problems, the matrix Equation (15) is over-determined and ill-conditioned. As an over-determined equation, a least square solution is sought, while

**A**is decomposed by the SVD technique (Heath 1997): where the superscript

*T*denotes the transpose,

**H**is an

*M*×

*M*orthogonal matrix,

**P**is an

*N*×

*N*orthogonal matrix, and

**is an**

*D**M*×

*N*diagonal matrix with

*N*non-negative diagonal elements

*φ*

_{1},

*φ*

_{2},…,

*φ*which satisfy the condition:

_{N}*φ*

_{1}>

*φ*

_{2}> …>

*φ*> 0. The least square solution of Equation (15) requires the minimization of the Euclidean norm ‖

_{N}**Aw**–

**b**‖

^{2}which gives: The small singular values

*φ*cause instability in the solution of

_{i}**w**, and this can be ameliorated by using the zero order Tikhonov regularization technique which requires the minimization of the Euclidean norm ‖

**Aw**–

**b**‖

^{2}+

*α*

^{2}‖

**w**‖

^{2}that yields the solution for

**w:**The factor serves to dampen the contribution of the small singular values. The choice of the regularization parameter,

*α*, is crucial to obtaining meaningful solutions. Various methods have been proposed for selecting the regularization parameter and they include the L-curve (Hansen 1992), generalized cross validation (Golub

*et al.*1979) and the discrepancy principle (Goncharskii

*et al.*1972). In this work, we utilize the L-curve method. It is a continuous curve consisting of all the points for . It shows the compromise between the perturbation error and the regularization error (over-smoothing). A suitable regularization is the one corresponding to a regularized solution near the corner of the L-curve, at which both errors are minimized. The values of

*α*used in the GEM simulations are obtained by a computer code which automatically identifies the corner of the L-curve.

## RESULTS AND DISCUSSION

*M*is the number of unknowns of the variable

*C*whose solution is sought and the superscripts

*e*and

*n*denote the exact and numerical solutions, respectively.

### Test case 1

*C*introduced along

_{s}*x*= 0 for infinite time duration in a uniform velocity field (

**V**=

**i**

*u*). The contaminant transport problem is in one spatial direction

*x*, and the exact solution for the conservative transport case with

*μ*= 0 is well known (Ogata & Banks 1961): Posed as an inverse problem, the strength of the pollution source

*C*is unknown. Although this is a 1-D problem, it is implemented with the GEM formulation in a 2-D rectangular domain with

_{s}*D*=

*D*=

_{x}*D*. The top and bottom boundaries of the domain are considered as no-flux boundaries, and the lateral extent in the

_{y}*x*-direction for the right boundary,

*L*, is set as a Dirichlet boundary with the concentration values obtained from Equation (21). To minimize the influence of the imposed no-flux condition on the top and bottom boundaries on the 1-D solution, the boundaries are kept as far apart from each other in the y direction by an arbitrary dimension of 10

^{3}km. Along the left boundary (

*x*= 0) the source concentration and flux are not known. The concentration everywhere in the domain at the initial time

*t*= 0 is zero.

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 km^{2}/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.

*γ*

_{1}is presented in Figure 2, which, as expected, indicates better prediction of the pollution source strength with smaller distance between the observation point and the source. The value of

*ε*is less than 2% when

*γ*

_{1}does not exceed 0.15 km, but for larger distances of the observation point, GEM is unable to give any meaningful prediction of the source strength. Figure 3 shows the plots of the GEM predicted pollution source strength for

*γ*

_{1}= 0.025, 0.1 and 0.15 km against the exact

*C*value of 1.0 g/L (gram/litre). To provide indicative values of the regularization parameter

_{s}*α*at each time step when

*γ*

_{1}= 0.025 km, these range between 1.6 × 10

^{−4}to 4.0 × 10

^{−5}with an average value of 6.5 × 10

^{−5}. This suggests a narrow band for the variation of the regularization parameter. The results in Figure 3 confirm our earlier finding of enhanced solution accuracy of the pollution source strength with closeness of the observation points to it. This is expected and has been confirmed by previous numerical investigations. The integrity of inverse pollution source solutions relies on the quality of the available information at observation points which depend on their number, spatial distribution, and positioning in relation to the flow field (Mahar & Datta 1997; Sun

*et al.*2006; Li & Mao 2011; Jha & Datta 2013). The CPU time for this simulation on a single CPU desktop 64-bit Intel

^{®}Core™ i7-2600, 3.40 GHz, 8 GB memory is 0.26 seconds which represents a CPU time of 1.63 × 10

^{−4}seconds/element/time step.

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.

Transport case . | Pe
. | D (km^{2}/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 (km^{2}/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 |

*C*=1.0 g/L for all times, its computed value for the three values of the Courant number for the advection dominant case (i) is presented in Figure 4. From Figure 4, it is observed that there is correct prediction of the pollution source strength by GEM for

_{s}*t*≥ 0.08 yr. At early times the numerical solutions oscillate about the exact value, with the simulation with

*Cr*

*=*1.0 exhibiting the least amount of oscillation compared to those with

*Cr*

*=*0.1 and 0.5. It should be pointed out that unlike in direct numerical modelling with GEM of contaminant transport where the value of the Courant number

*Cr*influences the numerical stability (Taigbenu 1999), its value in inverse modelling plays a dual role of influencing the numerical stability and determining the frequency at which data are available to the numerical scheme from the observation points. With respect to the latter role, because the transport process with

*Pe*

*=*50 is predominantly advection driven, the GEM simulation with

*Cr*

*=*0.1 takes about ten time steps before information of the concentration from the pollution source along

*x*

*=*0 becomes available at the two observation points along

*x*

*=*0.025 km, hence the wild oscillations of the numerical results at the early simulations times, but at later times information of the concentration is frequently available for the numerical calculations as the plume passes these points. Conversely, the GEM simulation with

*Cr*

*=*1.0 takes only one time step before concentration information from the pollution source is available at the two observation points along

*x*

*=*0.025 km, hence the more subdued oscillations at early times, but information of the concentration is less frequently available for the numerical calculations at later times. Thus, the GEM simulation with

*Cr*

*=*1.0 gives better prediction of

*C*at early times than that with

_{s}*Cr*

*=*0.1, and conversely at later times. However, it should be pointed out that it is difficult to ascertain whether these two influences of

*Cr*on the overall GEM solution are conflicting or complementing, and which influence is more significant for any contaminant transport mode. The GEM solutions for all

*Cr*values correctly predict the concentration strength

*C*at later times, and this indicates that GEM converges to the exact solution as the numerical calculations proceed with time. Another reason for the poor numerical solutions at early times is the discontinuity in the concentration along

_{s}*x*

*=*0 which is zero at

*t*

*=*0 and

*C*for

_{s}*t*> 0.

*Cr*

*=*1.0 and

*Cr*

*=*0.5 than with

*Cr*

*=*0.1, while the latter has trailing oscillations which are not observed for

*Cr*

*=*0.5 and

*Cr*

*=*1.0. The mean error

*ε*for the computed concentration for

*Cr*

*=*0.1,

*Cr*

*=*0.5 and

*Cr*

*=*1.0 is 0.48%, 0.52% and 0.53%, respectively.

*Pe*

*=*1), the GEM solutions for the pollution source strength

*C*and its exact value are presented in Figure 6. The GEM simulation with

_{s}*Cr*

*=*0.1 correctly predicts

*C*for times

_{s}*t*≥ 0.05 yr, while that with

*Cr*

*=*0.5 does so correctly for times

*t*≥ 0.075 yr, and that with

*Cr*

*=*1.0 correctly predicts

*C*for times

_{s}*t*≥ 0.1 yr. Because of the greater dispersion of the contaminant that occurs in this case than the previous, a wider spread of the plume takes place, indicating that concentration information at the observation points along

*x*

*=*0.025 km is more readily available at early times than with the earlier case. Furthermore, we do not observe the same level of oscillations of the numerical solutions for

*C*at early times as in the advection-dominant case (i), and this is most probably due to the better numerical stability characteristics observed with the direct modelling of this transport case (Taigbenu 1999). Figure 7 shows the concentration fronts at times

_{s}*t*

*=*0.1, 0.25 and 0.5 yr for Courant numbers values of

*Cr*

*=*0.1, 0.5 and 1.0. The GEM solution of the concentration fronts with

*Cr*

*=*0.1 is excellent and better than those with

*Cr*

*=*0.5 and

*Cr*

*=*1.0, and this is confirmed by the values of the mean error:

*ɛ*= 0.04% for

*Cr*

*=*0.1,

*ɛ*= 0.10% for

*Cr*

*=*0.5 and

*ɛ*= 0.17% for

*Cr*

*=*1.0.

*Pe*= 1) with noise levels of

*σ*= 2% and 5%. With

*Cr*

*=*1 used in the GEM simulation, the results for the pollution source strength are presented in Figure 8. It is observed that GEM predicts the pollution source strength

*C*reasonably well up to noise level of 2% but the solution oscillates more about the exact for noise level of 5%.

_{s}### Test case 2

*x*-dimension with

*R*

*=*1,

*μ*

*=*0 and

**V = i**

*u*. The velocity flow field is uniform, and the initial and boundary conditions for the direct problem, which had been solved by Mariño (1974), are given as:

*x*

*=*0, while the right boundary is specified as a Dirichlet boundary and the top and bottom boundaries are no-flux boundaries. A value of

*γ*= −1.0 yr

^{−1}is used in the exact solution. The parameters used for this example are:

*D*

*=*6.0 km

^{2}/yr,

*u*

*=*6.0 km/yr and Δ

*x*

*=*1.0 km which correspond to the case of marginal advection transport with

*Pe*

*=*1. There are two observation points where concentration measurements are available and these are on the top and bottom boundaries along

*x*

*=*1.0 km, and a uniform time step Δ

*t*= 3.33 × 10

^{−2}yr, corresponding to a value of 0.2 for the Courant number. The range of values of the regularization parameter

*α*at each time step when there is no noise in observation data is between 5.5 × 10

^{−6}and 5.7 × 10

^{−6}with an average value of 1.9 × 10

^{−5}, indicating small variations in regularization parameter over the simulation period. Figure 9 shows the exact and computed source strength, and the results indicate that GEM overestimates the source strength for the first ten time steps and gives correct estimates after

*t*≥ 0.4 yr. The mean error

*ɛ*for the computed source strength with and without errors randomly introduced into the observed data is 0.1566% for

*σ*= 0%, 0.1568% for

*σ*= 2% and 0.1660% for

*σ*= 5%. For the same noise levels, their influence on the concentration distribution is more muted, with

*ɛ*= 0.0508% for

*σ*= 0%,

*ɛ*= 0.0518% for

*σ*= 2% and

*ɛ*= 0.0563% for

*σ*= 5%. Figure 10 shows a comparison of the GEM predicted concentration profiles at times of 0.5 yr, 1.0 yr, 1.5 yrs and 2.0 yrs and their exact profiles. The results indicate that GEM correctly predicts the exact profiles.

### Test case 3

*x*-direction that is governed by Equation (1) with

*R*

*=*0,

*μ*

*=*0 and

**V = i**

*u*. For the direct problem, solved by Qiu

*et al.*(1998)), all the boundaries have Dirichlet boundary conditions, with the left boundary along

*x*

*=*0 having a sinusoidal concentration distribution

*C*= sin(πy) while all other boundaries have a zero concentration value. The analytical solution to this problem is given as: where Posed as an inverse problem, this problem is simulated considering the boundary along

_{s}*x*

*=*0 as a Γ

_{3}boundary where neither the concentration nor the flux is known, and boundary along

*x*

*=*1.0 km as a Neumann boundary. On the top and bottom boundaries,

*C*

*=*0. The [1 km × 1 km] domain is discretized uniformly in the

*x*and

*y*directions with Δ

*x*

*=*0.1 km and Δ

*y*

*=*0.1 km, the dispersion coefficient for the isotropic medium,

*D*

*=*1.0 km

^{2}/yr, and a uniform velocity

*u*= 10.0 km/yr, thus corresponding to the transport case with local Peclet number

*Pe*

*=*1. Using this 2-D example, the influence of the location of observation points is examined as with the first 1-D test case. The examined locations of observation points are nine internal points in the

*y*-direction along

*x*= 0.1, 0.2, 0.3, 0.4 and 0.5 km (

_{p}*p*= 1, 2, …, 9). The plot of the mean error

*ɛ*in predicting the source strength

*C*in relation to the location of the observation points is presented in Figure 11. The results exhibit a trend similar to that observed for the first 1-D test case which confirms that closeness of the observation points to the source of the pollution signal enhances the prediction capability of its strength by GEM.

_{s}*x*

*=*0.1 km, the GEM inverse solution for the pollution source signal, obtained with a regularization parameter value of 4.6 × 10

^{−4}, is presented in Figure 12 along with the exact solution for noise levels of 0, 2 and 5%. We observe that the recovered pollution source is correctly predicted by GEM for noise-free data and with noise level of 2%, but for noise level of 5% the numerical solution oscillates about the exact value. Figure 13 shows the concentration plume obtained by the exact solution and that by GEM with

*σ*= 0%, and both solutions are in good agreement. The mean error

*ε*for the concentration distribution obtained by GEM is 1.51 × 10

^{−2}% with

*σ*= 0%, 2.40 × 10

^{−2}% with noise

*σ*= 2% and 4.91 × 10

^{−2}% with noise

*σ*= 5%, while for the source strength prediction, the mean error values are 6.76 × 10

^{−2}% for

*σ*= 0%, 3.17 × 10

^{−1}% for

*σ*= 2% and 7.94 × 10

^{−1}% for

*σ*= 5%, indicating that the influence of noise on the concentration distribution is less significant compared to predicting the pollution source strength.

*α*on the numerical solutions for the spatial concentration distribution

*C*(

*x*,

*y*) and the source strength

*C*when there is no noise in the observation data. The variations of the mean error

_{s}*ɛ*for

*C*(

*x*,

*y*) and

*C*in relation to the regularization parameter are presented in Figure 14, and the results indicate that minimum value of

_{s}*ɛ*for the concentration is attained when

*α*= 4.1 × 10

^{−4}and for the source strength when

*α*= 6.2 × 10

^{−4}. The L-curve gave an optimum value of 4.6 × 10

^{−4}which lies between those two values.

### Test case 4

*et al.*(2008). The flow is uniform in the

*x*-direction with a velocity

*u*= 6 m/day and contaminant transport takes place in an orthotropic, homogeneous medium with thickness

*H*

*=*25 m and length

*L*

*=*1,000 m. The contaminant source is located at the top boundary within the section

*L*

_{2}–

*L*

_{1}, as shown in Figure 15. When the space and time dependent contaminant flux release function

*q*(

*x*,

*t*) is known, the solution procedure follows the direct approach. For the inverse model

*q*(

*x*,

*t*) is not known and additional information on the concentration is provided at the observation points shown in Figure 15. The problem is governed by Equation (1) with the initial condition:

*L*

_{1}= 200 m,

*L*

_{2}= 600 m,

*c*

_{0}= 0.2 g/L,

*D*= 2 × 10

_{x}^{5}m

^{2}/day,

*D*= 10

_{y}^{4}m

^{2}/day,

*μ*= 2.0 day

^{−1},

*R*= 1.0, Δ

*x*= 25 m, Δ

*y*= 1.25 m and Δ

*t*

*=*0.025 day. The computational domain is discretized into 800 rectangular uniform elements. The expression for the release history of the contaminant source is given by:

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.

*q*(

*x*,

*t*) and the concentration distribution

*C*(

*x*,

*y*,

*t*) are presented as error plots that are based on the mean error

*ɛ*that is given by Equation (20). The error plots for the pollution flux when

*σ*= 0%,

*σ*

*=*2% and

*σ*= 5% are presented in Figure 16, while the error plots for the concentration distribution are presented in Figure 17. Comparing Figures 16 and 17, the errors for the spatial concentration distribution are about three orders of magnitude lower than those from predicting the pollution source flux. Furthermore, introduction of randomized noise into observation data tends to produce a more muted influence on the inverse solution for the concentration distribution than for the pollution flux (Taigbenu 2015). The CPU time for running each of the simulations with a given noise level is 13.61 seconds, representing a CPU time of 4.3 × 10

^{−4}seconds/element/time step. The values of the regularization parameter that were used in the simulation of the noise free data range from 1.03 × 10

^{−4}to 1.23 × 10

^{−4}with an average value of 1.06 × 10

^{−4}.

### Test case 5

*L*

_{1}= 200 m,

*L*

_{2}= 600 m,

*c*

_{0}= 0.35 g/L,

*D*= 2 × 10

_{x}^{5}m

^{2}/day,

*D*= 4 × 10

_{y}^{4}m

^{2}/day,

*μ*= 2.0 day

^{−1},

*R*= 1.0, Δ

*x*= 25 m, Δ

*y*= 2.5 m and Δ

*t*

*=*0.025 day. The true flux of the pollution source occurs over an active period of 1 day, and it is given by the expression: This expression for the flux is used in the direct formulation to generate the concentrations over a period of 2 days at the 17 observation points located at a depth of 5 m from the top where the pollution source is located. These concentration values at the observation points are considered as the actual values. With these observed concentration values, the inverse GEM simulation is carried on the domain with the top boundary between

*x*= 200 m and

*x*= 600 m considered as a Γ

_{3}boundary where the concentration and the flux are not specified. The computational domain was discretized into 800 uniform rectangular elements for both the direct and inverse GEM simulations. The inverse GEM solution for the release history of the pollution flux

*q*(

*x*,

*t*) is presented in Figure 18. The results indicate that GEM correctly predicts the pollution source flux evolution, as well as the active period of 1 day over which the pollution occurs. The error plots of the inverse GEM solutions for the release history of pollution flux

*q*(

*x*,

*t*) and the concentration distribution

*C*(

*x*,

*y*,

*t*) are presented in Figures 19 and 20, respectively, for noise levels of

*σ*= 0%,

*σ*

*=*2% and

*σ*= 5%. It should be noted that the error plots for the pollution flux are presented only over the active period of the pollutant

*t*≤ 1 day. As with Test case 4, the errors for the spatial concentration distribution are about three orders of magnitude lower than those from predicting the pollution source flux, and introduction of randomized noise in the observation data produces a more muted influence on the inverse solution for the concentration distribution than for the pollution flux. The CPU time for the simulation run with no noise in the data is 24.36 seconds, and this represents a CPU time of 3.8 × 10

^{−4}seconds/element/time step. The values of the regularization parameter for the noise free data range from 1.03 × 10

^{−4}to 1.16 × 10

^{−4}with an average value of 1.05 × 10

^{−4}.

## 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.