ABSTRACT

The analytical solution of the convection diffusion equation is considered by two-dimensional Fourier transform and the inverse Fourier transform. To get the numerical solution, the Crank-Nicolson finite difference method is constructed, which is second-order accurate in time and space. Numerical simulation shows excellent agreement with the analytical solution. The dynamic visualization of the simulating results is realized on ArcGIS platform. This work provides a quick and intuitive decision-making basis for water resources protection, especially in dealing with water pollution emergencies.

Keywords:

Groundwater Pollution, Two-Dimensional Convection Diffusion Equation, Finite Difference Method, Visualization, Numerical Simulation

1. Introduction

Numerical simulations of groundwater pollution problems have 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, in particular, nitrate pollution of groundwater in rivers basins and agricultural watersheds are going from bad to worse all over the world [1] . Groundwater protection is an issue with both social and economic significant (e.g., Chen et al. [2] ), so to simulate the movement of the contaminated groundwater, lots of mathematical models have been applied extensively and the use of numerical simulations seems to be inevitable.

Generally, analytical solutions can not be derived for most classical models, consequently, the development of numerical solutions is required. S. Hasnain and M. Saqib originate results with finite difference schemes to approximate the solution of the classical Fisher Kolmogorov Petrovsky Piscounov (KPP) equation from population dynamics, their results show that Crank-Nicolson is very efficient and reliably numerical scheme for solving one-dimension fishers KPP equation [3] . During the last three decades, numerous transport problems have been solved numerical [4] , theoretical numerical models are necessary tools were presented [5] , the technique of intimating the movement of groundwater flow is improved greatly, see [6] [7] [8] and the references therein. Dillon gave many mathematical models and numerical methods for solving the groundwater problems [6] . Li and Jiao derived the analytical solutions of tidal groundwater flow in coastal two-aquifer system [7] . Sun applied a sort of numerical methods to simulate the movement of contaminants in groundwater [8] .

At present, the researches on groundwater pollution problems are mainly divided into two categories at home and abroad [9] [10] . On one hand, people generally study various discrete numerical schemes of mathematical models (i.e. the corresponding partial differential equations) for groundwater pollution, by doing lots of related works, the numerical solution is obtained and the convergence is analyzed, Lin, L. et al. [9] derived a simplified numerical model of groundwater and solute transport. On the other hand, people develop the simulation software of groundwater numerical [11] [12] [13] , particularly, numerical simulation software of groundwater system has been developed, which is widely used because of its modularity, visualization, interaction and diversification, such as, the popular GIS, which is used to provide the visualization methods and approaches of groundwater pollution diffusion simulation. Qin, R. et al. [11] presented a GIS-based software for forecasting pollutant drift on coastal water surfaces using fractional Brownian motion, it can be used to study on red tide drift. Valocchi, A. J. et al. [12] developed a series of interactive web simulation models to help students understand the coupled physical, chemical, and microbiological processes that affect the transport and fate of pollutants in groundwater. Li, J. et al. [13] mentioned that HYDRUS-1D software can simulate different concentrations of pollutants reaching the shallow aquifer under some vadose zone conditions, he presented a method for quantitative groundwater pollution assessment based on grey relational analysis (GRA).

A lot of research literatures about the former have been represented in China. In recent years, with the wide application of new technology and new method, many scholars have made innovations in theory and methodology, by means of the combination of theory and research direction of numerical model theory, they combine the theory of numerical model with the related acknowledges, so as to improve the reliability of simulation results [10] . However, most of the previous literatures only simulate the surface water. By studying the literature [11] [12] [13] , we found that the visual simulation of the actual problems of groundwater pollution is rare. If the above two were combined, it will not only in theory and application of important values, but also innovative, advanced and applied.

In book [14] , Kovarik, K. sets his sights on reviewing the whole group of numerical methods from the oldest (the finite differences method), and discusses the basic equations of a groundwater flow and of the transport of pollutants in a porous medium. Therefore, we would like to use the finite difference method to study the numerical simulation methods about mathematical models with seepage of groundwater pollution. In order to improve the accuracy in the temporal direction, we propose a second-order scheme which based on centered Crank- Nicolson finite difference scheme [15] - [20] . And we simulate the water pollution problems in a certain area and verify the validity and practicability of the model and its algorithm. Meanwhile, the dynamic visualization of simulation results is realized on ArcGIS platform. We hope that our work can provide an important basis for water pollution accident emergency response and decision-making, and can be used for environmental protection personnel to deal with water pollution emergencies.

The paper is organized as follows. In Section 2, we give the analytical solution of the Equations (1)-(5) by two-dimensional Fourier transform and the inverse Fourier transform. In Section 3, we introduce some notations, present the Crank-Nicolson finite difference discretization of the governing equation and derive the truncation error. Numerical experiment is given, exact solution comparisons with numerical solution are also discussed in Section 4. Visualization of simulation results based on GIS (ArcGIS figures) is presented in Section 5. Finally, conclusions and suggestions are drawn in Section 6.

2. Analytical Solution

In this paper, we consider the following physical problem. The leakage of the sewage pool in a paper mill causes the seepage of sewage, and the concentration of some substance in any point in the underground water is a function of space coordinate and time, i.e. $C\left(x\mathrm{,}y\mathrm{,}t\right)$ . We take an micro-body in the underground water, the concentration change of it is caused by two aspects: One is diffusion, including molecular diffusion and osmotic dispersion, another is the mass flux caused by the average liquid motion. In this problem, we assume that the seepage area is an infinite plane, and the groundwater flow is a one-dimensional one, the diffusion of pollutants is a two-dimensional dispersion, and the medium is a porous medium.

