# New Second-Order Finite Difference Scheme for the Problem of Contaminant in Groundwater Flow

Authors: Quanyong Zhu; Quanxiang Wang; Ju Fu; Zhiyue Zhang
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We develop a new efficient second-order finite difference scheme for
two-dimensional problem of contaminant in groundwater flow. Theoretical analysis shows that
the scheme is second-order convergence in theL2 norm and is unconditionally stable. Numerical
results demonstrate that the error measures of the second-order scheme are several times smaller
than those of the standard implicit finite difference scheme.


---

## Body

## 1. Introduction

Prediction of the movement of contaminants in groundwater flows by models has been given increased emphasis in recent years. Groundwater is often contaminated by, for example, the sewage out of factories or mines and the chemical fertilizer and pesticide in agriculture. Groundwater protection is an issue with both economic and social significant. Mathematical models and numerical methods have been used extensively to simulate the movement of the contaminated groundwater. Analytical solutions can be derived only for a few classical models that are not suitable for complex situations normally encountered in the field. Consequently, the development of numerical solutions is required. During the last two decades, a variety of numerical approaches were presented; the technique of intimating the movement of groundwater flow is improved greatly, see [1–4] and the references therein. Dillon concludes more completely concludes many mathematical models and numerical methods for solving the groundwater problem [1]. Sun applied a kind of numerical methods to simulate the movement of contaminants in groundwater [2]. Li and Jiao give the analytical solutions of tidal groundwater flow in coastal two-aquifer system [3]. Pelovska improved an explicit scheme for age-dependent population model with spatial diffusion [4].In this paper, we give preliminary results for the numerical approximation of the following equation:(1.1)∂C∂t+∂(vLC)∂x+∂(vTC)∂y-[∂∂x(DL∂C∂x)+∂∂y(DT∂C∂y)]=CsQsθb,∀(x,y,t)∈Ω×[0,T],C(x,y,0)=C0(x,y),(x,y)∈Ω,C=g1(x,y,t),(x,y,t)∈Γ1×[0,T],C=g2(x,y,t),(x,y,t)∈Γ2×[0,T],
where C is concentration of contaminants, vL and vT are the longitudinal and transversal velocity components, DL and DT are the longitudinal and transversal dispersion coefficients, Cs is the concentration of the solute or sink fluid, Qs is the volume flow of the source or sink fluid per unite aquifer area, θ is volumetric fraction of the fluid phase, b is thickness of the aquifer, and C0(x,y) is the initial condition. In addition, let region of interest be Ω=[0,a]×[0,b] and the boundary of Ω be ∂Ω, where ∂Ω=Γ1⋃‍Γ2, Γ1={(x,y)|x=0,y∈[0,b] or y=0, x∈[0,a]} and Γ2={(x,y)|x=1, y∈[0,b] or y=1, x∈[0,a]}. In this paper, we assume that g1(x,y), g2(x,y) are known smooth functions, and θ and b are constants.Equation (1.1) is employed widely in the problem of contaminant in groundwater flow, or the water flow with any chemical solute. In general, the analytical solution for the above problem is not available. Many numerical methods can be used to solve (1.1). Wang and Zhang used finite volume element method to discretize the air pollution model [5, 6]. Borggaard and coauthors examine two-level finite element approximation schemes applied to the Navier-Stokes equations with r-Laplacian subgrid scale viscosity [7]. Minãmbres and De La Sen apply numerical methods of acceleration of the convergence to the case of one parameter in adaptive control algorithms [8].High performance algorithms are generally associated with large stencils, which increase the band-width of the resulting matrix and lead to a large number of arithmetic operations, especially for higher-dimensional problems. To obtain satisfactory high-performance numerical results with reasonable computational cost, many efforts have been made to develop efficient schemes based on Crank-Nicolson scheme. Deng and Zhang propose a new high-order algorithm for a class of nonlinear evolution equations, and the algorithm is based on Crank-Nicolson finite difference scheme [9]. Dai and Nassar design an unconditionally stable finite difference algorithm which is also based on Crank-Nicolson finite difference scheme [10].The aim of the present paper is to improve the accuracy in the temporal direction. We propose a second-order scheme which based on centered Crank-Nicolson finite difference scheme [11–16]. For the advection term of (1.1), we employ the average of the time level of n+1 and n, which guarantee that the discrete scheme for (1.1) is unconditionaly stable and second-order accurate in space and time. To the best of the author's knowledge, the method to deal with the advection term in this paper has not been studied. The numerical experiments demonstrate that the method is efficient and robust with respect to mesh refinement and the time step size. In order to illustrate that the scheme in this paper is efficient and exercisable, we carry out the standard implicit center finite difference scheme (ICFDS) for the equation.The paper is organized as follows. In Section2, we introduce some notations and describe our finite difference schemes. Stability analysis for initial value is obtained in Section 3. L2 error estimates are derived in Section 4 and numerical experiments are given in Section 5. Finally, conclusions are drawn in Section 6.Throughout this paper,K and ε denote generic positive constants and may be different at different places.


