Soil-Skeleton and Soil-Water Heavy Metal Contamination by Finite Element Modelling with Freundlich Isotherm Adsorption Parameters

Abstract

World research results indicate that untreated leachate contains high contents of heavy metals (HM) that are likely to pollute the soil and groundwater (GW) environment and contribute to the increase of HMs in soil and GW. The Freundlich isotherm adsorption parameters are essential to soil input parameters for modelling of HMs’ transport to access the soil skeleton and soil pore water contamination by HMs. Finite element (FE) modelling of advection-dispersion transport of HMs by GW movement along with Freundlich isotherm adsorption parameters which continuously change with space in the model domain and with time is sophisticated to accurately evaluate the HMs’ concentrations in soil skeleton and pore water. The chapter describes the background of the existing isotherm adsorption theory, the adaptation of the Freundlich isotherm adsorption in the soil skeleton and soil pore water contamination by HMs, method of determination of the Freundlich isotherm adsorption parameters, the FE procedure of modelling of advection-dispersion transport of HMs by GW movement in general and along with Freundlich isotherm adsorption parameters in particular. A case study modelling has been demonstrated.

1. Introduction

In the world, in many study areas, most of the heavy metal (HM) content in the samples did not exceed the local standards, which can be mentioned as the study of [1] conducted a determination of soil samples of HM components in the North-Western area of Thessaloniki, North Greece near the insanitary landfill. Soil samples were taken at a depth of 2.5–17.5 m. Heavy metals Cd, Cr, Cu, Ni, Pb and Zn were analysed, however, although this area has a high degree of industrialisation, the soil is not contaminated by HMs. Research by Agamuthu and Fauziah [2] conducted sampling of typical soil at different locations in two landfills of Panchang Bedena and Kelana Jaya (Malaysia) to analyse HM content. Soil sampling depth is 2 m–35 m. The analysis results of samples at the Panchang Bedena landfill showed that all the analysed HMs have lower concentrations than the Dutch standard [3]. Opaluwa et al. [4] studied HMs in the soil at a depth of 0–15 cm and leaves on the campus of the Agriculture Faculty of National Polytechnic University, Nasarawa State, Nigeria and the areas near landfills and got the soil concentration of metals As, Cd, Co, Cu, Fe, Ni, Pb and Zn lower than the levels allowed by the World Health Organisation. Similarly, the translocation of HMs from the polluted soil to the aboveground parts of plants and lichens leads to a 1.5- to 5-fold increase in the content of HMs in all species, which fortunately neither exceed the toxicity threshold and nor prevent their growth in the experimental plots as by the work of Lyanguzova et al. [5]. According to Piyada and Suksaman [6], the landfill on the edge of Nai Muang Phichai district in Uttaradit Province, Thailand is one of the most polluted landfills in the world, leachate has a high content of HMs. Samples of soil and groundwater from the landfill and in the vicinity of the landfill are collected and analysed for HM content. Fortunately, the concentrations of Cd, Pb, Cu, Zn and Fe are all low, within the limits of soil quality standards. Research by Siti et al. [7], in Selangor state, Malaysia, there are 20 landfills, including the Ampar Tenang landfill closed in January 2010. However, the landfill is not covered with protective soil according to operational design standards, and before that the garbage is dumped directly onto the ground without any insulating material. Surface soil is relatively polluted by As, Pb, Fe, Cu and Al. As and Pb concentrations are greater than the allowable levels which are 5.90 mg/kg and 31.0 mg/kg, respectively. In addition, only Cu tends to decrease concentration with depth. Kamarudin et al. [8] studied the distribution of HMs in underground aquifers in the solid waste treatment area in Taiping, Perak, Malaysia. Soil samples were taken in 6 boreholes at a depth of 6 m–30 m, and a sampling distance with depth was 1 m. HMs Pb, Mn, Cr, Fe, Zn and Cd were analysed by ICP-MS. The concentrations of Pb, Mn, Fe and Zn are quite high, exceeding the allowable concentrations in the drinking water standard.

It can be seen that most of the HM contents in the soil environment of the landfill do not exceed the allowable levels. In cases where the content of HMs exceeds the standards, there will be an insanitary landfill that does not operate properly with waste burial techniques.

From the above review of several studies, it can be seen that untreated leachate, containing high levels of HMs, is a clear pollution source. Best. Leachate from landfills is capable of polluting the soil and groundwater environment if there are no measures for collection and treatment and to prevent the release of leachate to the surrounding environment. Most studies show that soil at the depth of 2 m and greater is not contaminated with HMs at a level greater than the allowable levels for agricultural land. However, the transport of HMs in the upper soil layer is extremely slow, especially thanks to the effect of adsorption.

2. Modelling of heavy metal transport in soil by groundwater movement

2.1 Theory of heavy metal advection-dispersion transport with soil adsorption

The general two-dimensional partial differential equation of the contaminant transport by advection-dispersion is as follows [9]:

$$D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}-v_x\frac{\partial C}{\partial x}-v_y\frac{\partial C}{\partial y}+Q=\frac{1-n}{n}{\rho }_s\frac{\partial q_e}{\partial t}+\frac{\partial C}{\partial t}$$E1

where C is the contaminant concentration (M/L³, e.g., mg/L); t is the time (d); D_x and D_y are the hydrodynamic dispersion coefficient in x and y direction, respectively (m²/d); x and y are the distances (m), v_x and v_y are the seepage velocity in x and y direction, respectively (m/d); Q is the distributed source of contaminant (mg/d); ρ_s is the solid particle density (note that ρ_d = ρ_s(1-n) in which ρ_d is the unit weight of the dry soil); n is the soil porosity and q_e is the contaminant mass adsorbed per adsorbent mass (mg/kg).

Solid particles are capable of adsorption of dissolved ions of HMs in the soil pore water. The two most common models used to represent the adsorption isotherm are Freundlich and Langmuir isotherms [10]. The Freundlich isotherm is the most common isotherm model, used to describe physical adsorption in a solid-liquids system [11]. Besides, the Langmuir isotherm includes the maximum adsorption capacity of the considered soil, which requires a further special study for the study site.

Freundlich’s adsorption isotherm is used in this study and is described as follows [12, 13] (refer to Patiha et al.):

$$q_e=K_F{C}^{\mathit{1}/\eta }$$E2

where C is the concentration in solution at equilibrium (mg/L); K_F and 1/η are fitting constants [13] and K_F is termed as the Freundlich coefficient (adsorption coefficient) (the unit for the Freundlich constant is ${\mathrm{mg}}^{\frac{1-\eta }{\eta }}{l}^{\frac{1}{\eta }}$/kg) and 1/η is the adsorption intensity (dimensionless). The value of K_F is obtained from the intercept and 1/η from the slope of the logarithmic plot of log q_e versus log C of the equation:

$$\mathit{log}{q}_e=\mathit{log}{K}_F+\frac{1}{\eta }\mathit{log}C$$E3

$$\frac{d{q}_e}{\mathit{dC}}=K_F\frac{1}{\eta }{C}^{\frac{1-\eta }{\eta }}=K_d$$E4

where K_d is the distribution coefficient.

K_F is the Freundlich constant and 1/η depends on the linearity of the isotherm and varies between 0 and 1. Only when 1/η = 1, the isotherm is linear and K_F = K_d.

From (4) the source term $\frac{{\rho }_s}{n}\frac{\partial q_e}{\partial t}$ in (1) is:

