Numerical Simulation of a Two-Dimensional Vertically Averaged Groundwater Quality Assessment in Homogeneous Aquifer Using Explicit Finite Difference Techniques

Leachate contamination in a landfill causes of pollution that flowing down to the groundwater. There are many methods to measure the groundwater quality. Mathematical models are often used to describe the groundwater flow. In this research, the affection of landfill construction to groundwater-quality around rural area are focused. Three mathematical models are combined. The first model is a two-dimensional groundwater flow model. It provides the hydraulic head of the groundwater. The second model is the velocity potential model. It provides the groundwater flow velocity. The third model is a two-dimensional vertically averaged groundwater pollution dispersion model. The groundwater pollutant concentration is provided. The forward time centered technique with the centered in space, the forward in space and the backward in space with all boundaries are used to obtain approximate hydraulic head, the flow velocity in xand ydirections, respectively. The approximated groundwater flow velocity is used to input into a two-dimensional vertically averaged groundwater pollution dispersion model. The forward time centered space technique with the centered in space, the forward in space and the backward in space with all boundaries are used to obtain approximate the groundwater pollutant concentration. The proposed explicit forward time centered spaced finite difference techniques to the groundwater flow model the velocity potential model and the groundwater pollution dispersion model give good agreement approximated solutions.


Introduction
Industrial development often causes pollution problems. Contamination of surface water causes pollution of chemicals in water into groundwater, such as the removal of toxins in landfills, solid waste into the ground. For example, groundwater gets contaminated by leaching of nitrate generated from fertilizer used on agricultural lands and waste dumps in rural and urban areas. The nitrate decrease ability of blood to carry oxygen, resulting in oxygen deficiency in different body parts linked to blue baby syndrome [1]. Salinity contamination in groundwater in a rice field near marine shrimp aquaculture farm leads to poor crop production [2]- [3]. As water percolated down through soil and rocks, bacteria, fungi, and other such biological pollutants were naturally filtered out or diluted. But in recent years, groundwater natural defence systems have been vastly overextended. The sheer volume of pollutants sent underground has escalated. Especially in the last century, when the global population increased. There is a growing demand of food and water. At the same time, the industry is quickly advancing, in this way contaminating the environment. Rivers and streams are also damaged by consumption and consumption. Whenever the water is on the surface, it is damaged, and it affects the quality of the groundwater. Since groundwater is the water found under of the land that is accumulated on the ground and sand, it is filtered from the ground, leaving the groundwater clean. Groundwater exists mainly in an aquifer. The aquifer is a body of saturated rock through which water can move. Since underground water is moving in the underground, we cannot tell the direction of groundwater movement where the groundwater is. If we do not dig and drill underground. Henry Darcy was interested in the description of groundwater flow. Until nowadays, Darcy's law is utilized describe in the mathematical formulations in scientific disciplines concerned with measuring water flow. Therefore, we will be using groundwater modeling of Darcy's law. The groundwater modeling [5] describes phenomena or to predict future behaviour. Groundwater systems can be modeled using partial differential equations. The equation can solve by finite difference method and finite element method. There have been studies of the groundwater modeling simulation algorithm with finite difference method [6] by using spreadsheets [7]- [8], MODFLOW [9].
In this research, a combination of three related groundwater quality measurement models is proposed. The first model is a two-dimensional vertically averaged groundwater flow model. It is providing the hydraulic head of groundwater. The hydraulic head is transformed to be the groundwater flow velocity in the second model that is the velocity potential model. The calculated flow velocity input into the third model. The last model is a two-dimensional vertically averaged groundwater pollution dispersion model. The dispersion model provides the pollutant concentration that is contaminated in the groundwater. Explicit finite difference method, forward time central space with a centered space technique, a forward space technique and a backward space technique on the boundary's solution, are employed to approximate their solutions of all models.

Governing equation
The groundwater flow through soil is governed by the Darcy's law that be described by partial differential equation.

