Numerical and experimental simulation of pollutants migration in porous media M. Zairi 7 M. J. Rouis

Abstract Quantitative prediction of pollutant migration is now of great interest in the field of waste management and storage. Numerical modelling of contaminant transport through porous media is an efficient tool to make this prediction possible. Site selection and design of protective liners for landfills and storage facilities are largely related to the behaviour of contaminants in such structures. This behaviour may be well defined by numerical simulations. This paper proposes a twodimensional finite element model for pollutant migration. The implementation of the advectiondispersion equation in the numerical model and validation tests by comparison with analytical and semianalytical solutions are presented. Experimental laboratory tests on cadmium and chloride migration through liner samples are also discussed and compared with the results from the numerical model. Résumé La prévision quantitative de la migration des polluants présente un grand intérêt dans le domaine de la gestion et du stockage des déchets. La modélisation numérique du transport des polluants à travers un milieu poreux est un outil efficace rendant possible cette prévision. La sélection de sites et la conception de barrières de protection pour des installations de stockage dépendent fortement du

comportement des polluants dans de telles structures. Ce comportement peut être bien défini par des simulations numériques. L’article présente un modèle 2D en éléments finis pour la migration de polluants. L’établissement des équations d’advection-dispersion dans le modèle numérique et les tests de validation par comparaison avec des solutions analytiques et semi-analytiques sont présentés. Des expériences en laboratoire sur la migration d’ions cadmium et chlorure à travers des barrières de protection reconstituées sont également discutées et les résultats comparés à ceux de la modélisation numérique. Keywords Environmental protection 7 Groundwater contamination 7 Laboratory studies 7 Models Mots clés Protection de l’environnement 7 Pollution des eaux souterraines 7 Études en laboratoire 7 Modèles numériques

Introduction

Numerical modelling of pollutant migration in porous media has recently received a great deal of attention due to an increased interest in the preservation of the quality of the environment and particularly the protection of groundwater from various pollutants. Mathematical modelling of contaminant behaviour in soils is considered to be a powerful tool for a wide range of problems concerned with Received: 5 July 1999 7 Accepted: 15 April 2000 groundwater quality preservation or rehabilitation. M. Zairi (Y) 7 M. J. Rouis Ecole Nationale d’Ingénieurs de Sfax, BP W, 3038 Sfax, Tunisia Enormous efforts have been made in recent years to understand contaminant behaviour in porous media and e-mail: Moncef.Zairi6enis.rnu.tn to represent it by mathematical models. Ogata (1970) was Tel.: c216-4-274088 the first to give a non-dimensional analytical solution to Fax: c216-4-275595

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

231

M. Zairi 7 M. J. Rouis

the mass transport equation in porous media. Others models were then proposed, e.g. those of Guerghian et al. (1980) and Rowe and Booker (1985). Although there are now many numerical methods that can be used to solve the mass transport equation, the finite element model is particularly appropriate for porous media with complex geometry and flows. It also permits a good representation of porous media heterogeneity.

Numerical model

librium sorption-desorption model of retardation, requires particularly careful use in modelling. At high pH values, removal of heavy metals from solutions is generally greater than the exchange capacity of the soil, due to precipitation as metal carbonate or hydroxide. To overcome this, an apparent greater Kd value is used to include precipitation and fixation in modelling retardation of heavy metals. In the finite element transformation of Eq. (1), the porous media (or Eq. 1 domain) is discretized in a number (m) of elements (e), which constitute the finite element mesh or grid. The approximated solution for this equation may be written for each element in the form: n

Pollutant behaviour in the soil is subject to many processes. The water–soil system is dynamic and many chemical, biological and physical reactions occur which may affect pollutant behaviour. Contaminant transport can be reasonably approximated by one-dimensional solutions, and the two-dimensional analysis will be adequate for the vast majority of practical applications. In the case of a clay layer acting as a liner below a landfill, for example, the one-dimensional analysis may be adequate when the width of the landfill is large compared to the thickness of the clay liner and the direction of contaminant transport is predominantly vertical. However, if the dimensions of the landfill are comparable to the clay liner thickness, horizontal and vertical transport may occur and the twodimensional analysis is more appropriate, especially when concentration at points outside the landfill is of concern. The governing equation of two-dimensional pollutant migration in saturated porous media, also named the advection-dispersion equation, is given by: f 2C f 2C fC fC fC p D*x 2 c D*y 2 P V*x P V*y ft fx fy fx fy