We take O as the coordinate origin (pollution source), take the infinite plane as the plane O-xy, the flow direction and the x axis direction are consistent. Then, our problem can be illustrated by the following two-dimensional parabolic equation with convection term

$\frac{\partial C}{\partial t}=D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}-v\frac{\partial C}{\partial x},\left(0<x<+\infty ,t>0\right),$ (1)

$C\left(x,y,0\right)=0,\left(x,y\right)\ne \left(0,0\right),$ (2)

${\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }nC\text{d}x\text{d}y=m,$ (3)

$\underset{x\to \pm \infty }{\lim }C\left(x,y,t\right)=0,\left(t>0\right),$ (4)

$\underset{y\to \pm \infty }{\lim }C\left(x,y,t\right)=0,\left(t>0\right).$ (5)

Problem (1)-(5) arises in the mathematical modeling of transport processes that exhibit diffusion, in which, the unknown C stands for the concentration of a solute, x and y are the horizontal coordinates, t is the time, $D_x$ and $D_y$ are the longitudinal and transversal dispersion coefficients respectively (namely, $D_x$ and $D_y$ are the aquifer transmissivity, subscripts x and y indicate the respective directions). m is the instantaneous injected solute mass per unit length of porous medium, v is the mean pore velocity, and n indicates the effective porosity. At the initial stage, we suppose, there is no contaminant in the river, and when $t\ge 0$ the concentration at $x=0$ remains at $C=C_0$ .

There are different approaches to solve two-dimensional parabolic equation, a series of analytical solutions derived from the basic physical principles have been presented which are mostly suitable under special boundary conditions [21] . In this part, we give the analytical solution of the Equations (1)-(5) by using two-dimensional Fourier transform and the inverse Fourier transform.

First of all, we give the definition of the two-dimensional Fourier transform and the inverse Fourier transform and some properties to be used in the following. The Fourier transform of $f\left(x\mathrm{,}y\right)$ is written to

$\widehat{f}\left({\xi }_1\mathrm{,}{\xi }_2\right)={\iint }_{R^2}f\left(x\mathrm{,}y\right)\cdot exp\left(-2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right)\text{d}x\text{d}y\mathrm{,}$ (6)

the definition of the inverse Fourier transform for $f\left(x\mathrm{,}y\right)$ as follows

$f\left(x\mathrm{,}y\right)={\iint }_{R^2}\widehat{f}\left({\xi }_1\mathrm{,}{\xi }_2\right)\cdot exp\left(2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right)\text{d}{\xi }_1\text{d}{\xi }_2\mathrm{,}$ (7)

and the derivative properties of the Fourier transform are

$\widehat{\left(\frac{\partial f}{\partial x}\right)}=2\text{π}i{\xi }_1\widehat{f},\widehat{\left(\frac{\partial f}{\partial y}\right)}=2\text{π}i{\xi }_2\widehat{f},$

$\widehat{\left(\frac{{\partial }^{2}f}{\partial x^2}\right)}={\left(2\text{π}i\right)}^2{\xi }_1^2\widehat{f},\widehat{\left(\frac{{\partial }^{2}f}{\partial y^2}\right)}={\left(2\text{π}i\right)}^2{\xi }_2^2\widehat{f}.$

