**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.

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)

#### **Figure 3.**

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

**Figure 4.** *The soil profile of Kieu Ky landfill site.*

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

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

**Figure 5.** *Typical model domains in the study site.*

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.

#### **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 secondorder 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.

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

#### *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.

The special and temporal concentrations are monitored, for which the observed (*Cobs*) and model (*Cmod*) 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.

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

*L*

**Figure 6.** *A tracer column test scheme.*

**Figure 7.** *Plots of observed vs. model tracer concentration.*

where *E*(**p**) - objective function; *L*- number of observed variables; *Cl* mod- model estimated values of the concentration; *Cl* 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:

$$\dim E(p), p \in p\_{ct} \tag{44}$$

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 *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{\partial E}{\partial p\_m} \right|\_{\hat{p}} = 0; (m = 1, 2, \dots, M) \tag{45}$$

where *M* is the number of parameter variables. Hessian matrix *G =* ∇<sup>2</sup> *E*(**p**):

$$G = \begin{vmatrix} \frac{\partial^2 E}{\partial p\_1^2} & \frac{\partial^2 E}{\partial p\_1 \partial p\_2} & \cdots & \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} & \cdots & \frac{\partial^2 E}{\partial p\_2 \partial p\_M} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \frac{\partial^2 E}{\partial p\_M \partial p\_1} & \frac{\partial^2 E}{\partial p\_M \partial p\_2} & \cdots & \frac{\partial^2 E}{\partial p\_M^2} \end{vmatrix} \tag{46}$$

is a positive definite matrix at *p*^.

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

Selection of the initial values of parameters **p**0.

Determination of the search sequence: **p**0, **p**1, **p**2, ..., **p***<sup>n</sup>* ... in such a way that *E* (**p***<sup>n</sup>* + 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:

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

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.

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

Suppose there is a set of initial values of parameters **p**0, it is required to find out the search sequence **p**0, **p**1, **p**2, ..., **p***<sup>n</sup>* so that *E*(**p***<sup>n</sup>* + 1) < *E* (**p***n*) for all *n*. Gradient vector **g**(**p***n+1*) at vicinity of **p***<sup>n</sup>* may be determined as follows:

$$\text{g}\,(p\_{n+1}) \approx \text{g}\_n + \text{G}\_n\\\Delta p; \Delta p = p - p\_n; \text{g}\_n = \text{g}\,(p\_n); \text{G}\_n = \text{G}\,(p\_n) \tag{48}$$

The necessary condition of extreme existence is **g**(**p***<sup>n</sup>* + 1) ≈ 0. This can be done if **p***n*+1 = **p***<sup>n</sup>* + Δ **p***n*, where Δ**p***<sup>n</sup>* are the solution of the following equation:

$$\mathbf{g}\_n + \mathbf{G}\_n \Delta p = \mathbf{0} \tag{49}$$

$$
\Delta p\_n = \Delta p = -\mathbf{g}\_n \mathbf{G}\_n^{-1} \Rightarrow p\_{n+1} = p\_n - \mathbf{g}\_n \mathbf{G}\_n^{-1} \tag{50}
$$

This process has to be repeated until the convergence criterion is reached. Thus, the displacement direction **d***<sup>n</sup>* is equal to -**G***<sup>n</sup>* �1 **g***<sup>n</sup>* and the optimal search step *λ<sup>n</sup>* is always equal to 1. This method is called the Newton method.

In Quasi-Newton methods, the matrix **G***<sup>n</sup>* �<sup>1</sup> 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 \mathbf{g}\_n \tag{51}$$

Definition of the next search point:

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

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***<sup>n</sup>* by **H***n*+1.

Initial Hessian matrix **H**<sup>1</sup> 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 p\_n \Delta p\_n^T}{\Delta p\_n^T \Delta \mathbf{g}\_n} - \frac{H\_n \Delta \mathbf{g}\_n \Delta \mathbf{g}\_n^T H\_n}{\Delta p\_n^T H\_n \Delta \mathbf{g}\_n} \tag{53}$$

where: Δ**p***<sup>n</sup>* = **p***n*+1-**p***n*, Δ**g***<sup>n</sup>* = **g***n*+1-**g***n*, and superscript *T* indicates transposed matrix. The parameters estimation finishes when either of the following criteria is observed:

$$|\mathbf{p}\_n - \mathbf{p}\_{n-1}| < \varepsilon\_1 \& | \; E(\mathbf{p}\_n)| < \varepsilon\_2 \tag{54}$$

$$|\mathbf{p}\_n - \mathbf{p}\_{n-1}| < \varepsilon\_1 \& |E(\mathbf{p}\_n) - E(\mathbf{p}\_{n-1})| < \varepsilon\_3 \tag{55}$$

where *ε*1, *ε*<sup>2</sup> and *ε*<sup>3</sup> are given small arbitrary positive values. The block scheme of the parameter estimation process is given in **Figure 8**.

#### *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

**Figure 8.** *Block-scheme of inverse analysis by Quasi-Newton method.*

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

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.

#### **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,


#### **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.*

#### **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.*

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

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

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

#### **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.*

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**).
