**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 dx dy = 0 \qquad \ell = 1, 2, \dots, M \tag{13}$$

where *W<sup>ℓ</sup>* is the weighting function. Using the Green lemma:

$$\int\_{\varDelta} \left( D\_x \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} + D\_y \frac{\partial^2 \mathbf{C}}{\partial \mathbf{y}^2} \right) W\_t dx dy = - \int\_{\varDelta} \left( D\_x \frac{\partial W\_\ell}{\partial \mathbf{x}} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + D\_y \frac{\partial W\_\ell}{\partial \mathbf{y}} \frac{\partial \mathbf{C}}{\partial \mathbf{y}} \right) dx dy + \int\_{\varGamma} \left( D\_x \frac{\partial \mathbf{C}}{\partial \mathbf{x}} n\_x + D\_y \frac{\partial \mathbf{C}}{\partial \mathbf{y}} n\_y \right) W\_t d\Gamma \tag{14}$$

The integral Ð <sup>Γ</sup> is present only for the elements having sides in boundary *Γqc* or *Γqυ<sup>c</sup>* which are generally denoted as *Γ<sup>q</sup>* ta có:

$$-\int\_{\Omega} \left( D\_x \frac{\partial W\_\ell}{\partial \mathbf{x}} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + D\_y \frac{\partial W\_\ell}{\partial \mathbf{y}} \frac{\partial \mathbf{C}}{\partial \mathbf{y}} \right) d\mathbf{x}d\mathbf{y} - \int\_{\Omega} \left( u\_x \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + u\_y \frac{\partial \mathbf{C}}{\partial \mathbf{y}} \right) W\_\ell d\mathbf{x}d\mathbf{y} = + \int\_{\Gamma} \left( D\_x \frac{\partial \mathbf{C}}{\partial \mathbf{x}} u\_x + D\_y \frac{\partial \mathbf{C}}{\partial \mathbf{y}} u\_y \right) W\_\ell d\Gamma = \mathbf{0} \tag{15}$$

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

$$\mathbf{C} \approx \hat{\mathbf{C}} = \sum\_{m=1}^{M} \mathbf{C}\_{m} \mathbf{N}\_{m} \tag{16}$$

where *Cm* is the approximation of the contaminant concentration at node *m* and *Nm* is the shape functions.

Equation (15) becomes:

$$\begin{aligned} -\int\_{\varOmega} \left( D\_{\mathbf{x}} \frac{\partial \mathcal{W}\_{\ell}}{\partial \mathbf{x}} \frac{\partial \mathcal{N}\_{\mathbf{m}}}{\partial \mathbf{x}} \mathcal{C}\_{\mathbf{m}} + D\_{\mathbf{y}} \frac{\partial \mathcal{W}\_{\ell}}{\partial \mathbf{y}} \frac{\partial \mathcal{N}\_{\mathbf{m}}}{\partial \mathbf{y}} \mathcal{C}\_{\mathbf{m}} \right) d\mathbf{x} d\mathbf{y} - \int\_{\varOmega} \left( v\_{\mathbf{x}} \frac{\partial \mathcal{N}\_{\mathbf{m}}}{\partial \mathbf{x}} \mathcal{W}\_{\ell} \mathcal{C}\_{\mathbf{m}} + v\_{\mathbf{y}} \frac{\partial \mathcal{N}\_{\mathbf{m}}}{\partial \mathbf{y}} \mathcal{W}\_{\ell} \mathcal{C}\_{\mathbf{m}} \right) d\mathbf{x} d\mathbf{y} \\ = \int\_{\varGamma\_{\mathbf{q}}} \left( \bar{q}\_{\mathbf{x}} \mathcal{W}\_{\ell} n\_{\mathbf{x}} + \bar{q}\_{\mathbf{x}} \mathcal{W} n\_{\mathbf{y}} \right) d\Gamma = \mathbf{0} \end{aligned} \tag{17}$$

