`Journal of Applied MathematicsVolume 2012, Article ID 575493, 13 pageshttp://dx.doi.org/10.1155/2012/575493`
Research Article

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

1Jiangsu Provincial Key Laboratory for Numerical Simulation of Large-Scale Complex Systems, School of Mathematical Sciences, Nanjing Normal University, Nanjing 210046, China
2School of Science, Lishui University, Lishui 323000, China

Received 5 April 2012; Revised 9 May 2012; Accepted 10 May 2012

Copyright © 2012 Quanyong Zhu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

#### 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 [14] 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: where is concentration of contaminants, and are the longitudinal and transversal velocity components, and are the longitudinal and transversal dispersion coefficients, is the concentration of the solute or sink fluid, is the volume flow of the source or sink fluid per unite aquifer area, is volumetric fraction of the fluid phase, is thickness of the aquifer, and is the initial condition. In addition, let region of interest be and the boundary of be , where , or , and , or , . In this paper, we assume that , are known smooth functions, and and 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 [1116]. For the advection term of (1.1), we employ the average of the time level of and , 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 Section 2, we introduce some notations and describe our finite difference schemes. Stability analysis for initial value is obtained in Section 3. 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, 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 . For the spatial approximation, we take the step sizes and , respectively. In this way, the spatial nodes can be denoted by , , ; , . Next, we denote the following: and similarly for the direction. We set .

We give the following assumptions: (H1) , , , , < . (H2) .

We denote as the solutions of (1.1), and in terms of Taylor formula, we can obtain where , truncation errors, and .

We denote as numerical solutions of (1.1), so the centered second-order finite difference scheme is:

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 let . We employ freezing coefficient method to show the stability of scheme (2.3). We let and constants be without loss of generality, we denote them by , and , so scheme (2.3) becomes

Let , the amplification factor , for stability, it has to satisfy the relation . By substituting the expressions of and into the equation above, we get the amplification factor as follows: where , . With a simple calculation we have . Hence, the following result can be obtained.

Theorem 3.1. Let be the solution of scheme (2.3) at time level 0, be the solution of scheme (2.3) at time level , , then there exists a positive constant , such that , for , and is independent of and . It follows that scheme (2.3) is unconditionally stable for initial value.

#### 4. Error Estimate

In this section, we give the error estimate for scheme (2.3). First, we define inner product in : For all grid functions ,

Set , the difference between (2.2) and (2.3) is

Employing the inner product with in both sides of error equation (4.2), it follows that

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 , and similarly for the direction.

Proof. From summation by parts formula, we can get the results easily.

Now, we talk about each item of (4.3),

By lemma 4.1 and (H1), with employing inequality, we can obtain Employing summation by parts formula, and noting that is bounded, we can get and similarly for the direction.

Denoting , then

In terms of Taylor formula and inequality, we can obtain

Substituting (4.4)–(4.8) into (4.3) we get

We let in (4.9), multiplying both sides by , and sum over ; using discrete Gronwall lemma, we can obtain

Hence, error estimate is derived as follows:

Theorem 4.2. Let and be the accurate solution of (1.1) and the solution of second-order finite difference scheme, respectively. If (H1)-(H2) hold, then , for all .

#### 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 both and directions. , , denote absolute error, relative error, and error, respectively.

Example 5.1. The equations to be solved are

The exact solutions of this test problem are and we take in the experiment. We take different spatial steps , and , respectively, in the numerical experiment. Table 1 demonstrates the maximum absolute error, maximum relative error, and error between approximation solution and exact solution. For all of the above cases, we take in the numerical experiment. Figures 1 and 2 are pictures of absolute errors by employing second-order finite difference scheme at , in Figure 1, , , and in Figure 2, , . Figure 3 presents the approximation order of in norm. It is shown that approximation orders of in different norms are almost 2 orders, respectively. This is the same as the theoretical results. Transverse axis denotes , and longitudinal axis denotes in Figure 3.

Table 1: The comparison of maximum absolute error, maximum relative error, and error when .
Figure 1: Absolute error at , , and .
Figure 2: Absolute error at , , and.
Figure 3: Approximation order of in norm, Example 5.1.

From Table 1, we can see that maximum absolute error, maximum relative error, and error of for each case in , and are second-order accuracy in space and time, respectively. With decreasing, the accuracies of maximum absolute error, maximum relative error, and error of are increasing. Roughly speaking, we can adopt larger time step to solve this class of equation if we need the error accuracy . This would save computational loads and satisfy our requirement. We can also see from Table 1 that the error denoted by shows a second-order accuracy. This is illustrated by the fact that the ratios of and 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:

The exact solutions of this test problem are ; we take . 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