(1)

where C is the concentration of the considered contaminant at a point of coordinates (x, y) at time t; D*pD/R, V*pV/R; R is the retardation factor, defined by Bear (1972, in Auriault and Lewandowska 1993) which is given by: Rp1crKd/n

(1a)

where Vx and Vy are pore water flow velocities in X and Y directions; Dx and Dy are hydrodynamic dispersion coefficients in X and Y directions; r is soil bulk density; Kd is distribution coefficient; and n is porosity. In the case of conservative species which have no interaction with porous medium particles, the retardation factor is equal to 1. For species adsorbed by soil particles, the value of this factor depends on the type of element studied and the soil adsorption capacity expressed by the distribution or partitioning coefficient. It is appreciated that this is a very simple description of the interaction between the contaminant and the porous medium particles, but whilst it is approximate, it has often been found to be an adequate approximation. Kd, the distribution coefficient which reflects the degree of retardation by reversible ion exchange, known as the equi232

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

C e (x, y)p A N ei Ci

(2)

ip1

with C e(x,y): approximate concentration of element (e) Nie: interpolation functions for each node of element (e) Ci: nodal concentration (unknown of the problem) n: number of nodes of the element (e). If C in Eq. (1) is substituted by the expression of Eq. (2), an approximate solution for Eq. (1) is obtained. The porous medium properties (e.g. hydrodynamic dispersion, flow velocities, porosity, distribution coefficient) are considered constant in each element but may vary from one element to another. The second derivatives of the pollutant-migration-governing equation are generally undefined for a variety of grid elements. They may be evaluated if written in the form of first derivatives, by using an integration by parts of the diffusive terms of the equation as given by Green’s theorem for integration in two dimensions (Zienkiewicz 1977). The integral on the boundary of the domain or contour integral resulting from the integration by parts is introduced in the form of a boundary condition. On the domain boundary, with an imposed ingressing mass flux, the global quantity QCi is substituted for the contour integral, where Q is the water flux and Ci its concentration with the compound i. If the previous developments are written for each element of the grid, the elementary equation may be written in the matrix form as: iC1 it C 1 c [M e] [K e] p0 Cn iCn it

3 4

34

(3)

with [K e]p[B] T[D][B]c[NN] T[V][B] [M e]p[N] T[N] where elements of matrices [B] and [B] T are composed of the nodal values of interpolation function first derivatives. The elements of matrices [N], [N] T and [NN] T are the nodal values of interpolation functions and those of matrices [D] and [V] are the elementary hydrodynamic dispersion coefficients and flow velocities respectively. The

Pollutant migration in porous media