$$K = -\int\_{\varDelta} \left( D\_x \frac{\partial W\_\ell}{\partial \mathbf{x}} \frac{\partial \mathbf{N}\_m}{\partial \mathbf{x}} + D\_y \frac{\partial W\_\ell}{\partial \mathbf{y}} \frac{\partial \mathbf{N}\_m}{\partial \mathbf{y}} \right) d\mathbf{x}d\mathbf{y} - \int\_{\varDelta} \left( v\_x \frac{\partial \mathbf{N}\_m}{\partial \mathbf{x}} \mathbf{W}\_\ell + v\_y \frac{\partial \mathbf{N}\_m}{\partial \mathbf{y}} \mathbf{W}\_\ell \right) d\mathbf{x}d\mathbf{y} \tag{18}$$

*Environmental Impact and Remediation of Heavy Metals*

$$F = -\int\_{\Gamma\_q} (\bar{q}\_c \mathcal{W}\_\ell n\_\mathfrak{x} + \bar{q}\_c \mathcal{W}\_\ell n\_\mathfrak{y}) d\Gamma \tag{19}$$

$$\mathbf{K}\mathbf{C} = \mathbf{F} \tag{20}$$

The shape functions Nm and weighting functions *W<sup>ℓ</sup>* can be linear or higher order functions. For unsteady problems, i.e., *<sup>R</sup>∂C/∂<sup>t</sup>* 6¼ 0 we have:

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

The square matrix E is:

$$E = \int\_{\varOmega} RW\_{\ell}N\_{m}d\mathbf{x}dy\tag{22}$$

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

$$\left(\theta \mathcal{K} + \frac{E}{\Delta t}\right) \mathcal{C}^{t + \Delta t} = \left(-(\mathbb{1} - \theta) \mathcal{K} + \frac{E}{\Delta t}\right) \mathcal{C}^{t} + \frac{\mathbb{1}}{2} \left(F^{t} + F^{t + \Delta t}\right) \tag{23}$$

The typical schemes are:

i. Forward difference (*θ* **=** 0):

$$R\frac{[E]}{\Delta t}\left\{\mathbf{C}^{t+\Delta t}\right\} + \left([K] - R\frac{[E]}{\Delta t}\right)\{\mathbf{C}^{t}\} = \{\mathbf{F}^{t}\}\tag{24}$$

ii. Backward difference (*θ* **=** 1):

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

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

$$
\left(\frac{\mathfrak{1}}{2}[K] + \frac{[E]}{\Delta t}\right)\{\mathcal{C}^{t+\Delta t}\} + \left(\frac{\mathfrak{1}}{2}[K] - \frac{[E]}{\Delta t}\right)\{\mathcal{C}^{t}\} = \frac{\mathfrak{1}}{2}\left(\{F^{t}\} + \{F^{t+\Delta t}\}\right) \tag{26}
$$

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 *hx* and *hy*.

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):

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

The four-row terms of column *i* are:

*Ke* ð Þ *<sup>i</sup>*,*j*,*k*,*<sup>l</sup> <sup>i</sup>* ¼ Ð *he y 0* Ð *he x <sup>0</sup> Wivx ∂Ni ∂x* þ *Wivy ∂Ni ∂y* þ *Dx ∂Wi ∂x ∂Ni ∂x* þ *Dy ∂Wi ∂y ∂Ni ∂y* � �*dxdy* Ð *he y 0* Ð *he x <sup>0</sup> Wjvx ∂Ni ∂x* þ *Wjvy ∂Ni ∂y* þ *Dx ∂Wj ∂x ∂Ni ∂x* þ *Dy ∂Wj ∂y ∂Ni ∂y* � �*dxdy* Ð *he y 0* Ð *he x <sup>0</sup> Wkvx ∂Ni ∂x* þ *Wkvy ∂Ni ∂y* þ *Dx ∂Wk ∂x ∂Ni ∂x* þ *Dy ∂Wk ∂y ∂Ni ∂y* � �*dxdy* Ð *he y 0* Ð *he x <sup>0</sup> Wlvx ∂Ni ∂x* þ *Wlvy ∂Ni ∂y* þ *Dx ∂Wl ∂x ∂Ni ∂x* þ *Dy ∂Wl ∂y ∂Ni ∂y* � �*dxdy* 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 (27)

