**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]:

*Soil-Skeleton and Soil-Water Heavy Metal Contamination by Finite Element Modelling… DOI: http://dx.doi.org/10.5772/intechopen.101828*

$$D\_x \frac{\partial^2 C}{\partial \mathbf{x}^2} + D\_y \frac{\partial^2 C}{\partial \mathbf{y}^2} - \nu\_x \frac{\partial C}{\partial \mathbf{x}} - \nu\_y \frac{\partial C}{\partial \mathbf{y}} + Q = \frac{\mathbf{1} - \mathbf{n}}{n} \rho\_s \frac{\partial q\_\varepsilon}{\partial t} + \frac{\partial C}{\partial t} \tag{1}$$

where *C* is the contaminant concentration (M/L3 , e.g., mg/L); *t* is the time (d); *Dx* and *Dy* are the hydrodynamic dispersion coefficient in *x* and *y* direction, respectively (m<sup>2</sup> /d); *x* and *y* are the distances (m), *vx* and *vy* are the seepage velocity in *x* and *y* direction, respectively (m/d); *Q* is the distributed source of contaminant (mg/d); *ρ<sup>s</sup>* is the solid particle density (note that *ρ<sup>d</sup>* = *ρs*(1*-n*) in which *ρ<sup>d</sup>* is the unit weight of the dry soil); *n* is the soil porosity and *qe* 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 = \mathbf{K}\_F \mathbf{C}^{1/\eta} \tag{2}$$

where *C* is the concentration in solution at equilibrium (mg/L); *KF* and 1/η are fitting constants [13] and *KF* is termed as the Freundlich coefficient (adsorption coefficient) (the unit for the Freundlich constant is mg<sup>1</sup>�*<sup>η</sup> <sup>η</sup> l* 1 *<sup>η</sup>*/kg) and 1/*η* is the adsorption intensity (dimensionless). The value of *KF* is obtained from the intercept and 1/*η* from the slope of the logarithmic plot of log *qe* versus log *C* of the equation:

$$
\log q\_{\varepsilon} = \log K\_F + \frac{1}{\eta} \log \mathcal{C} \tag{3}
$$

$$\frac{dq\_e}{d\mathcal{C}} = K\_F \frac{1}{\eta} \mathcal{C}^{\frac{1-\eta}{\eta}} = K\_d \tag{4}$$

where *Kd* is the distribution coefficient.

*KF* 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 *KF* = *Kd*. *∂qe*

From (4) the source term *<sup>ρ</sup><sup>s</sup> n <sup>∂</sup><sup>t</sup>* in (1) is:

$$\frac{1-n}{n}\rho\_s \frac{\partial q\_\varepsilon}{\partial t} = \frac{1-n}{n}\rho\_s K\_F \frac{1}{\eta} C^{\frac{1-\eta}{\eta}} \frac{\partial C}{\partial t} \tag{5}$$

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

$$D\_x \frac{\partial^2 C}{\partial \mathbf{x}^2} + D\_y \frac{\partial^2 C}{\partial \eta^2} + D\_z \frac{\partial^2 C}{\partial \mathbf{x}^2} - \nu\_x \frac{\partial C}{\partial \mathbf{x}} - \nu\_y \frac{\partial C}{\partial \eta} - \nu\_z \frac{\partial C}{\partial \mathbf{z}} = \left(\mathbf{1} + \frac{\mathbf{1} - \boldsymbol{n}}{\boldsymbol{n}} \rho\_i K\_F \frac{\mathbf{1}}{\eta} C^{\frac{1 - \boldsymbol{v}}{\boldsymbol{q}}}\right) \frac{\partial C}{\partial \mathbf{z}}\tag{6}$$

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

$$R = \mathbf{1} + \frac{\mathbf{1} - n}{n} \rho K\_F \frac{\mathbf{1}}{\eta} C^{\frac{1}{\eta} - 1} = \mathbf{1} + \frac{\mathbf{1} - n}{n} \rho K\_d \tag{7}$$

$$D\_{\mathbf{x}} \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} + D\_{\mathbf{y}} \frac{\partial^2 \mathbf{C}}{\partial \mathbf{y}^2} - v\_{\mathbf{x}} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} - v\_{\mathbf{y}} \frac{\partial \mathbf{C}}{\partial \mathbf{y}} + Q = R \frac{\partial \mathbf{C}}{\partial t} \tag{8}$$

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 advectiondispersion 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*0:

$$\mathbf{C} = \mathbf{C}\_{\theta}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \tag{9}$$

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:

$$\mathbf{C} = \mathbf{C}\_{\mathfrak{c}} \text{ on} \\ \Gamma\_{\mathfrak{c}} \tag{10}$$

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

$$\frac{\partial C}{\partial n} = f\_C \text{ on} \\ \Gamma\_{fC} \tag{11}$$

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

$$
v\_n \mathbf{C} - D \frac{\partial \mathbf{C}}{\partial n} = q\_c \text{ one} \\
\Gamma\_{qw} \tag{12}$$

where *vn* 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 *Kd,* 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

*Soil-Skeleton and Soil-Water Heavy Metal Contamination by Finite Element Modelling… DOI: http://dx.doi.org/10.5772/intechopen.101828*

about 10 at higher contaminant concentrations, it could be safer and more accurate to make use of the unified sorption variable *KF* to calculate *R*. Since *KF* changes constantly with *C* and using a constantly changing *KF* would be complicated, one solution is to select a small number of discrete values of *KF* that can be used to approximate and slightly underestimate the actual values of *KF* and each of these values can be used to calculate *R* over the range of contaminant concentrations that they are applicable.