elementary matrices [K e] and [M e] of Eq. (3) are assem2 (7) bled for all domain elements to give the following global Dtc p lP2a lmax system, in which m is the number of domain elements. where lmax is the greatest eigenvalue of the product of –1 iC (4) matrices [M] and [K] (Dhatt and Tozot 1981). In order to p0 [K] Cc[M] it solve Eq. (6), the Newton-Raphson method (as described by Dhatt and Tozot 1981) may solve linear or non-linear with non-stationary problems. This method was generally used m to solve non-linear problems. In the constructed algo[K]p A [K e] rithm, for each iteration Eq. (6) is solved by the Gauss ip1 elimination method. m The above developments were coded in a FORTRAN e [M]p A [M ] program named MIGPO (Zairi 1996), using eight nodes ip1 quadrilateral isoparametric elements for grid generation. The global system of Eq. (4) is a first-order differential More details about these elements, interpolation functions equation in time. and coordinate transformations are given in Zienkiewicz (1977). Definition of boundary conditions Two types of boundary conditions may exist on the Conditions for numerical solution stability domain limits. The first is a fixed concentration limit The discretization is the preliminary stage in transforming which is a Dirchelet-type boundary condition. The second a system of partial derivatives equations into an algebraic is a fixed concentration flux limit and a Neumann-type one. In this process derivative approximation generates boundary condition. The first condition is introduced truncation errors which cause the numerical dispersion of when solving the global equations system (Eq. 4) by a the solution. Solution stability studies are generally carried restructuring of the global matrix [K] to eliminate equa- out by a Fourier series development of the amplitude and tions corresponding to an imposed concentration node the phase of the numerical solution; they may then be and replacing the corresponding concentration variables in compared to those of the analytical solution. Guerghian et other equations by the specified values. The imposed flux al. (1980) and Fletcher (1988) have demonstrated that the condition is represented by the quantity QCi which is stability of the advection-dispersion equation by a added to the second member of the elementary equation at Galerkin finite element method formulation is related to the node i. Each node affected by this condition is consid- two parameters: the Peclet number (Pe) and the Courant ered as a source term. number (Cr). These two parameters are expressed as: (8) PepVDL/D Solution implementation The finite element spatial discretization of the governing (9) equation of pollutant migration in saturated homogenous CrpVDt/DL porous media has given a first-order differential equation with in time: V: flow velocity iC D L: characteristic length of the element (DX or DY) (5) p F for t 1 t0 [K] Cc[M] D: hydrodynamic dispersion coefficient it Dt: time step. C(to)pCo The solution stability studies relative to these parameters In the previous equation, the right-hand side vector F allow definition of the domain of applicability of a given corresponds to the boundary conditions. The resolution of numerical scheme. These parameters also allow the deterEq. (5) is in finding a C(t) function set that satisfies Eq. (5) mination of the optimum grid and time step for a stable at each time t and the initial conditions for tpt0. Replacing solution of the problem under consideration. In the literaC and its time derivatives by their approximate finite ture, a value of less than or equal to 2 is considered as a differences, writing the boundary conditions vector F for limit for the Peclet number; the Courant number generally times t and tcDt and using the concentration increment ranges between 0 and 1.6 for a Galerkin finite element DC gives: numerical scheme. [M]caDt [K] DCpDta FcDtc1Pa FtP[K] Ct (6) In order to determine the stability domain for the solution scheme used in the MIGPO model construction relative to with the Peclet and Courant numbers, a 40-element, 147-node grid is considered. Many simulations were executed for a DCpCtcDtPCt large range of these two factors in order to test the convera is the Euler constant. For a 1 0.5, Eq. (6) solution is gence and the stability of the MIGPO solution scheme. The unconditionally stable. If 0~a~0.5, the solution is stable simulation tests have shown that convergent and stable only for a time step Dt~Dtc (critical time step) given by: solutions are observed for a Peclet number ranging from

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

233

M. Zairi 7 M. J. Rouis

0.5 to 2 and a Courant number between 0.25 and 2 and for a null value of these two numbers. Non-convergent solutions are defined as those not reaching the convergence criterion for high iteration numbers. Numerical oscillation is expressed by negative or very high concentrations relative to initial concentration (contamination source). Practically, the stability domain is reached by grid adjustment, choosing an element size (Dx and Dy) that allows Peclet and Courant numbers in their optimal ranges. It is also possible to vary the time step and more accurate hydrodynamic dispersion coefficients in order to obtain a stable solution. Numerical model validation tests The efficiency of a mathematical model depends on its reliability. This is generally checked by comparisons with experimental results or analytical solutions or with numerical models whose efficiency is proved. For efficiency testing of the “MIGPO”, three validations were undertaken: (1) comparison with Ogata’s analytical solution (Ogata 1970); (2) comparison with Rowe’s semi-analytical model (Rowe and Booker 1985); (3) experimental validation by comparison with the results of laboratory tests. Validation by comparison with Ogata’s analytical solution Ogata (1970) gave a one-dimensional analytical solution of pollutant migration in porous media. In order to compare the MIGPO results with those obtained by this analytical method, a conservative solute with a concentration of 2 mg/l flowing through a 1 m thick clay layer having a porosity of 0.43 was considered. The flow velocities were of 0.05 m/year in the flow direction and zero laterally. The hydrodynamic dispersion coefficient was 1.25 and 0.46 m 2/ year in the parallel and lateral flow directions respectively. The grid used for this simulation incorporated four elements with a time step of 5 years. At xpo, a concentration-type boundary condition of the form C(o,t)pCo is assumed for both Ogata and the finite element model. For the lower boundary of the layer, at xpL, the following condition is applied to the two solutions: Fig. 1 Comparative concentration breakthrough curves calculated from MIGPO and Ogata solutions for time periods of a 5 years and b 50 years

234

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

fC L, tp0 fx

(10)

