### 1. Introduction

### 2. Materials and Methods

### 2.1. Data Collection and Analysis

**Step 1**: Determining the arsenic concentration in the nearby areas of the Kieu Ky landfill by the empirical research based on the statistical methods.**Step 2**: The result in**Step 1**is used as the initial value to implement the proposed mathematical algorithm to study the arsenic propagation along the leachate in the ground.

### 2.2. Mathematical Modelling of Pollutant Transport in the Leachate Immigration from a Landfill

*C*(

*x, y, z, t*) as a scalar function of concentration (mg/kg) of a pollutant in the leachate. The function

*C*(

*x, y, z, t*) is represented in a Cartesian coordinate system

*Oxyz*and changes over time

*t*. As usual, the origin of the coordinate system

*Oxyz*places on the ground surface of soil media, the axis

*Oz*points vertically downward the depth of the soil layers. The concentration function

*C*(

*x, y, z, t*) depends on four variables

*x, y, z*and

*t*.

*x, y, z*) in the soil space, the change of the pollutant concentration

*C*(

*x, y, z, t*) with velocities $\mathit{v}={[\begin{array}{ccc}u& v& w\end{array}]}^{T}$ [m/s] in three directions x, y and z, respectively can be represented effectively with Eq. (1), by using the concept of the concentration gradient vector

*∇C*(

*x, y, z, t*) in the coordinate system

*Oxyz*as demonstrated in Fig. 4.

*∇C*and other relevant mathematical operators can be used as well to formulate the dispersion process model of a pollutant dispersion in the leachate.

*C*(

*x, y, z, t*) of a pollutant at a given point (

*x, y, z*) is low somewhere compared to the surrounding points, physically the substance diffuses into the point (

*x, y, z*) from the surroundings, so the value of the concentration

*C*(

*x, y, z, t*) is increased. Conversely, if the concentration

*C*(

*x, y, z, t*) is high as compared with the surroundings, then the substance moves away and the concentration of the pollutant at the point (

*x, y, z*) is decreased. Finally, the net dispersion is proportional to the second derivative of the concentration function if the dispersion coefficient $\mathit{D}={[\begin{array}{ccc}{D}_{x}& {D}_{y}& {D}_{z}\end{array}]}^{T}$ is a constant vector.

##### (2)

$$\frac{\partial C}{\partial t}=-\mathit{\nabla}C\xb7v+(\mathit{\nabla}\xb7)(D\odot \mathit{\nabla}C)$$*∇*·)(

**⊙**

*D**∇C*), which represents the pollutant transport by dispersion, the notation (

*∇*·) represents the divergence of the gradient field

*∇C*, and the operator ⊙ is the element-wise product (the Hadamard product) of two vectors.

*R*which represents the change of the concentration

*C*(

*x, y, z, t*) in the soil media space due to other chemical reactions and physicochemical interactions, the governing equation Eq. (2) can be rewritten in a generic form as follows:

##### (3)

$$\frac{\partial C}{\partial t}=-\mathit{\nabla}C\xb7v+(\mathit{\nabla}\xb7)(D\odot \mathit{\nabla}C)-R$$*R*≃ 0. In other words, expressing Eq. (3) in terms of the partial derivatives of the concentration function

*C*(

*x, y, z, t*) yields

##### (4)

$$\begin{array}{c}\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+w\frac{\partial C}{\partial y}+v\frac{\partial C}{\partial z}\\ =\frac{\partial}{\partial x}\left({D}_{x}\frac{\partial C}{\partial x}\right)+\frac{\partial}{\partial y}\left({D}_{y}\frac{\partial C}{\partial y}\right)+\frac{\partial}{\partial z}\left({D}_{z}\frac{\partial C}{\partial z}\right)\end{array}$$### 2.3. Finite Element Algorithm to Solve Numerically the Governing Equation over Time in 3D Space for Effective Analysis of Variation of Pollutants Concentrations

*C*(

*x, y, z, t*) simultaneously depends on 4 variables:

*x, y, z*and

*t*. Mathematically, the variation of

*C*(

*x, y, z, t*) over time t and x, y, z is explicitly described by the specific derivative relations as showed in Eq. (4). Moreover, the variation of

*C*(

*x, y, z, t*) depends on 6 parameters as u, w, v, Dx, Dy and Dz. In practice, these parameters can be changed in soil domain even in time. Due to such the spatial and time-dependent complexity of the mathematical model, in this study, it is proposed that a FEM technique to solve the governing equation over time

*t*in 3D space x, y, z with an assumption that the permeability and dispersion of the pollutant concentration are similar in 2 horizontal directions x and y. To discrete the equation with respect to a 3D space of finite elements, the finite difference method is used for both time variable and space variables x and z. This method can also be called as the finite element method with evenly spaced grid nodes at the distance

*Δx*and

*Δz*. Accordingly, in a time that is small enough

*Δt*, at a distance

*Δx*and

*Δz*small enough, the governing equation Eq. (4) at time t = n is discrete as follows:

##### (5)