The four-row terms of column *j* are:

*Ke* ð Þ *<sup>i</sup>*,*j*,*k*,*<sup>l</sup> <sup>j</sup>* ¼ ð*he y 0* ð*he x 0 Wivx ∂Nj ∂x* þ *Wivy ∂Nj ∂y* þ *Dx ∂Wi ∂x ∂Nj ∂x* þ *Dy ∂Wi ∂y ∂Nj ∂y* � �*dxdy* ð*he y 0* ð*he x 0 Wjvx ∂Nj ∂x* þ *Wjvy ∂Nj ∂y* þ *Dx ∂Wj ∂x ∂Nj ∂x* þ *Dy ∂Wj ∂y ∂Nj ∂y* � �*dxdy* ð*he y 0* ð*he x 0 Wkvx ∂Nj ∂x* þ *Wkvy ∂Nj ∂y* þ *Dx ∂Wk ∂x ∂Nj ∂x* þ *Dy ∂Wk ∂y ∂Nj ∂y* � �*dxdy* ð*he y 0* ð*he x 0 Wlvx ∂Nj ∂x* þ *Wlvy ∂Nj ∂y* þ *Dx ∂Wl ∂x ∂Nj ∂x* þ *Dy ∂Wl ∂y ∂Nj ∂y* � �*dxdy* (28)

**89**

The four-row terms of column *k* are:

*Ke* ð Þ *<sup>i</sup>*,*j*,*k*,*<sup>l</sup> <sup>k</sup>* ¼ ð*he y* 0 ð*he x* 0 *Wivx ∂Nk ∂x* þ *Wivy ∂Nk ∂y* þ *Dx ∂Wi ∂x ∂Nk ∂x* þ *Dy ∂Wi ∂y ∂Nk ∂y* � �*dxdy* ð*he y* 0 ð*he x* 0 *W jvx ∂Nk ∂x* þ *W jvy ∂Nk ∂y* þ *Dx ∂W <sup>j</sup> ∂x ∂Nk ∂x* þ *Dy ∂W <sup>j</sup> ∂y ∂Nk ∂y* � �*dxdy* ð*he y* 0 ð*he x* 0 *Wkvx ∂Nk ∂x* þ *Wkvy ∂Nk ∂y* þ *Dx ∂Wk ∂x ∂Nk ∂x* þ *Dy ∂Wk ∂y ∂Nk ∂y* � �*dxdy* ð*he y* 0 ð*he x* 0 *Wlvx ∂Nk ∂x* þ *Wlvy ∂Nk ∂y* þ *Dx ∂Wl ∂x ∂Nk ∂x* þ *Dy ∂Wl ∂y ∂Nk ∂y* � �*dxdy* (29)

The four-row terms of column *l* are:

*Ke* ð Þ *<sup>i</sup>*,*j*,*k*,*<sup>l</sup> <sup>l</sup>* ¼ ð*he y 0* ð*he x 0 Wivx ∂Nl ∂x* þ *Wivy ∂Nl ∂y* þ *Dx ∂Wi ∂x ∂Nl ∂x* þ *Dy ∂Wi ∂y ∂Nl ∂y* � �*dxdy* ð*he y 0* ð*he x 0 Wjvx ∂Nl ∂x* þ *Wjvy ∂Nl ∂y* þ *Dx ∂Wj ∂x ∂Nl ∂x* þ *Dy ∂Wj ∂y ∂Nl ∂y* � �*dxdy* ð*he y 0* ð*he x 0 Wkvx ∂Nl ∂x* þ *Wkvy ∂Nl ∂y* þ *Dx ∂Wk ∂x ∂Nl ∂x* þ *Dy ∂Wk ∂y ∂Nl ∂y* � �*dxdy* ð*he y 0* ð*he x 0 Wlvx ∂Nl ∂x* þ *Wlvy ∂Nl ∂y* þ *Dx ∂Wl ∂x ∂Nl ∂x* þ *Dy ∂Wl ∂y ∂Nl ∂y* � �*dxdy* 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (30)