$$\frac{1-n}{n}{\rho }_s\frac{\partial q_e}{\partial t}=\frac{1-n}{n}{\rho }_s{K}_F\frac{1}{\eta }{C}^{\frac{1-\eta }{\eta }}\frac{\partial C}{\partial t}$$E5

Therefore, the Eq. (1) may be written in the following form:

$$D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}+D_z\frac{{\partial }^{2}C}{\partial x^2}-v_x\frac{\partial C}{\partial x}-v_y\frac{\partial C}{\partial y}-v_z\frac{\partial C}{\partial z}=\left(1+\frac{1-n}{n}{\rho }_s{K}_F\frac{1}{\eta }{C}^{\frac{1-\eta }{\eta }}\right)\frac{\partial C}{\partial t}$$E6

The so-called coefficient of retardation (retardation factor) R is also used:

$$R=1+\frac{1-n}{n}\rho K_F\frac{1}{\eta }{C}^{\frac{1}{\eta }-1}=1+\frac{1-n}{n}\rho K_d$$E7

$$D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}-v_x\frac{\partial C}{\partial x}-v_y\frac{\partial C}{\partial y}+Q=R\frac{\partial C}{\partial t}$$E8

The retardation factor R is always greater or equal to 1. It is equal to 1 when 1/η = 1.

The partial differential equation of the contaminant transport by advection-dispersion equation is subject to initial and boundary conditions for a particular problem in reality over a certain domain. The initial condition defines the known contaminant concentration over the whole domain at the initial time t = t₀:

$$C=C_o\left(x,y,z\right)$$E9

The boundary condition (BC) would be one of the following kinds:

• The first kind BC (the Dirichlet BC) defines a known concentration on the boundary:

$$C=C_c\mathrm{on}{\Gamma }_c$$E10

• The second kind BC (the Neumann BC) defines a known gradient of contaminant concentration across the boundary:

$$\frac{\partial C}{\partial n}=J_C\mathrm{on}{\Gamma }_{\mathit{JC}}$$E11

• The third kind BC (the Cauchy BC) defines a known rate of contaminant flux through the boundary:

$$v_{n}C-D\frac{\partial C}{\partial n}=q_c\mathrm{one}{\Gamma }_{\mathit{q\upsilon c}}$$E12

where v_n is the velocity normal to the boundary and C is the contaminant concentration at the boundary.

Eq. (8) has an analytical solution only for simple domain configurations, unchanged boundary conditions and constant spatial and temporal values of parameters, i.e., hydrodynamic dispersion coefficient, seepage velocity and retardation factor. Among the transport parameters, the retardation factor is the most sensitive and variable value in time and space as it is a non-linear function of the HM concentrations. This issue always needs to be kept in mind in numerical simulation of solute transport in groundwater in a soil medium with adsorption ability. Numerical methods, e.g., the finite element method (FEM), are capable of solving the equation for any domain configuration, spatial and temporal changing boundary conditions and parameters’ values.

Due to the adsorption isotherm, to more accurately estimate retardation and contaminant transport other than the use of a single value is required in accordance with the relationship in Eq. (7). However, defining transport in terms of a retardation coefficient based on nonlinear adsorption could be complicated. Therefore, Coles [14] examined how the Freundlich model can be used to predict retardation by presenting a simpler way of accounting for nonlinear adsorption and by employing a more appropriate parameter than the Freundlich constant. The linear distribution coefficient K_d, was used by the author to calculate the retardation factor R. Based on the results, the author concluded that the actual ratio of adsorbate to adsorbent may be smaller by a factor of about 10 at higher contaminant concentrations, it could be safer and more accurate to make use of the unified sorption variable K_F to calculate R. Since K_F changes constantly with C and using a constantly changing K_F would be complicated, one solution is to select a small number of discrete values of K_F that can be used to approximate and slightly underestimate the actual values of K_F and each of these values can be used to calculate R over the range of contaminant concentrations that they are applicable.

2.2 FEM of the heavy metal advection-dispersion transport with adsorption by the soil

Let the domain Ω bed be divided into a number of elements E with a total number of nodes M. Let us temporarily not consider the right-hand side term R∂C/∂t of Eq. (8), take the weighting of Eq. (8) and let it be zero [15]:

$${\int }_{\Omega }\left(D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}D-v_x\frac{\partial C}{\partial x}-v_y\frac{\partial C}{\partial y}+Q\right)W_{\ell }\mathit{\text{dxdy}}=0\ell =1,2,\dots ,M$$E13

where $W_{\ell }$ is the weighting function.

Using the Green lemma:

$${\int }_{\Omega }\left(D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}\right)W_{\ell }\mathit{\text{dxdy}}=-{\int }_{\Omega }\left(D_x\frac{\partial W_{\ell }}{\partial x}\frac{\partial C}{\partial x}+D_y\frac{\partial W_{\ell }}{\partial y}\frac{\partial C}{\partial y}\right)\mathit{\text{dxdy}}+{\int }_{\Gamma }\left(D_x\frac{\partial C}{\partial x}{n}_x+D_y\frac{\partial C}{\partial y}{n}_y\right)W_{\ell }\mathit{d\Gamma }$$E14

The integral ∫_Γ is present only for the elements having sides in boundary Γ_qc or Γ_qυc which are generally denoted as Γ_q ta có:

$$-{\int }_{\Omega }\left(D_x\frac{\partial W_{\ell }}{\partial x}\frac{\partial C}{\partial x}+D_y\frac{\partial W_{\ell }}{\partial y}\frac{\partial C}{\partial y}\right)\mathit{\text{dxdy}}-{\int }_{\Omega }\left({\upsilon }_x\frac{\partial C}{\partial x}+{\upsilon }_y\frac{\partial C}{\partial y}\right)W_{\ell }\mathit{\text{dxdy}}=+{\int }_{\Gamma }\left(D_x\frac{\partial C}{\partial x}{n}_x+D_y\frac{\partial C}{\partial y}{n}_y\right)W_{\ell }\mathit{d\Gamma }=0$$E15

With the approximation function of the contaminant concentration is as follows [15]:

$$C\approx \widehat{C}=\sum _{m=1}^M{C}_m{N}_m$$E16

where C_m is the approximation of the contaminant concentration at node m and N_m is the shape functions.

Equation (15) becomes:

$$\begin{array}{c}-{\int }_{\Omega }\left(D_x\frac{\partial W_{\ell }}{\partial x}\frac{\partial N_m}{\partial x}{C}_m+D_y\frac{\partial W_{\ell }}{\partial y}\frac{\partial N_m}{\partial y}{C}_m\right)\mathit{\text{dxdy}}-{\int }_{\Omega }\left(v_x\frac{\partial N_m}{\partial x}{W}_{\ell }{C}_m+v_y\frac{\partial N_m}{\partial y}{W}_{\ell }{C}_m\right)\mathit{\text{dxdy}}\\ ={\int }_{{\Gamma }_q}\left({\overline{q}}_x{W}_{\ell }{n}_x+{\overline{q}}_{x}W{n}_y\right)\mathit{d\Gamma }=0\end{array}$$E17

$$\mathit{K}=-{\int }_{\Omega }\left(D_x\frac{\partial W_{\ell }}{\partial x}\frac{\partial N_m}{\partial x}+D_y\frac{\partial W_{\ell }}{\partial y}\frac{\partial N_m}{\partial y}\right)\mathit{\text{dxdy}}-{\int }_{\Omega }\left(v_x\frac{\partial N_m}{\partial x}{W}_{\ell }+v_y\frac{\partial N_m}{\partial y}{W}_{\ell }\right)\mathit{\text{dxdy}}$$E18

