#### Abstract

We present a model of a polluted groundwater site. The model consists of a coupled system of advection-diffusion-reaction equations for the groundwater level and the concentration of the pollutant. We use the complete flux scheme for the space discretization in combination with the -method for time integration and we prove a new stability result for the scheme. Numerical results are computed for the Guarani Aquifer in South America and they show good agreement with results in literature.

#### 1. Introduction

There is a great interest in developing groundwater models in the last decades. Groundwater is an important component of water resource systems. For example, more than half of the United States population depend on groundwater for drinking water [1]. Unfortunately, since groundwater is susceptible to pollutants, the quality of groundwater can be poor and is sometimes even deteriorating. Groundwater pollution occurs because of, among other things, disposal of industrial wastes such as gasoline, oil, and chemicals and agricultural activities such as the use of fertilizers. Also, uncontrolled hazardous waste, chemicals, road salts, and contaminants from dumps can make their way down into the groundwater [2]. Moreover, the groundwater is currently withdrawing which makes the quality even worse. This all indicates that groundwater might be unsafe for human use.

To solve these problems, clean water can be injected into aquifers and polluted groundwater can be pumped out. For this, the location of injecting and pumping as well as their rates is of great importance. Decisions have to be made for the proper locations and the total volume of groundwater that should be injected and/or pumped. Numerical simulations of groundwater flow enable us to make predictions about groundwater level and the amount of pollution and are therefore very useful. Therefore, mathematical models are needed for groundwater flow.

In groundwater modeling many different types of models are considered, for example, flow models, which describe the hydraulic head, or solute transport models, which describe the concentration level [3]. Our model is a combination of both a flow model and a solute transport model. Thus, the model includes both the concentration of the pollutant and the groundwater level. New in our model is that we include the effect of pumping and injecting water.

Our model consists of two equations, a pure diffusion equation for the groundwater level and an equation of advection-diffusion-reaction type that describes the concentration of the pollutants. For space discretization of the advection-diffusion-reaction equation we use the finite volume-complete flux scheme developed in [4, 5]. This scheme has proven to be a robust method to discretize advection-diffusion-reaction equations, which is uniformly second-order accurate, even for dominant advection, yet it does not produce any spurious oscillations, which makes the scheme very suitable for groundwater simulations. For the diffusion equation, the complete flux scheme reduces to the central difference scheme, which we use to discretize the groundwater level equation. We have applied the complete flux scheme to various problems from continuum physics, such as the flow of multicomponent mixtures and plasmas [6, 7], incompressible fluid flow [8], or the numerical simulation of plankton populations [9].

For the temporal integration the -method is used. New stability conditions for the combination of the -method with the finite volume-complete flux scheme for spatial discretization are derived in this paper.

We have organized our paper as follows. First, in Section 2, we describe our groundwater model. Next, in Section 3, we describe two versions of the complete flux scheme for our model. In Section 4 we derive stability conditions for the fully discrete scheme. The numerical results are discussed in Section 5. Finally, we end with conclusions in Section 6.

#### 2. Mathematical Model

The Guarani Aquifer in South America is one of the largest aquifers in the world. It is a very important source of drinking water in Argentina, Brazil, Paraguay, and Uruguay and should obviously not be polluted. However, the aquifer is very vulnerable to pollution due to land-surface activities [10]. Therefore it is interesting to define a model problem for the transport of polluted water based on the parameters of this aquifer [11–14].

We consider a groundwater site next to an aquifer. The groundwater in this basin is polluted and flowing towards the aquifer. Since the polluted water should not enter the basin we assume a well at the boundary between basin and aquifer, where all polluted water is removed. We restrict ourselves to a one-dimensional domain , with the width of the basin, adjacent to the aquifer at . The soil in the basin consists mainly of clay and clay-like material, so we choose a uniform one-layer model consisting of clay only. We measure the groundwater level in the layer relative to a flat impermeable bottom. Furthermore, we are not particularly interested in the detailed composition of the polluted water, and therefore we only consider one typical polluting species, namely, lodide. The unknowns of our model are therefore the groundwater level in the layer and the concentration of lodide. A schematic picture of the groundwater site is shown in Figure 1.

Polluted groundwater can be removed from the basin by pumping out polluted water and/or injecting clean water, which advances the polluted water towards the well, where it is taken out. We investigate the impact of pumping and injecting on the (rate of) removal of polluted water.

