A COUPLED THERMO-HYDRO-MECHANICAL MODEL OF EXPANSIVE CLAYS SUBJECTED TO HEATING AND HYDRATATION nadia laredj, hanifi missoum, karim bendani and mustapha maliki about the authors corresponding author Nadia Laredj Université Abdelhamid Ibn Badis de Mostaganem, Construction, Transport and Protection of the Environment laboratory D28 . A5. Rue Benanteur Charef, CP°27000, Mostaganem, Algeria E-mail: nad27000@yahoo.fr Hanifi Missoum Université Abdelhamid Ibn Badis de Mostaganem, Construction, Transport and Protection of the Environment laboratory BP°227, Route Belhacel, CP°27000, Mostaganem, Algeria E-mail: hanifimissoum@yahoo.fr Karim Bendani Université Abdelhamid Ibn Badis de Mostaganem, Construction, Transport and Protection of the Environment laboratory BP°227, Route Belhacel, CP°27000, Mostaganem, Algeria E-mail: kbendani@yahoo.fr Mustapha Maliki Université Abdelhamid Ibn Badis de Mostaganem, Construction, Transport and Protection of the Environment laboratory 12. Bâtiment F, Cité Yassmine, zone CIA, CP°27000, Mostaganem, Algeria E-mail: mus27000@yahoo.fr Abstract The focus of this work is to provide a numerical formulation for coupled thermo-hydro-mechanical processes in unsatu-rated expansive clays, especially in compacted bentonite, with a multiphase fluid flow. The model is characterized by the presence of a deformable solid matrix filled with two fluid phases (liquid water and air). In the proposed model, both pore-water and air transfers are assumed to be governed by the generalized Darcy's law. Fully coupled, nonlinear partial differential equations are established and then solved by using a Galerkin weighted residual approach in the space domain and an implicit integrating scheme in the time domain. The model has been validated against an experimental test from the literature, which involves bentonite under laboratory conditions. The calculated relative errors between the experimental and numerical results are 3% for the temperature and 7% for the stresses. Consequently, the developed numerical model predicts satisfactory results, when compared to the experimental test measures. The model is applicable to two-dimensional problems with various initial and boundary conditions; non linear soil parameters can be easily included in this model. Keywords thermo-hydro-mechanical process, unsaturated bentonite, finite element, numerical modelling, expansive clays 1 INTRODUCTION The thermo-hydro-mechanical modelling of a partially saturated soil considered as a multiphase porous medium composed of a deformable solid skeleton and fluid phases filling the pore spaces of the soil is of great interest in widely different fields of engineering. This article specifically focuses on unsaturated swelling clays, which are widely distributed in nature. In agriculture, water adsorption by the clay determines the ability of the soils to transport and supply the water and nutrients. Compacted bentonites play a critical role in various high-level, nuclear-waste isolation scenarios and in barriers for commercial landfills [1, 2]. In engineering and construction the swelling and compaction of clayey soils induce stresses that are very troublesome in the foundations and structures of buildings. The engineering behaviour of unsaturated soil has been the subject of numerous experimental and theoretical investigations [3-8]. In recent years, the coupled thermo-hydro-mechanical behaviour in porous media has been a subject of great importance in many engineering disciplines. Many mathematical models in the literature were proposed for fully saturated conditions [9-12], whereas few numerical approaches have been reported regarding the partially saturated porous materials [13-18]. These models are generally based on typical simplifying assumptions, such as the rigid soil skeleton (no solid deformation) [19], ACTA GEOTECHNICA SLOVENICA, 2011/1 29. n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation the static gas phase (no gas flow) [14, 16, 17, 20] and the quasi-static condition [14, 16, 17, 20], which are in contrast to the physics of the problem in many situations. A more commonly used formulation consists of a one-phase flow model where it is assumed that the flow of gas is negligible, and the gas phase remains constantly at the ambient external atmospheric pressure through the partially saturated soil [14-17]. Therefore, the continuity equation of the gas flow is ignored. Another problem is that the porosity is often considered as a constant or a simple function of bulk strain [21, 22]. In the present paper, the above-mentioned simplifying assumptions are excluded from the numerical model, and a fully coupled formulation is presented to simulate the behaviour of the partially saturated clays. In the model, the pore-water and pore-air transfers are assumed to be governed by the generalized Darcy's law. The fully coupled, nonlinear, partial differential equations are established and then solved by using a Galerkin weighted residual approach in the space domain and an implicit integrating scheme in the time domain. The obtained model was finally validated by means of some case tests for the prediction of the thermo-hydro-mechanical behaviour of unsaturated swelling soils. q = -i vt ~(vvPv + vr )l + (CpiViPi + CpvVvPi + CpvVaPv + CpdaVaPda )(T - Tr ) f= Hc (T -T) + LnSaPv (3) where Hc is the specific heat capacity of the soil, T is the temperature, Tr is the reference temperature, L is the latent heat of vaporization of the soil water, Cpl , Cpv and Cpa are the specific heat capacities of the soil water, the soil vapour and the soil dry air, respectively, and xt is the coefficient of the thermal conductivity of the soil. Three modes of heat transfer are included: thermal conduction, sensible heat transfer associated with the liquid, vapour and air flow and latent heat flow with vapour. 2.2 MOISTURE TRANSFER For the moisture transfer, the mass-transfer balance equation, accommodating both liquid and vapour, can be expressed as: d(pnS,) d(pvn(S, -1)) 2 THEORETICAL FORMULATION In this work a three-phase porous material consisting of solid, liquid and air requires is considered. A set of coupled governing differential equations are presented below to describe the coupled multiphase flow in the soil. The model is based on combinations of equations or derivations from conservation principles and the classical laws of known physical phenomena for coupled flow. The governing differential equations for pore water pore air and heat transfer in unsaturated soil are derived as follows: 2.1 H6RT TRANSFER Considering heat transfer by means of conduction, convection and the latent heat of vaporisation effects, and applying the principle of the conservation of energy, the following equation is derived: df = -V.Q dt (1) where 0 is the heat content of the soil and Q is the total heat flux, defined as: dt dt = -PiVvi -P,Vvv -PVV.v„ (4) where n is the porosity, p is the density, St is the degree of saturation, t is the time and v is the velocity. The subscripts l, a and v refer to the liquid, air and water vapour, respectively. In this simulation, a generalized Darcy's law is used to describe the velocities of the pore water and air: K vl =--l- (Vu, + gVz ) g (5) v = -K„ V (6) where Kt and Ka are the hydraulic conductivities of the liquid and air, respectively, yl is the unit weight of the liquid, ya is the unit weight of the air, ui is the pore water pressure, ua is the pore air pressure, z is the elevation and V is the gradient operator. The hydraulic conductivities of the water and air through the soil may be expressed in terms of the saturation degree or water content, as follows: K = K (Si) (7) 30. ACTA GeûTeCHNICA SLOVENICA, 2011/1 n. laredl et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation K = Ka (S, ) (8) where qa is the dynamic viscosity of the air. The water-vapour density pv is evaluated from the thermodynamic assumptions, and when the liquid and vapour phases are in equilibrium, it can be evaluated with the following relationship [23]: Pv =P0.ht (9) where ht is the total relative humidity calculated by the following expression: ht = exp PRT (10) and p0 is the total saturated water vapour defined as [23]: r0 = [l94.4exp(-0.06374(T - 2 73) + 0.1634.103 (T - 273)2) (11) where Rv is the gas constant for water vapour and T is the temperature. 2.3 PORE AIR MASS TRANSFER Using Henry's law to take account of the dissolved air in the pore water, the following equation is derived for the dry air phase from the principle of mass conservation: d[nPa (Sa + HS, )] dt = -V.[Pa (Va + HV, )] (12) where, H is Henry's volumetric coefficient of solubility and pa is the dry air density. The dry air density pa can be evaluated from Dalton's law as: u R Pa = — -—Pv (13) R and Rv are the gas constants for the dry air and the water vapour respectively, and T is the temperature. where the subscripts a, s and T refer to the net stress, suction and temperature contributions. The stress-strain relationship can, therefore, be expressed as: da" = D(de — dss — deT ) (15) where: SJ = K s S ^ t^ j (16) where a" is the net stress and D is the elastic matrix. A number of constitutive relationships can be employed, for example, an elasto-plastic constitutive relationship [24]. 2.5 COUPLED EQUATIONS This leads to a set of coupled, non-linear differential equations, which can be expressed in terms of the primary variables T, u, ua and u of the model as the energy balance: C Oh.+C 0L+C 0uL+C — - CTl n. ' CTT n. ' CTa n. ' CTu n. _ at at at at (17) - v[ktívu1 ] + [ktt VT ] + [KTaVua ] + jt the mass balance: ^ + ClT ^ + Cla ^ + Cu — = dt ,T dt ,a dt ,u dt (18) = V[K„Vu, ] + [KlT VT ] + [KuVua ] + J, du¡ ~dt Ca, —'- + CaT — + C C du= dt ' aa dt ' au dt (19) dT + C dt aa dt = V[Kd Vu \ + [KmVua ] + J the stress equilibrium: Culdur + CuTduT + Cuadua + Cuudu + db = 0 (20) where K^ and C^ represent the corresponding terms of the governing equations (i,j = l, T, a, u) 2.4 CONSTITUTIVE STRESS-STRAIN RELATIONSHIP For problems in unsaturated swelling porous media the total strain e is assumed to consist of components due to the suction, temperature and stress changes. This can be given in an incremental form as: de = des+ des + deT (14) 3 DISCRETIZATION TECHNIQUES The numerical solution of the theoretical models commonly used in geo-environmental problems is often achieved by a combination of numerical discretization techniques. For the example presented in this paper, the finite-element method was employed for the spatial discretization and a finite-difference time-stepping scheme for temporal discretihzation. ACTA G£OT£CHNICfl SLOVENICA, 201l/l 31. n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation In particular, the Galerkin weighted residual method [25] is used to formulate the finite-element discretization. An implicit, mid-interval backward difference algorithm is implemented to achieve the temporal discretization, since it has been found to provide a stable solution for highly non-linear problems [26]. With the appropriate initial and boundary conditions, the set of typically non-linear coupled partial differential equations can be solved. Applying a Galerkin formulation for the finite-element method we obtain a system of matrix equations, as follows: [ K ]{j} + [C ]{j} + {/ } = 0 (21) where: [ K ] = Ktt KT KTa 0 CTT c T C Ta C Tu KlT K¡¡ Ka 0úúú , [C] = c ^ T c Ca Cu 0 Kal K aa 0úú CaT CaL C aa C au 0 0 0 0 c ^ uT Cur C ua Cuu j} = (T u Ua u) , (22) j } = dT du, du du dt dt dt dt {J } = (( h Ja Ju) (23) To solve equation (21), a general form of fully implicit, mid-interval, backward difference, time-stepping algorithm is used to discretise the governing equation temporally. Therefore, equation (21) can be rewritten as: k (j )[1-q]{"+1 }+ [{j }-j} -c (j At (24) -J (j ) = {0} where (jn) is the level of time at which the matrices K, C and J are to be evaluated and is given by: (j )=w{jn+1} + (1-w){j } (25) where w is the integration factor, which defines the required time interval [0,1]) and 9 = 0, 0.5, 1 for the backward, central and forward difference schemes, respectively. For a mid-interval backward difference scheme, w = 0.5 and 9 = 0. Therefore, equation (24) reduces to: K + c Kj"} j} + IK1 } (K1 } + J (K1 } 2 At 2 (26) = {0} Eq. (26) can be rewritten as: [V+1 }-{j} K"+a5 {j} + C" At + J"+05 ={0} (27) A solution for {j"+1} can be obtained, provided the matrices K, C and J at the time interval (n + 0.5) can be determined. This is achieved by the use of a predictor-corrector iterative solution procedure. 4 APPLICATIONS AND RESULTS 4.1 EXAMPLE 1 4.1.1 PROBLEM DEFINITION The following example is investigated to demonstrate the swelling pressure calculation by using the back hydro-mechanical model. The simulation begins with an isothermal two-phase flow coupled with deformation. A free extension is allowed at the first stage to calculate the free-extension displacement on the boundary. The geometric set-up and boundary conditions are as shown in Fig. 1. The example is a compacted bentonite block, 0.025 m in length and 0.024 m in height. The element discretization is Ax = 0.00125m and Ay = 0.00124m , eight noded composed elements. The initial conditions of the system are: atmospheric air pressure and liquid saturation Sl = 0.357. A water solution enters the sample from the bottom under a pressure described by the curve in Fig. 2. The corresponding part on the boundary is fully saturated. 0.02S m IC : S1 = 0.357 u,= 1.01e5T| a Water intrusion Figure 1. Model set-up of the example. 32. ACTA GeûTeCHNICA SLOVENICA, 2011/1 n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation 1.2E+07 O.OE+OO 2.5E+05 4.0E+05 Time (sec) 6.0E+05 Figure 2. Curve of the liquid pressure on the boundary. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Liquid saturation [-] Figure 3. Computed profiles of liquid saturation along the vertical symmetric axis. The material properties for this example are based on data in the literature [27, 28] and summarized in table 1. The deformation under free extension conditions is assumed to be non-linear elastic. Table 1. Porous medium properties [27, 28]. Symbols Functions and constants Units Pi 1000 kg/m3 Pg 1.26 kg/m3 Ps 1600 kg/m3 St 1-0.85 1-exp(-1.20x10-7 (ua - u,)) ft 1.20 x10-3 Pa s ft 1.80 x10-5 Pa s Ki 1.2x10-12 M, [1 + 1.3x10-10 (ua -ui)17] m/s Ka 1.3x10-19 ^ [e(1- S,))30 Ma m/s n 0.37 S0 31.80 m2/g E 3.5 MPa V 0.3 Tr 293 °K 4.1.2 RESULTS AND DISCUSSION The simulation results of the free extension processes after 5.8 x 105 s (6.7 days) are shown in Fig. 3 and Fig. 4. With the intrusion of water from the bottom, the saturation process starts. This phenomenon can clearly be seen from the saturation evolution profile along the vertical symmetric axis (Fig. 3). At the early stage, the value of the liquid saturation increases quite quickly. After 4.9 x 105 s (5.7 days) the liquid pressure on the bottom sinks at 5.0 x105s (5.8 days) and reaches zero. Because of the low permeability of the bentonite, the liquid pressure at the centre of the specimen sinks more slowly, which results in the higher-pressure region in the specimen after 6.7 days, as shown in Fig. 4 (a). a) b) 0.01 X(m) X(m) Figure 4. Simulation results of the free extension process a) Distribution of liquid pressure, b) Distribution of liquid saturation ACTA GeOTeCHNICA SLOVeNICA, 2011/1 33. n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation With the intrusion of water, the sample begins to expand. After 6.7 days, the shape of the specimen should be as shown in fig. 5. The maximum width of the specimen increases by about 20% (Fig. 5). In this example, the experimental data from free swelling tests for compacted bentonites were used [27]. 0.02 0.015 0.01 0.005 ^ X(m) Figure 5. Simulated shape of the sample, and distribution of the liquid saturation for the material at t = 5.8x105 s. 4.2 EXAMPLE 2 4.2.1 PROBLEM DEFINITION The model was verified through simulating a test problem of bentonite tested in laboratory conditions [29]. The compacted bentonite is planned to be used to limit the flow velocity of groundwater in the near field of a geological repository for radioactive waste. The sample of bentonite has a height of 203 mm. To minimize the heat losses, the cell was insulated with a heat-proof envelope. Heat was applied at the bottom plate of the cylinder, while the temperature at the other end was kept constant and equal to 20 °C. A maximum temperature of 150 °C was applied. A constant water pressure was applies to the end opposite the end where the temperature variation was prescribed. Constant volume conditions were ensured in the test. The main variables that were measured during the test include the temperature, the relative humidity and the total axial stress. The model geometry is the same as the sample with a height equal to 203 mm. The temperature variation at the bottom end of the model is as specified during the test; it was raised in steps until it reached 150 °C. The temperature at the top end of the specimen was kept constant at 20 °C. The other surfaces are adiabatic. According to the test procedure, the hydraulic boundary condition is impermeable for all surfaces of the model. The gas pressure is equal to the atmospheric pressure. The initial value of the porosity is 0.3242, the initial temperature is 20 °C, the initial gas pressure is 0.10132 MPa and the initial stress of the sample is 0.5 MPa. The initial values were obtained according to the measured data from experiments [28]. According to the measured data from the bentonite [29, 30, 31], the intrinsic permeability is 1.0E-21 m2 and the specific heat capacity of the solids is 920 J/Kg. 4.2.2 RESULTS AND DISCUSSION Fig. 6 shows a comparison between the measured temperatures as a function of time. At the bottom end of the sample, the measured and simulated temperatures agree well. At the top end of the sample, the simulation underestimates the temperature slightly at the initial heating stage, but trends generally agree well. In this validation example, the calculated relative error percentage with respect to the temperature is 3%. The good agreement between the simulated and measured temperature indicates that the model can simulate the thermal response of processes having high temperature gradients. The numerical simulation shows that the thermal conductivity of bentonite, treated as a multiphase material, plays an important role in the thermal response of the whole medium. The numerical results can be further improved if the effects of the material heterogeneity can be quantified during testing. The simulated results of the vapour pressure are shown in Fig. 7. Although there are no measured results with which to compare, it shows that the variation of the vapour pressures is realistic. The evolution of the axial stress was also reasonably well simulated (Fig. 8), in the trend, and the obtained relative error percentage is 7% for this validation. The numerical calculation shows that the results of stress strongly depend on the constitutive relationship between stress and strain. A more comprehensive development of the material properties is needed to further improve the numerical capability of the code. Complex boundary conditions for further work. 34. ACTA GeûTeCHNICA SLOVENICA, 2011/1 n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation 140 80 60 40 20 0 J r Mlf 1/ T12 T12 simul T11 -T11simul T10 -T10 simul T9 T9 simul T8 T8 simul T7 -T7 simul T6 T6 simul T5 T5 simul T4 T4 simul T3 T3 simul T2 T2 simul T1 T1 simul T0 T0 simul 09/05/2003 28/06/2003 17/08/2003 06/10/2003 25/11/2003 14/01/2004 04/03/2004 23/04/2004 12/06/2004 01/08/2004 Time (days) Figure 6. Comparison of results between the simulated and measured temperatures at different locations (heights) of the sample as a function of time. Time (days) Figure 7. The simulated results of vapour pressure. ACTA GeOTeCHNICA SLOVeNICA, 2011/1 35. n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation 12 10 * 8 Q. ■-------- —1 — y y /— y / / / / Kf / / // / / V -Simulated ■ Measured 0 09/05/2003 28/06/2003 17/08/2003 06/10/2003 25/11/2003 14/01/2004 04/03/2004 23/04/2004 12/06/2004 01/08/2004 Time (days) Figure 8. Comparison of the results between the simulated and measured axial stress. 5 CONCLUSION A coupled thermo-hydro-mechanical model for a description of the behaviour of swelling porous media has been presented in this paper. A set of fully coupled, nonlinear, partial differential equations was established and then solved by using a Galerkin weighted residual approach in the space domain and using an implicit integrating scheme in the time domain. A range of simulation results was presented, detailing with the hydraulic and thermal behaviour of the expansive soil. The simulation results were compared with the experimentally measured results and it was shown that a good correlation was found in the hydraulic regime and a reasonable correlation in the thermal field. The results of the validations indicate that the model is general and suitable for the analyses of many different problems in unsaturated soils. REFERENCES [1] Gens, A. and Olivella, S. (2001). Chemo-mechan-ical modelling of expansive materials. 6th International Workshop on Key Issues in Waste Isolation Research. November, 2001. Paris (France). 463-495. [2] Andra (2005). Dossier argile. Evaluation de la faisabilité du stockage géologique en formation argileuse, Rapport Andra C.RP.ADP.04.0002. [3] Fredlund, D.G. and Morgenstern, N.R. (1977). Stress state variables for unsaturated soils. J. Geotch. Engng. Div., ASCE 103, N°. 5, 447-466. [4] Alonso, E.E., Gens, A. and Josa, A. (1990). A constitutive model for partially saturated soils. Geotechnique 40, N° .3, 405-430. [5] Cui, J. and Delage, P. (1996). Yielding and plastic behaviour of an unsaturated compacted silt. Geotechnique 46, N°. 2,291-311. [6] Cui, Y.J., Yahia-Aissa, M. and Delage, P. (2002). A model for the volume change behaviour of heavily compacted swelling clays. Engineering Geology 64, N°. 2, 233-250. [7] Fleureau, J.M., Verbrugge, J.C., Huergo, P.J., Correia, A.G. and Kheirbek Saoud, S. (2002). Aspects of the behaviour of compacted clayey soils on drying and wetting paths. Can. Geotech. J. 39, 1341-1357. [8] Saiyouri, N., Tessier D. and Hicher P.Y. (2004). Experimental study of swelling in unsaturated compacted clay, Clay minerals, 39, N° 4, 469-479. [9] Zienkiewicz, O.C., Chan, A.H.C., Pastor, M. and Paul, D.K. (1999). Static and dynamic behavior of soils: a rational approach to quantitative solution, I. Fully saturated problems, Proc. Roy. Soc. London; 429, 285-309. [10] Schrefler, B.A., Zhang, H.W., Pastor, M. and Zienkiewicz, O.C. (1998). Strain localization modeling and pore pressure in saturated and samples. Compt. Mech., 22, 266-280. 6 4 2 36. ACTA GeûTeCHNICA SLOVENICA, 2011/1 n. laredj et al.: a coupled thermo-hydro-mechanical model of expansive clays subjected to heating and hydratation [11] Li, C., Borja, R.I. and Reguciro, R.A. (2004). Dynamic of porous media at finite strain. Comput. Method Appl. Mech. Eng., 193, 3837-3870. [12] Zhang, H.W., Fu, Z.D. and Wu. JK. (2009). Coupling multiscale finite element method for consolidation analysis of heterogeneous saturated porous media. Adv. Water Resour., 32, 268-279. [13] Li, X., Thomas, H.R. and Fan, Y. (1999). Finite element method and constitutive modeling and computation for unsaturated soils. Comput. Geotech., 196, 135-159. [14] Sheng, D., Sloan, S.W., Gens, A. and Smith, D.W. (2003). Finite element formulation and algorithms for unsaturated soils. Int. J. Numer. Anal. Methods Geomech., 27, 745-765. [15] Stelzer, R. and Hofsletter, G. (2005). Adaptive finite element analysis of multi-phase problems in geotechnics. Comput. Geotech., 32, 458-481. [16] Sheng, D., Gens, A., Fredlund, DG. and Sloan, SW. (2008). Unsaturated soils: from constitutive modelling to numerical algorithms. Comput. Geotech, 35, 810-824. [17] Wang, W., Kosakowski, G. and Kolditz, O. (2009). A parallel finite element scheme for thermo-hydro-mechanical coupled problems in porous media. Comput. Geotech., 35, 1631-1641. [18] Di Rado H.A., Beneyto, PA., Mroginski, JL. and Awruch, A.M. (2009). Influence of the saturation suction relationship in the formulation of non saturated soil consolidation models. Math. Compt. Model, 49, 1058-1070. [19] Wu, YS and Forsyth PA. (2001). On the selection of primary variables in numerical formulation for modeling multiphase flow in porous media. J. Contam. Hydrol., 48, 277-304. [20] Callari C. and Abati, A. (2009). Finite element methods for unsaturated porous solids and their application to dam engineering problems. Comput. Struct., 87, 458-501. [21] Nuth, M. and Laloui, L. (2008). Advances in modeling hysteretic water retention curve in deformable soils. Soil Tillage Res., 100, 7-72. [22] Thang, A.M., Cui, Y.J. and Le, T.T. (2008). A study on the thermal conductivity of compacted benton-ites. Appl. Clay Sci., 41,179-181. [23] Edlefsen, N.E. and Anderson, A.B.C. (1943). The thermodynamics of soil moisture. Hilgardia 16, 31-299 [24] Thomas, H.R. and He, H. (1998). Modelling the behaviour of unsaturated soil using an elastoplastic constitutive relationship. Geotechnique 48, N°. 5, 589-603. [25] Zienkiewicz, O.C. and Taylor, R.L. (2000). The finite element method. Butterworth-Heinemann, 5th ed., Oxford. [26] Thomas, H.R and King, S.D. (1991).Simulation of fluid flow and energy processes associated with high level radioactive waste disposal in unsaturated alluvium. Water Ressour. Res. 22, N°. 5, 765-775. [27] Agus, S. and Schanz, T. (2005). Swelling pressures and wetting-drying curves of a highly compacted bentonite-sand mixture, in: Schanz, T. (2003), Unsaturated Soils: Experimental Studies, Proceedings of the International Conference. From Experimental Evidence Towards Numerical Modelling of Unsaturated Soils. Weimar, Germany, September 18-19, 2003, Volume 1 Series: Springer Proceedings in Physics, Vol. 93 Volume Package: Unsaturated Soils: Experimental Studies 2005, XIV, 533 p., Springer. [28] Poling, B. E., Prausnitz, J. M. and O'Connell, J. P. (2001). The Properties of Gases and Liquids, Mc Graw Hill, 5th ed., New York. [29] Gatabin, C. and Billaud, P. (2005). Bentonite THM mock up experiments. Sensor data report. CEA, Report NT-DPC/SCCME 05-300-A. [30] Loret, B. and Khalili, N. (2000). A three phase model for unsaturated soils. Int. J. Numer. Anal. Meth. Geomech. 24, 893-927. [31] Marshall, T.J. and Holmes, J.W. (1988). Soil physics. Bristol. Cambridge University Press, 2nd ed., New York. ACTA GeOTeCHNICA SLOVeNICA, 2011/1 37.