The concentration evolution relative to the initial concentration C/C0 with depth is presented in Fig. 1 for time periods of 5 and 50 years. Each of the curves plotted for the MIGPO corresponds to the concentration evolution at grid nodes which are lined in the flow direction. Figure 1 shows that the three profiles obtained from the MIGPO results have the same trend and evolution in time as that of the analytical solution. Concentration values obtained by the two methods are of the same magnitude, the greatest difference between the two results being some 0.01 (for a time of 5 years). Validation by comparison with the Rowe semi-analytical model The model POLUTE, established by Rowe and Booker (1985), was also used to verify the validity of the MIGPO results. This method involves a transformation of all the governing equations by introducing a Laplace transform; the transformed equations are solved analytically and the solutions are numerically inverted. The flow of a 2 mg/l cadmium solute through a 1 m thick saturated layer having a permeability of 10 –10 m/s and a porosity of 0.4 was considered to compare the results of the two models. In this theoretical example, only molecular diffusion was considered as a transport mechanism. The dispersion coefficient used was cadmium in an aqueous solution (7.17!10 –8 m 2/s; in Quigley et al. 1987) multiplied by a tortuosity coefficient of 0.4. The cadmium adsorption by soil particles was also modelled considering a retardation factor of 3. The evolution with depth of the ratio of the initial cadmium concentration and that after 5 and 50 years is presented in Fig. 2. The results from the two models differ only slightly for the time period of 5 years (B0.02 mg/l) with no difference discernible for the 50-year period. Concentration profiles through the layer appear similar for both models and show an increasing concentration with time.

Pollutant migration in porous media

Fig. 2a,b Comparative concentration breakthrough curves calculated from MIGPO and Rowe solutions for time periods of a 5 years and b 50 years

Table 1 Chemical composition of kiln dust and fly ash Oxide (%)

CaO

MgO

Na2O

K2O

Fe2O3

MnO

SiO2

Al2O3

SO4

CO3

Cl

Kiln dust Fly ash

42.06 19.74

1.56 3.66

1.45 7.71

2.7 1.0

1.83 6.06

0.1 0.13

13.67 31.12

4.3 18.2

6.8 7.76

24.8 4.3

0.73 0.3

Many simulation tests have demonstrated a good agreement between MIGPO and POLUTE results both quantitatively and qualitatively. This is more evident for long time periods and narrow grids. Experimental study of pollutant migration through liners The experimental study of pollutant migration through liner samples is part of a research programme on the use of fly ash and stabilised kiln dust as major components for the construction of liners. It is hoped that this may provide a potential solution to the environmental problems caused by this industrial refuse storage, particularly in five Tunisian cities having cement works and a projected landfill in their vicinity. The aim of the experimental work was to outline the mechanical and hydrogeological characteristics of the liners being developed. Numerous experimental tests were undertaken on a variety of liners, e.g. permeability, leaching and longevity tests. In this paper the pollutant migration tests for two liners are presented. The mineralogical and chemical compositions of the kiln dust and fly ash used are presented in Tables 1 and 2. As can be seen, both the alkali (Na, K) and sulphate concentrations were high. The kiln dust is characterised by high calcite and lime concentrations, while the fly ash was a CASTM type with a moderate lime concentration compared to that of kiln dust and particularly high concentrations of alkalis, iron and alumina. A number of liners were constructed using a variety of kiln dust and fly ash mixtures. The compositions of liners 1 and 2 discussed in this work are presented in Table 3. Silica steam was added to liners 1 and 2 to ensure a silica source

and to stabilise the mixtures (Rouis 1991). The prepared mixtures were put into a PVC cylinder, 70 mm in diameter and 20 mm in height. The mixture was slightly compacted in three layers and the liners thus formed were matured for 1 week in a wet room prior to being left saturated in distilled water for 3 weeks. Four weeks after preparation of the mixture, the distilled water was replaced by a solution with a known concentration of cadmium chloride (CdCl2, 2H2O). Samples of the floating solution were taken weekly over the 67 days of the test. Each time a sample was taken, it was replaced by an equivalent volume of water so that a 150-mm head was maintained above the liner. A sufficient quantity of cadmium chloride was dissolved in this additional water in order to maintain the chloride (Cl –) Table 2 Mineral composition of kiln dust and fly ash Mineral phase

Main phase