A two-dimensional vertically averaged groundwater flow model
Groundwater flow can be described by a transient hydraulic head flow model. The aquifer is assumed by a porous medium. If we consider the homogeneous aquifer, the soil transmitting will be uniform. The groundwater flow model is used to represent the behavior of hydraulic head. The transmission property of geologic is hydraulic conductivity including the considered aquifer and their space. A primary tool for illustration of groundwater flow is Darcy's law. Assuming that the hydraulic head in the considered area is vertically averaged. The domain considered is assumed by Ω = {(x, y) : 0 ≤ x ≤ I, 0 ≤ y ≤ J}. The governing equation of Darcy's flow is [6], where h(x, y, t)(m) is the hydraulic head, K x (m/day) is the hydraulic conductivity in x-direction, K y (m/day) is the hydraulic conductivity in y-direction, W (day −1 ) is the source or sink function, S is the specific storage coefficient for all (x, y, t) ∈ Ω × [0, T ] as shown in Fig. 1. If we consider homogeneous aquifer system, the hydraulic conductivity will be constant K x = K y = K. Then ( 1 ) becomes

Initial condition
If the potential hydraulic head in the area is static, the initial condition is assumed by where h 0 (m) is a given averaged potential hydraulic head in the considered area.

Boundary condition
The rates of change of the hydraulic head along the domain boundaries are given by where B hN , B hS , B hW and B hE are the rates of change of hydraulic head on the north, south, west and east domain boundaries as shown in Fig. 2.  Velocity potential is a powerful tool in analyzing irrotational flows. The groundwater flows model provides the hydraulic head. The velocity potential model is also provided by the hydraulic head. The groundwater flow velocity in xand ydirections can be obtained by [6], 154 Numerical Simulation of a Two-Dimensional Vertically Averaged Groundwater Quality Assessment in Homogeneous Aquifer Using Explicit Finite Difference Techniques where u, v (m/day) are groundwater flow velocity in xand ydirections, respectively.

A two-dimensional vertically averaged groundwater pollution dispersion model
An advection-diffusion model provides a continuous description of groundwater pollutant transport in the groundwater. A two-dimensional vertically averaged groundwater pollution dispersion model is [10] ∂C ∂t where C(x, y, t) (kg/m 3 ) is groundwater pollutant concentration, u (m/day) is groundwater flow velocity in x-direction, v (m/day) is groundwater flow velocity in y-direction, D x (m 2 /day) is the diffusion coefficient of groundwater pollution through specified soil in x-axis, D y (m 2 /day) is the diffusion coefficient of groundwater pollution through specified soil in yaxis and W s (day −1 ) is groundwater pollutant sources or sink function by contaminators.

Initial condition
If the potential groundwater pollutant concentration in the consider area is described by where C 0 (kg/m 3 ) is a averaged potential groundwater pollutant concentration in the considered area.

Boundary condition
The rates of change of the pollutant concentration along the domain boundaries are given by where B CN , B CS , B CW and B CE are the rates of change of the pollutant concentration on the north, south, west and east domain boundary as shown in Fig. 3.

Numerical techniques
In this research, approximated solutions are obtained by finite difference techniques. A numerical solutions are based on the concept that the partial differential equation can be replaced We now discretize ( 2 ) by dividing the interval [0, I] into M subintervals such that M ∆x = I, the interval [0, J] into N subintervals such that N ∆y = J, and the interval [0, T ] into P subintervals such that P ∆t = T . We can then approximate h(x i , y j , t n ) by h n i,j , value of the difference approximation of h(x, y, t) at point x i = i∆x, y j = j∆y and t n = n∆t, when 0 ≤ i ≤ M , 0 ≤ j ≤ N and 0 ≤ n ≤ P , in which M, N and P are positive integers.

Explicit finite difference technique for a twodimensional vertically averaged groundwater flow model
The approximated hydraulic head, h(x, y, t) ≈ h(i ∆x , j ∆y , n ∆t ) = h n i,j .