We do Fourier transformation of x and y (this transformation has no effect on the independent variables, so there is no Fourier transformation for t, so it has nothing to do with the t to the Equation (1) firstly.

$\frac{\partial \widehat{C}}{\partial t}=-4{\text{π}}^2{D}_x{\xi }_1^2\widehat{C}-4{\text{π}}^2{D}_y{\xi }_2^2\widehat{C}-2{\text{π}}^{2}iv{\xi }_1\widehat{C},$ (8)

deal with the conditional Eqution (3), we have

${\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }C\text{d}x\text{d}y=\frac{m}{n}.$ (9)

Then, we do the Fourier transformation to the initial condition Eqution (2), by Eqution (9), we have

$\widehat{C}\left({\xi }_1,{\xi }_2,0\right)={\iint }_{R^2}\frac{m}{n}\delta \left(x,y\right)\cdot \exp \left(-2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right)\text{d}x\text{d}y,$ (10)

in which

$\delta \left(x,y\right)=\{\begin{array}{l}\infty ,\left(x,y\right)=\left(0,0\right)\\ 0,\text{other}\text{points}\end{array},{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\delta \left(x,y\right)\text{d}x\text{d}y=1,$ (11)

it is called the Dirac function, also called the generalized function. We denote

$g\left(x,y\right)=\delta \left(x,y\right)\cdot \exp \left(-2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right),$ (12)

it can be easily to proved

$\frac{m}{n}\cdot g\left(0,0\right)=\frac{m}{n}.$ (13)

Notice Eqution (8) is a separable equation, so we have

$\frac{\text{d}\widehat{C}}{\widehat{C}}=\left(-\left(4{\text{π}}^2{D}_x{\xi }_1^2+4{\text{π}}^2{D}_y{\xi }_2^2\right)-2\text{π}iv{\xi }_1\right)\text{d}t,$ (14)

integrate the both sides of the Eqution (14) about variable t, we have

$\ln \widehat{C}=\left(-\left(4{\text{π}}^2{D}_x{\xi }_1^2+4{\text{π}}^2{D}_y{\xi }_2^2\right)-2\text{π}iv{\xi }_1\right)t+C_0.$ (15)

where $C_0$ is a constant and it has nothing with other variables (the same below).

So, we get the general solution of Eqution (8)

$\widehat{C}=C_0\cdot \exp \left(\left(-\left(4{\text{π}}^2{D}_x{\xi }_1^2+4{\text{π}}^2{D}_y{\xi }_2^2\right)-2\text{π}iv{\xi }_1\right)t\right),$ (16)

substitute 0 for t in formula (16), using Eqution (10), we obtain

$\widehat{C}\left({\xi }_1,{\xi }_2,0\right)=C_0\cdot \exp \left(0\right)=\frac{m}{n},C_0=\frac{m}{n},$ (17)

in fact, by Equtions ((16) and (17)) we have the solution of Eqution (8)

$\widehat{C}=\frac{m}{n}\cdot \exp \left(\left(-\left(4{\text{π}}^2{D}_x{\xi }_1^2+4{\text{π}}^2{D}_y{\xi }_2^2\right)-2\text{π}iv{\xi }_1\right)t\right).$ (18)

Next, for Eqution (16), we do the inverse Fourier transform

$\begin{array}{c}C\left(x,y,t\right)={\iint }_{R^2}\widehat{C}\left({\xi }_1,{\xi }_2,t\right)\cdot \exp \left(2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right)\text{d}{\xi }_1\text{d}{\xi }_2\\ =\frac{m}{n}{\iint }_{R^2}\left(\exp \left(-\left(4{\text{π}}^2\left(D_x{\xi }_1^2+D_y{\xi }_2^2\right)\right)-2\text{π}iv{\xi }_1\right)t\cdot \exp \left(2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right)\right)\text{d}{\xi }_1\text{d}{\xi }_2\\ =\frac{m}{n}{\iint }_{R^2}\exp \left(\left(-4{\text{π}}^2\left(D_x{\xi }_1^2+D_y{\xi }_2^2\right)-2\text{π}iv{\xi }_1\right)t+\left(2\text{π}i\left(x{\xi }_1+y{\xi }_2\right)\right)\right)\text{d}{\xi }_1\text{d}{\xi }_2\\ =\frac{m}{n}{\iint }_{R^2}\exp \left(-4{\text{π}}^2{D}_x{\xi }_1^{2}t-4{\text{π}}^2{D}_x{\xi }_2^{2}t+2\text{π}i{\xi }_1\left(x-vt\right)+2\text{π}iy{\xi }_2\right)\text{d}{\xi }_1\text{d}{\xi }_2\\ =\frac{m}{n}{\iint }_{R^2}\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1\right)}^2-{\left(2\text{π}\sqrt{D_{y}t}{\xi }_2\right)}^2+2\text{π}i{\xi }_1\left(x-vt\right)+2\text{π}iy{\xi }_2\right)\text{d}{\xi }_1\text{d}{\xi }_2\\ =\frac{m}{n}{\iint }_{R^2}\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right)}^2-{\left(2\text{π}\sqrt{D_{y}t}{\xi }_2-\frac{iy}{2\sqrt{D_{y}t}}\right)}^2-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right)\text{d}{\xi }_1\text{d}{\xi }_2\end{array}$

$\begin{array}{c}=\frac{m}{n}\cdot \exp \left(-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right)\cdot {\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right)}^2\right)\text{d}{\xi }_1\\ \cdot {\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{y}t}{\xi }_2-\frac{iy}{2\sqrt{D_{y}t}}\right)}^2\right)\text{d}{\xi }_2.\end{array}$ (19)

And because $exp\left(-t^2\right)$ is an even function, it is obviously to get

${\int }_{-\infty }^{+\infty }exp\left(-t^2\right)\text{d}t=2{\int }_0^{+\infty }exp\left(-t^2\right)\text{d}t=\sqrt{\text{π}}\mathrm{,}$ (20)

therefore, take Equation (19) into account, we have

$\begin{array}{l}{\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right)}^2\right)\text{d}{\xi }_1\\ =\frac{1}{2\text{π}\sqrt{D_{x}t}}{\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right)}^2\right)\text{d}\left(2\text{π}\sqrt{D_{x}t}{\xi }_1\right)\\ =\frac{1}{2\text{π}\sqrt{D_{x}t}}{\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right)}^2\right)\text{d}\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right),\end{array}$ (21)

where $\frac{i\left(x-vt\right)}{2\text{π}\sqrt{D_{x}t}}$ is a constant, using Equaiton (20), then Equaiton (21) can be written to

${\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{x}t}{\xi }_1-\frac{i\left(x-vt\right)}{2\sqrt{D_{x}t}}\right)}^2\right)\text{d}{\xi }_1=\frac{\sqrt{\text{π}}}{2\text{π}\sqrt{D_{x}t}}=\frac{1}{2\sqrt{\text{π}}\sqrt{D_{x}t}},$ (22)

similarly, we have

${\int }_{-\infty }^{+\infty }\exp \left(-{\left(2\text{π}\sqrt{D_{y}t}{\xi }_2-\frac{iy}{2\sqrt{D_{y}t}}\right)}^2\right)\text{d}{\xi }_2=\frac{1}{2\sqrt{\text{π}}\sqrt{D_{y}t}}.$ (23)

If we plug Equation (19) and Equation (22) back into Equation (23), we obtain

$\begin{array}{c}C\left(x,y,t\right)=\frac{m}{n}\cdot \frac{1}{2\sqrt{\text{π}}\sqrt{D_{x}t}}\cdot \frac{1}{2\sqrt{\text{π}}\sqrt{D_{y}t}}\cdot \exp \left(-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right)\\ =\frac{m}{n}\cdot \frac{1}{4\text{π}t\sqrt{D_x{D}_y}}\cdot \exp \left(-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right),\end{array}$ (24)

it follows from the above equations that

$C\left(x,y,t\right)=\frac{\frac{m}{n}}{4\text{π}t\sqrt{D_x{D}_y}}\cdot \exp \left(-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right).$ (25)

As mentioned above, for Equation (19), taking two-dimensional Fourier transform and the inverse Fourier transform of both sides, its analytical solution