Secondary phase

Kiln dust

Calcite, lime, quartz, anhydrite, arcanite, alite Amorphous

Magnesite

Fly ash

Quartz, magnesite

Table 3 Composition of liners 1 and 2 Liner

Kiln dust (%)

Fly ash (%)

Silica steam (%)

Liner 1 Liner 2

95.24 76.19

0 19.05

4.76 4.76

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

235

M. Zairi 7 M. J. Rouis

concentration at 1000 mg/l and the cadmium (Cd 2c) concentration at 1500 mg/l. At the end of the test, the 100mm liner sample was cut into five 20-mm-thick slices. Each slice was dried for 24 h at 105 7C, then crushed and the chloride concentrations determined on a solid sample dissolved by acid. During the 67 days of testing, no water flow took place in the apparatus, which had a leachate collecting device. The pH of the liner varied from 7 on the first day to 11.5 on the seventh day and then stabilised at between 11.5 and 12. Chloride migration The chloride concentration of the interstitial water in the liners was obtained by reducing the quantities given by the chemical analyses to those of the water in the initial Fig. 4 sample prior to drying. This water content is considered to equal effective porosity of the liner, measured by the Variation in cadmium concentration in solute over liner 1 mercury porosimeter. In this way, chloride concentration variation in the solute over the liner can be determined together with the ion concentration profile in the liner. The chloride concentration of the source solution decreased with this ion migration into the liner, from 1000 mg/l at the beginning of each week (after the adjustment to the initial concentration) to 750 mg/l at the end, just before water sampling. A chloride concentration of 450 mg/l in the interstitial water at the liner base was reached after 67 days of flow (Fig. 3). This is related to the absence of any interaction between the chloride ion and the liner constituents. Cadmium migration The variation in cadmium concentration of the source solution with time for liner 1 is presented in Fig. 4. A sharp decrease in cadmium concentration was recorded from the first day of the experiment, from 1500 to 11.5 mg/l. It then stabilised at about 0.05 mg/l, even though a cadmium-rich

Fig. 5 Cadmium concentration in liners 1 and 2 at end of diffusion test

solute was added after each weekly sampling. A similar behaviour was recorded for all compositions of liner. This decrease is related to the precipitation of cadmium carbonate as a result of the increase in pH due to portlandite dissolution from the liner into the solute. This is confirmed by the cadmium concentration profiles through liners 1 and 2 (Fig. 5) which show a sharp decrease in the first 200 mm of sample. In their study of the influence of pH on selectivity and retention of heavy metals in clay soils, Yong and Phadungchewit (1993) found that the presence of carbonates in the soil contributes significantly to their retention capability. At high pH values, precipitation mechanisms, (for example as hydroxides and/or as carbonates) dominate the cation exchange capacity. It remains difficult to be confident of the real cadmium concentration in the interstitial water in the liner. It is Fig. 3 Chloride concentration in liner interstitial water at end of diffu- likely that most of the analysed cadmium in the liners is from the precipitates which cover their pore surfaces. This sion test 236

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

Pollutant migration in porous media

