# A Finite Volume-Complete Flux Scheme for a Polluted Groundwater Site

Authors: M. F. P. ten Eikelder; J. H. M. ten Thije Boonkkamp; M. P. T. Moonen; B. V. Rathish Kumar
Journal: Advances in Numerical Analysis
 (2016)

---

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


---

## Body

## 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 Section2, 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 domain0
≤
x
≤
L, with L the width of the basin, adjacent to the aquifer at x
>
L. 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.Figure 1
Layer of a polluted groundwater site adjacent to an aquifer.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 read(1a)
∂
h
∂
t
=
a
∂
2
h
∂
x
2
+
s
h
,
(1b)
∂
c
∂
t
+
∂
∂
x
q
c
=
D
∂
2
c
∂
x
2
+
s
c
,
(1c)
q
=
-
κ
∂
h
∂
x
,where h
(
m
) is the groundwater level, c
(
k
g
m
-
3
) the concentration of the pollutant, q
(
m
s
-
1
) the velocity of the lodide pollution, and D
(
m
2
s
-
1
) the diffusion coefficient of lodide. Furthermore, a
(
m
2
s
-
1
) is defined as κ
/
S
S with κ
(
m
s
-
1
), the hydraulic conductivity, and S
S
(
m
-
1
), 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 D
,
κ, and S
S to be constant coefficients. The constituent relation in (1c) is Darcy’s law. The source terms s
h
(
m
s
-
1
) and s
c
(
k
g
m
-
3
s
-
1
) are given by(2)
s
h
=
-
p
+
m
,
s
c
=
m
f
m
,with p and m being the pumping and injection function, respectively, and with m
f
[
k
g
m
-
4
] 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 p
,
m by the following Gaussian profiles:(3a)
p
x
=
p
0
exp
⁡
-
p
w
x
-
x
p
2
,
(3b)
m
x
=
m
0
exp
⁡
-
m
w
x
-
x
m
2
,where p
0
(
m
s
-
1
) is the pumping rate, x
p
(
m
) is the pumping position, m
0
(
m
s
-
1
) is the injecting rate, and x
m
(
m
) is the injecting position. The parameters p
w
(
m
-
2
) and m
w
(
m
-
2
) 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,(4a)
h
x
,
0
=
h
0
,
(4b)
c
x
,
0
=
c
0
exp
⁡
-
c
w
x
-
x
c
2
,where h
0
(
m
s
-
1
) is the initial groundwater level, c
0
(
k
g
m
-
3
) is the maximum initial concentration, x
c
(
m
) is the location of pollution, and c
w
(
m
-
2
) determines the width of the initial pollution. To complete our model, we prescribe the following boundary conditions:(5a)
∂
h
∂
x
0
,
t
=
0
,
∂
c
∂
x
0
,
t
=
0
,
(5b)
h
L
,
t
=
H
,
c
L
,
t
=
0
,where H is the depth of the aquifer. The boundary conditions (5a) imply that there is no influx of polluted water at x
=
0 since q
=
0 and ∂
c
/
∂
x
=
0. The boundary condition c
(
L
,
t
)
=
0 means that there is a well located at x
=
L where all polluted water is immediately removed. We refer to (1a)–(5b) as thegroundwater model.


## 3. Numerical Scheme for the Groundwater Model