$C\left(x\mathrm{,}y\mathrm{,}t\right)=\frac{\frac{m}{n}}{4\text{π}t\sqrt{D_x{D}_y}}\cdot exp\left(-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right)\mathrm{.}$ (26)

is obtained.

Using MATLAB software to program, under the condition of $L_x=L_y=0:30$ , $T=1$ or $T=5$ , $\frac{m}{n}=5$ , $v=0.1$ , we give the figures of analytical solution

when the diffusion coefficient is $D_x=D_y=1$ or $D_x=D_y=2$ , respectively. The the height value represents the size of concentration in Figure 1 and Figure 2.

As shown in Figure 1, we can see that the concentration becomes smaller when the diffusion coefficient becomes larger at the same time. From the Figure 2, we can see that when the diffusion coefficient is constant, the concentration decreases as time increases. The results as what we have anticipated.

3. Numerical Methods

In this section, we study the structures and the properties of the numerical methods. As stated in the Section 1 and 2, Equation (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, so many numerical methods can be used to solve Equation (1), this is one of the most significant problem. Numerical modeling of the groundwater flow in an aquifer is adopted from the detailed study of Prickett and Lonnquist [22] , similarly, we develop finite difference equations for the advective- dispersive contaminant transport. Here to simulate the law of movement about pollutant in the medium, we present a second-order scheme to discretize the governing equation, which is based on centered Crank-Nicolson finite difference scheme [15] - [20] , moveover, the discretization of the physical domain for contaminant transport and the groundwater flow is given in figures [23] .

For the presentation of our finite difference method, we first introduce some notations which will be used later. Let region of interest be $\Omega =\left[\mathrm{0,}{L}_x\right]\times \left[\mathrm{0,}{L}_y\right]$ and the boundary of $\Omega$ be $\partial \Omega$ , we denote temporal increment by tau. For the spatial approximation, take two positive integers $M_1$ and $M_2$ , and take the step sizes $h_1=L_x/M_1$ , $h_2=L_y/M_2$ , respectively. In this way, the spatial nodes can be denoted by $\left(x_i\mathrm{,}{y}_i\right)$ , $x_i=i{h}_1\mathrm{<mi>\; </mi>}\left(0\le i\le M_1\right)$ , $y_j=j{h}_2\mathrm{<mi>\; </mi>}\left(0\le j\le M_2\right)$ . In addition, we define ${\overline{\Omega }}_h=\left\{\left(x_i,y_j\right)|0\le i\le M_1,0\le j\le M_2\right\}$ , ${\Omega }_h={\overline{\Omega }}_h\cap \Omega$ , ${\Gamma }_h={\overline{\Omega }}_h\cap \partial \Omega$ .

3.1. Derivation of the Difference Scheme

In this part, we mainly consider the difference scheme and give preliminary results for the numerical approximation of the following equations

$\frac{\partial C}{\partial t}=D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}-v\frac{\partial C}{\partial x},\left(x,y\right)\in \Omega ,t>0,$ (27)

$C\left(x,y,0\right)=0,\left(x,y\right)\in \Omega ,$ (28)

$C\left(x\mathrm{,}y\mathrm{,}t\right)=\varphi \left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}\left(x\mathrm{,}y\right)\in \partial \Omega \mathrm{.}$ (29)

where $\Omega =\left(\mathrm{0,}{L}_x\right)\times \left(\mathrm{0,}{L}_y\right)$ . We assume that $\varphi \left(x\mathrm{,}y\mathrm{,}t\right)$ is known smooth functions, and the diffusion coefficients $D_x\mathrm{,}{D}_y$ and convection coefficient v are constant. For simplicity, introduce

$\omega =\left\{\left(i,j\right)|\left(x_i,y_j\right)\in {\Omega }_h\right\},\gamma =\left\{\left(i,j\right)|\left(x_i,y_j\right)\in {\Gamma }_h\right\},\overline{\omega }=\omega \cup \gamma .$

Define ${\mathcal{U}}_h=\left\{u|u=\left\{u_{ij}|\left(i,j\right)\in \overline{\omega }\right\}\right\}.$ For any $u\mathrm{,}v\in {\mathcal{U}}_h$ , introduce the following notations of difference quotients

${\delta }_x{u}_{i+\frac{1}{2},j}=\frac{1}{h_1}\left(u_{i+1,j}-u_{ij}\right),{\delta }_x{u}_{i-\frac{1}{2},j}=\frac{1}{h_1}\left(u_{ij}-u_{i-1,j}\right),$

${\Delta }_x{u}_{ij}=\frac{1}{2h_1}\left(u_{i+1,j}-u_{i-1,j}\right),{\delta }_x^2{u}_{ij}=\frac{1}{h_1^2}\left(u_{i-1,j}-2u_{ij}+u_{i+1,j}\right),$

and similarly for the y direction $\left({\delta }_y{u}_{i\mathrm{,}j-\frac{1}{2}}\mathrm{,}{\delta }_y{u}_{i\mathrm{,}j+\frac{1}{2}}\mathrm{,}{\Delta }_y{u}_{ij}\mathrm{,}{\delta }_y^2{u}_{ij}\right)$ . It is easy to observe that

${\Delta }_x{u}_{ij}=\frac{1}{2}\left({\delta }_x{u}_{i+\frac{1}{2},j}+{\delta }_x{u}_{i-\frac{1}{2},j}\right),{\delta }_x^2{u}_{ij}=\frac{1}{h_1}\left({\delta }_x{u}_{i+\frac{1}{2},j}-{\delta }_x{u}_{i-\frac{1}{2},j}\right),$