## 2. Second-Order Difference Scheme

For the presentation of our finite difference method, we first introduce some notations which will be used later. We denote temporal increment byΔt. For the spatial approximation, we take the step sizes hx=a/M and hy=b/N, respectively. In this way, the spatial nodes can be denoted by (xi,yj), xi=ihx, i=0,1,…,M; yj=jhy, j=0,1,…,N. Next, we denote the following:
(2.1)δxui,jn=ui+1,jn-ui,jnhx,δx¯ui,jn=δxui-1,jn=ui,jn-ui-1,jnhx,δx¯(DLδx)ui=1hx[(DL)i+1/2δxui-(DL)i-1/2δxui-1],
and similarly for the y direction. We set F(x,y,t)=CsQs/θb.We give the following assumptions:(H1)
forall(x,y,t)∈Ω×[0,T], 0<vL≤K1, 0<vT≤K2, 0<K3≤DL≤K4, 0 < K5 ≤ DT ≤ K6.(H2)
||C||L∞(0,T,L∞(Ω))≤K.We denote{Cijn} as the solutions of (1.1), and in terms of Taylor formula, we can obtain
(2.2)Cijn+1-CijnΔt+12[(δx+δx¯)(vLC)ijn+12+(δx+δx¯)(vLC)ijn2]+12[(δy+δy¯)(vTC)ijn+12+(δy+δy¯)(vTC)ijn2]-δx¯(DL)n+1/2δxCijn+1+Cijn2-δy¯(DT)n+1/2δyCijn+1+Cijn2=Fijn+1/2+Rijn+1/2,i=1,2,…,M-1,j=1,2,…,N-1,C(x,y,0)=C0(x,y),(x,y)∈Ω,C=g1(x,y,t),(x,y,t)∈Γ1×[0,T],C=g2(x,y,t),(x,y,t)∈Γ2×[0,T],
where Fn+(1/2)=F(xi,yj,tn+1/2), Rijn+1/2 truncation errors, and Rijn+1/2=O(Δt2+hx2+hy2).We denote{Uijn} as numerical solutions of (1.1), so the centered second-order finite difference scheme is:
(2.3)Uijn+1-UijnΔt+12[(δx+δx¯)(vLU)ijn+12+(δx+δx¯)(vLU)ijn2]+12[(δy+δy¯)(vTU)ijn+12+(δy+δy¯)(vTU)ijn2]-δx¯(DL)n+1/2δxUijn+1+Uijn2-δy¯(DT)n+1/2δyUijn+1+Uijn2=12(Fijn+Fijn+1),i=1,2,…,M-1,j=1,2,…,N-1,U(x,y,0)=C0(x,y),(x,y)∈Ω,U=g1(x,y,t),(x,y,t)∈Γ1×[0,T],U=g2(x,y,t),(x,y,t)∈Γ2×[0,T].The existence and uniqueness of the solution of scheme (2.3) are easily known by the positive-definite property.


## 3. Stability Analysis