The governing equations are the conservation laws of volume, or rather mass, of the ground water and the mass balance equation for lodide. Pumping and injecting are included as source terms in these equations. Our model is inspired by [15] and its governing equations readwhere is the groundwater level, the concentration of the pollutant, the velocity of the lodide pollution, and the diffusion coefficient of lodide. Furthermore, is defined as with , the hydraulic conductivity, and , the specific storage (due to the presence of the well of the right side, the pollutant cannot enter the aquifer. Therefore, the considered groundwater site is classified as confined. The effect of phreatic storage usually occurs in unconfined sites and is therefore not included). The specific storage describes the amount of water per unit volume of a saturated aquifer that is absorbed or expelled due to changes in the compression of the fluid and medium caused by a change in hydraulic head [16]. Since we consider a uniform site, we assume , and to be constant coefficients. The constituent relation in (1c) is Darcy’s law. The source terms and are given bywith and being the pumping and injection function, respectively, and with being the injecting factor. From (2) it is evident that both pumping and injecting have an impact on the water level; however, the concentration of lodide is only influenced by injection. We have chosen to describe the pumping and injecting functions by the following Gaussian profiles:where is the pumping rate, is the pumping position, is the injecting rate, and is the injecting position. The parameters and determine the pumping width and injecting width, respectively.

We choose a uniform initial groundwater level and describe the initial lodide concentration as a Gaussian profile, that is,where is the initial groundwater level, is the maximum initial concentration, is the location of pollution, and determines the width of the initial pollution. To complete our model, we prescribe the following boundary conditions:where is the depth of the aquifer. The boundary conditions (5a) imply that there is no influx of polluted water at since and . The boundary condition means that there is a well located at where all polluted water is immediately removed. We refer to (1a)–(5b) as the* groundwater model*.

#### 3. Numerical Scheme for the Groundwater Model

In this section we outline the* finite volume-complete flux* (FV-CF) scheme for the discretization of the groundwater model; for a detailed derivation and validation of the scheme, see [4].

First, we present the stationary flux approximation and subsequently its extension to time-dependent problems. Consider the stationary version of (1b) formulated aswhere is the flux corresponding to . Let us introduce an equidistant grid with being the grid size and being the number of grid points. Furthermore, we cover the domain with control volumes , where . Integrating the first equation in (6) over the control volume and using the midpoint rule for the integral of the source term , we obtain the discrete conservation lawwhere is the numerical flux at and .

The complete flux approximation for is derived from the local BVP for (6) on , subject to the boundary conditions and . The derivation proceeds in the following steps. First, we integrate the conservation law in (6) from the interface to to obtain the integral balance. ConsiderSecond, introducing the auxiliary function , defined bywe can rewrite the flux in terms of its integrating factor, that is,Next, substituting this expression in (8), isolating the derivative, integrating the resulting equation from to , and applying the boundary conditions, we obtain the following expressions for the flux:where and are the homogeneous and inhomogeneous part of the flux, corresponding to the advection-diffusion operator and the source term, respectively. Finally, to relate the inhomogeneous flux to the source term, we substitute in (11c) and change the order of integration. This way we findwhere is referred to as Green’s function for the flux; see [4] for more details.

To derive expressions for the numerical flux, we have to apply appropriate quadrature rules to all integrals involved. Introducing the Péclet function , with an abuse of notation, defined aswe obtain for the numerical flux the following expressions:where the functions and are defined as see Figure 2. Furthermore, (14a)–(14c) contain the average Péclet number and the upwind value of the source term , defined aswhere is the central difference approximation of , that is, and where is the average velocity. Clearly, also the numerical flux contains two parts, namely, the homogeneous flux , corresponding to the advection-diffusion operator, and the inhomogeneous flux , taking into account the effect of the source term. We refer to the flux approximations in (14a)–(14c) as the* complete flux* (CF) scheme.

**(a)**

**(b)**

Note that for or the homogeneous flux reduces to the central difference flux or the upwind flux , respectively. Indeed, rearranging terms in (14b) we can write the expression for the homogeneous flux as a weighted average of and , that is,where is the upwind value of at the interface . Consequently, the numerical flux does not produce spurious oscillations for , contrary to the central difference scheme which generates these oscillations for . This allows us to use much coarser grid for the CF scheme than would be required for the central difference scheme. Moreover, we can show that the inclusion of the inhomogeneous term guarantees second-order convergence in the limit ; see [17].