${\Delta }_y{u}_{ij}=\frac{1}{2}\left({\delta }_y{u}_{i,j+\frac{1}{2}}+{\delta }_y{u}_{i,j-\frac{1}{2}}\right),{\delta }_y^2{u}_{ij}=\frac{1}{h_2}\left({\delta }_y{u}_{i,j+\frac{1}{2}}-{\delta }_y{u}_{i,j-\frac{1}{2}}\right).$

The aim of the present part is to improve the accuracy in the temporal direction, so for the advection term of (27), we employ the average of the time level of $n+1$ and $n$ , which guarantees that the discrete scheme for (27) is unconditionally stable and second-order accurate in space and time.

For the time approximation, we take a positive integer N, let $\tau =T/N$ , partition the interval $\left[\mathrm{0,}T\right]$ into N equal parts of width $\tau$ . Also we assume $t_n=n\tau$ , $n=0,1,\cdots ,N$ , where $\tau$ is the time grid stepsize. Let

$t_{n+\frac{1}{2}}=\left(t_n+t_{n+1}\right)/2,{\Omega }_{\tau }=\left\{t_n|0\le n\le N\right\},$

${\mathcal{U}}_{\tau }=\left\{w|w=\left(w^0\mathrm{,}{w}^1\mathrm{,}\cdots \mathrm{,}{w}^N\right)\right\}\mathrm{,}$

for any $w\in {\mathcal{U}}_{\tau }$ , introduce some notations as follows

$w^{n+\frac{1}{2}}=\frac{1}{2}\left(w^{n+1}+w^n\right),{\delta }_t{w}^{n+\frac{1}{2}}=\frac{1}{\tau }\left(w^{n+1}-w^n\right).$

We define grid functions on ${\overline{\Omega }}_h\times {\Omega }_{\tau }$

$C_{ij}^n=u\left(x_i,y_j,t_n\right),\left(i,j\right)\in \overline{\omega },0\le n\le N,$

considering (27) at the point $\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+\frac{1}{2}}\right)$ , it holds that

$\begin{array}{c}\frac{\partial C}{\partial t}\left(x_i,y_j,t_{n+\frac{1}{2}}\right)=D_x\frac{{\partial }^{2}C}{\partial x^2}\left(x_i,y_j,t_{n+\frac{1}{2}}\right)+D_y\frac{{\partial }^{2}C}{\partial y^2}\left(x_i,y_j,t_{n+\frac{1}{2}}\right)\\ -v\frac{\partial C}{\partial x}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+\frac{1}{2}}\right)\mathrm{.}\end{array}$ (30)

where $\left(i\mathrm{,}j\right)\in \omega \mathrm{,\; <mi>\; </mi>0}\le n\le N-1$ .

Thus, we obtain the discretized form of Equaiton (27).

3.2. Local Truncation Error

A numerical method used in the derivation of the finite difference equations, for groundwater contaminant transport, can be described as a flux conserving scheme which also proved that the equations don't contain spatial error with respect to the model it mentioned [24] . However truncation errors exist, due to the approximations used to describe the flow of groundwater contaminant. Truncation errors might cause significant solution inaccuracies for the numerical models, especially for the advective dispersive contaminant transport case. We know the fact that when advective terms in the advection-dispersion equation become dominant, the equation behaves more like a hyperbolic than a parabolic [25] . But in our paper, the advective terms in our model is very small and not dominant, so it can be seem as a classic two-dimensional parabolic equation.

We denote $C_{ij}^n$ as the solutions of (27), and take Taylor formula in to account, applying the Taylor expansion for (30), it generates

${\delta }_t{C}_{ij}^{n+\frac{1}{2}}=D_x{\delta }_x^2{C}_{ij}^{n+\frac{1}{2}}+D_y{\delta }_y^2{C}_{ij}^{n+\frac{1}{2}}-v{\Delta }_x{C}_{ij}^{n+\frac{1}{2}}+R_{ij}^{n+\frac{1}{2}},$ (31)

where $\left(i\mathrm{,}j\right)\in \omega \mathrm{,\; <mi>\; </mi>0}\le n\le N-1$ .

$R_{ij}^{n+\frac{1}{2}}$ is truncation errors, and it is easy to obtain

$\left|R_{ij}^{n+\frac{1}{2}}\right|=O\left({\tau }^2+h_1^2+h_2^2\right),\left(i,j\right)\in \omega ,\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\le n\le N-1.$ (32)

which show that this scheme is second order accurate in time and space.

3.3. Crank-Nicolson Scheme

From the initial and boundary condition Equations (27)-(29), we have

$C_{ij}^0=0,\left(i,j\right)\in \omega ,$ (33)

$C_{ij}^n=\varphi \left(x_i,y_j,t_n\right),\left(i,j\right)\in \gamma ,0\le n\le N,$ (34)

Ignoring the higher order terms $R_{ij}^{n+\frac{1}{2}}$ in (31), and denote $C_{ij}^n$ as numerical

solutions of (27). Replacing $C_{ij}^n$ with its approximation $c_{ij}^n$ , we construct, for Equations (27)-(29), a Crank-Nicolson finite difference scheme

${\delta }_t{c}_{ij}^{n+\frac{1}{2}}=D_x{\delta }_x^2{c}_{ij}^{n+\frac{1}{2}}+D_y{\delta }_y^2{c}_{ij}^{n+\frac{1}{2}}-v{\Delta }_x{c}_{ij}^{n+\frac{1}{2}},\left(i,j\right)\in \omega ,0\le n\le N-1,$ (35)

$c_{ij}^0=0,\left(i,j\right)\in \omega ,$ (36)

$c_{ij}^n=\varphi \left(x_i,y_j,t_n\right),\left(i,j\right)\in \gamma ,0\le n\le N,$ (37)

we take the preceding sign into Equation (31), we have