Forward time centered space
We introduced the forward time centered space technique(FTCS) to a two-dimensional vertically averaged groundwater flow model ( 2 ), we obtained that For i = 0, j = 0 and n = 0, the approximated fictitious points on the boundaries, are obtained by For i = 0 ,1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = 0 ,j = J and n = 0, the approximated fictitious points on the boundaries, are obtained by h n 0,J+1 = h n 0,J−1 + 2B hN ∆y.
Substituting ( 24 ) and ( 25 ) into ( 18 ), it is obtained that For 1 < i < I ,j = 0 and n = 0, the approximated fictitious point on the boundary, is obtained by For 1 < i < I ,j = J and n = 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 29 ) into ( 18 ), it is obtained that For i = I ,j = 0 and n = 0, the approximated fictitious points on the boundaries, For i = I ,1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = I , j = J and n = 0, the approximated fictitious points on the boundaries,, are obtained by Substituting ( 36 ) and ( 37 ) into ( 18 ), it is obtained that

Forward time centered space technique with a forward space technique to approximate the boundaries solution
For i = 0 ,j = 0 and n = 0, the approximated fictitious points on the boundaries, are obtained by 156 Numerical Simulation of a Two-Dimensional Vertically Averaged Groundwater Quality Assessment in Homogeneous Aquifer Using Explicit Finite Difference Techniques For i = 0 ,1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = 0, j = J and t > 0, the approximated fictitious points on the boundaries, are obtained by For 1 < i < I, j = 0 and n = 0, the approximated fictitious point on the boundary, is obtained by For 1 < i < I, j = J and n = 0, the approximated fictitious point on the boundary is obtained by For i = I, j = 0 and n = 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 51 ) and ( 52 ) into ( 18 ), it is obtained that For i = I, 1 < j < J and t > 0, the approximated fictitious point on the boundary, is obtained by For i = I, j = J and n = 0, the approximated fictitious points on the boundaries, are obtained by

Forward time centered space technique with a backward space technique to approximate the boundaries solution
For i = 0, j = 0 and n = 0, the approximated fictitious points on the boundaries, are obtained by For i = 0, 1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = 0, j = J and t > 0, the approximated fictitious points on the boundaries are obtained by For 1 < i < I, j = 0 and n = 0, the approximated fictitious point on the boundary, is obtained by For 1 < i < I, j = J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = I, j = 0 and t > 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 71 ) and ( 72 ) into ( 18 ), it is obtained that For i = I, 1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = I, j = J and n = 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 76 ) and ( 77 ) into ( 18 ), it is obtained that

Finite difference technique for groundwater flow velocity approximation
We introduced the centered space method to the velocity potential model in two-dimension velocity field ( 7 )-( 8 ), we have where h is the approximated hydraulic head, 3.3 Explicit schemes for a two-dimensional vertically averaged groundwater pollution dispersion model

Forward time central space
Taking the forward time centered space technique in a twodimensional vertically averaged groundwater pollution dispersion model ( 9 ), we get the following discretization, Rearranging ( 86 ) for 0 ≤ i ≤ M , 0 ≤ j ≤ N and 0 ≤ n ≤ P , we have 158 For i = 0, j = 0 and n = 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 88 )-( 89 ) into ( 87 ), it is obtained that For i = 0 ,1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = 0, j = J and n = 0, the approximated fictitious points on the boundary, are obtained by Substituting ( 93 )and ( 94 ) into ( 87 ), it is obtained that For 1 < i < I, j = 0 and t > 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 96 ) into ( 87 ), it is obtained that For 1 < i < I, j = J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = I, j = 0 and n = 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 100 ) and ( 101 ) into ( 87 ), it is obtained that For i = I, 1 < j < J and t > 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 103 ) into ( 87 ), it is obtained that For i = I, j = J and n = 0, the approximated fictitious points on the boundaries, are obtained by C n I,J+1 = C n I,J−1 + 2B CN ∆y.