$$\mathit{F}=-{\int }_{{\Gamma }_q}\left({\overline{q}}_c{W}_{\ell }{n}_x+{\overline{q}}_c{W}_{\ell }{n}_y\right)\mathit{d\Gamma }$$E19

$$\mathrm{KC}=\mathrm{F}$$E20

The shape functions Nm and weighting functions $W_{\ell }$ can be linear or higher order functions. For unsteady problems, i.e., R∂C/∂t ≠ 0 we have:

$$\mathit{E}\left(R\frac{\mathit{dC}}{\mathit{dt}}\right)+\mathit{KC}=\mathit{F}$$E21

The square matrix E is:

$$\mathit{E}={\int }_{\Omega }R{W}_{\ell }{N}_m\mathit{\text{dxdy}}$$E22

Eq. (21) has the following general form in regard to temporal derivative:

$$\left(\mathit{\theta K}\mathit{+}\frac{\mathit{E}}{\mathit{\Delta t}}\right){\mathit{C}}^{\mathit{t}\mathit{+}\mathit{\Delta t}}\mathit{=}\left(\mathit{-}\left(\mathit{1}\mathit{-}\mathit{\theta }\right)\mathit{K}\mathit{+}\frac{\mathit{E}}{\mathit{\Delta t}}\right){\mathit{C}}^{\mathit{t}}\mathit{+}\frac{\mathit{1}}{\mathit{2}}\left({\mathit{F}}^{\mathit{t}}\mathit{+}{\mathit{F}}^{\mathit{t}\mathit{+}\mathit{\Delta t}}\right)$$E23

The typical schemes are:

1. Forward difference (θ = 0):

   $$R\frac{\left[E\right]}{\mathit{\Delta t}}\left\{C^{t+\mathit{\Delta t}}\right\}+\left(⌈K⌉-R\frac{\left[E\right]}{\mathit{\Delta t}}\right)\left\{C^t\right\}=\left\{F^t\right\}$$E24

2. Backward difference (θ = 1):

   $$\left(\left[\mathit{K}\right]\mathit{+}\frac{\left[\mathit{E}\right]}{\mathit{\Delta t}}\right)\left\{{\mathit{C}}^{\mathit{t}\mathit{+}\mathit{\Delta t}}\right\}\mathit{-}\frac{\left[\mathit{E}\right]}{\mathit{\Delta t}}\left\{{\mathit{C}}^{\mathit{t}}\right\}\mathit{=}\left\{{\mathit{F}}^{\mathit{t}\mathit{+}\mathit{\Delta t}}\right\}$$E25

3. Central difference (the Crank–Nicolson scheme) (θ = 0.5):

$$\left(\frac{\mathit{1}}{\mathit{2}}\left[\mathit{K}\right]\mathit{+}\frac{\left[\mathit{E}\right]}{\mathit{\Delta t}}\right)\left\{{\mathit{C}}^{\mathit{t}\mathit{+}\mathit{\Delta t}}\right\}\mathit{+}\left(\frac{\mathit{1}}{\mathit{2}}\left[\mathit{K}\right]\mathit{-}\frac{\left[\mathit{E}\right]}{\mathit{\Delta t}}\right)\left\{{\mathit{C}}^{\mathit{t}}\right\}\mathit{=}\frac{\mathit{1}}{\mathit{2}}\left(\left\{{\mathit{F}}^{\mathit{t}}\right\}\mathit{+}\left\{{\mathit{F}}^{\mathit{t}\mathit{+}\mathit{\Delta t}}\right\}\right)$$E26

Let us consider two-dimensional in xy direction problems. The domain is divided into a mesh of triangular or quadrangular finite elements. For the mesh of rectangular elements (Figure 1) which have the side of h_x and h_y.

Figure 1.

A mesh of rectangular elements.

In the above equations, the matrix K at the element level is a square matrix Ne × Ne in which Ne is the number of the vertices of the elements (square matrix 3 × 3 or 4 × 4 for triangular or quadrangular elements, respectively). As an illustration, Galerkin FEM with linear shape functions and for a rectangular element with nodes i, j, k and l numbered counter-clockwise (Figure 1) which has sides of hxe and hye the matrices K, E and F are determined as follows. Since each term of the matrix is very long, each column containing Ne rows is to be written (columns 1, 2, 3, 4 are denoting nodes i, j, k and l, respectively, and rows 1, 2, 3, 4 are denoting nodes i, j, k and l, respectively):

The four-row terms of column i are:

$${\mathit{K}}_{\left(\mathit{i},\mathit{j},\mathit{k},\mathit{l}\right)\mathit{i}}^{\mathit{e}}=\left[\begin{array}{c}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{i}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{i}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{i}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{i}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{j}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{j}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{j}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{j}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{k}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{k}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{k}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{k}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{l}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{l}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{l}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{l}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{i}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\end{array}\right.$$E27

The four-row terms of column j are:

$${\mathit{K}}_{\left(\mathit{i},\mathit{j},\mathit{k},\mathit{l}\right)\mathit{j}}^{\mathit{e}}\mathit{=}\begin{array}{c}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{i}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{W}}_{\mathit{i}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\mathit{+}{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{i}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{i}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{j}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{W}}_{\mathit{j}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\mathit{+}{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{j}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{j}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{k}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{W}}_{\mathit{k}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\mathit{+}{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{k}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{k}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{l}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{W}}_{\mathit{l}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\mathit{+}{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{l}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial x}}\mathit{+}{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{l}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{j}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\end{array}$$E28

The four-row terms of column k are:

$$K_{\left(i,j,k,l\right)k}^e=\begin{array}{c}{\int }_0^{h_y^e}{\int }_0^{h_x^e}\left[W_i{v}_x\frac{\partial N_k}{\partial x}+W_i{v}_y\frac{\partial N_k}{\partial y}+D_x\frac{\partial W_i}{\partial x}\frac{\partial N_k}{\partial x}+D_y\frac{\partial W_i}{\partial y}\frac{\partial N_k}{\partial y}\right]\mathit{\text{dxdy}}\\ {\int }_0^{h_y^e}{\int }_0^{h_x^e}\left[W_j{v}_x\frac{\partial N_k}{\partial x}+W_j{v}_y\frac{\partial N_k}{\partial y}+D_x\frac{\partial W_j}{\partial x}\frac{\partial N_k}{\partial x}+D_y\frac{\partial W_j}{\partial y}\frac{\partial N_k}{\partial y}\right]\mathit{\text{dxdy}}\\ {\int }_0^{h_y^e}{\int }_0^{h_x^e}\left[W_k{v}_x\frac{\partial N_k}{\partial x}+W_k{v}_y\frac{\partial N_k}{\partial y}+D_x\frac{\partial W_k}{\partial x}\frac{\partial N_k}{\partial x}+D_y\frac{\partial W_k}{\partial y}\frac{\partial N_k}{\partial y}\right]\mathit{\text{dxdy}}\\ {\int }_0^{h_y^e}{\int }_0^{h_x^e}\left[W_l{v}_x\frac{\partial N_k}{\partial x}+W_l{v}_y\frac{\partial N_k}{\partial y}+D_x\frac{\partial W_l}{\partial x}\frac{\partial N_k}{\partial x}+D_y\frac{\partial W_l}{\partial y}\frac{\partial N_k}{\partial y}\right]\mathit{\text{dxdy}}\end{array}$$E29

The four-row terms of column l are:

$${\mathit{K}}_{\left(\mathit{i},\mathit{j},\mathit{k},\mathit{l}\right)\mathit{l}}^{\mathit{e}}=\left.\begin{array}{c}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{i}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{i}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{i}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{i}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{j}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{j}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{j}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{j}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{k}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{k}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{k}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{k}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\\ {\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left[{\mathit{W}}_{\mathit{l}}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{W}}_{\mathit{l}}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}+{\mathit{D}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{l}}}{\mathit{\partial x}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial x}}+{\mathit{D}}_{\mathit{y}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{l}}}{\mathit{\partial y}}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial y}}\right]\mathit{\text{dxdy}}\end{array}\right]$$E30