In this section, for convenience, we letF=0. We employ freezing coefficient method to show the stability of scheme (2.3). We let vL,vT,DL,DT and constants be without loss of generality, we denote them by b1,b2,d1, and d2, so scheme (2.3) becomes
(3.1)Uijn+1-UijnΔt+b12(Ui+1,jn+1-Ui-1,jn+1+Ui+1,jn-Ui-1,jn2hx)+b22(Ui,j+1n+1-Ui,j-1n+1+Ui,j+1n-Ui,j-1n2hy)-d12hx2(Ui+1,jn+1-2Ui,jn+1+Ui-1,jn+1+Ui+1,jn-2Ui,jn+Ui-1,jn)-d22hy2(Ui,j+1n+1-2Ui,jn+1+Ui,j-1n+1+Ui,j+1n-2Ui,jn+Ui,j-1n)=0.LetUijn=wne-1(iξxhx+jξyhy), the amplification factor G(ξx,ξy,Δt)=wn+1/wn, for stability, it has to satisfy the relation |G(ξx,ξy,Δt)|≤1. By substituting the expressions of Uijn+1 and Uijn into the equation above, we get the amplification factor as follows:
(3.2)G=1-4λ1d1sin2(ξxhx/2)-4λ2d2sin2(ξyhy/2)--1(λ1b1hxsinξxhx+λ2b2hysinξyhy)1+4λ1d1sin2(ξxhx/2)+4λ2d2sin2(ξyhy/2)+-1(λ1b1hxsinξxhx+λ2b2hysinξyhy),
where λ1=Δt/hx2, λ2=Δt/hy2. With a simple calculation we have |G|<1. Hence, the following result can be obtained.Theorem 3.1.
LetU0 be the solution of scheme (2.3) at time level 0,  Un+1 be the solution of scheme (2.3) at time level n+1, e=C-U, then there exists a positive constant K, such that ||en+1||2≤K||e0||2,  for (n+1)Δt≤T, and K is independent of Δt and hx,hy. It follows that scheme (2.3) is unconditionally stable for initial value.


## 4. Error Estimate