$\begin{array}{c}\frac{c_{ij}^{n+1}-c_{ij}^n}{\tau }=\frac{1}{2}{D}_x\left[\frac{c_{i+\mathrm{1,}j}^{n+1}-2c_{ij}^{n+1}+c_{i-\mathrm{1,}j}^{n+1}}{h_1^2}+\frac{c_{i+\mathrm{1,}j}^n-2c_{ij}^n+c_{i-\mathrm{1,}j}^n}{h_1^2}\right]\\ +\frac{1}{2}{D}_y\left[\frac{c_{i\mathrm{,}j+1}^{n+1}-2c_{ij}^{n+1}+c_{i\mathrm{,}j-1}^{n+1}}{h_2^2}+\frac{c_{i\mathrm{,}j+1}^n-2c_{ij}^n+c_{i\mathrm{,}j-1}^n}{h_2^2}\right]\\ -\frac{1}{4}v\left[\frac{c_{i+\mathrm{1,}j}^{n+1}-c_{i-\mathrm{1,}j}^{n+1}}{h_1}+\frac{c_{i+\mathrm{1,}j}^n-c_{i-\mathrm{1,}j}^n}{h_1}\right]\end{array}$ (38)

Based on the formulation (38) we can get the concentration distribution of pollutants and obtain the important data of velocity, pressure and so on. In order to solve the above Crank-Nicolson difference equation, we introduce

$S_x=\frac{\tau D_x}{2h_1^2}$ , $S_y=\frac{\tau D_y}{2h_2^2}$ , $S_v=\frac{\tau v}{4h_1}$ , Equaiton (38) can be written as

$\begin{array}{c}{c}_{ij}^{n+1}-c_{ij}^n=S_x\left[\left(c_{i+1,j}^{n+1}-2c_{ij}^{n+1}+c_{i-1,j}^{n+1}\right)+\left(c_{i+1,j}^n-2c_{ij}^n+c_{i-1,j}^n\right)\right]\\ +S_y\left[\left(c_{i,j+1}^{n+1}-2c_{ij}^{n+1}+c_{i,j-1}^{n+1}\right)+\left(c_{i,j+1}^n-2c_{ij}^n+c_{i,j-1}^n\right)\right]\\ -S_v\left[\left(c_{i+1,j}^{n+1}-c_{i-1,j}^{n+1}\right)+\left(c_{i+1,j}^n-c_{i-1,j}^n\right)\right],\end{array}$ (39)

Extracting the coefficient of $c_{ij}^{n+1}$ , $c_{i+\mathrm{1,}j}^{n+1}$ , $c_{i-\mathrm{1,}j}^{n+1}$ , $c_{i\mathrm{,}j+1}^{n+1}$ , $c_{i\mathrm{,}j-1}^{n+1}$ , $c_{ij}^n$ , $c_{i+\mathrm{1,}j}^n$ , $c_{i-\mathrm{1,}j}^n$ , $c_{i\mathrm{,}j+1}^n$ , $c_{i\mathrm{,}j-1}^n$ from both sides of the above Equation (39), collection the like terms, we get

$\begin{array}{l}\left[1+2\left(S_x+S_y\right)\right]c_{ij}^{n+1}+\left(S_v-S_x\right)c_{i+1,j}^{n+1}-\left(S_x+S_v\right)c_{i-1,j}^{n+1}-S_y\left(c_{i,j+1}^{n+1}+c_{i,j-1}^{n+1}\right)\\ =\left[1-2\left(S_x+S_y\right)\right]c_{ij}^n+\left(S_x-S_v\right)c_{i+1,j}^n+\left(S_x+S_v\right)c_{i-1,j}^n+S_y\left(c_{i,j+1}^n+c_{i,j-1}^n\right).\end{array}$ (40)

As shown in Figure 3, Equaiton (40) involves unknown concentrations of five grid nodes center around $\left(i\mathrm{,}j\right)$ at the moment $t_{n+1}$ and $t_n$ , we may easily give Equaiton (40) according to all of the the interior grids. Using the known initial conditions and boundary conditions, we get a five diagonal equations, and for each advancing time step, a set of such equations must be solved.

We assumption that the pollutants follow the flow faithfully, we can get $c_{ij}^0$ from the exact solution, at the start of each time step, predictions for the principal unknowns, which are based on the response of the model at previous time steps, are performed. Combined with boundary condition Equation (4), $c_{ij}^n$ can be obtained by layer by layer. However, when the n layer is calculated by the $n+1$ layer, it is necessary to solve the system of linear algebraic equations (because the coefficient matrix is strictly diagonally dominant, the equations can be uniquely solved), the distribution of the network as shown in Figure 4.

With a simple calculation, we rewrite the Crank-Nicolson scheme (40) into the following matrix form

$A{c}^{n+1}=B{c}^n+c^{*}.$ (41)

where

$c^{n+1}={\left[c_{11}^{n+1},c_{21}^{n+1},\cdots ,c_{M1}^{n+1},c_{12}^{n+1},c_{22}^{n+1},\cdots ,c_{M2}^{n+1},\cdots ,c_{1M}^{n+1},c_{2M}^{n+1},\cdots ,c_{MM}^{n+1}\right]}_{M^2\times 1}^{\text{T}}$

$c^n={\left[c_{11}^n,c_{21}^n,\cdots ,c_{M1}^n,c_{12}^n,c_{22}^n,\cdots ,c_{M2}^n,\cdots ,c_{1M}^n,c_{2M}^n,\cdots ,c_{MM}^n\right]}_{M^2\times 1}^{\text{T}}$