In Table 2, we compare second-order finite difference scheme and implicit center finite difference Scheme (5.3) in maximum absolute error of solutions and error at . We let in second-order finite difference scheme, and in the implicit center finite difference scheme. In Table 3, we give maximum absolute error, maximum relative error of solutions, and error at . Figures 4 and 5 are the pictures of absolute error by employing second-order finite difference scheme, in Figure 4, , , and in Figure 5, , . Figure 6 presents the approximation order of in norm.

Table 2: The comparison of SOFDS and ICFDS in different partitions.
Table 3: The comparison of the maximum absolute error, the maximum relative error, and error at .
Figure 4: Absolute error when , , and .
Figure 5: Absolute error when , , and .
Figure 6: Approximation order of in norm, Example 5.2.

We use two constants: convergence rate and bound constant , that is, 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 norm, respectively. Meanwhile, the bound constants for the ICFDS scheme are 0.0517 and 0.0810 in the maximum norm and in the 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, and error of in Table 3. When grid size is decreasing, maximum absolute error, maximum relative error, and error of 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 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 error of for each case are second-order accurate in spacial and temporal directions, respectively. This is illustrated by the fact that the ratios of 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 Section 5 confirm the theoretical results. Compared with standard implicit center finite difference scheme, the second-order finite difference scheme is more accurate.

#### Acknowledgments

This project is supported by the National Natural Science Foundation of China (no.11071123), the Project of Graduate Education Innovation of Jiangsu Province (no. CXLX12_0388, CXZZ12_0382) and Scientific Research Fund of Zhejiang Provincial Education Department (no. Y201223607).

#### References

1. P. J. Dillon, “An analytical model of contaminant transport from diffuse sources in saturated porous media,” Water Resources Research, vol. 25, no. 6, pp. 1208–1218, 1989.
2. N.-Z. Sun, “Applications of numerical methods to simulate the movement of contaminants in groundwater,” Environmental Health Perspectives, vol. 83, pp. 97–115, 1989.
3. H. L. Li and J. J. Jiao, “Analytical solutions of tidal groundwater flow in coastal two-aquifer system,” Advances in Water Resources, vol. 25, no. 4, pp. 417–426, 2002.
4. G. Pelovska, “An improved explicit scheme for age-dependent population models with spatial diffusion,” International Journal of Numerical Analysis and Modeling, vol. 5, no. 3, pp. 466–490, 2008.
5. P. Wang and Z. Zhang, “Quadratic finite volume element method for the air pollution model,” International Journal of Computer Mathematics, vol. 87, no. 13, pp. 2925–2944, 2010.
6. Z. Zhang, “Error estimates for finite volume element method for the pollution in groundwater flow,” Numerical Methods for Partial Differential Equations, vol. 25, no. 2, pp. 259–274, 2009.
7. J. Borggaard, T. Iliescu, and J. P. Roop, “Two-level discretization of the Navier-Stokes equations with r-Laplacian subgridscale viscosity,” Numerical Methods for Partial Differential Equations, vol. 28, no. 3, pp. 1056–1078, 2012.
8. J. J. Miñambres and M. de la Sen, “Application of numerical methods to the acceleration of the convergence of the adaptive control algorithms. The one-dimensional case,” Computers and Mathematics with Applications Part A, vol. 12, no. 10, pp. 1049–1056, 1986.
9. D. Deng and Z. Zhang, “A new high-order algorithm for a class of nonlinear evolution equation,” Journal of Physics A, vol. 41, pp. 1–17, 2008.
10. W. Z. Dai and R. Nassar, “An unconditionally stable finite difference scheme for solving a 3D heat transport equation in a sub-microscale thin film,” Journal of Computational and Applied Mathematics, vol. 145, no. 1, pp. 247–260, 2002.
11. W. Q. Wang, “The alternating segment Crank-Nicolson method for solving convection-diffusion equation with variable coefficient,” Applied Mathematics and Mechanics, vol. 24, no. 1, pp. 29–38, 2003.
12. M. Ganesh and K. Mustapha, “A Crank-Nicolson and ADI Galerkin method with quadrature for hyperbolic problems,” Numerical Methods for Partial Differential Equations. An International Journal, vol. 21, no. 1, pp. 57–79, 2005.
13. Z. Zhang, “An economical difference scheme for heat transport equation at the microscale,” Numerical Methods for Partial Differential Equations, vol. 20, no. 6, pp. 855–863, 2004.
14. Z. Z. Sun, Numerical Methods for Partial Differential Equations, Science Press, 2005.
15. J. F. Lu and Z. Guan, Numerical Methods for Partial Differential Equations, Qinghua University Press, 2003.
16. W. Hundsdorfer and J. G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, vol. 33 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 2003.
17. D. Liang and W. D. Zhao, “An optimal weighted upwinding covolume method on non-standard grids for convection-diffusion problems in 2D,” International Journal for Numerical Methods in Engineering, vol. 67, no. 4, pp. 553–577, 2006.