In this section, we give theL2 error estimate for scheme (2.3). First, we define inner product in Ω: (u,v)=∑i,juijvijhxhy For all grid functions eij,
(4.1)∥e∥={∑i=0M‍∑j=0N‍eij2hxhy}1/2,|e|1=[(δx¯e,δx¯e)+(δy¯e,δy¯e)]1/2.Sete=C-U, the difference between (2.2) and (2.3) is
(4.2)eijn+1-eijnΔt+12[(δx+δx¯)(vLe)ijn+12+(δx+δx¯)(vLe)ijn2]+12[(δy+δy¯)(vTe)ijn+12+(δy+δy¯)(vTe)ijn2]-δx¯(DL)n+1/2δxeijn+1+eijn2-δy¯(DT)n+1/2δyeijn+1+eijn2=Fijn+1/2-12(Fijn+Fijn+1)+Rijn+1/2.i=1,2,…,M-1,j=1,2,…,N-1.Employing the inner product with(en+1+en)/2 in both sides of error equation (4.2), it follows that
(4.3)(en+1-enΔt,en+1+en2)+12((δx+δx¯)(vLe)n+12+(δx+δx¯)(vLe)n2,en+1+en2)+12((δy+δy¯)(vTe)n+12+(δy+δy¯)(vTe)n2,en+1+en2)-12(δx¯(DL)n+1/2δxen+1+en2,en+1+en2)-12(δy¯(DT)n+1/2δyen+1+en2,en+1+en2)=(Fn+1/2-12(Fn+Fn+1)),(en+1+en2)+(Rn+1/2,en+1+en2).In the following, we will carry on the estimation to each item in (4.3); before estimating (4.3), we give a lemma.Lemma 4.1.
Consider(δxun,vn)=-(δx¯vn,un)+hxhy∑j=0N((uM,jnvM,jn-u0,jnv0,jn)/hx), and similarly for the y direction.Proof.
From summation by parts formula, we can get the results easily.Now, we talk about each item of (4.3),
(4.4)(en+1-enΔt,en+1+en2)=12Δt(∥en+1∥2-∥en∥2).By lemma4.1 and (H1), with employing ε inequality, we can obtain
(4.5)|((δx+δx¯)(vLe)n+12+(δx+δx¯)(vLe)n2,en+1+en2)|+|((δy+δy¯)(vTe)n+12+(δy+δy¯)(vTe)n2,en+1+en2)|≤K[|(en+1+en,δx¯(en+1+en))|+|(en+1+en,δx(en+1+en))|+|(en+1+en,δy¯(en+1+en))|+|(en+1+en,δy(en+1+en))|]≤K(∥en+1∥2+∥en∥2)+ε|en+1+en2|12.
Employing summation by parts formula, and noting that DL is bounded, we can get
(4.6)-12(δx¯(DL)n+1/2δxen+1+en2,en+1+en2)≥K3(δx¯en+1+en2,δx¯en+1+en2),
and similarly for the y direction.DenotingL=Dmin=min{K3,K5}, then
(4.7)-12(δx¯(DL)n+1/2δxen+1+en2,en+1+en2)-12(δy¯(DT)n+1/2δyen+1+en2,en+1+en2)≥L|en+1+en2|12.In terms of Taylor formula andε inequality, we can obtain
(4.8)|(Fn+1/2-12(Fn+Fn+1),en+1+en2)|≤K(∥en+1∥2+∥en∥2+Δt4).(Rn+1/2,en+1+en2)≤K(Δt2+hx2+hy2)2+ε(∥en+1∥2+∥en∥2).Substituting (4.4)–(4.8) into (4.3) we get
(4.9)12Δt(∥en+1∥2-∥en∥2)+L|en+1+en2|12≤K(∥en+1∥2+∥en∥2)+ε|en+1+en2|12+K(Δt2+hx2+hy2)2.We letε<L in (4.9), multiplying both sides by 2Δt, and sum over n; using discrete Gronwall lemma, we can obtain
(4.10)max0≤n≤[T/Δt]∥Cn-Un∥=max0≤n≤[T/Δt]∥en∥≤K(Δt2+hx2+hy2).Hence,L2 error estimate is derived as follows:Theorem 4.2.
LetC and U be the accurate solution of (1.1) and the solution of second-order finite difference scheme, respectively. If (H1)-(H2) hold, then ||Cn-Un||=O(Δt2+hx2+hy2), for all n(0≤n≤[T/Δt]).


## 5. Numerical Experiments