Next, we consider the time-dependent version of the conservation law, which reads with corresponding semidiscrete conservation law. Considerwhere . For the numerical flux we have two options. First, we can simply take the flux approximation (14a)–(14c), henceforth referred to as the* stationary complete flux* (SCF) scheme. Combining the semidiscrete conservation law (20) with the complete flux approximation (14a)–(14c) we obtainwhere the coefficients are defined byThe alternative is to include the time derivative in the source term; that is, we introduce the modified source term as and we replace in (14c) by . The inhomogeneous numerical flux then changes towhere denotes the upwind value of at the interface . Obviously, the homogeneous flux remains the same. We refer to this flux approximation as the* transient complete flux* (TCF) scheme. The corresponding semidiscretization then reads

To discretize (1a) we also use the FV-CF scheme. In this case the corresponding Péclet function is identically , and since , the inhomogeneous flux vanishes. Consequently, the scheme reduces to the standard central difference scheme. Consider

For time integration of either (21) or (25) coupled with (26) we use the -method [18]. However, since the Péclet numbers depend on the numerical approximation , we apply a predictor-corrector approach where we first integrate (26), update the Péclet numbers, and subsequently integrate either (21) or (25).

#### 4. Stability Analysis

In this section we investigate stability of the -method applied to both (21) and (25). To that purpose we assume that and ignore boundary conditions and source term. For the latter ODE-system we moreover assume that .

Both semidiscretizations give rise to an ODE system of the formwhere the discretization matrices and represent the homogeneous and inhomogeneous flux differences, respectively. Furthermore, for the SCF scheme and for the TCF scheme. The -method applied to (27) readswhere denotes the approximation of at time level with being the time step, and so forth. In the stability analysis which follows, we ignore the source term and the boundary term . We first consider (25); the stability result for (21) follows readily as a special case.

Theorem 1. *The -method applied to the TCF semidiscretisation (25) is stable if the following conditions hold:with and and where the function is defined by *

*Proof. *Substituting the parameters , and so forth, defined in (22a) and (22b), the -method for (25) reads where we have used since . To investigate stability, we substitute the discrete Fourier mode. Consider and we obtain the following expressions for the amplification factor , that is, with . Elaborating the stability requirement, we find with and . Straightforward evaluation gives And, consequently, reduces to Next, straightforward differentiation gives for , from which we conclude that and . It is obvious that, for , the inequality unconditionally holds. Consider next the case . A sufficient condition for stability is then from which the stability result readily follows.

Corollary 2. *The -method applied to (21) is stable ifwith , and defined in Theorem 1.*

*Proof. *Take and repeat the proof of Theorem 1.

To conclude, the standard stability result for the -method applied to (26) follows, if we substitute in the stability condition (37); see [18].

#### 5. Numerical Results

In this section we present numerical solutions for three different scenarios of the groundwater model: first, no injecting and no pumping and consequently also no flow, second, injecting and no pumping, and third, both injecting and pumping. We refer to these scenarios as no flow, injecting, and injecting-pumping, respectively. The corresponding parameter values are given in Table 1 and are based on [11–14]. All numerical solutions are computed in dimensionless form using the TCF scheme in combination with the trapezoidal rule for time integration , on a grid with grid points and with (dimensionless) time step , and subsequently converted to full dimensional form. Numerical solutions for the SCF scheme are similar and are not included.

The choice of and is justified by the grid validation of the numerical solution presented in Table 2. In this table we present the numerical approximations of and at and at time computed with (dimensionless) time step on a sequence of grids, to verify the dependence of the numerical solution on the grid. Clearly, for the numerical solution is practically independent of the grid.

To verify the order of convergence of the SCF and TCF semidiscretizations, we repeatedly compute by successively halving the grid size and apply Richardson extrapolation [19]. For the (dimensionless) time step we take . Let denote the numerical approximation of computed with grid size . If the semidiscretization is th order convergent, we have We will present values for the injecting and injecting-pumping scenario.

##### 5.1. No-Flow Scenario