The Galerkin FEM with linear shape functions results in: The four-row terms of column *I,* the Eq. (27) become:

$$K\_{(j,k,l)}^{\epsilon} = \begin{bmatrix} \left[D\_x^{\epsilon} \left(\frac{(k\_j'-\gamma)}{k\_j k\_j'}\right)^2 + D\_y^{\epsilon} \left(\frac{(k\_i'-\gamma)}{k\_j' k\_j'}\right)^2\right] + \left[\frac{(k\_x'-\chi)\left(k\_j'-\chi\right)}{k\_j' k\_j'} - \left(\frac{(k\_j'-\chi)}{k\_j' k\_j'} - \nu\_{y\_j}' \frac{(k\_x'-\chi)}{k\_j' k\_j'}\right)\right] \\\\ -D\_x^{\epsilon} \frac{(k\_j'-\chi)\left(k\_j'-\chi\right)}{k\_j' k\_j'} - D\_y^{\epsilon} \frac{x}{k\_j' k\_j'} - D\_y^{\epsilon} \frac{x}{k\_j' k\_j'} - \left[\frac{x\left(k\_j'-\chi\right)}{k\_j' k\_j'} - \left(\frac{(k\_j'-\chi)}{k\_j' k\_j'} - \nu\_{y\_j}' \frac{(k\_x'-x)}{k\_j' k\_j'}\right)\right] \\\\ \left[-D\_x^{\epsilon} \frac{y}{k\_j' k\_j'} - \nu\_{y\_j}' \frac{(k\_x'-\chi)}{k\_j' k\_j'} - D\_y^{\epsilon} \frac{x}{k\_j' k\_j'} \frac{(k\_x'-\chi)}{k\_j' k\_j'}\right] + \left[\frac{xy}{k\_j' k\_j'} \left(-\nu\_x' \frac{(k\_j'-\chi)}{k\_j' k\_j'} - \nu\_{y\_j}' \frac{(k\_x'-x)}{k\_j' k\_j'}\right)\right] \\\\ \left[D\_x^{\epsilon} \frac{y}{k\_j' k\_j'} - D\_y^{\epsilon} \frac{(k\_x'-x)}{k\_j' k\_j'} - D\_y^{\epsilon} \frac{(k\_x'-x)}{k\_j' k\_j'} \frac{(k\_x'-x)}{k\_j'$$

$$\begin{split} F^{\epsilon}\_{i} &= \int\_{0}^{h^{\epsilon}\_{x}} \int\_{0}^{h^{\epsilon}\_{y}} \frac{\mathcal{Q}^{\epsilon}}{h^{\epsilon}\_{x}h^{\epsilon}\_{y}} \left(h^{\epsilon}\_{x} - x\right) \left(h^{\epsilon}\_{y} - y\right) dx dy - \underbrace{\int\_{0}^{h^{\epsilon}\_{x}} \frac{q^{\epsilon}}{h^{\epsilon}\_{x}} \left(h^{\epsilon}\_{x} - x\right) dx}\_{} - \underbrace{\int\_{0}^{h^{\epsilon}\_{y}} \frac{q^{\epsilon}}{h^{\epsilon}\_{y}} \left(h^{\epsilon}\_{y} - y\right) dy}\_{} \\ &= \frac{1}{4} \mathcal{Q}^{\epsilon} h^{\epsilon}\_{x} h^{\epsilon}\_{y} - \underbrace{\frac{1}{2} q^{\epsilon} h^{\epsilon}\_{x} - \frac{1}{2} q^{\epsilon} h^{\epsilon}\_{y}}\_{} \end{split} \tag{35}$$

The matrix E is:

$$\begin{aligned} \text{node}: \qquad &i \qquad j \qquad k \qquad l \\ \text{E}\_{i}^{1} = R\_{\epsilon} \begin{bmatrix} \int\_{\Omega} W\_{i} N\_{i} d\mathbf{x} dy & \int\_{\Omega} W\_{i} N\_{j} d\mathbf{x} dy & \int\_{\Omega} W\_{i} N\_{k} d\mathbf{x} dy & \int\_{\Omega} W\_{i} N\_{\ell} d\mathbf{x} dy \\ \int\_{\Omega} W\_{j} N\_{i} d\mathbf{x} dy & \int\_{\Omega} W\_{j} N\_{j} d\mathbf{x} dy & \int\_{\Omega} W\_{j} N\_{k} d\mathbf{x} dy & \int\_{\Omega} W\_{j} N\_{\ell} d\mathbf{x} dy \\ \int\_{\Omega} W\_{k} N\_{i} d\mathbf{x} dy & \int\_{\Omega} W\_{k} N\_{j} d\mathbf{x} dy & \int\_{\Omega} W\_{j} N\_{k} d\mathbf{x} dy & \int\_{\Omega} W\_{j} N\_{\ell} d\mathbf{x} dy \\ \int\_{\Omega} W\_{\ell} N\_{i} d\mathbf{x} dy & \int\_{\Omega} W\_{\ell} N\_{j} d\mathbf{x} dy & \int\_{\Omega} W\_{\ell} N\_{k} d\mathbf{x} dy & \int\_{\Omega} W\_{\ell} N\_{\ell} N\_{\ell} d\mathbf{x} dy \end{bmatrix} \begin{bmatrix} C\_{i} \\ C\_{j} \\ C\_{k} \\ C\_{\ell} \end{bmatrix} \end{aligned} \tag{36}$$

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

node : *i j <sup>E</sup>*<sup>1</sup> <sup>¼</sup> *Re* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � � *he <sup>y</sup>* � *y* � � *he xhe y he <sup>x</sup>* � *<sup>x</sup>* � � *<sup>h</sup><sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � � *he <sup>y</sup>* � *y* � � *he xhe y x h<sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω x h<sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y he <sup>x</sup>* � *<sup>x</sup>* � � *<sup>h</sup><sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω x h<sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y x he <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω xy he xhe y he <sup>x</sup>* � *<sup>x</sup>* � � *<sup>h</sup><sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω xy he xhe y x he <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y he <sup>x</sup>* � *<sup>x</sup>* � � *<sup>h</sup><sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y dxdy* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y x h<sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y dxdy* 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 (37)

node : *k l*

Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � � *he <sup>y</sup>* � *y* � � *he xhe y xy he xhe y dxdy* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � � *<sup>h</sup><sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y dxdy* Ð *Ω x h<sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y xy he xhe y dxdy* Ð *Ω x h<sup>e</sup> <sup>y</sup>* � *y* � � *he xhe y he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y dxdy* Ð *Ω xy he xhe y xy he xhe y dxdy* Ð *Ω xy he xhe y he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y dxdy* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y xy he xhe y dxdy* Ð *Ω he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y he <sup>x</sup>* � *<sup>x</sup>* � �*<sup>y</sup> he xhe y dxdy* 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 *Ci Cj Ck C<sup>ℓ</sup>* 2 6 6 6 6 6 4 3 7 7 7 7 7 5 (38)

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

$$Pe = \frac{v\_{\text{x}}h\_{\text{x}}^{\epsilon}}{D\_{\text{x}}} \le 2; Pe = \frac{v\_{\text{y}}h\_{\text{y}}^{\epsilon}}{D\_{\text{y}}} \le 2; Cr = \frac{v\_{\text{x}}\Delta t}{h\_{\text{x}}^{\epsilon}} \le 1; Cr = \frac{v\_{\text{y}}\Delta t}{h\_{\text{y}}^{\epsilon}} \le 1\tag{39}$$

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


#### **Table 1.**

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

#### **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<sup>+</sup> in the solution. In solution with pH greater than 7, all the ions H<sup>+</sup> 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 *KF* 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 *KF* and 1/*η* of the fitting curves (**Figure 2**).

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 *KF* and 1/*η* of the soils: *KF* = 0.348 (mg/g) and 1/*η* = 0.235 (σ = 0.059) for Zn, *KF* = 0.462 and 1/*η* = 0.320 for Cu (σ = 0.054), *KF* = 0.279 and 1/*η* = 0.437 (σ = 0.059) for Ni (with the *qe* in mg/g and *C* unit in mg/L).

#### **Figure 2.**

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

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: *KF* = 1.560 (mg/g) and 1/*η* = 0.327 for Cu, *KF* = 0.363 and 1/*η* = 0.441 for Ni, *KF* = 1.363 and 1/ *η* = 0.351 for Pb, *KF* = 0.463 and 1/*η* = 0.426 for Zn, *KF* = 0.540 and 1/*η* = 0.293 for Cd (with the *qe* 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: *KF* = 2.540 (mg/g) and 1/*η* = 0.91 for Cu, *KF* = 0.702 and 1/*η* = 0.510 for Ni, *KF* = 0.998 and 1/*η* = 0.440 for Pb, *KF* = 1.016 and 1/*η* = 0.440 for Zn (with the *qe* in mg/g and *C* unit in mg/L).

#### **2.4 Dispersion parameters**

The coefficient of longitudinal hydrodynamic dispersion *DL* 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

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

*D*<sup>0</sup> and the coefficient of molecular diffusion in a porous medium *D\*d*, i.e., *DL = D*<sup>0</sup> *+ D\*d* [9].

The coefficient of mechanical dispersion *D*<sup>0</sup> depends upon the microstructure of the soil, the seepage velocity and the molecular dispersion in water as follows by Jacob and Arnold [9]:

$$^jD'\_{\vec{\eta}} = a\_L v \mathbf{f}(\text{Pe}, \xi); \text{Pe} = \frac{Lv}{D\_d} \tag{40}$$

in which: *v* is the seepage velocity (m/d); *Pe* is the Peclet number; *L* is the characteristic length of the pores (m); *Dd* is the molecular dispersion in water; *ξ* is the ratio between the pores'size and the characteristic length through the pores; *f*(Pe*,ξ*)=Pe/ (Pe+2 + 4*ξ 2* ) 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*,ξ*) ffi 1; *aL* is the longitudinal dispersivity.

For a one-directional groundwater flow, the coefficient of mechanical dispersion *D*<sup>0</sup> is the multiplication of longitudinal dispersivity (*aL*) and seepage velocity [9]. The longitudinal dispersivity is of the order of the average soil particle [9], e.g., particle size *d*50. The coefficient of molecular diffusion in a porous medium *D\*d* is as follows [23]:

$$D\_d^\* = \frac{D\_d}{nF\_R} \tag{41}$$

in which: *F*<sup>R</sup> is the formation factor which is specified by the geophysicists as the ratio between soil resistivity and water resistivity.

The formation factor *F*<sup>R</sup> ranges from 0.1 (for clay) to 0.7 (for sand) [23], and always less than 1.

The coefficient of molecular diffusion in water *Dd* is:

$$D\_d = \frac{RT}{N} \frac{1}{6\pi r\mu} \tag{42}$$

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 *Dd* 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*<sup>0</sup> ) is negligible, i.e., *DL* ≈ *D\* <sup>d</sup>*. 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*.