In this part, we will give two experiments to test our second-order finite difference scheme. The numerical results will be presented to illustrate the efficiency and order of accuracy of the algorithm. In the following numerical experiments, we use the same number of uniform subintervals in bothx and y directions. ECij=|Cij-Uij|, RCij=|Cij-Uij/Cij|, ||C-U||=[∑ijN(Cij-Uij)2h2]1/2 denote absolute error, relative error, and L2 error, respectively.Example 5.1.
The equations to be solved are(5.1)∂C∂t+∂C∂x+∂C∂y-12∂2C∂x2-12∂2C∂y2=2etcos(x+y)-2etsin(x+y),(x,y)∈[0,1]×[0,1],t∈[0,T],C(x,y,0)=cos(x+y),(x,y)∈[0,1]×[0,1],C(0,y,t)=etcosy,C(x,0,t)=etcosx,C(x,1,t)=etcos(x+1),C(1,y,t)=etcos(1+y),t∈[0,T].The exact solutions of this test problem areC(x,t)=etcos(x+y) and we take T=1 in the experiment. We take different spatial steps h=1/4,1/8,1/16, and 1/32, respectively, in the numerical experiment. Table 1 demonstrates the maximum absolute error, maximum relative error, and L2 error between approximation solution and exact solution. For all of the above cases, we take Δt=h in the numerical experiment. Figures 1 and 2 are pictures of absolute errors by employing second-order finite difference scheme at T=1, in Figure 1, h=0.1, Δt=0.1, and in Figure 2, h=0.05, Δt=0.05. Figure 3 presents the approximation order of C in L2 norm. It is shown that approximation orders of C in different norms are almost 2 orders, respectively. This is the same as the theoretical results. Transverse axis denotes lgh, and longitudinal axis denotes lg||C-U|| in Figure 3.Table 1
The comparison of maximum absolute error, maximum relative error, andL2 error when Δt=h.
The step size in space
h
=
1
/
4
h
=
1
/
8
h
=
1
/
16
h
=
1
/
32
max
E
C
3.9203
e
-
3
1.0783
e
-
3
2.7107
e
-
4
6.8318
e
-
5
max
R
C
2.0388
e
-
2
5.1581
e
-
3
9.8663
e
-
3
2.5216
e
-
3
|
|
C
-
U
|
|
2.1257
e
-
3
5.4288
e
-
4
1.3603
e
-
4
3.4020
e
-
5
max
E
C
/
h
2
0.06273
0.06902
0.06939
0.06996
|
|
C
-
U
|
|
/
h
2
0.03401
0.03474
0.03483
0.03484Figure 1
Absolute error atT=1, h=0.1, and Δt=0.1.Figure 2
Absolute error atT=1, h=0.05, andΔt=0.05.Figure 3
Approximation order ofC in L2 norm, Example 5.1.From Table1, we can see that maximum absolute error, maximum relative error, and L2 error of C for each case in h=1/4,1/8,1/16, and 1/32 are second-order accuracy in space and time, respectively. With h decreasing, the accuracies of maximum absolute error, maximum relative error, and L2 error of C are increasing. Roughly speaking, we can adopt larger time step to solve this class of equation if we need the error accuracy 10-2. This would save computational loads and satisfy our requirement. We can also see from Table 1 that the error denoted by maxEC shows a second-order accuracy. This is illustrated by the fact that the ratios of maxEC/h2 and ||C-U||/h2 keep roughly a constant as the computational grid is being refined. From Figures 1 and 2, we can intuitively see that the degree of numerical solutions is approximating the exact solutions in different grid points. The slope of the line in Figure 3 is nearly 2, which shows second-order accuracy.Example 5.2.
We consider the following problem:(5.2)∂C∂t+∂∂x(Csinxsiny)+∂∂y(Csinxsiny)-2∂2C∂x2-2∂2C∂y2=2e-2tsinxsiny(1+cosxsiny+sinxcosy),(x,y,t)∈[0,π]×[0,π]×[0,T],C(x,y,0)=sinxsiny,(x,y)∈[0,π]×[0,π],C(x,0,t)=C(0,y,t)=0,(x,π,t)=C(π,y,t)=0,t∈[0,T].The exact solutions of this test problem areC(x,y,t)=e-2tsinxsiny; we take T=2. In this example, we use implicit center finite difference scheme (ICFDS) to compare with second-order finite difference scheme(SOFDS), and the implicit center finite difference scheme in this example is
(5.3)Uijn+1-UijnΔt+sinxi+1sinyjUi+1,jn+1-sinxi-1sinyjUi-1,jn+12h+sinxisinyj+1Ui,j+1n+1-sinxisinyj-1Ui,j-1n+12h-2Ui+1,jn+1-2Uijn+1+Ui-1,jn+1h2-2Ui,j+1n+1-2Uijn+1+Ui,j-1n+1h2=2e-2nΔtsinxisinyj(1+cosxisinyj+sinxicosyj).In Table2, we compare second-order finite difference scheme and implicit center finite difference Scheme (5.3) in maximum absolute error of solutions and L2 error at T=2. We let Δt=h in second-order finite difference scheme, and Δt=h2 in the implicit center finite difference scheme. In Table 3, we give maximum absolute error, maximum relative error of solutions, and L2 error at T=2. Figures 4 and 5 are the pictures of absolute error by employing second-order finite difference scheme, in Figure 4, h=π/20, Δt=0.1, and in Figure 5, h=π/40, Δt=0.1. Figure 6 presents the approximation order of C in L2 norm.Table 2
The comparison of SOFDS and ICFDS in different partitions.
The step size in space
max
E
C (SOFDS)
max
E
C (ICFDS)
|
|
C
-
U
|
| (SOFDS)
|
|
C
-
U
|
| (ICFDS)
h
=
π
/
10
1.2324
e
-
3
5.2935
e
-
3
1.8081
e
-
3
8.3026
e
-
3
h
=
π
/
20
2.8791
e
-
4
1.2583
e
-
3
4.2181
e
-
4
1.9726
e
-
3
h
=
π
/
40
7.1142
e
-
5
3.1205
e
-
4
1.0366
e
-
4
4.8671
e
-
4
Rater
2.04
2.03
2.04
2.03
BoundM0
0.0119
0.0517
0.0174
0.0810Table 3
The comparison of the maximum absolute error, the maximum relative error, andL2 error at T=2.
The step sizes in space and time
max
E
C
max
R
C
|
|
C
-
U
|
|
|
|
C
-
U
|
|
/
h
2
h
=
π
/
10
,
Δ
t
=
0.4
1.8944
e
-
3
1.2227
e
-
1
2.9256
e
-
3
0.02964
h
=
π
/
20
,
Δ
t
=
0.2
4.8450
e
-
4
3.4301
e
-
2
7.4876
e
-
4
0.03035
h
=
π
/
40
,
Δ
t
=
0.1
1.2189
e
-
4
8.5131
e
-
3
1.8837
e
-
4
0.03053Figure 4
Absolute error whenT=2, h=π/20, and Δt=0.1.Figure 5
Absolute error whenT=2, h=π/40, and Δt=0.1.Figure 6
Approximation order ofC in L2 norm, Example 5.2.We use two constants: convergence rater and bound constant M0, that is,
(5.4)∥C-U∥=M0*hr,
which will be obtained by the least square method with error in norm ||·|| [17]. The bound constants in error estimates for our SOFDS scheme are 0.0119 and 0.0174 in the maximum norm and in the L2 norm, respectively. Meanwhile, the bound constants for the ICFDS scheme are 0.0517 and 0.0810 in the maximum norm and in the L2 norm, respectively. This shows that the bound constants in error estimates for our SOFDS scheme are much smaller than the ones for the ICFDS scheme even when they have the same order.In this example, we obtain numerical results of maximum absolute error, maximum relative error, andL2 error of C in Table 3. When grid size is decreasing, maximum absolute error, maximum relative error, and L2 error of C are also decreasing. Those numerical results confirm the theoretical analysis presented in the previous section. This is the example with exact solutions against which we can compare the numerical solution to prove the efficiency and order of accuracy of our algorithm in both the spatial and temporal directions. From Table 3, it is clear that grid size h reduces a half, and the computational accuracy increases one magnitude. In order to illustrate that the scheme in this paper is efficient and exercisable, we carry out the implicit center finite difference Scheme (5.3) for the equation. Numerical results are shown in Table 2. From Table 2, we can see that errors of our scheme in different norms for the same spacial step are smaller than the ones of implicit center finite difference Scheme (5.3) in different norms. It is shown from Table 3 that maximum absolute error, maximum relative error, and L2 error of C for each case are second-order accurate in spacial and temporal directions, respectively. This is illustrated by the fact that the ratios of ||C-U||/h2 keep roughly a constant as the computational grid is being refined. From Figures 4 and 5, we can intuitively see that the degree of numerical solutions is approximating the exact solutions in different grid points. The slope of the line in Figure 6 is nearly 2, which shows second-order accuracy. This is consistent with the theoretical results.


## 6.   Conclusions

In this paper, we design an efficient second-order finite difference scheme for solving a class of groundwater problem. Theoretically analysis shows that the proposed scheme is unconditionally stable and second-order accurat in time and space. Numerical results given in Section5 confirm the theoretical results. Compared with standard implicit center finite difference scheme, the second-order finite difference scheme is more accurate.



---

*Source: 575493-2012-09-25.xml*