In this case we do not include pumping nor injecting. Since the initial groundwater level is constant, there is no flow of the pollutant and the concentration equation is just a diffusion equation. Consequently, the SCF and TCF schemes reduce to the central difference scheme. We show the result in Figure 3. The figure shows a little spreading of the initial concentration, with the maximum value reduced to approximately 90% of its initial value. No pollutant is removed.

**(a)**

**(b)**

##### 5.2. Injecting Scenario

In this scenario we inject clean water at the left of the domain and do not pump any polluted water. The clean water forces the polluted water to the well on the right where it is removed. The solution is shown in Figure 4. The source term is virtually outside a very small layer near and this implies that in the steady limit for the groundwater level tends to a linear profile, and, consequently, the velocity approaches a constant. Since the diffusion coefficient is very small, transport of the pollutant is dominated by advection, as is evident from Figure 4. A zoom-in of the lodide concentration near is presented in Figure 5(a), confirming that all polluted water is removed at the right end of the domain.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

After about 25 years the total pollution, integrated over the entire domain, has decreased to less than 5% of its original value, which is shown in Figure 5(b). This means that the pollution moves with a speed of less than 1 meter a year. At the same time, the groundwater level increases.

In Table 3 we present the values for both the SCF and TCF schemes as a function of the number of grid points . Clearly, both schemes show second-order convergence for . In addition, we have included the (in absolute value) maximum Péclet number , indicating that advection is dominant on the coarser grids.

##### 5.3. Injecting-Pumping Scenario

In this scenario we add pumping to the model. We inject water just left of the concentration plume at and pump polluted water out at the right of the plume at . We show the results in Figure 6. From this figure we can distinguish different regimes of the solution. First, left of the injection peak, the groundwater level increases in time. Second, between the injection and pumping peaks both source terms are virtually leading to a linear decreasing groundwater level and a concentration profile propagating at constant velocity, just as in the injecting scenario. Third, near , the source term resulting in the flattening of the groundwater level and a decrease of the velocity . Thus, left of the pumping peak, the pollutant moves faster than at the right, and since , the pollutant piles up. Finally, beyond the pumping peak the groundwater level is again linear and the pollutant is transported with constant velocity towards the well. A zoom-in of the lodide concentration near is given in Figure 7(a).

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

Note that the groundwater level profile has not increased that much compared to the previous situation. This results in a longer period for the pollutant to be removed, and after about 33 years the pollution has decreased to less than 1% of the original amount. We have shown this in Figure 7(b).

The values are presented in Table 4, from which we conclude that also in this scenario the SCF and TCF schemes demonstrate second-order convergence for the grid size . We see that, for removing the pollutant from the site, the steepness of the groundwater level is of great importance, that is, the steeper the faster it is removed. In the situations we have considered, it takes approximately 20 to 30 years to remove the lodide. At first sight, this seems to be very long. However, one has to keep in mind that the soil type is clay, which has a low hydraulic conductivity. Our results correspond to some of the results in the literature, which also find that pollutants in groundwater move with only a few feet a year.

#### 6. Summary and Conclusions

In this section we present our main conclusions and propose some extensions to the model. Our model describes a polluted groundwater site next to the Guarani Aquifer in South America. We have included both injecting of clean and pumping of polluted water.

Since the governing equation for the concentration of pollutant is of advection-diffusion-reaction type, we employed the complete flux scheme for space discretization. The scheme has proven to be very suitable, since it is uniformly second-order accurate, even for dominant advection, and it does not produce any spurious oscillations. For time integration we use the -scheme. A stability condition for the fully discrete scheme is derived.

We have investigated the effect of injecting and pumping on the (rate of) removal of pollutants and conclude that the injection of clean water upstream of the pollution peak together with pumping out polluted water downstream of the pollution peak enhances the removal of polluted water.

Our model can be extended and improved in many different ways of which we mention a few. First, the source terms could be made more realistic by letting them depend on the concentration of the pollutant as well as the groundwater level. Also a possibility is to make the source terms time dependent. In this way one could choose not to pump anymore if it would not affect the removal any further. Another possibility is to extend the model to two dimensions. Also here, one could still use the finite volume-complete flux scheme to solve the advection-diffusion-reaction equation. A two-dimensional extension of the scheme is included in [4].

#### Competing Interests

The authors declare that they have no competing interests.