$\begin{array}{c}{c}^{*}=[\left(S_x+S_v\right)c_{01}^{n+1}+S_y{c}_{10}^{n+1},\left(S_x+S_v\right)c_{02}^{n+1},\cdots ,\left(S_x+S_v\right)c_{0M}^{n+1}+S_y{c}_{1,M+1}^{n+1},\\ S_y{c}_{20}^{n+1},0,\cdots ,S_y{c}_{2,M+1}^{n+1},\cdots ,\left(S_x-S_v\right)c_{M+1,1}^{n+1}+S_y{c}_{M0}^{n+1},\left(S_x-S_v\right)c_{M+1,2}^{n+1},\\ {\cdots ,\left(S_x-S_v\right)c_{M+1,M}^{n+1}+S_y{c}_{M,M+1}^{n+1}]}_{M^2\times 1}^{\text{T}}\\ +[\left(S_x+S_v\right)c_{01}^n+S_y{c}_{10}^n,\left(S_x+S_v\right)c_{02}^n,\cdots ,\left(S_x+S_v\right)c_{0M}^n+S_y{c}_{1,M+1}^n,\\ S_y{c}_{20}^n,0,\cdots ,S_y{c}_{2,M+1}^n,\cdots ,\left(S_x-S_v\right)c_{M+1,1}^n+S_y{c}_{M0}^n,\left(S_x-S_v\right)c_{M+1,2}^n,\\ {\cdots ,\left(S_x-S_v\right)c_{M+1,M}^n+S_y{c}_{M,M+1}^n]}_{M^2\times 1}^{\text{T}}\end{array}$

the M²-by-M² matrix $A$ is expressed as

$A=\left(\begin{array}{ccccc}{P}_1& Q_1& 0& \cdots & 0\\ Q_1& P_1& Q_1& \cdots & 0\\ 0& Q_1& Q_1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & P_1\end{array}\right)$

the M²-by- M² matrix $B$ is given by

$B=\left(\begin{array}{ccccc}{P}_2& Q_2& 0& \cdots & 0\\ Q_2& P_2& Q_2& \cdots & 0\\ 0& Q_2& Q_2& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & P_2\end{array}\right)$

For simplicity of presentation, the entries of matrix $A$ and matrix $B$ are given by matrix $P_1$ , matrix $Q_1$ , matrix $P_2$ and matrix $Q_2$ .

the M-by-M matrix $P_1$ is

$P_1=\left(\begin{array}{ccccc}1+2\left(S_x+S_y\right)& S_v-S_x& 0& \cdots & 0\\ -\left(S_x+S_v\right)& 1+2\left(S_x+S_y\right)& S_v-S_x& \cdots & 0\\ 0& -\left(S_x+S_v\right)& 1+2\left(S_x+S_y\right)& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1+2\left(S_x+S_y\right)\end{array}\right)$

the M-by-M matrix $Q_1$ is

$Q_1=\left(\begin{array}{ccccc}-S_y& & & & \\ & -S_y& & & \\ & & \ddots & & \\ & & & -S_y& \\ & & & & -S_y\end{array}\right)$

the M-by-M matrix $P_2$ is

$P_2=\left(\begin{array}{ccccc}1-2\left(S_x+S_y\right)& S_x-S_v& 0& \cdots & 0\\ S_x+S_v& 1-2\left(S_x+S_y\right)& S_x-S_v& \cdots & 0\\ 0& S_x+S_v& 1-2\left(S_x+S_y\right)& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1-2\left(S_x+S_y\right)\end{array}\right)$

the M-by-M matrix $Q_2$ is

$Q_2=\left(\begin{array}{ccccc}{S}_y& & & & \\ & S_y& & & \\ & & \ddots & & \\ & & & S_y& \\ & & & & S_y\end{array}\right).$

4. Numerical Simulation

In this part, we will give an experiments to test our Crank-Nicolson difference scheme (40). 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 x and y directions.

Before conducting numerical experiments, we first make three hypotheses. Firstly, the main direction of the dispersion of the material in the groundwater is in the same direction as the coordinate axis. Secondly, the main seepage direction of groundwater is also consistent with the coordinate axis. Thirdly, the diffusion coefficient and permeability coefficient are constant in the axis direction.

So, we can choose the following parameters in the numerical experiment

$L_x=L_y=20,T=10,M_1=M_2=20\text{or}40,$

$N=100,D_x=D_y=1,m=n=1,v=\mathrm{0.1.}$

Example: The equations to be solved are

$\frac{\partial C}{\partial t}=D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}-v\frac{\partial C}{\partial x},\mathrm{<mi>\; </mi>}\left(x,y\right)\in \left(0,L_x\right)\times \left(0,L_y\right),t\in \left(0,T\right],$ (42)

$C\left(x,y,0\right)=0,\left(x,y\right)\in \Omega ,$ (43)

$C\left(x\mathrm{,}y\mathrm{,}t\right)=\varphi \left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}\left(x\mathrm{,}y\right)\in \Gamma \mathrm{,}$ (44)

where $\Omega =\left(\mathrm{0,}{L}_x\right)\times \left(\mathrm{0,}{L}_y\right)$ , $\Gamma =\partial \Omega$ . Let the exact solution is