The Galerkin FEM with linear shape functions results in:

The four-row terms of column I, the Eq. (27) become:

$${\mathit{K}}_{\left(\mathit{i},\mathit{j},\mathit{k},\mathit{l}\right)\mathit{i}}^{\mathit{e}}\mathit{=}\left[\begin{array}{l}\left[{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left\{\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right\}}^{\mathit{2}}\mathit{+}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left\{\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right\}}^{\mathit{2}}\right]\mathit{+}\left[\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left(\mathit{-}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right)\right]\\ \left[\mathit{-}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}\frac{\mathit{x}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right]\mathit{+}\left[\frac{\mathit{x}\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left(\mathit{-}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right)\right]\\ \left[\mathit{-}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}\frac{\mathit{y}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}\frac{\mathit{x}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right]\mathit{+}\left[\frac{\mathit{xy}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left(\mathit{-}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right)\right]\\ \left[{\mathit{D}}_{\mathit{x}}^{\mathit{e}}\frac{\mathit{y}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right]\mathit{+}\left[\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)\mathit{y}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left(\mathit{-}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{-}\mathit{y}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\mathit{-}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}\frac{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\mathit{-}\mathit{x}\right)}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\right)\right]\end{array}\right.$$E31

The four-row terms of column j, the Eq. (28) become:

$${\mathit{K}}_{\left(\mathit{i},\mathit{j},\mathit{k},\mathit{l}\right)\mathit{j}}^{\mathit{e}}\mathit{=}\begin{array}{c}\frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[\mathit{-}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}\mathit{+}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]\mathit{+}\frac{\mathit{1}}{\mathit{12}}\left(\mathit{-}\mathit{2}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{+}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\\ \frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}\mathit{+}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]\mathit{+}\frac{\mathit{1}}{\mathit{12}}\left(\mathit{-}\mathit{2}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{+}\mathit{2}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\\ \frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}\mathit{-}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]\mathit{+}\frac{\mathit{1}}{\mathit{12}}\left(\mathit{-}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{+}\mathit{2}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\\ \frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[\mathit{-}\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}\mathit{-}\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]\mathit{+}\frac{\mathit{1}}{\mathit{12}}\left(\mathit{-}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}\mathit{+}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\end{array}$$E32

The four-row terms of column k, the Eq. (29) become:

$$K_{\left(i,j,k,l\right)k}^e=\begin{array}{c}\frac{1}{3h_x^e{h}_y^e}\left[-\frac{1}{2}{D}_x^e{\left(h_x^e\right)}^2-\frac{1}{2}{D}_y^e{\left(h_y^e\right)}^2\right]+\frac{1}{12}\left(-v_x^e{h}_y^e-v_y^e{h}_x^e\right)\\ \frac{1}{3h_x^e{h}_y^e}\left[\frac{1}{2}{D}_x^e{\left(h_x^e\right)}^2-D_y^e{\left(h_y^e\right)}^2\right]+\frac{1}{12}\left(-v_x^e{h}_y^e-2v_y^e{h}_x^e\right)\\ \frac{1}{3h_x^e{h}_y^e}\left[D_x^e{\left(h_x^e\right)}^2+D_y^e{\left(h_y^e\right)}^2\right]+\frac{1}{12}\left(-2v_x^e{h}_y^e-2v_y^e{h}_x^e\right)\\ \frac{1}{3h_x^e{h}_y^e}\left[-D_x^e{\left(h_x^e\right)}^2+\frac{1}{2}{D}_y^e{\left(h_y^e\right)}^2\right]+\frac{1}{12}\left(-2v_x^e{h}_y^e-v_y^e{h}_x^e\right)\end{array}$$E33

The four-row terms of column l, the Eq. (30) become:

$${\mathit{K}}_{\left(\mathit{i},\mathit{j},\mathit{k},\mathit{l}\right)\mathit{l}}^{\mathit{e}}=\left.\begin{array}{c}\frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}-{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]+\frac{\mathit{1}}{\mathit{12}}\left({\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}-\mathit{2}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\\ \frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[-\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}-\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]+\frac{\mathit{1}}{\mathit{12}}\left({\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}-{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\\ \frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[-{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}+\frac{\mathit{1}}{\mathit{2}}{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]+\frac{\mathit{1}}{\mathit{12}}\left(\mathit{2}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}-{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\\ \frac{\mathit{1}}{\mathit{3}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left[{\mathit{D}}_{\mathit{x}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)}^{\mathit{2}}+{\mathit{D}}_{\mathit{y}}^{\mathit{e}}{\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}\right)}^{\mathit{2}}\right]+\frac{\mathit{1}}{\mathit{12}}\left(\mathit{2}{\mathit{v}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}-\mathit{2}{\mathit{v}}_{\mathit{y}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}\right)\end{array}\right]$$E34

The contribution of the loading vector f is determined by taking the integral over the element. For example, for node i at a side along the boundary:

$$\begin{array}{c}{\mathit{F}}_{\mathit{i}}^{\mathit{e}}={\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{{\mathit{Q}}^{\mathit{e}}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}-\mathit{x}\right)\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}-\mathit{y}\right)\mathit{\text{dxdy}}-\underset{¯}{{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\frac{{\bar{\mathit{q}}}^{\mathit{e}}}{{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}\left({\mathit{h}}_{\mathit{x}}^{\mathit{e}}-\mathit{x}\right)\mathit{dx}}-\underset{¯}{{\int }_{\mathit{0}}^{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\frac{{\bar{\mathit{q}}}^{\mathit{e}}}{{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\left({\mathit{h}}_{\mathit{y}}^{\mathit{e}}-\mathit{y}\right)\mathit{dy}}\\ =\frac{\mathit{1}}{\mathit{4}}{\mathit{Q}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}-\underset{¯}{\frac{\mathit{1}}{\mathit{2}}{\bar{\mathit{q}}}^{\mathit{e}}{\mathit{h}}_{\mathit{x}}^{\mathit{e}}}-\underset{¯}{\frac{\mathit{1}}{\mathit{2}}{\bar{\mathit{q}}}^{\mathit{e}}{\mathit{h}}_{\mathit{y}}^{\mathit{e}}}\end{array}$$E35

The underlined terms are existing only if the side i, j and l are lying in boundary Γ _q, and therefore they are not existing for the inside elements:

The matrix E is:

$$\text{node}:ijkl$$

$${\mathit{E}}^1=R_e\left[\begin{array}{cccc}{\int }_{\Omega }{W}_i{N}_i\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_i{N}_j\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_i{N}_k\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_i{N}_{\ell }\mathit{\text{dxdy}}\\ {\int }_{\Omega }{W}_j{N}_i\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_j{N}_j\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_j{N}_k\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_j{N}_{\ell }\mathit{\text{dxdy}}\\ {\int }_{\Omega }{W}_k{N}_i\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_k{N}_j\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_j{N}_k\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_j{N}_{\ell }\mathit{\text{dxdy}}\\ {\int }_{\Omega }{W}_{\ell }{N}_i\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_{\ell }{N}_j\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_{\ell }{N}_k\mathit{\text{dxdy}}& {\int }_{\Omega }{W}_{\ell }{N}_{\ell }\mathit{\text{dxdy}}\end{array}\right]\left[\begin{array}{c}{C}_i\\ C_j\\ C_k\\ C_{\ell }\end{array}\right]$$E36

Putting the weighting function W and shape function N into Eq. (36) results in:

$$\text{node}:ij$$

$${\mathit{E}}^1=R_e\left[\begin{array}{ll}{\int }_{\Omega }\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\\ {\int }_{\Omega }\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\\ {\int }_{\Omega }\frac{\mathit{xy}}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{\mathit{xy}}{h_x^e{h}_y^e}\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\\ {\int }_{\Omega }\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\end{array}\right.$$E37

$$\text{node}:kl$$

$$\left.\begin{array}{cc}{\int }_{\Omega }\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{\mathit{xy}}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{\left(h_x^e-x\right)\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\\ {\int }_{\Omega }\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{\mathit{xy}}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{x\left(h_y^e-y\right)}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\\ {\int }_{\Omega }\frac{\mathit{xy}}{h_x^e{h}_y^e}\frac{\mathit{xy}}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{\mathit{xy}}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\\ {\int }_{\Omega }\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\frac{\mathit{xy}}{h_x^e{h}_y^e}\mathit{\text{dxdy}}& {\int }_{\Omega }\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\frac{\left(h_x^e-x\right)y}{h_x^e{h}_y^e}\mathit{\text{dxdy}}\end{array}\right]\left[\begin{array}{c}{C}_i\\ C_j\\ C_k\\ C_{\ell }\end{array}\right]$$E38

By assembling all the element matrices K, E and F a global system of equations can be obtained the solutions of which are the approximated contaminant concentrations at all nodes.

For linear elements, the element sizes and time steps need to be selected based on the following criteria [16]:

$$\mathit{Pe}=\frac{v_x{h}_x^e}{D_x}\le 2;\mathit{Pe}=\frac{v_y{h}_y^e}{D_y}\le 2;\mathit{Cr}=\frac{v_x\mathit{\Delta t}}{h_x^e}\le 1;\mathit{Cr}=\frac{v_y\mathit{\Delta t}}{h_y^e}\le 1$$E39

2.3 Soil heavy metal adsorption parameters

The adsorption capacities of HMs change with physical parameters such as pH, temperature etc. The adsorption of heavy metals As, Cu, Zn, Pb, Cr, Cd and Hg on the soil at different pH was experimentally investigated by He et al. [17]. The adsorption capacities of Cr decreased with increasing pH, which may be caused by the unique physical properties of Cr. The adsorption capacities of the remaining HMs are increased with increasing pH. One of the reasons is that the increase in pH effectively reduces the concentration of H⁺ in the solution. In solution with pH greater than 7, all the ions H⁺ are in par with ion OH⁻. Therefore, HM ions with a positive charge can more effectively be absorbed by the soil colloids. It means the adsorption capacities of HM ions increased with increasing pH value.

One of the aspects of the influence of pH on the adsorption of HMs by soil particles is that pH has an influence on the solubility of HMs in solution [18] and controls various adsorption reactions on the surface of solid particles, and the increase in pH, which promoted an increase in the adsorption point of the soil colloid since soil colloids generally have a negative charge [19]. The chapter will deal with the HM adsorption capacity at pH around 7.

To investigate the effect of temperature on the adsorption of As, Cu, Zn, Pb, Cr, Cd and Hg, He et al. [17] used soil samples at a different temperature from 30–50°C. The data obtained by the authors show the increase of adsorption capacities of HMs in the soil material with increasing temperature.

The experiment data for Cr, Cu, Pb and Zn by He et al. [17] at the temperature of 25°C have been extracted from the authors’ publication’s figures. The least squared error method was used by the authors of this study for determining the Freundlich coefficient K_F and adsorption intensity 1/η, the results of which are also presented in Table 1 which is showing a very high correlation of the experiment data point and the Freundlich isotherm coefficients K_F and 1/η of the fitting curves (Figure 2).

KF	1/η	R^2	KF	1/η	R^2	KF	1/η	R^2	KF	1/η	R^2
With qe in mg/g
Cr			Cu			Pb			Zn
0.264	0.260	0.999	0.131	0.450	0.991	0.300	0.290	0.996	0.050	0.569	0.999
As			Cd			Hg
0.404	0.340	0.989	0.080	0.420	0.964	0.510	0.300	0.985

Table 1.

Freundlich coefficient KF and adsorption intensity 1/η for the soil of the study of He et al. [17].

Figure 2.

Adsorption capacities of heavy metals adsorbed on soil material at equilibrium concentration.

The Freundlich isotherm coefficients KF and 1/η of silty soils were also studied by some other authors. The study of Noppadoland [20] investigates the adsorption of the most common HMs (Cu, Ni, and Zn) by various soils. Fifteen soil samples were collected from various areas of North-Eastern Thailand.

They were excavated from different depths, ranging from 20 cm to 50 cm below the soil surface. The average soil pH is about 6.5. The areas near watercourses, communities and industries were selected as sites from which the soil samples were taken. The authors have received the following average values of the Freundlich isotherm coefficients K_F and 1/η of the soils: K_F = 0.348 (mg/g) and 1/η = 0.235 (σ = 0.059) for Zn, K_F = 0.462 and 1/η = 0.320 for Cu (σ = 0.054), K_F = 0.279 and 1/η = 0.437 (σ = 0.059) for Ni (with the q_e in mg/g and C unit in mg/L).

Soils in some regions of North-Western Spain have been the subject of agricultural management practices involving the use of fertilisers and various types of organic waste containing HMs. Although such practices have facilitated crop growth, they have also raised the natural contents in HMs of the soils. Therefore, Emma et al. [21] researched the ability of the soils with high concentrations of Cr and Ni to adsorb and retain Cd, Cu, Ni, Pb and Zn. The soil pH is about 6.5 and the experiments’ temperature is 25°C. They have obtained the following Freundlich coefficients: K_F = 1.560 (mg/g) and 1/η = 0.327 for Cu, K_F = 0.363 and 1/η = 0.441 for Ni, K_F = 1.363 and 1/η = 0.351 for Pb, K_F = 0.463 and 1/η = 0.426 for Zn, K_F = 0.540 and 1/η = 0.293 for Cd (with the q_e in mg/g and C unit in mg/L).

Claudia et al. [22] carried out a specific adsorption evaluation through the amounts of adsorbed Cu, Pb, Cr, Ni and Zn after desorption experiments in ten different soils. The HM adsorption isotherm Freundlich parameters at temperature 25°C and for the neutral pH soils are as follows: K_F = 2.540 (mg/g) and 1/η = 0.91 for Cu, K_F = 0.702 and 1/η = 0.510 for Ni, K_F = 0.998 and 1/η = 0.440 for Pb, K_F = 1.016 and 1/η = 0.440 for Zn (with the q_e in mg/g and C unit in mg/L).

2.4 Dispersion parameters

The coefficient of longitudinal hydrodynamic dispersion D_L in the water flow direction which is the coefficient of hydrodynamic dispersion coefficient D in Eq. (1) consists of two components: the coefficient of mechanical dispersion D′ and the coefficient of molecular diffusion in a porous medium D*_d, i.e., D_L = D′ + D*_d [9].

The coefficient of mechanical dispersion D′ depends upon the microstructure of the soil, the seepage velocity and the molecular dispersion in water as follows by Jacob and Arnold [9]:

$$D{\prime }_{\mathit{ij}}=a_L\mathit{vf}\left(\mathit{Pe},\xi \right);\mathit{Pe}=\frac{\mathit{Lv}}{D_d}$$E40

in which: v is the seepage velocity (m/d); Pe is the Peclet number; L is the characteristic length of the pores (m); D_d is the molecular dispersion in water; ξ is the ratio between the pores’ size and the characteristic length through the pores; f(P_e,ξ) = P_e/(P_e+2 + 4ξ²) is a function which is expressing the transport of the HMs or solutes via molecular dispersion between the neighbouring flow streams at the macro scale, and in most cases f(Pe,ξ) ≅ 1; a_L is the longitudinal dispersivity.

For a one-directional groundwater flow, the coefficient of mechanical dispersion D′ is the multiplication of longitudinal dispersivity (a_L) and seepage velocity [9]. The longitudinal dispersivity is of the order of the average soil particle [9], e.g., particle size d₅₀. The coefficient of molecular diffusion in a porous medium D*_d is as follows [23]:

$$D_d^{\ast }=\frac{D_d}{n{F}_R}$$E41

in which: F_R is the formation factor which is specified by the geophysicists as the ratio between soil resistivity and water resistivity.

The formation factor F_R ranges from 0.1 (for clay) to 0.7 (for sand) [23], and always less than 1.

The coefficient of molecular diffusion in water D_d is:

$${\mathit{D}}_{\mathit{d}}\mathit{=}\frac{\mathit{RT}}{\mathit{N}}\frac{{1}}{{6}\mathit{\pi r\mu }}$$E42

In which: R is the gas constant; K is temperature unit Kelvin; N is the Avogadro number; T is the temperature (K); μ is the water viscosity; r is the average radius of the HM or solute.

The coefficients of molecular diffusion coefficients of inorganic cations and anions in water D_d may be found in the book by Henry [24].

Jacob and Arnold [9] have divided dispersion and diffusion into five zones (Figure 7 in [9]) in accordance to the Pectlet number, for which the roles of the molecular diffusion and the hydrodynamic dispersion are described. Zone I is corresponding to very slow water movement with the Pectlet number less than 0.4 so that the molecular diffusion predominates and the mechanical dispersion (D′) is negligible, i.e., D_L ≈ D^*_d. In our case, the Pectlet number is equal to 0.0011, for which the hydrodynamic dispersion D is approximately equal to the molecular diffusion in saturated porous medium D^*_d.

3. The FE modelling application to Kieu Ky landfill, Hanoi, Vietnam

Kieu Ky waste landfill is located in Gia Lam district in the South-East of the Center of Hanoi in the Bac Bo plain, the second largest plain in Vietnam. The waste landfill facility has an area of 13 ha consisting of composting compartments, a leachate reservoir and five landfill cells (Figure 3). The landfill cells have bottoms at the depth of 4.5 m and the thickness of dumped waste of 5 m–15 m (Figure 4). The facility handled 175 tons of solid waste in a day. It is operated from 2002 to 2020. The area is covered by Holocene formation, under which a rich and with good quality Pleistocene aquifer is underlying.

Figure 3.

The layout of landfills, geotechnical boreholes and boreholes for soil quality sampling.

Figure 4.

The soil profile of Kieu Ky landfill site.

Kieu Ky landfill area has a natural ground surface of elevation around 4.5 m above mean sea level. The local shallow geological and hydrogeological conditions are as follows (Figure 4): (1) Surface cultivated soil of about 0.8 m in thickness, which consists of grey-yellow silt with some small construction solid waste pieces, and (2) Layer of Upper Holocene silt of grey-yellow, grey-green and grey-brown colours, the thickness is around 6 m. The silt’s porosity (n) and hydraulic conductivity (K) have been determined and are 0.455 m/d and 0.0045 m/d, respectively.

Two model domains (MD) have been selected: one is the natural soil profile next to the landfill (from the ground surface to the depth of 6 m, i.e. to the groundwater aquifer surface, with the length of 6 m) (MD1) and the second one is the soil profile beneath the bottom of the landfill (from the depth 4.5 m to the surface of the groundwater aquifer with the length of 1.5 m) (MD2). The three characteristic values (minimum, average and maximum) of the Freundlich isotherm adsorption parameters are considered in the two selected model domains.

The hydraulic conditions of the two model domains are determined based on Figure 5 and on that the water level of the Upper Holocene aquifer is 2 m below the ground surface, the level of leachate and the water level of the leachate pond are the same. Domain 1 is a former rice field and almost is constantly wet. This creates a saturated soil profile. Direct leakage of leachate from the landfills to the land slot to supply HMs to penetrate the soil profile. Domain 2 is underneath the bottoms of the landfills, which are lower than the groundwater level of the beneath aquifer. Similar to domain 1, this also creates a saturated soil profile. Besides, it is most likely that some landfill leachate may leak into the domain. The seepage velocity of which is determined by Darcy law through the hydraulic gradient and soil hydraulic conductivity. The soil hydraulic conductivity in the vertical direction was determined by the laboratory permeameter. Subsurface soil samples have been collected for laboratory experiments for the determination of saturated hydraulic conductivity.

Figure 5.

Typical model domains in the study site.

3.1 Parameter calibration

The problem of aquifer parameter calibration has been studied extensively. In modelling of groundwater flow and transport, besides the specification of the aquifer geometry and its boundary conditions, the determination of aquifer’s geohydraulic parameters, e.g., hydraulic conductivity, porosity, dispersivity, source or sink and prescribed boundary fluxes is necessary. The inverse problem of parameter calibration can be defined as the optimal determination of the parameters based on the observation data of the dependent variables, such as hydraulic heads or solute concentrations, collected in space and time. The inverse methods have been classified by Neuman [25] into two groups: indirect and direct. The indirect approach is based on the output error criterion, where the accuracy of the parameters is improved by an iterative process until the model response is close enough to the observed one. The direct approach is based on the use of super-determinate equations derived from rearranging the discretisation equations in such a way that the parameters are considered as unknown variables and their optimal values are such that minimise the residuals of the equations in a certain sense. The modern inverse techniques are often imbedded with the numerical models, usually finite difference and finite element models. All the soil HMs relevant transport parameters may be calibrated simultaneously. However, this would result in a high uncertainty of the obtained calibrated parameters as the overall modelling results may have a good optimisation error while the calibrated parameters are not reliable as they are beyond the physical limits. Therefore, some parameters are better determined by experimental tests and the remaining parameters are calibrated by inverse analyses. This procedure is particularly suitable for the soil adsorption parameters and dispersion parameters of low permeable soils.

Generally, the objective function (E(k)) to be minimised in the inverse analysis can be expressed as the sum of weighted squares of the differences between the model responses and the observation ones and the sum weighted squares of the difference between the estimated model parameter and prior parameter. The indirect method using this kind of objective function is called regularised Output Least Squares (OLS). If the second term of sum weighted squares of the difference between the estimated model parameter and the prior parameter has vanished, e.g., the regulation coefficient is equal to zero, the method is called generalised OLS. In the latter, if the optimal weighting coefficients all are equal to the unit, the method becomes original OLS.

The numerical methods in the solution of OLS problems are unconstrained nonlinear optimisation, which includes search method, gradient method and second-order method (Quasi-Newton methods). Within the chapter, one-dimensional dispersion testing for the determination of soil dispersivity by Quasi-Newton methods is described for demonstration.

3.1.1 One-D dispersion testing for determination of soil dispersivity by Quasi-Newton methods

A tracer column test is illustrated in Figure 6 in which a constant tracer concentration is maintained in the left boundary (a relative concentration of 1 is usually used) and a constant zero-concentration in the other boundary.

Figure 6.

A tracer column test scheme.

The special and temporal concentrations are monitored, for which the observed (C^obs) and model (C^mod) concentration at time t are presented in Figure 7. One-D dispersion determination of the soil dispersivity by Quasi-Newton methods are described as follows.

Figure 7.

Plots of observed vs. model tracer concentration.

The most common criterion in the evaluation of the difference between the model estimated and observed variables is the least squared root given as:

$$\mathit{Min}E\left(\mathit{p}\right)=\sum _{l=1}^L{\left[C_l^{\mathit{mod}}\left(\mathit{p}\right)-C_l^{\mathit{obs}}\right]}^2$$E43

where E(p) - objective function; L- number of observed variables; C_l^mod- model estimated values of the concentration; C_l^obs- observed measured values of the concentration; p- parameters of the physical medium (hydrodynamic dispersion coefficient as a function of pore velocity, porosity and dispersivity).

Let us consider the following multi-dimensional optimisation problem:

$$\mathit{Min}\mathit{E}\left(\mathit{p}\right),{p}\in {{p}}_{\mathit{ct}}$$E44

where p_ct - a set of possible values of parameter variables p. This set of parameter values may be selected based on the existing data of the parameters of similar physical materials, of materials at the same locations, statistical data etc.

If the objective function E(p) has a second-order derivative then the necessary and sufficient conditions for $\widehat{\mathit{p}}$ to be the stationary point, i.e., E(p) has a local extreme value [26]:

Gradient g = ∇E(p) = 0 at p, i.e.:

$${\left|\frac{\mathit{\partial E}}{\mathit{\partial }{\mathit{p}}_{\mathit{m}}}\right|}_{\widehat{\mathit{p}}}\mathit{=}{0}\mathit{;}\left(\mathit{m}\mathit{=}{1},{2},\dots ,\mathit{M}\right)$$E45

where M is the number of parameter variables.

Hessian matrix G = ∇²E(p):

$$G=\left|\begin{array}{cccc}\frac{{\partial }^{2}E}{\partial p_1^2}& \frac{{\partial }^{2}E}{\partial p_1\partial p_2}& \dots & \frac{{\partial }^{2}E}{\partial p_1\partial p_M}\\ \frac{{\partial }^{2}E}{\partial p_2\partial p_1}& \frac{{\partial }^{2}E}{\partial p_2^2}& \dots & \frac{{\partial }^{2}E}{\partial p_2\partial p_M}\\ \dots \dots & \dots \dots & \dots \dots & \dots \dots \\ \frac{{\partial }^{2}E}{\partial p_M\partial p_1}& \frac{{\partial }^{2}E}{\partial p_M\partial p_2}& \dots & \frac{{\partial }^{2}E}{\partial p_M^2}\end{array}\right|$$E46

is a positive definite matrix at $\widehat{\mathit{p}}$.

The optimisation algorithms in the determination of parameters consist of the following steps:

Selection of the initial values of parameters p₀.

Determination of the search sequence: p₀, p₁, p₂, …, p_n … in such a way that E(p_n + 1) < E (p_n) for all n.

Checking the convergence criterion. If the convergence is observed, then the local extremes have been reached and the parameter values are considered to be estimated.

Commonly, the search sequence has the following general form:

$${\mathit{p}}_{\mathit{n}\mathit{+}\mathbf{1}}\mathit{=}{\mathit{p}}_{\mathit{n}}\mathit{+}{\mathit{\lambda }}_{\mathit{n}}{\mathit{d}}_{\mathit{n}}$$E47

where d_n- vector of displacement directions; λ_n- step size (that must be most optimally selected).

Three main groups of optimisation algorithms may be classified for solving optimisation problems: (1) Search method, when only the values of the objective function are considered, (2) Gradient method, when the gradients of the objective function are utilised and (3) Second-order method, if the second derivatives of the objective functions are utilised. The Quasi-Newton method belongs to the third group.

Suppose there is a set of initial values of parameters p₀, it is required to find out the search sequence p₀, p₁, p₂, …, p_n so that E(p_n + 1) < E (p_n) for all n. Gradient vector g(p_n+1) at vicinity of p_n may be determined as follows:

$$\mathit{g}\left({\mathit{p}}_{n+1}\right)\approx {\mathit{g}}_n+{\mathit{G}}_n\Delta \mathit{p};\mathrm{\Delta }{p}=\mathit{p}-{\mathit{p}}_n;{{g}}_n=\mathit{g}\left({\mathit{p}}_n\right);{{G}}_n=\mathit{G}\left({\mathit{p}}_n\right)$$E48

The necessary condition of extreme existence is g(p_n + 1) ≈ 0. This can be done if p_n+1 = p_n + Δ p_n, where Δp_n are the solution of the following equation:

$${\mathit{g}}_{\mathit{n}}\mathit{+}{\mathit{G}}_{\mathit{n}}\mathit{\Delta p}\mathit{=}\mathbf{0}$$E49

$$\mathit{\Delta }{\mathit{p}}_{\mathit{n}}\mathit{=}\mathit{\Delta p}\mathit{=}\mathit{-}{\mathit{g}}_{\mathit{n}}{\mathit{G}}_{\mathit{n}}^{\mathit{-}\mathbf{1}}\mathit{\Rightarrow }{\mathit{p}}_{\mathit{n}\mathit{+}\mathbf{1}}\mathit{=}{\mathit{p}}_{\mathit{n}}\mathit{-}{\mathit{g}}_{\mathit{n}}{\mathit{G}}_{\mathit{n}}^{\mathit{-}\mathbf{1}}$$E50

This process has to be repeated until the convergence criterion is reached. Thus, the displacement direction d_n is equal to -G_n⁻¹g_n and the optimal search step λ_n is always equal to 1. This method is called the Newton method.

In Quasi-Newton methods, the matrix G_n⁻¹ is replaced by symmetric positive H_n, which is updated from iteration to iteration. The following steps are included in the n iteration:

Initiate search direction:

$${{d}}_n=-{{H}}_n{{g}}_n$$E51

Definition of the next search point:

$${{p}}_{n+1}={{p}}_n+{\lambda }_n{{d}}_n$$E52

This may be done by any line search method such as blanket method, golden section search, Fibonacci section search, quadratic interpolation method.

Replacement of matrix H_n by H_n+1.

Initial Hessian matrix H₁ can be any symmetric positive definite and the simplest one is a unit matrix I. Matrices H_n+1 have been proposed by different authors. Davidson, Fletcher and Powell have proposed the following [26]:

$$H_{n+1}=H_n+\frac{\Delta {\mathit{p}}_n\Delta {\mathit{p}}_n^T}{\Delta {\mathit{p}}_n^T\Delta {\mathit{g}}_n}-\frac{{\mathit{H}}_n\Delta {\mathit{g}}_n\Delta {\mathit{g}}_n^T{\mathit{H}}_n}{\Delta {\mathit{p}}_n^T{\mathit{H}}_n\Delta {\mathit{g}}_n}$$E53

where: Δp_n = p_n+1-p_n, Δg_n = g_n+1-g_n, and superscript T indicates transposed matrix.

The parameters estimation finishes when either of the following criteria is observed:

$$\mid {\mathbf{p}}_n-{\mathbf{p}}_{n-1}\mid <{\epsilon }_1\&\mid E\left({\mathbf{p}}_n\right)\mid <{\epsilon }_2$$E54

$$\mid {\mathbf{p}}_n-{\mathbf{p}}_{n-1}\mid <{\epsilon }_1\&\mid E\left({\mathbf{p}}_n\right)-E\left({\mathbf{p}}_{n-1}\right)\mid <{\epsilon }_3$$E55

where ε₁, ε₂ and ε₃ are given small arbitrary positive values.

The block scheme of the parameter estimation process is given in Figure 8.

Figure 8.

Block-scheme of inverse analysis by Quasi-Newton method.

3.1.2 Determination of Freundlich’s adsorption isotherm parameters

The soil samples were taken from borehole BH5 in April 2016, which is 15 years from the operation of landfill cell No. 5 (Figure 3). The hydrodynamic dispersion coefficient was determined based on the above-described values of the coefficient of molecular diffusion, the soil porosity and the formation factor in paragraph 2.4 which were used as the input parameters. The element size and time step need to be not greater than 0.63 m and 422 days, respectively. Element size of 0.01 m and a time step of 1 day were used in this modelling for having sufficient data points along with a short distance of the concentration breakthrough curve.

Freundlich’s adsorption isotherm parameters are calibrated with the obtained HMs’ contents in soil taken from BH5. The trial and error method of calibration is used. Table 2 summarised the calibration results, i.e., the values of Freundlich’s adsorption isotherm parameters and mean error between the analysed and model HMs’ contents in the soil. Figure 9 illustrates the calibrated model HMs’ concentrations versus the analysed HMs’ concentrations. In general, the calibration models have a good fitting with a relative error of less than 7%, except the zinc.

Metal	KF	1/η	Mean error (mg/g)	Relative mean error (%)
Cr	0.264	0.260	0.0051	2.76
Cu	0.131	0.450	0.0033	5.03
Pb	0.073	0.850	0.0069	6.81
Zn	0.144	0.279	0.0312	19.41

Table 2.

The calibrated Freundlich’s adsorption isotherm parameters.

Figure 9.

Analysed and modelled results with the calibrated Freundlich’s adsorption isotherm parameters. (a) Lead in soil. (b) Chromium in soil water. (c) Copper in soil. (d) Zinc in soil water.

3.2 The FE model results

As the analysis results presented in Figure 5, four heavy metals Cr, Cu, Pb and Zn expose high concentrations on the surface 1–2 m of the soil layer. Modelling the transport of those four HMs was carried out. The breakthrough curves of concentrations of the four HMs in soil and pore water in MD1 are presented in Figure 10, where the allowable limits [27, 28] are also indicated. Thanks to the HMs’ adsorbability of the soil, only the upper layer of the soil horizon would be contaminated with HMs at levels higher than the allowable limits for agricultural land. For a period of 30 years, the soil would be contaminated in the upper 1 m, 2 m and 3 m by Cr, Zn and Pb, respectively (Figure 10: a1, c1 and d1). The concentrations of Cr, Pb and Zn in the soil pore water are higher than allowable limits in the upper 1.5 m, 6 m (i.e., the whole soil layer) and 2.2 m, respectively (Figure 10: a2, c2 and d2).

Figure 10.

Heavy metal concentrations prediction by FEM for 30 years - MD1. (a1) Chromium in soil. (a2) Chromium in pore water. (b1) Copper in soil. (b2) Copper in pore water. (c1) Lead in soil. (c2) Lead in pore water. (d1) Zinc in soil. (d2) Zinc in pore water.

Since the soil layer is under the landfill cells and leachate pond, only HMs in the soil pore water in MD2 are described here. MD2 with a very short length (1.5 m) presents a more problematic contamination situation. The MD2’s results are described here. Since the 27th year from the beginning of the landfill operation, the pore water with a concentration of Cr greater than the allowable limit begins to discharge into the upper Holocene aquifer (Figure 11a). The situation is more severe regarding Pb: the pore water with Pb concentration greater than the allowable limit begins to discharge into the upper Holocene aquifer from the 9th year (Figure 11b). The Arsenic concentration greater than the allowable limit in pore exists only in the upper 0.4 m after 30 years of the landfill operation (Figure 11c).

Figure 11.

Heavy metal concentrations in pore water prediction by FEM for 30 years: MD2. (a) Chromium in pore water - MD2. (b) Lead in pore water - MD2. (c) Arsenic in pore water - MD2.

4. Conclusions

FE modelling of advection-dispersion transport of HMs by GW movement along with Freundlich isotherm adsorption parameters which continuously change with space in the model domain and with time is sophisticated, but is capable of accurately evaluating the HMs’ concentrations in soil skeleton and pore water. The chapter describes the background of the existing isotherm adsorption theory. The chapter has provided a detailed mathematical formulations of the FEM in solving the advection-dispersion contaminant transport in soil water. It also demonstrates that the Freundlich isotherm adsorption parameters are essential to soil input parameters for modelling of HMs’ transport to access the soil skeleton and soil pore water contamination by HMs. In designing the experiments for the determination of the Freundlich isotherm adsorption parameters, the range of the HMs’ concentrations in water is suggested to be corresponding to the actual HMs’ concentrations under study. Besides, the background of the existing isotherm adsorption theory, the adaptation of the Freundlich isotherm adsorption in the soil skeleton and soil pore water contamination by HMs has been introduced.

The methodology has been applied to a case study of Kieu Ky waste landfill in Hanoi, Vietnam. The transport of HMs in soil water is determined not only by hydrodynamic dispersion but also largely by the adsorption of the metals by the soil. With the use of the collected interpreted isotherm adsorption parameters, the magnitudes of soil and soil water contamination by HMs from the waste leachate are very much different from each other due to both the HMs’ concentrations in leachate and the soil isotherm adsorption parameters. Unlike the pollutant transport in aquifers with coarse grain size particles like sand and gravel without clay materials with nearly-zero adsorption, the transport of pollutants in silty soils essentially requires adsorption parameters to have reliable modelling results.

The application modelling results show that HMs Cr, Pb and Zn present soil, soil pore water and groundwater contamination vulnerability, specifically as follows.

• Soil contamination with Cr, Pb and Zn by the direct spreading of the metals with dust and leachate from the waste landfills. For a period of 30 years, the soil would be contaminated in the upper 1 m, 2 m and 3 m by Cr, Zn and Pb, respectively

• The concentrations of Cr, Pb and Zn in the pore water in the silt layer are higher than allowable limits in the upper 1.5 m, 6 m (i.e., the whole soil layer) and 2.2 m, respectively.

• Since the 9th and 27th year from the beginning of the landfill operation, the pore water with Pb and Cr concentrations greater than the allowable limits begins to discharge into the upper Holocene aquifer, respectively.

• From a quarter of a century from the landfill operation start, Cr and Zn in the soil water would reach the Upper Holocene aquifer to pollute the aquifer. The waste leachate would cause the Upper Holocene aquifer polluted with Cr.

Conflict of interest

The author declares no conflict of interest.