$$\begin{array}{c}\frac{{C}_{i,j}^{n+1}-{C}_{i,j}^{n-1}}{2\mathit{\Delta}t}+\left(\frac{{u}_{i+1,j}-{u}_{i-1,j}}{2\mathit{\Delta}x}\right){C}_{i,j}^{n}+\frac{u}{2\mathit{\Delta}x}{C}_{i+1,j}^{n}-\\ \frac{u}{2\mathit{\Delta}x}{C}_{i-1,j}^{n}+\frac{{v}_{i,j+1-}{v}_{i,j-1}}{2\mathit{\Delta}z}{C}_{i,j}^{n}+\frac{v}{2\mathit{\Delta}z}{C}_{i,j+1}^{n}-\frac{v}{2\mathit{\Delta}z}{C}_{i,j-1}^{n}=\\ \frac{Dx}{\mathit{\Delta}{x}^{2}}{C}_{i+1,j}^{n}-\frac{Dx}{\mathit{\Delta}{x}^{2}}{C}_{i-1,j}^{n}+\frac{Dz}{\mathit{\Delta}{z}^{2}}{C}_{i+1,j}^{n}-\\ \frac{Dz}{\mathit{\Delta}{z}^{2}}{C}_{i-1,j}^{n}-\left(\frac{Dx}{\mathit{\Delta}{x}^{2}}+\frac{Dz}{\mathit{\Delta}{z}^{2}}\right)\hspace{0.17em}{C}_{i,j}^{n+1}-\left(\frac{Dx}{\mathit{\Delta}{x}^{2}}+\frac{Dz}{\mathit{\Delta}{z}^{2}}\right)\hspace{0.17em}{C}_{i,j}^{n-2}\end{array}$$*u*

*and are*

_{i,j}*v*

*the permeability velocities in 2 directions x and z at the grid node i, j.*

_{i,j}##### (6)

$$\begin{array}{c}\left(\frac{1}{2\mathit{\Delta}t}+\frac{Dx}{\mathit{\Delta}{x}^{2}}+\frac{Dz}{\mathit{\Delta}{z}^{2}}\right)\hspace{0.17em}{C}_{i,j}^{n+1}=\\ \left(\frac{u}{2\mathit{\Delta}x}+\frac{Dx}{\mathit{\Delta}{x}^{2}}\right)\hspace{0.17em}{C}_{i-1,j}^{n}+\left(\frac{v}{2\mathit{\Delta}z}+\frac{Dz}{\mathit{\Delta}{z}^{2}}\right)\hspace{0.17em}{C}_{i,j-1}^{n}+\\ \left(\frac{1}{2\mathit{\Delta}t}-\frac{{u}_{i+1,j-}{u}_{i-1,j}}{2\mathit{\Delta}x}-\frac{{v}_{i,j+1-}{v}_{i,j-1}}{2\mathit{\Delta}z}\right)\hspace{0.17em}{C}_{i,j}^{n}-\\ \left(\frac{Dx}{\mathit{\Delta}{x}^{2}}+\frac{Dz}{\mathit{\Delta}{z}^{2}}\right)\hspace{0.17em}{C}_{i,j}^{n-1}+\left(\frac{Dx}{\mathit{\Delta}{x}^{2}}+\frac{u}{2\mathit{\Delta}x}\right)\hspace{0.17em}{C}_{i+1,j}^{n}+\\ \left(\frac{Dz}{\mathit{\Delta}{z}^{2}}-\frac{v}{2\mathit{\Delta}z}\right)\hspace{0.17em}{C}_{i,j+1}^{n}\end{array}$$##### (7)

$$A{C}_{i,j}^{n+1}=B{C}_{i-1,j}^{n}+C{C}_{i,j-1}^{n}+D{C}_{i,j}^{n}+E{C}_{i,j}^{n-1}+F{C}_{i+1,j}^{n}+G{C}_{i,j+1}^{n}$$$A=\frac{1}{2\mathit{\Delta}t}+\frac{Dx}{\mathit{\Delta}{x}^{2}}+\frac{Dz}{\mathit{\Delta}{z}^{2}}$;

$B=\frac{v}{2\mathit{\Delta}z}+\frac{Dz}{\mathit{\Delta}{z}^{2}}$;

$C=\frac{u}{2\mathit{\Delta}x}+\frac{Dx}{\mathit{\Delta}{x}^{2}}$;

$D=\frac{1}{2\mathit{\Delta}t}-{u}_{x}-{v}_{z}$;

$E=-\left(\frac{Dz}{\mathit{\Delta}{z}^{2}}+\frac{Dx}{\mathit{\Delta}{x}^{2}}\right)$;

$F=\frac{Dz}{\mathit{\Delta}{z}^{2}}+\frac{v}{2\mathit{\Delta}t}$;

$G=\frac{Dx}{\mathit{\Delta}{x}^{2}}-\frac{u}{2\mathit{\Delta}x}$

*C*(

*x, y, z, t*) in the x, z directions over time t. In this study, the following FEM algorithm is constructed and the Matlab programming syntaxes are employed to implement the proposed algorithm.

**FEM Algorithm**: Calculating ${C}_{i,j}^{t}$ at a grid node i, j at a time n

**Inputs**: Initial concentration

*C*

_{0}; Parameters

*Dx*=

*Dy, Dz, u*=

*v, w*; Time step

*Δt*; Increments of the grid

*Δx*=

*Δy, Δz*; Upper limitations of the increments

*n*

_{t}*, n*

_{x}*, n*

*;*

_{z}**Output**: ${C}_{i,j}^{t}$

### 3. Results and Discussion

*∇C*(

*x, y, z, t*), the divergence of the gradient field ( ) and the element-wise product operator ⊙, the advection–diffusion mechanism of a heavy metal contaminant is formulated explitly and effectively, which briefly and clearly describe the physical nature of the disperation process. In addition, the use of FEM for constructing the numerical algorithm is robust and efficient in the aspects of mathematical computation and simulation. The theoretical investigation and the implemented case study show that, the proposed mathematical modelling procedure and the FEM-based algorithms can be potentially and effectively applied to model and predict the spread of different heavy metals in the MSW landfill sites.

### 4. Conclusions

*∇C*(

*x, y, z, t*), the divergence of the gradient field ( ) and the element-wise product operator ⊙, the advection – diffusion model of the leachate migration is formulated in a compact and simplified manner.