Forward time centered space technique with a forward space technique to approximate the boundaries solution
For i = 0, j = 0 at t > 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 108 )-( 109 ) into ( 87 ), it is obtained that For i = 0, 1 < j < J and t > 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 111 ) into ( 87 ), it is obtained that For i = 0, j = J and t > 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 113 ) and ( 114 ) into ( 87 ), it is obtained that For 1 < i < I, j = 0 and t > 0, the approximated fictitious point on the boundary, is obtained by For 1 < i < I, j = J and n = 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 118 ) into ( 87 ), it is obtained that For i = I, j = 0 and t > 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 120 ) and ( 121 ) into ( 87 ), it is obtained that For i = I, 1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 123 ) into ( 87 ), it is obtained that For i = I, j = J and n = 0, the approximated fictitious points on the boundaries, are obtained by C n I+1,J = C n I,J + B CE ∆x, C n I,J+1 = C n I,J + B CN ∆y.

Forward time centered space technique with a backward space technique to approximate the boundaries solution
For i = 0, j = 0 and t > 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 128 ) and ( 129 ) into ( 87 ), it is obtained that 160 Numerical Simulation of a Two-Dimensional Vertically Averaged Groundwater Quality Assessment in Homogeneous Aquifer Using Explicit Finite Difference Techniques For i = 0, 1 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by For i = 0, j = J and t > 0, the approximated fictitious points on the boundaries, are obtained by Substituting ( 133 ) and ( 134 ) into ( 87 ), it is obtained that For 1 < i < I, j = 0 and n = 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 136 ) into ( 87 ), it is obtained that For 1 < i < I, j = J and n = 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 138 ) into ( 87 ), it is obtained that For i = I, j = 0 and n = 0, the approximated fictitious points on the boundaries,are obtained by Substituting ( 140 ) and ( 141 ) into ( 87 ), it is obtained that For i = I, 0 < j < J and n = 0, the approximated fictitious point on the boundary, is obtained by Substituting ( 143 ) into ( 87 ), it is obtained that For i = I, j = J and n = 0, the approximated fictitious points on the boundary, are obtained by Substituting ( 145 ) and ( 146 ) into ( 87 ), it is obtained that

Numerical simulations
The demonstration analyzed the behavior of groundwater flow, hydraulic head, vector velocity and pollution dispersion of groundwater. The experimented area has homogeneous aquifer.

Simulation 1: Groundwater pollution measurement in an area close to 3 land fill sites
The considered area is surrounded by a landfill site as shown in Fig. 4. The area has dimension, 10 km. × 10 km. The physical parameters are the rate of change of hydraulic head along the north boundary, the west boundary, the south boundary and the east boundary is -5 m., -5 m., -5 m. and 5 m. respectively. The potential hydraulic head is 10 m. The hydraulic conductivity is 15 m./day and the storage capacity is 100. The simulation parameters are shown in Table 1.
The grid spacing is taked by ∆x = ∆y = 50 m and ∆t = 1 day. The stationary simulation time is 10 years. In this simulation, the forward time centered space with a centered space  We get the approximated hydraulic head by using 3 different techniques that are shown in Tables 3 and Fig. 5. Next, the groundwater flow velocity is also simulated. The groundwater velocity is approximated by using ( 79 )-( 80 ). The approximated flow velocity is shown in Tables 4, 5 and Fig. 6. Finally, the approximated groundwater flow velocity will be input into the groundwater pollution dispersion model as field data. The rate of change of groundwater pollutant matter along the north boundary, the west boundary, the south boundary and the east boundary is -0.05 m., -0.05 m., -0.05 m. and 0.05 m. respectively. The potential groundwater pollutant is 0.01 m., The diffusion coefficient of groundwater pollution in x-, y-axis is 1(m. 2 /day), 1(m. 2 /day) respectively. We can summarize their groundwater pollutant matter parameters as shown in Table 2 by using the forward time centered space with a centered space technique on the boundaries solution ( 88 )-( 107 ), a forward space technique on the boundaries solution ( 108 )-( 127 ) and a backward space technique on the boundaries solution ( 128 )-( 147 ). We get the approximated groundwater pollutant concentration as shown in Table 6 and Fig. 7.