leads to difficulties in calculating distribution and hydro- The inability of the model to reproduce the Cd data is dynamic dispersion coefficients. related to the fact that speciation of Cd with chloride is important in evaluating its mobility. Results from theoreNumerical simulation of experimental tests tical and experimental studies (Yong and Sheremata 1991) The results of the pollutant migration tests on liners have indicated that a better insight into the characteristics of been used for the experimental validation of the MIGPO heavy metal adsorption modelling and prediction of model. For the numerical simulation of these tests the contaminant transport in soils can be obtained if the effect following procedure was adopted: of the other constituents in the solutions is taken into 1. The domain of the liner sample was divided into ten account. Thermodynamic calculations may be used in elements. predicting the effect of the speciation and complexation by 2. Throughout the test, chloride concentration was main- inorganic and organic ligands. The model presented here tained at 1000 mg/l at nodes on the limit between the should be completed with a thermodynamic calculation cadmium chloride solute and the liner. sequence to allow the prediction of such contaminant 3. Chloride diffusion coefficients were taken as transport. DxpDyp0.02 m 2/year. 4. The flow velocities Vx and Vy were taken as zero by assuming that the transport mechanisms are limited to Conclusions molecular diffusion only. 5. Conservative values of chloride ions were considered – This paper has outlined a two-dimensional finite element the retardation factor equalling 1. 6. Simulations were carried out for ten time steps of solution for the advection-dispersion equation. The techniques discussed have been implemented in the MIGPO 6.7 days each. program. In order to optimise the use of this model, a limited study of the numerical dispersion and oscillation Simulation results The simulated and measured interstitial water chloride criteria was undertaken by an examination of a large range concentrations through the liner sample over a time period of Peclet and Courant numbers. of 67 days are shown in Fig. 6a. The curves obtained with Validation of the MIGPO results was carried out by the MIGPO model represent the average concentration for comparing the results with those obtained using two other each 20-mm-thick slice of liner sample. As can be seen techniques: the Ogata analytical solution and Rowe semifrom the figure, the simulated concentrations are close to analytical solution, together with experimental results. The those obtained experimentally, with both profiles having a numerical simulation results were in good agreement with those of analytical and semi-analytical solutions, while the similar appearance. Experiments were also undertaken using the same simula- concentrations determined were of the same order of tion conditions but with a finer grid of 20 elements and 79 magnitude and with a similar profile. nodes (Dxp1.75 cm and Dyp2 cm). The time interval in The experimental study of pollutant migration through this case was 2.68 days. Chloride concentration profiles for liners indicates that the behaviour of contaminants in such the simulated liners are shown in Fig. 6b. As can be seen structures depends on the type of species considered, the from the figure, there is a closer correlation with the composition of the liner and the physical and chemical experimental results than in the first tests, suggesting that conditions and properties of the medium. The numerical the grid refinement reduces the difference between the simulation test for chloride migration gave results in good MIGPO and experimental concentration profiles. agreement with those measured during the laboratory Fig. 6a,b Comparison of liner chloride concentration profiles from MIGPO and experimental tests for a 10 elements grid and b 20 elements grid

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

237

M. Zairi 7 M. J. Rouis

tests, suggesting that this is an appropriate, conservative Ogata A (1970) Theory of dispersion in a granular medium. Fluid movement in earth materials. US Geol Surv Prof Pap model of contaminant migration. Laboratory measurement 411-I of the distribution of cadmium in the liner material indiQuigley RM, Yanful EK, Fernandez F (1987) Ion transfer by cated that more complex models are required to represent diffusion trough clayey barrieres. Proc Conf on Geotechnical the chemical and physical processes controlling its migraPractice for Waste Disposal, ASCE, Ann Arbor, 15–17 June, tion. pp 137–158

References Auriault JL, Lewandowska J (1993) Homogenisation analysis of diffusion and adsorption macro transport in porous media: macro transport in the absence of advection. Géotechnique 43(3) : 457–469 Bear J (1972) Dynamics of fluids in porous media. Elsevier, New York, 764 pp Dhatt G, Tozot G (1981) Une présentation de la méthode des éléments finis. Maloine, Paris Fletcher CAJ (1988) Computational techniques for fluid dynamics – fundamental and general techniques. Springer, Berlin Heidelberg New York Guerghian AB, Ward DS, Cleary RW (1980) A finite element model for the migration of leachate from a sanitary landfill in Long Island, New York. Water Resour Bull 16(5) : 900–906

238

Bull Eng Geol Env (2000) 59 : 231–238 7 Q Springer-Verlag

Rouis MJ (1991) Conception et caractérisation de barrières hydrogéologiques à partir de poussières de four de cimenterie. Thèse de doctorat et sciences appliquées, University of Sherbrooke, Canada Rowe RK, Booker JR (1985) 1-D pollutant migration in soils of finite depth. ASCE J Geotech Eng 111(4) : 137–158 Yong RN, Phadungchewit Y (1993) pH influence on selectivity and retention of heavy metals in some clay soils. Rev Can Géotech 30 : 821–833 Yong RN, Sheremata TW (1991) Effect of chloride ion on adsorption of cadmium from a landfill leachate Can Geotech J 28 : 378–387 Zairi M (1996) Modélisation de la migration des polluants en milieux poreux. Cas des sites de stockage de phosphogypse de Sfax. Thèse de doctorat de spécialité, Ecole Nationale d’Ingénieurs de Sfax, Sfax Zienkiewicz OC (1977) The finite element method in engineering science. McGraw-Hill, New York, 787 pp