$C\left(x,y,t\right)=\frac{\frac{m}{n}}{4\text{π}t\sqrt{D_x{D}_y}}\cdot \exp \left(-\frac{{\left(x-vt\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right).$ (45)

Using MATLAB software, we can get the function figures of the numerical solution and analytical solution, as indicated in Figure 5 and Figure 6. It proves that the degree of numerical solutions is approximating the exact solutions in different grid points.

Numerical experiments are performed to show the efficiency of the Crank- Nicolson schemes, which is reliably for solving two-dimensional parabolic

equation with convection term, it is apparent from the Figure 5 and Figure 6 that the exact solution and numerical solution are better fitted with mesh generation encrypting, it turns out that the difference scheme we used can be a good approximate to analytical solution.

5. Visualization of Simulation Results Based on GIS (ArcGIS Figures)

The result of water pollution simulation is a large amount of data stored in the form of documents. It is the basic purpose of water pollution diffusion simulation to analyze and study these data, so as to master the law of water quality change, using the traditional method to deal with the simulation results, although the method is feasible, but it is extremely time-consuming. We are going to use programs to carry out numerical simulation of groundwater pollution problems and realize the time and concentration prediction of pollutants transport.

In this section, a new system concluding the design of the visual area and the dispersion of simulation interval has been developed within the framework of ArcGIS based on two dimensional groundwater pollution simulations, it also provides a perfect interface which allows construction of user defined criteria, such as running computations and visualization of the results. At the same time, we can make a conclude that the use of GIS feature layer demonstrate that the GIS-based two dimensional groundwater pollution simulations system can provide the user better decision making aid.

5.1. The Design of the Visual Area and the Dispersion of the Simulation Interval

Geographic information system (GIS) has the ability to express strong geographic data, is a powerful tool for water pollution simulation results visualization. In the simulation of groundwater pollution, there are various barrier and recharge boundaries, a river, and a contaminant source. The aquifer domain is uniformly discretized with a grid interval of certain a feet we defined according our needed. The initial heads in both the aquifer and groundwater contamination are set to zero, so we view port to rectangle region, x axis to plane, y axis to formation height, point $\left(\mathrm{0,0}\right)$ to source of pollution, we can get the simulation of the discrete interval, as shown in Figure 7.

5.2. ArcGIS Visualization of Density Images at Different Time

The calculation results of MATLAB along the river equidistant distribution grid, each grid point according to the time sequence of stored information of the coordinates of the concentration of pollutants. According to our information and grid coordinate information, create the space layer, and into the ArcMap. To the data at every moment in the MATLAB, constantly modify the value of t, we import the concentration matrix of each time t into the EXCEL file every 20 time points. From the time point corresponding to the EXCEL file derived concentration coordinates table, the concentrations of each column were $x\mathrm{,}y$ . Import the final five EXCEL tables into ArcMap, the concentration image at each moment can be obtained.

Because the distribution of pollutant concentration is a continuous field of three-dimensional space, and the calculation results for the spatial discretization of the grid point distribution, in order to achieve the pollutant concentration diffusion field visualization continuous two-dimensional representation, with a certain height h as the observation height, the leakage of space in parallel to the O-xy plane in the direction of cutting, with a height of $h+\delta$ and $h-\delta$ ( $\delta$ is error tolerance) for planes in space, between the two cut point set as h height of the data source, extract t (diffusion time) in highly h data source point deduction collection time. Taking Galerkin spatial interpolation, the reduced leakage of discrete space the expression for a continuous space, we can get h height, pollutant concentration distribution of the t moment.

Based on the numerical experiments in Section 4, we import the pollution concentration spatiotemporal distribution data (1 - 5 hours) into the ArcMap platform to realize two-dimensional visualization of pollution diffusion. The figures of ArcMap display intuitively the concentration distribution and dynamic changes of pollutants at different times and locations. Figures 8-12 indicate diffusion changes after pollution flows into the river 1 h, 2 h, 3 h, 4 h and 5 h, respectively. In Figures, the pollutant concentration is divided into grade 1 - grade 5 according to the degree of risk from high to low, pollutant concentration range were >900, 601 - 900, 301 - 600, 101 - 300 and <100 mg/m³, we render the different areas in red, orange, yellow, green and blue, respectively.

We set the level of concentration range in terms of the needed, and divide the contours into different regions, such as lethal areas, dangerous areas, warning area and safety area etc. We define different color values in different regions, which can represent the spatial distribution of pollutant concentration at each moment, thus we achieve two-dimensional visualization of the simulation results of river water pollution.

The present study shows risk and uncertainty analysis based on two- dimensional numerical simulation results, GIS can significantly improve the accuracy of groundwater pollution hazard assessment. This approach efficiently assists in evaluation and ranking of groundwater pollution control management strategies, and future design of groundwater pollution proofing works. The

accuracy of predictions can probably be improved by taking into account other detailed information provided by two-dimensional groundwater pollution modeling, such as groundwater pollution velocity and duration. Further research is needed to develop such improved damage relationships.

6. Conclusion and Suggestions

Numerical simulation of water quality is a powerful tool to study groundwater pollution. In this paper, we give analytical solution of the Equations (1)-(5) by two-dimensional Fourier transform and the inverse Fourier transform firstly. Next, we design an efficient second-order finite difference scheme for solving a class of groundwater problem and simulate the migration of pollutants and get the concentration distribution of pollutants, theoretically analysis shows that the proposed scheme is second-order accurate in time and space. Then, exact solution comparisons with numerical solution are also discussed, it observes that the computed results show excellent agreement with the analytical solution, and the MATLAB software is used to simulate the numerical simulation of groundwater flow problem, which can be used to predict the time and concentration of contaminant transport. Finally, we present the development process and technical methods of visual expression of water pollution simulation based on GIS technology, we visualize the simulation results (ArcGIS figures) to study the temporal and spatial variation of water pollution problems. The simulation system is stable, fast, no jumping, and the results are intuitive and reasonable. It provides technical support and guarantee for the comprehensive management and management of water environment.

Acknowledgements

This work is supported in part by the National Natural Science Foundation of China under Grant nos. 11501335, 11371229 and by a Project of Shandong Province Higher Educational Science and Technology Program under Grant nos. J14LI03.

Cite this paper

Li, L.Y. and Yin, Z. (2017) Numerical Simulation of Groundwater Pollution Problems Based on Convection Diffusion Equation. American Journal of Computational Mathematics, 7, 350-370. https://doi.org/10.4236/ajcm.2017.73025

References