Simulation 2: Groundwater pollution measurement in an area close to 3 land fill sites and pumping term
We consider an area in simulation 1 that have system of water pumping, sink is −50 m/day and 2 sources are 25 m/day as shown in Fig. 8. We can summarize their groundwater pollutant matter parameters as show in Table 7. We get the approximated hydraulic head by using centered space techniques are   -0.14 -0.14 -0.14 -0.14 -0.14 0.00 10 -0. -0.14 -0.14 -0.14 -0.14 -0.14 0.00 7 -0.      shown in Table 9 and Fig. 9. Next, The groundwater velocity as shown in Table 10, 11 and Fig. 10. Finally, the approximated groundwater flow velocity will be input into the groundwater pollution dispersion model as field data. The system of water pumping, We can summarize their groundwater pollutant pollutant matter parameters as show in Table 8. The approximated groundwater pollutant as shown in Table 12 and Fig. 11.

Discussion
In simulation 1, the computed hydraulic head moves from the higher hydraulic head to the lower hydraulic head when the    time has been passed for 10 years. The velocity vector gives the fields groundwater flow directions in the considered area. The groundwater pollutant along the landfill has been flowing into the considered area. The groundwater pollutant in the area will be growing up as well. The pollutant concentration in around landfill vicinity will be higher than another area. This means that the groundwater pollutant is introduced by groundwater flow and its hydraulic head. The pollutant will be flowing out through the east boundary.
In simulation 2, the system of water pumping and injection are considered. There are 1 water pumping well and 2 water injection well. The water injection well can be maintained the hydraulic head level entire the consider area. The water pumping well will be reduce the hydraulic head around them along 10 years. This means that the water injection well system can be used to reserve groundwater levels as nature groundwater resources. The groundwater pollutant around the injection wells will be increased due to the groundwater flow velocity. On the other hand, the groundwater pollutant around the pumping well will be reduced due to the removal mechanism. These means that the groundwater quality can be induced by the pumping or injecting system.

Conclusion
A two-dimensional groundwater flow model that gives the hydraulic head level is introduced. The techniques of the initial condition and boundary conditions of the groundwater flow model are proposed. The groundwater flow velocity model is also introduced. It gives groundwater flow directions. The forward time centered technique with the centered in space, the forward in space and the backward in space with all boundaries are used to obtain the approximated hydraulic head, the flow velocity in x-and y-directions, respectively.
The approximated groundwater flow velocity are used to input into a two-dimensional vertically averaged groundwater pollution dispersion model. The groundwater pollution due to a landfill around a considered area is focused. The techniques of the initial condition and boundary condition setting is proposed. The forward time centered space technique with the centered in space, the forward in space and the backward in space with all boundaries are used to approximate the groundwater pollutant concentration. The proposed explicit forward time centered spaced finite difference techniques to the groundwater flow model the velocity potential model and the groundwater pollution dispersion model give good agreement approximated solutions. They are economical finite difference techniques due to our proposed techniques are explicit methods. It is not produced any large system of linear equations. In both simulations, the hydraulic head is depends on the hydraulic head level on their boundaries. The groundwater flow velocity is also induced by the hydraulic head. The water pumping and injecting system a can be change the hydraulic head and its groundwater flow direction. The groundwater pollutant concentration can be induced by the water pumping and injecting system. Furthermore, If a landfill project is planned to construct in an area the environmental impact assessment should be considered to the simulated groundwater quality assessment in the future as well.