In this section we outline thefinite 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 as(6)
d
f
d
x
=
s
c
,
f
=
q
c
-
D
d
c
d
x
,where f is the flux corresponding to c. Let us introduce an equidistant grid x
j
=
(
j
-
1
)
Δ
x
(
j
=
1,2
,
…
,
N
) with Δ
x
=
L
/
(
N
-
1
) being the grid size and N being the number of grid points. Furthermore, we cover the domain with control volumes Ω
j
=
(
x
j
-
1
/
2
,
x
j
+
1
/
2
)
(
j
=
2,3
,
…
,
N
-
1
), where x
j
+
1
/
2
=
x
j
+
Δ
x
/
2. Integrating the first equation in (6) over the control volume Ω
j and using the midpoint rule for the integral of the source term s
c, we obtain the discrete conservation law(7)
F
j
+
1
/
2
-
F
j
-
1
/
2
=
s
c
,
j
Δ
x
,where F
j
+
1
/
2 is the numerical flux at x
j
+
1
/
2 and s
c
,
j
=
s
c
(
x
j
).The complete flux approximation forF
j
+
1
/
2 is derived from the local BVP for (6) on (
x
j
,
x
j
+
1
), subject to the boundary conditions c
(
x
j
)
=
c
j and c
(
x
j
+
1
)
=
c
j
+
1. The derivation proceeds in the following steps. First, we integrate the conservation law in (6) from the interface x
j
+
1
/
2 to x
∈
(
x
j
,
x
j
+
1
) to obtain the integral balance. Consider(8)
f
x
-
f
x
j
+
1
/
2
=
S
x
,
S
x
≔
∫
x
j
+
1
/
2
x
s
c
ξ
d
ξ
.Second, introducing the auxiliary function Λ, defined by(9)
Λ
x
≔
1
D
∫
x
j
+
1
/
2
x
q
ξ
d
ξ
=
-
κ
D
h
x
-
h
x
j
+
1
/
2
,we can rewrite the flux in terms of its integrating factor, that is,(10)
f
=
-
D
d
d
x
c
e
-
Λ
e
Λ
.Next, substituting this expression in (8), isolating the derivative, integrating the resulting equation from x
j to x
j
+
1, and applying the boundary conditions, we obtain the following expressions for the flux:(11a)
f
x
j
+
1
/
2
=
f
h
x
j
+
1
/
2
+
f
i
x
j
+
1
/
2
,
(11b)
f
h
x
j
+
1
/
2
=
-
e
-
Λ
x
j
+
1
c
j
+
1
-
e
-
Λ
x
j
c
j
∫
x
j
x
j
+
1
D
-
1
e
-
Λ
d
x
,
(11c)
f
i
x
j
+
1
/
2
=
-
∫
x
j
x
j
+
1
D
-
1
e
-
Λ
S
d
x
∫
x
j
x
j
+
1
D
-
1
e
-
Λ
d
x
,where f
h
(
x
j
+
1
/
2
) and f
i
(
x
j
+
1
/
2
) 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 S
(
x
) in (11c) and change the order of integration. This way we find(12)
f
i
x
j
+
1
/
2
=
Δ
x
∫
0
1
G
σ
s
c
x
σ
d
σ
,
σ
x
=
x
-
x
j
Δ
x
,where G
σ 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 functionP, with an abuse of notation, defined as(13)
P
≔
q
Δ
x
D
=
-
κ
Δ
x
D
d
h
d
x
,we obtain for the numerical flux F
j
+
1
/
2 the following expressions:(14a)
F
j
+
1
/
2
=
F
j
+
1
/
2
h
+
F
j
+
1
/
2
i
,
(14b)
F
j
+
1
/
2
h
=
-
D
Δ
x
B
P
¯
j
+
1
/
2
c
j
+
1
-
B
-
P
¯
j
+
1
/
2
c
j
,
(14c)
F
j
+
1
/
2
i
=
1
2
-
W
P
¯
j
+
1
/
2
s
c
,
u
,
j
+
1
/
2
Δ
x
,where the functions B and W are defined as(15)
B
z
≔
z
e
z
-
1
;
W
z
≔
e
z
-
1
-
z
z
e
z
-
1
, see Figure 2. Furthermore, (14a)–(14c) contain the average Péclet number P
-
j
+
1
/
2 and the upwind value of the source term s
c
,
u
,
j
+
1
/
2, defined as(16a)
P
-
j
+
1
/
2
≔
1
2
P
j
+
P
j
+
1
,
(16b)
s
c
,
u
,
j
+
1
/
2
≔
s
c
,
j
i
f
q
-
j
+
1
/
2
≥
0
,
s
c
,
j
+
1
i
f
q
-
j
+
1
/
2
<
0
,where P
j is the central difference approximation of P
(
x
j
), that is,(17)
P
j
=
-
κ
2
D
h
j
+
1
-
h
j
-
1
, and where q
-
j
+
1
/
2
=
1
/
2
q
j
+
q
j
+
1 is the average velocity. Clearly, also the numerical flux contains two parts, namely, the homogeneous flux F
j
+
1
/
2
h, corresponding to the advection-diffusion operator, and the inhomogeneous flux F
j
+
1
/
2
i, taking into account the effect of the source term. We refer to the flux approximations in (14a)–(14c) as thecomplete flux (CF) scheme.Figure 2
The functionsB (a) and W (b).
(a)
(b)Note that forP
=
0 or |
P
|
→
∞ the homogeneous flux F
j
+
1
/
2
h reduces to the central difference flux F
j
+
1
/
2
c
d or the upwind flux F
j
+
1
/
2
u
w, respectively. Indeed, rearranging terms in (14b) we can write the expression for the homogeneous flux as a weighted average of F
j
+
1
/
2
c
d and F
j
+
1
/
2
u
w, that is,(18)
F
j
+
1
/
2
h
=
q
¯
j
+
1
/
2
2
W
P
¯
j
+
1
/
2
1
2
c
j
+
c
j
+
1
+
1
-
W
P
¯
j
+
1
/
2
c
u
,
j
+
1
/
2
-
D
Δ
x
c
j
+
1
-
c
j
=
2
W
P
¯
j
+
1
/
2
F
j
+
1
/
2
c
d
+
1
-
2
W
P
¯
j
+
1
/
2
F
j
+
1
/
2
u
w
,where c
u
,
j
+
1
/
2 is the upwind value of c at the interface x
j
+
1
/
2. Consequently, the numerical flux F
j
+
1
/
2 does not produce spurious oscillations for P
→
∞, contrary to the central difference scheme which generates these oscillations for P
>
2. 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 P
→
∞; see [17].Next, we consider the time-dependent version of the conservation law, which reads(19)
∂
c
∂
t
+
∂
f
∂
x
=
s
c
,
f
=
q
c
-
D
∂
c
∂
x
, with corresponding semidiscrete conservation law. Consider(20)
c
˙
j
Δ
x
+
F
j
+
1
/
2
-
F
j
-
1
/
2
=
s
c
,
j
Δ
x
,where c
˙
j
t
≈
∂
c
/
∂
t
x
j
,
t. For the numerical flux F
j
+
1
/
2 we have two options. First, we can simply take the flux approximation (14a)–(14c), henceforth referred to as thestationary complete flux (SCF) scheme. Combining the semidiscrete conservation law (20) with the complete flux approximation (14a)–(14c) we obtain(21)
c
˙
j
Δ
x
-
α
j
-
1
/
2
c
j
-
1
+
α
j
+
1
/
2
+
β
j
-
1
/
2
c
j
-
β
j
+
1
/
2
c
j
+
1
=
Δ
x
γ
j
-
1
/
2
s
c
,
j
-
1
+
1
-
γ
j
+
1
/
2
+
δ
j
-
1
/
2
s
c
,
j
-
δ
j
+
1
/
2
s
c
,
j
+
1
,where the coefficients are defined by(22a)
α
j
+
1
/
2
≔
D
Δ
x
B
-
P
¯
j
+
1
/
2
,
β
j
+
1
/
2
≔
D
Δ
x
B
P
¯
j
+
1
/
2
,
(22b)
γ
j
+
1
/
2
≔
max
⁡
1
2
-
W
P
¯
j
+
1
/
2
,
0
,
δ
j
+
1
/
2
≔
min
⁡
1
2
-
W
P
¯
j
+
1
/
2
,
0
.The alternative is to include the time derivative in the source term; that is, we introduce the modified source term s
^
c as(23)
s
^
c
≔
s
c
-
∂
c
∂
t
, and we replace s
c in (14c) by s
^
c. The inhomogeneous numerical flux then changes to(24)
F
j
+
1
/
2
i
=
1
2
-
W
P
¯
j
+
1
/
2
s
c
,
u
,
j
+
1
/
2
-
c
˙
u
,
j
+
1
/
2
Δ
x
,where c
˙
u
,
j
+
1
/
2 denotes the upwind value of ∂
c
/
∂
t at the interface x
j
+
1
/
2. Obviously, the homogeneous flux remains the same. We refer to this flux approximation as thetransient complete flux (TCF) scheme. The corresponding semidiscretization then reads(25)
Δ
x
γ
j
-
1
/
2
c
˙
j
-
1
+
1
-
γ
j
+
1
/
2
+
δ
j
-
1
/
2
c
˙
j
-
δ
j
+
1
/
2
c
˙
j
+
1
-
α
j
-
1
/
2
c
j
-
1
+
α
j
+
1
/
2
+
β
j
-
1
/
2
c
j
-
β
j
+
1
/
2
c
j
+
1
=
Δ
x
γ
j
-
1
/
2
s
c
,
j
-
1
+
1
-
γ
j
+
1
/
2
+
δ
j
-
1
/
2
s
c
,
j
-
δ
j
+
1
/
2
s
c
,
j
+
1
.To discretize (1a) we also use the FV-CF scheme. In this case the corresponding Péclet function is identically 0, and since W
(
0
)
=
1
/
2, the inhomogeneous flux vanishes. Consequently, the scheme reduces to the standard central difference scheme. Consider(26)
h
˙
j
Δ
x
+
a
Δ
x
-
h
j
-
1
+
2
h
j
-
h
j
+
1
=
Δ
x
s
h
,
j
.For time integration of either (21) or (25) coupled with (26) we use the ϑ-method [18]. However, since the Péclet numbers P
¯
j
+
1
/
2 depend on the numerical approximation h
j, 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 q
x
=
C
o
n
s
t and ignore boundary conditions and source term. For the latter ODE-system we moreover assume that q
>
0.Both semidiscretizations give rise to an ODE system of the form(27)
M
c
˙
+
A
c
=
B
s
c
+
b
,where the discretization matrices A and B represent the homogeneous and inhomogeneous flux differences, respectively. Furthermore, M
=
Δ
x
I for the SCF scheme and M
=
B for the TCF scheme. The ϑ-method applied to (27) reads(28)
1
Δ
t
M
c
n
+
1
-
c
n
+
1
-
ϑ
A
c
n
+
ϑ
A
c
n
+
1
=
1
-
ϑ
B
s
c
n
+
b
n
+
ϑ
B
s
c
n
+
1
+
b
n
+
1
,
0
≤
ϑ
≤
1
,where c
n denotes the approximation of c
(
t
) at time level t
n
=
n
Δ
t
n
=
0,1
,
2
,
… with Δ
t
>
0 being the time step, and so forth. In the stability analysis which follows, we ignore the source term s
c and the boundary term b. 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:(29)
0
≤
ϑ
<
1
2
,
d
≤
1
1
-
2
ϑ
W
P
C
1
2
P
,
1
2
≤
ϑ
≤
1
,with d
=
D
Δ
t
/
Δ
x
2 and P
=
q
Δ
x
/
D and where the function C is defined by (30)
C
z
=
tanh
⁡
z
z
.Proof.
Substituting the parametersα
j
±
1
/
2, and so forth, defined in (22a) and (22b), the ϑ-method for (25) reads(31)
Δ
x
1
2
-
W
P
1
Δ
t
c
j
-
1
n
+
1
-
c
j
-
1
n
+
Δ
x
1
2
+
W
P
1
Δ
t
c
j
n
+
1
-
c
j
n
+
1
-
ϑ
D
Δ
x
-
B
P
c
j
+
1
n
+
B
-
P
+
B
P
c
j
n
-
B
-
P
c
j
-
1
n
+
ϑ
D
Δ
x
-
B
P
c
j
+
1
n
+
1
+
B
-
P
+
B
P
c
j
n
+
1
-
B
-
P
c
j
-
1
n
+
1
=
0
, where we have used c
˙
u
,
j
+
1
/
2
=
c
˙
j since q
>
0. To investigate stability, we substitute the discrete Fourier mode. Consider(32)
c
j
n
=
g
n
e
i
k
x
j
, and we obtain the following expressions for the amplification factor g, that is, (33)
g
=
a
1
+
i
a
2
-
1
-
ϑ
d
b
1
+
i
b
2
a
1
+
i
a
2
+
ϑ
d
b
1
+
i
b
2
,
a
1
=
1
2
-
W
P
cos
⁡
φ
+
1
2
+
W
P
,
a
2
=
-
1
2
-
W
P
sin
⁡
φ
,
b
1
=
2
B
+
+
B
-
sin
2
⁡
1
2
φ
,
b
2
=
P
sin
⁡
φ
,
φ
=
k
Δ
x
, with 0
≤
φ
<
π. Elaborating the stability requirementg
2
≤
1, we find∗
a
·
b
+
ϑ
-
1
2
d
b
2
≥
0
, with a
=
a
1
a
2
T and b
=
b
1
b
2
T. Straightforward evaluation gives (34)
a
·
b
=
4
p
1
φ
sin
2
⁡
1
2
φ
,
b
2
=
4
p
2
φ
sin
2
⁡
1
2
φ
,
p
1
φ
=
1
2
cos
2
⁡
1
2
φ
+
2
W
P
sin
2
⁡
1
2
φ
B
P
+
B
-
P
+
P
W
P
-
1
2
cos
2
⁡
1
2
φ
,
p
2
φ
=
B
P
+
B
-
P
2
sin
2
⁡
1
2
φ
+
P
2
cos
2
⁡
1
2
φ
, And, consequently, ∗ reduces to∗
∗
p
1
φ
+
ϑ
-
1
2
d
p
2
φ
≥
0
. Next, straightforward differentiation gives (35)
p
1
′
φ
=
B
P
W
P
-
1
2
sin
⁡
φ
≤
0
,
p
2
′
φ
=
2
B
P
B
-
P
sin
⁡
φ
≥
0
, for 0
≤
φ
<
π, from which we conclude that 1
=
p
1
0
≥
p
1
π
=
W
P
B
P
+
B
-
P
>
0 and 0
<
P
2
=
p
2
0
≤
p
2
(
π
)
=
B
P
+
B
P
2. It is obvious that, for 1
/
2
≤
ϑ
≤
1, the inequality ∗
∗ unconditionally holds. Consider next the case 0
≤
ϑ
<
1
/
2. A sufficient condition for stability is then(36)
p
1
π
+
ϑ
-
1
2
d
p
2
π
≥
0
, from which the stability result readily follows.Corollary 2.
Theϑ-method applied to (21) is stable if(37)
0
≤
ϑ
<
1
2
,
d
≤
1
1
-
2
ϑ
1
2
C
1
2
P
,
1
2
≤
ϑ
≤
1
,with d
,
P, and C defined in Theorem 1.Proof.
TakeW
(
P
)
=
1
/
2 and repeat the proof of Theorem 1.To conclude, the standard stability result for theϑ-method applied to (26) follows, if we substitute P
=
0 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 Table1 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 ϑ
=
1
/
2, on a grid with N
=
127 grid points and with (dimensionless) time step Δ
t
=
10
-
3, and subsequently converted to full dimensional form. Numerical solutions for the SCF scheme are similar and are not included.Table 1
Parameters values of the groundwater model.
Parameter
Value
Unit
a
1.5
×
1
0
-
1
m2 year−1
D
1.6
×
1
0
-
2
m2 year−1
κ
4.0
×
1
0
-
3
m year−1
S
S
2.6
×
1
0
-
3
m-
1
L
3.0
×
1
0
1
m
H
1.5
×
1
0
2
m
h
0
1.5
×
1
0
2
m
m
f
1
0
-
12
kgm
-
4
T
5
year
p
0
0,0
,
7.5
×
1
0
2
m year−1
p
w
—,1
0
2
,
1
0
2
m
x
p
—, —,21
m
m
0
0,7.5
×
1
0
3
,
3.0
×
1
0
3
m year−1
m
w
—,1
0
2
,
1
0
2
m
x
m
—,0,9
m
c
0
3.7
×
1
0
-
8
kg m−3
c
w
0.14
m−2
x
c
15
mThe choice ofN and Δ
t is justified by the grid validation of the numerical solution presented in Table 2. In this table we present the numerical approximations of c and h at x
=
(
1
/
2
)
L and at time t
=
T computed with (dimensionless) time step Δ
t
=
10
-
3 on a sequence of grids, to verify the dependence of the numerical solution on the grid. Clearly, for N
≥
127 the numerical solution is practically independent of the grid.Table 2
A grid validation for bothc and h. The values are computed for x
j
=
(
1
/
2
)
L and t
n
=
T.
N
No flow
Injecting
Injecting-pumping
c
(
(
1
/
2
)
L
,
T
)
h
(
(
1
/
2
)
L
,
T
)
c
(
(
1
/
2
)
L
,
T
)
h
(
(
1
/
2
)
L
,
T
)
c
(
(
1
/
2
)
L
,
T
)
h
(
(
1
/
2
)
L
,
T
)
2
3
+
1
3.64
E
-
08
1.50
E
02
8.30
E
-
09
2.56
E
02
3.30
E
-
08
1.57
E
02
2
4
+
1
3.62
E
-
08
1.50
E
02
2.07
E
-
08
2.05
E
02
1.94
E
-
08
1.88
E
02
2
5
+
1
3.60
E
-
08
1.50
E
02
2.63
E
-
08
1.98
E
02
2.28
E
-
08
1.86
E
02
2
6
+
1
3.60
E
-
08
1.50
E
02
2.72
E
-
08
1.98
E
02
2.36
E
-
08
1.86
E
02
2
7
+
1
3.60
E
-
08
1.50
E
02
2.74
E
-
08
1.98
E
02
2.38
E
-
08
1.86
E
02
2
8
+
1
3.60
E
-
08
1.50
E
02
2.74
E
-
08
1.98
E
02
2.38
E
-
08
1.86
E
02To verify the order of convergence of the SCF and TCF semidiscretizations, we repeatedly computec
(
1
/
2
)
L
,
T by successively halving the grid size and apply Richardson extrapolation [19]. For the (dimensionless) time step we take Δ
t
=
10
-
6. Let C
(
Δ
x
) denote the numerical approximation of c
(
(
1
/
2
)
L
,
T
) computed with grid size Δ
x. If the semidiscretization is pth order convergent, we have(38)
2
p
≈
C
Δ
x
/
2
-
C
Δ
x
C
Δ
x
/
4
-
C
Δ
x
/
2
≕
r
Δ
x
. We will present r
(
Δ
x
) 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 Figure3. 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.Figure 3
No-flow scenario: lodide concentration.
(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 Figure4. The source term s
h is virtually 0 outside a very small layer near x
=
0 and this implies that in the steady limit for t
→
∞ the groundwater level tends to a linear profile, and, consequently, the velocity q approaches a constant. Since the diffusion coefficient D 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 x
=
L is presented in Figure 5(a), confirming that all polluted water is removed at the right end of the domain.Figure 4
Injecting scenario: lodide concentration and groundwater level. The red profiles in (b) and (d) denote the initial conditions.
(a)
(b)
(c)
(d)Figure 5
Injecting scenario: (a) a zoom of the lodide concentration nearx
=
L and (b) total pollution relative to the initial pollution.
(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 Figure5(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 Table3 we present the r
(
Δ
x
) values for both the SCF and TCF schemes as a function of the number of grid points N
=
L
/
Δ
x
+
1. Clearly, both schemes show second-order convergence for Δ
x
→
0. In addition, we have included the (in absolute value) maximum Péclet number P
m
a
x
=
max
j
,
n
⁡
P
¯
j
+
1
/
2
n, indicating that advection is dominant on the coarser grids.Table 3
Injecting scenario:r
(
Δ
x
) values as a function of N.
N
SCF
TCF
P
m
a
x
2
5
+
1
1.18
11.77
51.0
2
6
+
1
1.52
-
1.86
25.6
2
7
+
1
2.24
1.40
12.8
2
8
+
1
3.18
2.60
6.4
2
9
+
1
3.74
3.31
3.2
2
10
+
1
3.93
3.67
1.6


### 5.3. Injecting-Pumping Scenario

In this scenario we add pumping to the model. We inject water just left of the concentration plume atx
=
x
m and pump polluted water out at the right of the plume at x
=
x
p. 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 0 leading to a linear decreasing groundwater level h and a concentration profile c propagating at constant velocity, just as in the injecting scenario. Third, near x
=
x
p, the source term s
h
<
0 resulting in the flattening of the groundwater level h and a decrease of the velocity q. Thus, left of the pumping peak, the pollutant moves faster than at the right, and since m
>
p, 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 x
=
L is given in Figure 7(a).Figure 6
Injecting-pumping scenario: lodide concentration and groundwater level. The red profiles in (b) and (d) denote the initial conditions.
(a)
(b)
(c)
(d)Figure 7
Injecting-pumping scenario: (a) a zoom-in of the lodide concentration nearx
=
L. (b) Total pollution relative to the initial pollution.
(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 Figure7(b).Ther
Δ
x 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 Δ
x
→
0. 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.Table 4
Injecting-pumping scenario:r
(
Δ
x
) values as a function of N.
N
SCF
TCF
P
m
a
x
2
5
+
1
-
5.55
6.61
32.4
2
6
+
1
0.29
-
7.58
16.3
2
7
+
1
1.51
0.79
8.2
2
8
+
1
2.85
2.38
4.1
2
9
+
1
3.63
3.21
2.0
2
10
+
1
3.94
3.62
1.0


## 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 Figure3. 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.Figure 3
No-flow scenario: lodide concentration.
(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 Figure4. The source term s
h is virtually 0 outside a very small layer near x
=
0 and this implies that in the steady limit for t
→
∞ the groundwater level tends to a linear profile, and, consequently, the velocity q approaches a constant. Since the diffusion coefficient D 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 x
=
L is presented in Figure 5(a), confirming that all polluted water is removed at the right end of the domain.Figure 4
Injecting scenario: lodide concentration and groundwater level. The red profiles in (b) and (d) denote the initial conditions.
(a)
(b)
(c)
(d)Figure 5
Injecting scenario: (a) a zoom of the lodide concentration nearx
=
L and (b) total pollution relative to the initial pollution.
(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 Figure5(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 Table3 we present the r
(
Δ
x
) values for both the SCF and TCF schemes as a function of the number of grid points N
=
L
/
Δ
x
+
1. Clearly, both schemes show second-order convergence for Δ
x
→
0. In addition, we have included the (in absolute value) maximum Péclet number P
m
a
x
=
max
j
,
n
⁡
P
¯
j
+
1
/
2
n, indicating that advection is dominant on the coarser grids.Table 3
Injecting scenario:r
(
Δ
x
) values as a function of N.
N
SCF
TCF
P
m
a
x
2
5
+
1
1.18
11.77
51.0
2
6
+
1
1.52
-
1.86
25.6
2
7
+
1
2.24
1.40
12.8
2
8
+
1
3.18
2.60
6.4
2
9
+
1
3.74
3.31
3.2
2
10
+
1
3.93
3.67
1.6


## 5.3. Injecting-Pumping Scenario

In this scenario we add pumping to the model. We inject water just left of the concentration plume atx
=
x
m and pump polluted water out at the right of the plume at x
=
x
p. 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 0 leading to a linear decreasing groundwater level h and a concentration profile c propagating at constant velocity, just as in the injecting scenario. Third, near x
=
x
p, the source term s
h
<
0 resulting in the flattening of the groundwater level h and a decrease of the velocity q. Thus, left of the pumping peak, the pollutant moves faster than at the right, and since m
>
p, 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 x
=
L is given in Figure 7(a).Figure 6
Injecting-pumping scenario: lodide concentration and groundwater level. The red profiles in (b) and (d) denote the initial conditions.
(a)
(b)
(c)
(d)Figure 7
Injecting-pumping scenario: (a) a zoom-in of the lodide concentration nearx
=
L. (b) Total pollution relative to the initial pollution.
(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 Figure7(b).Ther
Δ
x 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 Δ
x
→
0. 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.Table 4
Injecting-pumping scenario:r
(
Δ
x
) values as a function of N.
N
SCF
TCF
P
m
a
x
2
5
+
1
-
5.55
6.61
32.4
2
6
+
1
0.29
-
7.58
16.3
2
7
+
1
1.51
0.79
8.2
2
8
+
1
2.85
2.38
4.1
2
9
+
1
3.63
3.21
2.0
2
10
+
1
3.94
3.62
1.0


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



---

*Source: 1945958-2016-08-24.xml*