Backanalysis of thermohydraulic bentonite properties from laboratory tests

Engineering Geology 64 (2002) 91 – 115 www.elsevier.com/locate/enggeo

Backanalysis of thermohydraulic bentonite properties from laboratory tests X. Pintado*, A. Ledesma, A. Lloret Geotechnical Engineering and Geo-sciences Department, Universitat Polite`cnica de Catalunya, s/n 08034 Barcelona, Spain Received 10 March 2000; accepted 16 March 2001

Abstract The paper presents a general methodology to perform backanalysis of laboratory tests involving thermohydraulic behaviour of bentonite in a systematic manner. The procedure is based on a maximum likelihood approach that defines a probabilistic framework in which error measurements and the reliability of the parameters identified can be estimated. The method is applied to the identification of some thermal and hydraulic properties of a bentonite specimen, using temperature and water content measurements as input data. Three basic cases have been considered: (a) thermal case, identifying the global thermal conductivity, ke, and the global specific heat, ce; (b) hydraulic case, identifying the tortuosity factor, s, and the exponent of the unsaturated permeability law, n; and (c) thermohydraulic case, in which s, n, and the saturated thermal conductivity, ksat, have been estimated. A numerical code that solves the coupled thermohydraulic equations has been used as main computational tool. A detailed description of the experimental device and a discussion on the tests results and on the parameters identified are also included. The values obtained are within the normal range of these parameters, but the method provides a systematic and consistent procedure to find the best parameters that reproduce the measurements for the selected models. Also the method gives an insight into the model structure, and allows detecting dependence and coupling between parameters. D 2002 Published by Elsevier Science B.V. Keywords: Laboratory test; Backanalysis; Unsaturated soil; Unsaturated flow; Temperature; Swelling clay

1. Introduction Predictive modelling of termohydromechanic (THM) behaviour of clay barriers in radioactive waste repositories needs to use complex conceptual and numerical models where coupling between different processes that take place in the barrier must be taken into account. For instance, the thermal gradient originated by heat generation process in waste canisters induces flux of water in vapour phase that has influence

*

Corresponding author.

on the water movement in liquid phase and induces changes in water saturation and soil deformation. Also these changes alter notably water and gas permeabilities and the thermal conductivity of the barrier. During the recent years, numerical models have been developed to analyse the fully coupled THM behaviour of clay barriers (Olivella et al., 1994; Gawin et al., 1995; Selvadurai, 1996; Thomas and Li, 1997; Gens et al., 1998; CATSIUS CLAY project, 1998, among others). The number of parameters that is necessary to know in order to carry out these numerical modelizations is large. In a typical thermohydromechanical analysis the number of parameters needed to

0013-7952/02/$ - see front matter D 2002 Published by Elsevier Science B.V. PII: S 0 0 1 3 - 7 9 5 2 ( 0 1 ) 0 0 11 0 - 7

92

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

characterize each material may be about 15. Some of these parameters could be evaluated directly from specially designed laboratory tests, i.e. saturated water conductivity and their variation with dry density or the water retention curve and their variation with temperature. However, there are other parameters for which experimental methods to determine directly their value in highly compacted clay do not exist, i.e. variation of water permeability with saturation degree or tortuosity factor. In these cases, it is necessary to perform backanalysis of complex laboratory experiments designed to highlight the processes associated to the desired parameters. This type of tests must be considered as a boundary value problem and the backanalysis must be carried out with the same type of numerical model that is used to analyse the behaviour of the clay barrier. So far, in the context of nuclear waste disposal, laboratory test backanalysis is carried out by a trial and error technique (CATSIUS CLAY project, 1998). However, in other related disciplines as geophysics, groundwater engineering, or geotechnical engineering, it is usual to perform backanalysis in an automatic manner by means of optimization techniques (Backus and Gilbert, 1970; Carrera and Neuman, 1986; Ledesma et al., 1996a). The applicability of these techniques is clear when complex constitutive models are used and a large number of measurements of different physical magnitudes are available. In the present work a general methodology to obtain thermohydraulic parameters of bentonite from backanalysis of laboratory tests is described. The methodology is based on a maximum likelihood approach and thus it is defined within a probabilistic framework. The formulation is quite general so as to consider error measurements and prior information, and can be adopted with any general model. Its use in a thermo-hydromechanical context is theoretically possible if the laboratory measurements and the required computational tools are available. In order to obtain the laboratory data required for the identification of some thermohydraulic parameters, a new device has been developed. In the tests, a controlled heat flow is applied on a specimen with no deformation restrictions and with global water content remaining constant. These type of tests were already performed in a preliminary version back in the 1940s (Kersten, 1949). More recently Mohamed et al. (1993, 1996) and Kanno et al. (1996) have developed similar

experiments in a clay barrier context. In the tests, temperature evolution at some points of the clay specimen is monitored and final water content distribution along the sample is measured. The aim of the present work is to show the capabilities of the parameter identification techniques to analyse the results of thermohydraulic laboratory tests. First, the governing equations used to model the coupled heat and water flows are shortly described. The basis of the identification methodology is also briefly described. Then the experimental device developed and the tests results are presented. Finally, the identification techniques are applied to estimate some of the thermohydraulic parameters involved in the modelling of the test.

2. Numerical modelling of thermohydraulic behaviour of clays Finite element code ‘‘CODE BRIGHT’’ (Olivella et al., 1996) has been used to model thermohydraulic behaviour of clay. Although the code allows the study of the mechanical behaviour of soils in a coupled form, only the thermal and water flow capacities of the code have been considered. Initially, the code was developed for non-isothermal multiphase flow of brine and gas through porous deformable saline media (Olivella et al., 1994). Only a brief description of relevant parts of these equations is included here. In CODE BRIGHT, equations for mass balance were established following the compositional approach. That is, mass balance is performed for water, air and salt species instead of using solid, liquid and gas phases. Equation for balance of energy is established for the medium as a whole. The final objective is to find the unknowns (temperature T and liquid pressure, Pl) from water mass and energy balance equations. Therefore, the dependent variables will have to be related to the unknowns by means of equilibrium restrictions and constitutive equations. For example, degree of saturation will be computed using a retention curve, depending on temperature and liquid pressure. Gas pressure will be assumed constant in this analysis. The following notation will be used in writing governing equations: superscripts w refer to water species and a refer to air species; subscripts s, l and g refer to

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

solid, liquid and gas phase, respectively. P: pressure, /: porosity, q: density, j: total mass flux, i: non-advective mass flux, q: advective flux, x: mass fraction of species, h: mass content per unit volume of phase, i.e. h = xq; Sl, f Sg: degree of saturation of liquid and gaseous phases, i.e. fraction of pore volume occupied by each phase. E: specific internal energy, ic: conductive heat flux, jE: energy fluxes due to mass motion.

93

where ic is energy flux due to conduction through the porous medium, the other fluxes ( jEs, jEl, jEg) are advective fluxes of energy caused by mass motions and f Q is an internal/external energy supply. The fluxes in the divergence term include conduction of heat and advection of heat caused by the motion of every species in the medium. If soil is considered nondeformable, advective flux due to motion of solid particles (jEs) is zero.

2.1. Mass balance of water 2.3. Equilibrium restrictions Water is present in liquid and gas phases. The total mass balance of water is expressed as: @ w w w w ðh Sl / þ hw g Sg /Þ þ r  ðjl þ jg Þ ¼ f @t l

ð1Þ

where f w is an external supply of water. Volumetric mass of water in a phase (e.g. water in gas phase hgw) is the product of the mass fraction of water (xgw) and the bulk density of the phase (qg), i.e. hgw = xgwqg. The total mass flux of a specie in a phase (e.g. flux of water present in gas phase jgw) is, in general, the sum of three terms: the nonadvective flux: igw, i.e. diffusive/dispersive,  the advective flux caused by gas motion: hgwqg, where qg is the Darcy’s flux,  the advective flux caused by solid motion: /Sghgwdu/dt, where du/dt is the vector of solid velocities. 

If soil deformation is neglected and gas pressure is considered as constant, then advective fluxes of water vapour are not taken into account. Also, air balance equation is not required in this analysis. 2.2. Internal energy balance for the medium The equation for internal energy balance for the porous medium is established taking into account the internal energy in each phase (Es, El, Eg): @ ðEs qs ð1  /Þ þ El ql Sl / þ Eg qg Sg /Þ þ r @t ðic þ jEs þ jEl þ jEg Þ ¼ f Q

ð2Þ

Mass of water vapour per unit volume of gas phase (hgw) is related to soil suction, ( Pg  Pl), by psychrometric law: w 0 hw g ¼ ðhg Þ e

ðP P ÞM

l w  RðTgþ273:15Þq

l

ð3Þ

where (hgw)0 is the vapor density at null suction (in contact with a planar surface) and is dependent of temperature. Mw is the molecular mass of water and R the gas constant. The air dissolved in liquid phase is related to partial pressure of air ( Pa) through Henry’s law: ql Ma Pa ð4Þ HMw where Ma is the molecular mass of air and H the Henry’s constant.

hal ¼

2.4. Constitutive equations 2.4.1. Hydraulic Liquid advective flux follows Darcy’s law with a relative permeability dependent on liquid degree of saturation: kkrl ðrPl  ql gÞ ð5Þ ll where k is the intrinsic permeability tensor, ll is the dynamic viscosity of liquid phase, g the gravity vector and krl is the relative permeability. Dependence of relative permeability with degree of saturation (Sl) may be considered by different ways in CODE BRIGHT. In the present case, a potential law with an exponent ‘‘n’’ has been chosen: ql ¼ 

krl ¼ ðSe Þn ; Se ¼

Sl  Srl Sls  Srl

ð6Þ

94

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

where Sls and Srl are the maximum and residual degree of soil saturation. Retention curve may be characterised by a modified Van Genuchten model: Sl  Srl Sls  Srl  b   1
s  1b s bs r ¼ 1þ 1 ; P ¼ P0 P Ps r0

Se ¼

ð7Þ

w w iw g ¼ ð/qg Sg Dm IÞrxg

Dw m

ð8Þ

¼ sD ¼ s D5:9 10

12

ð273:15 þ T Þ2:3 Pg

ð9Þ

2.4.2. Thermal Thermal conductivity is used in Fourier’s law to compute conductive heat flux, i.e.: ð10Þ

Global thermal conductivity (k) depends on porosity and saturation by a geometric mean approximation: ð1Sl Þ

k ¼ kSsatl kdry

a a Eg ðJ=kgÞ ¼ Egw xw g þ Eg xg a ¼ ð2:5 106 þ 1900T Þxw g þ 1006T xg

Es ðJ=kgÞ ¼ cs T

ð12Þ

where cs is the solid phase specific heat and T the temperature expressed in
C. Alternatively, a simpler approach may be used by considering an equivalent global specific heat, ce, and an equivalent global thermal conductivity, ke, independent of water content.

!

where D is the molecular diffusion of vapour in open air in m2/s, s is the dimensionless tortuosity factor, T is the temperature in
C, Pg is the gas pressure in MPa and I the identity matrix.

ic ¼ krT

a a El ðJ=kgÞ ¼ Elw xw l þ El xl a ¼ 4180T xw l þ 1006T xl

where b, bs, Ps and P0 are material parameters, s is the soil suction, r0 the surface tension at temperature for which P0 was determined and r the surface tension at temperature T. Molecular diffusion and mechanical dispersion compose the nonadvective flux of a species in a phase. In this case, mechanical dispersion is neglected. Fick’s law governs molecular diffusion of vapour in air:

with (Pollock, 1986):

Internal energies per unit mass of liquid (El), gas (Eg) and solid phase can be expressed as (Batchelor, 1983; Gens et al., 1998):

ð11Þ

where ksat and kdry are the thermal conductivities of soil in saturated and dry conditions respectively. Other laws are also available in CODE BRIGHT modelling approach.

3. Methodology for backanalysis: maximum likelihood formulation 3.1. Basic formulation During the last 20 years, a few procedures to perform inverse analysis in a systematic manner in geo-technical and geo-science problems have been proposed (Maier and Gioda, 1981; Menke, 1984; Carrera and Neuman, 1986; Gioda and Sakurai, 1987; Ledesma et al., 1996b, among others). Most of them adopt a probabilistic point of view to cope with the problem of identifying a set of parameters from field or laboratory measurements, assumed a fixed model. This is convenient as reliability of the parameters identified and previous information on those parameters can be included in the formulation. All these methods can be considered as included in the general theory of system identification or inverse analysis (Eykhoff, 1974; Tarantola, 1987; Bui, 1993). In this work, a maximum likelihood approach has been adopted. The measurements, x*, are considered random quantities, and a relation between the corresponding computed values, x, and a set of parameters, p, is assumed: x = M(p). In this case, M represents the

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

model which is generally non-linear. Then, the best estimation of the parameters is found by maximizing the likelihood of an hypothesis, L (Edwards, 1972), defined as proportional to the joint probability of errors in measuring state variables and in the prior information of the parameters: L ¼ kPðx; pÞ ¼ kPðxÞPðpÞ

ð13Þ

where k is a proportionality constant, and the errors of state variables and parameters are assumed independent. In this approach, the model is considered correct, and therefore, the differences between computed and observed variables are assigned to a measurement error. Accordingly, differences between estimated parameters and previous information on them come from the error on that prior information. Maximizing the likelihood in Eq. (13) is equivalent to maximize the probability of reproducing the errors obtained in the measurement process and in the prior information. If those errors are assumed Gaussian, it is possible to write: 1 PðxÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ð2pÞmV ACx A  1

t 1

 ðx  xÞ ðCx Þ ðx  xÞ 2 1 PðpÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ð2pÞnV AC0p A  1 0 t 0 1 0  ðp  pÞ ðCp Þ ðp  pÞ 2

Maximizing the likelihood (Eq. (13)) is equivalent to minimize the function J =  2 ln L(p), that is:

J ¼ðx  MðpÞÞt C1 x ðx  MðpÞÞ

þ ðp0  pÞt ðC0p Þ1 ðp0  pÞ þ lnACx A þ lnAC0p A þ nln2p þ mln2p  2lnk

where (x *  x) is the vector of differences between measured and computed values using the model, (p0  p) is the vector of differences between prior information and parameters to be estimated, Cx is the measurements’ covariance matrix, Cp0 is the parameters’ covariance matrix of the prior information, mV is the number of measurements and nV the number of parameters. (.)t indicates transposed matrix and A.A is the determinant symbol. In general, field data and laboratory test data could be treated as measurement and prior information, respectively.

ð15Þ

where the last three terms can be disregarded in the minimization process, as they are constant. Expression (15) is called ‘‘objective function’’, and its minimum gives the parameters that solve the identification problem. It should be pointed out that in most cases the error structure is assumed fixed. Indeed the covariance matrices can be defined in advance according to the problem. For instance, if all measurements have independent error with the same variance, rx2, then: Cx = rx2I, where I is the identity matrix. Moreover, if prior information is not available, only the first term in Eq. (15) must be used in the minimization process. Note that in this case the objective function becomes proportional to rx2 and therefore the function to be minimized is finally: J ¼ ðx  MðpÞÞt ðx  MðpÞÞ

ð14Þ

95

ð16Þ

which is equivalent to the classical ‘‘least-squares’’ criterion. If the covariance matrix is kept for the sake of consistency, it can be seen as a weighting matrix, giving more importance to the measurements with less variance. The procedure can be generalized for the case in which the error structure is difficult to know in advance. In that case, the variances can be included as additional parameters to be identified (Ledesma et al., 1996b). 3.2. Optimization procedure There is a wide range of algorithms available to find the minimum of a function. Generally, methods that use the computed gradients of the objective function seem to be more efficient, and particularly, the method of Gauss –Newton is appropriate for this quadratic function (Fletcher, 1981). In case of convergency problems, the extension of the algorithm due to Levenberg– Marquardt (Marquardt, 1963) has been used. Essentially the procedure consists of iterating in

96

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

the parameters’ space starting from an initial value, and using an increment defined by (Tarantola, 1987):

0 1 Dp ¼ Dp0 þ At C1 x A þ ðCp Þ 1 0 þ lI At C1 ð17Þ x ðDx  Ap Þ

decided in advance from the experience with the model, as the numerical errors are very sensitive to that.

where Dp0 = p0  p and p0 is the vector of parameters’ prior information, l is a scalar controlling the convergence of the algorithm and A is the sensitivity matrix (mV nV), defined as:

The maximum likelihood approach provides a statistical framework that also gives information on the reliability of the parameters identified. On one hand, an analysis of the sensitivity matrix shows the importance of each parameter on the state variables. On the other hand, the covariance of the parameters identified can be estimated by means of:

A ¼ @x=@p

ð18Þ

which is in fact a part of the gradient of the objective function. It can be shown that the hessian of the objective function is a sum of two terms, one depending on AtA and the other one depending on second derivatives of the state variables. Thus the terms between brackets in expression (17) can be interpreted as an approximation to the hessian matrix involving only first derivatives of the objective function (this is to avoid computational difficulties). Close to the minimum l should be zero, and the original Gauss – Newton algorithm is recovered. It can be shown that a large value of l implies that the vector Dp follows the gradient direction, which gives robustness to the method although it is very slow in terms of convergence. Usually a large value is adopted at the beginning, and it is decreased while the objective function decreases, until reaching zero or a very small value close to the minimum. The computation of the sensitivity matrix is one of the main aspects of the method. In some cases, an analytical expression or an approximation based on a finite element formulation can be proposed (Ledesma et al., 1996a). However, a simple numerical derivation technique can be used when the model is highly nonlinear and there is not another choice: forward difference (Ledesma and Gens, 1997) @x xðp þ DpÞ  xðpÞ  @pi Dpi

ð19Þ

where Dp is a vector which only includes an increment of one of the parameters (i.e. pi). This requires computing two direct problems for each parameter, which has proven to be successful and simple when there are not too many parameters involved in the analysis. The increment Dpi must be

3.3. Reliability of the parameters identified

1 Cp ¼ ðC0p Þ1 þ At ðCx Þ1 A

ð20Þ

which is the ‘‘a posteriori’’ covariance matrix of the parameters, computed at the minimum of the objective function. Expression (20) gives in fact a lower bound of the variances (Ledesma et al., 1996a), but it provides a reasonable insight in the evaluation of the quality of the parameters identified. It can be seen that the ‘‘a posteriori’’ variances are related to the shape of the objective function near the minimum. A ‘‘flat’’ objective function will correspond to large variations on the parameters identified. The analysis of sensitivity matrix also provides information on the quality of the measurement points used. For instance, a large variation of a measurement with respect to a particular parameter will be more useful in the identification process. Thus the analysis of those derivatives helps in the selection of a measurement strategy. However, the information from sensitivity analysis should be combined with the error structure of the measurements. It can be achieved by means of the information density matrix, Q, (Tarantola, 1987; Ledesma et al., 1996a) computed as: 1 t 1 Q ¼ AðAt C1 x AÞ A Cx

ð21Þ

This is a mV by mV matrix that becomes the identity matrix when the system is complete (equal number of measurements and parameters). Each row of Q corresponds to a measurement and provides information on the interdependency of the observations established by the model and the redundancy of the data.

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

4. Experimental work 4.1. Description of the equipment and test procedure The objective of the test is to apply a controlled flux of heat on one of the ends of a cylindrical specimen (38 mm diameter, 76 mm height) and to maintain the other end at constant temperature. A latex membrane that allows the deformation and keeps constant the overall water content, and a layer (5.5 cm thick) of heat insulating materials (deformable foam, expanded polystyrene and glass fibre) surround the specimen. All parts are located into a perspex tube that gives stiffness to the whole equipment. In order to assure the knowledge of heat flux that crosses the sample, two specimens symmetrically placed with respect to the heater are used in the tests. Fig. 1 shows a scheme of the equipment developed. The heater is a copper cylinder (38 mm diameter, 50 mm height) with five small electrical resistances inside. The resistances are connected to an adjustable feed source of direct current that allows controlling

97

the input power from 0 to 5 W. In the tests, a constant power of 2.17 W has been used, reaching steady temperatures in the range of 70 – 80
C in the hotter end of the specimen. In the cold end, a constant temperature of 30
C is maintained by flowing water in a stainless steel head in contact with the soil. A temperature regulation system keeps the temperature of the contact between the head and the soil, constant with variations smaller than 0.5
C. In order to improve the sealing of the soil and the durability of the latex membrane, in the contact with the heater, the membrane is surrounded by liquid silicone that solidifies in a few hours. In the whole specimen, the water loss during the test due to diffusion trough the membrane is about 0.1 g/day. The specimens are located in vertical position. In order to assure a good contact between the heads and the sample, a light stress (about 0.05 MPa) is applied on the upper zone of the equipment. Thermal conductivity of insulating layer has been obtained from backanalysis of temperature measurements. A value of 0.039 W/m
C was obtained and was used in further analyses as a

Fig. 1. Scheme of the experimental device.

98

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

fixed and known value. Fig. 2 shows the temperature distribution in the device in steady state conditions, computed using the formulation presented above. Axisymmetrical analyses performed with CODE BRIGHT allowed evaluating the effect of lateral loss of heat, and it has been estimated as  60% of the total heater power. This indicates the importance of performing a 2D analysis of the experiment. During the tests, temperatures in both ends and in three internal points of the specimen, located at regular intervals, are monitored by means of a data acquisition system controlled by a personal computer. The system is also used to impose a constant temperature at the cold end of the sample. Temperature measurements are concentrated in one of the two specimens of each test, whereas in the other specimen the temperature is only

measured at the central point, just to check the symmetry of heat flux. At the end of the tests, diameter change is measured in some points of the specimen with an accuracy of 0.01 mm. Finally, the soil samples are cut in six small cylinders and the water content of each one is determined. The material tested is a bentonite from the SouthEast of Spain (FEBEX bentonite). Its geomechanical and geochemical properties are described extensively by Villar (1994) and ENRESA (1998). Bentonite has been compacted at dry density of 1.63 g/cm3 and with a water content of 15.33% (degree of saturation of 0.63). Initial temperature of specimen was 22
C and during 7 days a power of 2.17 W in the heater and a temperature of 30
C in the cold end of the specimen was maintained.

Fig. 2. Steady temperature distribution in
C computed in the specimen, heater and insulation materials.

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

99

Fig. 3. Isochrones of temperature measured in the specimen.

4.2. Experimental results Temperatures measured during the heating period are shown in Figs. 3 and 4. Temperature reaches a quasi-steady regime at 10 h after starting the test. Fig. 5 shows water content measured at the end of the test in two samples. The amount of water content changes is an evidence of the magnitude of water flows in

both liquid and vapour phases. The differences observed (lower than 1%) between both samples can be caused by the destructive method used to measure the final water content. Moreover, average final water content is slightly lower than the initial. This discrepancy is due to a small water loss through the elastic membrane and through the membrane-heater contact.

Fig. 4. Evolution of temperature in some points of the sample.

100

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

Fig. 5. Water content measured at the end of the test (7 days of heating).

The change of diameter measured in the soil sample after the end of test is shown in Fig. 6. In the most part of the specimen the shrinkage due to the drying process originated by flow of water to the cold end of the specimen prevails over the dilation originated by temperature increment. Near the zone where the temperature increment has been lower, a significant diameter increment can be observed due to swelling of bentonite.

5. Backanalysis of the experiments Using the information provided by a few heating experiments, an identification of some selected thermohydraulic properties of the bentonite used in FEBEX project (ENRESA, 1998) has been attempted. To do that, the procedure described above has been followed. When identifying parameters using optimization techniques, it is usual to perform first simple analyses, and

Fig. 6. Change of specimen diameter after 7 days of heating.

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

ed to be done in the future. However, it is not possible to cope with the full identification without acquiring experience from previous analyses considering only a few parameters, which are the ones presented here. Three cases have been considered: (a) Thermal case. Using temperatures from the specimen at different times as measured quantities, the thermal properties of the bentonite (equivalent global conductivity, ke, and the equivalent global specific heat, ce, are identified. Only the thermal problem is considered in this case. (b) Hydraulic case. Using water content measurements obtained at the end of the test, at different points of the sample, two hydraulic parameters of the bentonite (tortuosity factor, s, and the exponent of the unsaturated permeability law, n) are identified. In this case the thermo-hydraulic problem is being computed, but thermal parameters are assumed known. (c) Thermohydraulic case. Water content measurements at the end of the test and temperatures during the whole test are considered as measured variables. Two hydraulic parameters (the tortuosity factor, s, and the exponent of the unsaturated permeability law, n) and one thermal parameter (the saturated conductivity, ksat) are identified. The coupled thermo-hydraulic problem is solved in this case. The finite element mesh used to solve the direct problem is presented in Fig. 7. Note that a twodimensional analysis has been considered (with axi-

Fig. 7. Finite element mesh and thermal boundary conditions used in the analysis.

to assume some of the parameters involved already known. For instance, if only two parameters are to be identified (nV= 2), then it is possible to plot the objective functions and to follow the path of the iterative algorithm to the minimum. This is very convenient in order to understand the coupling between those two parameters as well. In this paper, a few examples considering only two or three parameters in each case have been identified. This is part of an ongoing research work, and a final identification of all the parameters involved is expect-

Table 1 Parameters and constitutive relationships used in the numerical analyses 
 1 b h ibs Retention curve Pg Pl 1b rl 1  Pss ; P ¼ P0 rr0 Se ¼ SSlsl S Srl ¼ 1 þ P

Liquid flow

kkrl ðrPl  ql gÞ ll Sl  Srl krl ¼ ðSe Þn ; Se ¼ Sls  Srl ð1Sl Þ

Sl kdry ic ¼ krT ; k ¼ ksat

Solid phase specific heat

cs ¼ 1:38T þ 732:52

Diffusion of vapour General

12 Dw m ¼ sD ¼ s 5:9 10

b = 0.30 bs = 1.5 P0 = 35 MPa Ps = 4000 MPa k = 2.77 10  21 I (m2) Sls = 1.0 Srl = 0.01 n = 3.0

ql ¼

Conductive heat flow

101

ksat = 1.15 W/m
C kdry = 0.47 W/m
C cs (J/kg
C), T (
C)




ð273:15þTÞ Pg

 2:3

w s = 0.8, Dm (m2) U = 0.402

Pg = 0.1 MPa

102

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

symmetry), and that has proven to be fundamental to take into account the lateral heat loss. The outer elements represent the insulation material, which has been discretized as well. Thermal boundary conditions are indicated in the figure. Hydraulic boundary conditions correspond to impervious limits. Calibration tests were used to identify the thermal properties of the insulation material, and therefore, all the conditions involving the simulation, apart from the bentonite properties, are defined. The default parameters used in all cases are indicated in Table 1. If any of them is to be identified, then its value is assumed unknown and obtained from the described procedure. These default values in Table 1 correspond to the ones used in analyses of the FEBEX project (Gens et al., 1998). In this case, measurements of temperature and measurements of water contents can be considered as statistically independent. Thus the corresponding covariance matrix Cx is diagonal. On the other hand, if the model (including geometry, constitutive laws, etc.) is assumed to be correct, the error can be assigned to the measurement process. For instance, regarding the thermocouples, the standard deviation of the measurements that can be assigned to those devices (rterm) is about

0.4
C. On the other hand, the water content measurements have an estimated standard error of about 0.1%, obtained from the error produced in the weighting procedures (rweight). Apart from that error, there are other sources of uncertainty in the measurement process, i.e. the error in the position of the thermocouples, the error when controlling the time, etc. The location error is considered the most important of that kind and has been estimated directly, according to the experience. For this test, it has been assumed that a standard deviation of rxT = 0.5 mm for the position of the thermocouples and of rxx = 1 mm for the location of the water content measurement in the layers in which the sample is divided. The variance of a temperature (rT2) and of a water content measurement (rx2) will have to include both the error assigned to the measurement device itself and the location error. That is: r2T ¼ r2term þ ðdT =dxÞ2 r2xT ; r2x ¼ r2weight þ ðdx=dxÞ2 r2xx ;

ð22Þ

The derivatives included in Eq. (22) take into account the variation of temperature and water content with the position. For instance, in a place where

Fig. 8. Variances of temperature measurements used in backanalysis.

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

temperature changes substantially with respective to the location ‘‘x’’, the error in temperature due to an error in the location of the measuring device will be more important. As those derivatives are not constant (they change along time), the measurement variance will also change for each situation. As a first estimation, derivatives computed numerically from the temperature and water content measured in the tests have been used in the estimation of the variances. Thus, each measurement of temperature is performed at a particular point and at a particular time, and a specific associated variance can be estimated using expression (22). The same applies for water content measurements, but they are only performed at the end of the test, so time is not a variable in this case. Fig. 8 shows the computed values of the variances of temperature measurements to be used in the identification analyses. A total of 16 values have been plotted. They correspond to four thermocouples at dif-

103

ferent locations and measured at four different times arbitrarily selected. Note that time tends to increase the variance estimated for the measurement, because the derivative (dT/dx) increases also. The minimum variance corresponds to the thermocouple located in the heater, for which the location error is zero, and only rterm2=(0.4)2 applies. The variance assigned for the water content measurements computed using expression (22) changes slightly between 0.04 and 0.08 (water contents in %). The number of measurements considered in the analyses has been 16 for temperatures (as indicated above) and 6 for water contents. Although the thermocouples allowed a large number of measurements at different times, only a selection of 16 values has been incorporated in the identification process. This is to avoid combining two sources of information with a very different number of measurements. Finally, no prior information has been considered in the analyses.

Fig. 9. Contours of objective function for thermal case (a), and path followed by the iterative procedure.

104

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

Fig. 10. Values of coefficients of sensitivity matrix associated at two temperature measurements.

5.1. Case (a): thermal analysis In this case, only the thermal problem is considered. Therefore, only the 16 temperature measurements above mentioned have been used as input data. Global thermal conductivity, ke, and global specific heat, ce, are considered as parameters to be identified. This problem becomes a classical case of solving the heat flow equation with the boundary conditions indicated, and it is a simple analysis suitable for checking purposes.

As only two parameters are identified, contours of the objective function have been plotted in Fig. 9, including the path followed by the identification procedure. Note that the minimum is well defined and that the objective function is oriented along the axes, which implies some uncoupling between both parameters. In fact, it is well known that for a case where heat flow is imposed, thermal conductivity controls the steady state situation, whereas the specific heat controls the transient state, and they are quite independent.

Fig. 11. Global thermal conductivity and global heat identified in case (a).

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

105

Fig. 12. Measured and computed temperatures using the parameters identified in case (a).

Fig. 10 shows some terms of the sensitivity matrix computed with the model for different times. Those elements are derivatives of temperatures with respect to thermal conductivity and with respect to specific heat computed at two specific locations. Note that the sensitivity corresponding to the thermal conductivity is larger at the end of the test, that is, for steady state conditions. On the contrary, the sensitivity associated to the specific heat is maximum after a few hours of starting the test. Also, note that measurements closer to the heater give more sensitivity for both parameters. Fig. 11 presents the parameters identified and its evolution with the iteration procedure. Note that in just six iterations the final values are obtained. These values are in agreement with the results obtained in other experiments performed with the same bentonite (ENRESA, 1998).

The parameters identified have been used to solve the direct problem and to simulate the temperatures measured in the test. Fig. 12 shows the comparison between the observed and the computed values and it can be seen that the agreement is very good. Significant results related to this analysis are presented in Table 2. It should be pointed out that the standard deviation estimated from the differences between measured and computed temperatures is 0.67
C, which is about the error assumed for the temperature data. The variation coefficients of the parameters indicate a high reliability of the parameters identified. Also, the correlation coefficient between both parameters is negative, because during the transient phase, a high thermal conductivity would lead to lower temperatures (as heat flow is imposed), and that should be compensated to adjust the measured temperatures with a lower specific heat in order to permit a faster temperature increase.

Table 2 Results of thermal case Parameters

Differences between measurements and calculated values (d)

Values

Variation coefficient from Eq. (20)

Correlation coefficient from Eq. (20)

Objective function

Estimation of variance of d

ke, W/m
C

ce, J/kg
C

rke/ke

rce/ce

qkece = rkece/rkerce

J

rd ¼

0.879

1464

0.0098

0.034

 0.290

18.478

0.67

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 RðTmeas Tcalc Þ ,
C mV1

106

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

Fig. 13. Parameters and objective function for each iteration in case (b).

5.2. Case (b): hydraulic analysis In this case, the thermohydraulic problem has been solved, but only two hydraulic parameters have been identified: the tortuosity, s, and the exponent ‘‘n’’ of the

unsaturated permeability law, as defined in Table 1, for diffusion of water vapour and for liquid flow relationships. The rest of the parameters have been assumed known, with the values indicated in that table. Only the measurements corresponding to the water contents

Fig. 14. Measured and computed water contents using the parameters identified in cases (b) and (c).

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

measured at the end of the experiment at different locations of both specimens have been used as input data. Fig. 13 presents the evolution of the parameters during the iteration process. The value of the objective function for each iteration is also depicted. Fig. 14 shows the comparison between measured and computed water contents. As it could be expected, the model with the parameters identified gives an ‘‘averaged’’ solution; however, the agreement obtained can be considered as good. Fig. 15 shows the contours of the objective function. It is interesting to point out that small values of tortuosity are associated to bigger exponents. Therefore, it seems that there is a slight coupling between two parameters if only final water contents are used as measurements. This is reasonable, as both parameters control the amount of water at each point, and the procedure cannot distinguish between the two effects. Nevertheless, an analysis of the water flow produced

107

in the sample (in liquid and vapour form) can help in the understanding of this coupling. Fig. 16 presents these computed flows, from the heater to the cold end, in vapour form, and from the cold end to the heater, in liquid form. Flow in vapour form is more important at the beginning of the test (t < 24 h). However, under steady state conditions, the magnitude of both flows has to coincide (same magnitude but opposite direction), as water content does not change anymore. That situation is almost achieved at 61 days, for which both curves in Fig. 16a and b are similar (small differences are still found in the hotter region, but they tend to decrease with time). That could explain the correlation between ‘‘n’’ and tortuosity. For instance, a high value of ‘‘n’’ gives in Eq. (6) low values of permeability, and therefore a small water flow. That should be compensated (at steady state conditions) with a low vapour flow, which in turn requires a low tortuosity factor.

Fig. 15. Contours of objective function for case (b), and path followed by the iterative procedure.

108

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

Fig. 16. Mass flow of water in liquid (a) and vapour (b) phases at various times.

Table 3 summarizes the results of this hydraulic case. Note the negative correlation coefficient of the parameters which corresponds to the above mentioned coupling. Note also that the variation coefficients for both parameters are higher than in the thermal case: 0.16 and 0.19 compared to 0.01 and 0.03.

5.3. Case (c): thermohydraulic analysis In this case, the fully coupled thermo-hydraulic model has been used to simulate the experiment. Both temperatures obtained at different times (the same 16 measurements used in thermal case) and water con-

Table 3 Results of hydraulic case Parameters

Differences between measurements and calculated values (d)

Values

Variation coefficient from Eq. (20)

Correlation coefficient from Eq. (20)

Objective function

Estimation of variance of d

s

n

rs/s

rn/n

qsn = rsn/rsrn

J

rd ¼

0.589

2.84

0.160

0.191

 0.851

4.412

0.22

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi Rðxmeas xcalc Þ ,% mV1

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

109

Fig. 17. Parameters and objective function for each iteration in case (c).

tents measured at the end of the test (the same six measurements as that in the hydraulic case) have been considered as input data. Three parameters have been identified: one corresponding to the thermal problem (the saturated thermal conductivity, ksat), and two parameters corresponding to the hydraulic problem

(the tortuosity factor, s, and the exponent in the relative permeability law, n). Fig. 17 shows the evolution of the parameters obtained in the iteration process, and the evolution of the objective function. Note that only in three iterations the final values are almost obtained. The

Fig. 18. Evolution of computed and measured temperatures at some points of the specimen.

110

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

Fig. 19. Evolution of computed degree of saturation in some points of the specimen case (c).

direct problem has been solved using the identified values, in order to compare computed water contents and temperatures with the measured quantities employed in the identification procedure. Figs. 18 and 19 present the evolution of temperature and degree of saturation computed using the model and the parameters identified. Computed temperatures are also compared with the measured ones. This comparison was not possible for degree of saturation as it was not measured during the test. Note that temperatures reached steady conditions in about 10 h whereas around 20 days were needed to reach steady conditions in the hydraulic regime. Thus the transient thermal problem and the transient hydraulic problem

are uncoupled. However, some thermohydraulic interdependence can be appreciated in these figures: for large times, water content is clearly bigger close to the cold end than close to the heater. The bigger the water content (or degree of saturation), the higher the thermal conductivity. Therefore, a non-constant k will develop in the sample, and a non-linear distribution of temperature will be obtained: temperature close to the heater will increase, and close to the cold end will decrease. This effect can be detected in the computed values of temperatures, in Fig. 18, for times greater than 100 h. However, it is difficult to detect that in the measured temperature values because of the oscillations and the error measurements.

Table 4 Results of thermo-hydraulic case Parameters

Objective function

Values

Variation coefficient from Eq. (20)

ksat, W/m
C s n 1.286 0.610 2.89

rk/ksat 0.016

rs/s 0.053

rn/n 0.121

Objective function (global)

Objective function Objective function (contribution of (contribution of water temperature measurements) content measurements)

J 53.54

JT 48.9

qsn = rsn/rsrn

qkn = rkn/rkrn

qsk = rsk/rsrk

Estimation of variance of d qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi Tcalc Þ xcalc Þ ,
C rd ¼ Rðxmeas ,% rd ¼ RðTmeas mT V1 mx V1

 0.54

0.427

 0.510

0.975

Correlation coefficient from Eq. (20)

0.223

Jx 4.63

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

On the other hand, Fig. 19 allows the analysis of the saturation process in the sample. Note that from an average initial value of 63%, a change of about F 20%

111

is produced from this value towards the both ends of the sample. A particular behaviour at around t = 20 h can be observed, when saturation degree increases and

Fig. 20. Coefficients of the information density matrix in the minimum corresponding to different temperature measurements.

112

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

then decreases in the central part of the specimen. This is due to the fact that flow in vapour form is very important at the beginning, moving from the heater to the cold end. However, when this vapour reaches cooler zones, a condensation is produced before the water flow in the opposite direction becomes important. Table 4 presents the main results of case (c). It is important to point out that temperature measurements provide the most important contribution to the objective function (i.e. JT J Jx). Note that the error corresponding to the temperature information, represented by the objective function JT, is now larger than the objective function in the thermal analysis—case (a). This is due to the fact that now water content measurements have been included. However, the value of Jx in case (c) is similar to the objective function at the minimum in the hydraulic case. As in previous cases, the variation coefficients are higher for the hydraulic parameters than for the thermal one. The correlation coefficients indicate some interdependency between parameters. The correlation between tortuosity factor and n is due to the coupling indicated above between the vapour and the liquid flows. The negative correlation between tortuosity factor and thermal conduc-

tivity indicates that an increment of tortuosity, which increases the vapour flow, also increases the advective heat flow, and that requires a decrease of k to transport the same heat. On the other hand, k and n have a positive correlation, as increasing n reduces the liquid flow and therefore the advective heat flow carried out by the liquid water, which in turn requires increasing thermal conductivity to transport the amount of heat prescribed. A comparison between measured and computed water contents using the parameters identified in this case was presented in Fig. 14, where computed water contents from case (b) were also depicted. It can be seen that differences between water contents estimated in cases (b) and (c) are small because the value of thermal conductivity imposed in the case (b) is similar to the value identified in case (c). Finally, an analysis of the density information matrix at the minimum has been performed for the thermohydraulic case. Fig. 20 represents some rows of this matrix that correspond to temperature measurements at different times and locations. The values depicted incorporate the sensitivity matrix (that is the values of the derivatives of temperature with respect to the parameters) and the error structure included in

Fig. 21. Global thermal conductivity estimated in cases (a) and (c) and directly measured by Villar (1994).

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

the model. Those values indicate the amount of information provided by each measurement, and the interdependency of the data. Note that measurements at x = 0 (close to the heater) provide the maximum values of these coefficients, which indicates that this is the best point to measure temperatures, either in transient or in steady state conditions. That position is convenient because of various reasons: there is an important change in temperature values, and on top of that, the location error is zero. That figure also indicates that temperatures measured at about 20 h and at around 60 mm from the heater do not provide significant information in the identification process. The parameters obtained in case (c) are consistent with values presented by other authors. For instance, regarding the global thermal conductivity, Fig. 21 presents a comparison of results from case (a) and (c) with laboratory results from Villar (1994) and ENRESA (1998), which have been obtained by direct measurement using a surficial probe. Note that case (a) gives a global value for ke, independent of degree of saturation, and thus an average value has been obtained. Case (c) provided the value of ksat and expression (11) has been used to plot the corresponding curve. It can be seen that the agreement with laboratory results is good in both cases, in the range of 0.4 to 0.8 of degree of saturation, which corresponds to the values obtained in the test. The hydraulic parameters identified in cases (b) and (c), s  0.6 and n  2.9, are similar to other values reported in the literature obtained in an heuristic manner (CATSIUS CLAY project, 1998). However, what it is important to emphasize is the coupling between both parameters. That is, there are combinations of parameters that give similar results in terms of water contents. This should be kept in mind when performing predictive modelling of larger boundary value problems involving thermohydraulic behaviour.

6. Conclusions A general procedure to perform backanalysis of laboratory tests involving thermohydraulic behaviour of bentonite has been presented. The problem has been formulated using a maximum likelihood approach, which has proven to be very useful in other geo-science problems. This framework has some

113

advantages, i.e. allows the introduction of the measurement error and of the prior information available in a consistent way. An estimation of the reliability of the parameters identified can be obtained as well. The method requires to minimize an objective function which depends on the differences between measured and computed variables. The method has been applied to the identification of the parameters controlling the thermohydraulic problem in an experiment involving the heating of two bentonite specimens. To solve the direct problem a numerical code (CODE BRIGHT) has been used, which solves the couple thermo-hydraulic equations within a finite element framework. The sensitivity matrix (derivatives of the measurements with respect to the parameters), which is required in the identification procedure, has been computed using a numerical differentiation technique. The Gauss –Newton algorithm has been used for minimization purposes. Details of the experimental set up have been described. The calibration phase of the experiment showed that in spite of the insulation of the sample, a lateral heat loss of 60% was produced. Therefore, a 2D analysis was required in order to simulate properly the problem. The calibration stage allowed the determination of the thermal properties of the insulation as well, and they have been used in further analyses. Temperatures at different locations of the samples at different times, and water contents at different sections at the end of the test have been used as measurements in the identification process. Three cases have been considered for identification purposes. The first case, (a), involved only the thermal problem, for which global conductivity and global specific heat were identified. Only temperatures were assumed as input data, and the rest of the parameters were assumed known. Values identified were consistent with the usual ones obtained by trial and error procedures. Case (b) involved a hydraulic problem in which tortuosity, s, and the exponent, ‘‘n’’, of the unsaturated permeability law were identified. The full thermohydraulic problem was computed, but thermal parameters were assumed known. Water contents from both samples at the end of the test were used as input measurements. The procedure showed that there are a few combinations of both parameters that give similar result in terms of objective function. This is reason-

114

X. Pintado et al. / Engineering Geology 64 (2002) 91–115

able, as measured water content is a global quantity, and it is difficult to distinguish between water transported by liquid flow (controlled by ‘‘n’’) and by vapour diffusion (controlled by s). Also, coupling between parameters can be explained considering that when steady state conditions apply, flow in vapour form (controlled by s) and liquid form (controlled by ‘‘n’’) must be equal in magnitude but opposite in direction, to keep water content constant. Finally, case (c) involved the identification of saturated thermal conductivity, tortuosity factor, s, and the exponent, ‘‘n’’, using the coupled thermohydraulic simulation, and considering temperatures and water contents, both as input data. The analyses indicate that the reliability of the identification is higher for the thermal parameters than for the hydraulic parameters. Also, the quality of the model employed seems to be reasonable, as the error between measured and computed values is similar to the measurement error considered. The values obtained are within the normal range of these parameters, but the method provides a systematic and consistent procedure to find the best parameters that reproduce the measurements for the selected models. Also the method gives an insight into the model structure, and allows detecting dependence and coupling between parameters. Note that, according to Fig. 6, a change of volume lower than 3% has been measured in the cool end of the specimen. That value is not negligible, but may introduce only small changes on the thermohydraulic parameters considered in this work. However, it is difficult to foresee the magnitude of these changes, because of the highly nonlinear coupling of the processes involved. Therefore, future work could include the analysis of the mechanical aspects of the test, extending the identification procedure for the THM problem.

Acknowledgements The work described has been supported by the Spanish National Agency for Nuclear Waste Disposal (ENRESA) and by the EU through Project PL950071. X. Pintado is the recipient of a fellowship from CIMNE. Basic research has been supported by the Spanish Scientific and Technical Research Agency (CICYT)

through research grant PB93-0964. The authors acknowledge the support of S. Olivella from UPC, mainly in the computational aspects of the model.

References Backus, G.E., Gilbert, F., 1970. Uniqueness in the inversion of inaccurate gross Earth data. Philos. Trans. R. Soc. A 226, 123 – 192. Batchelor, G.K., 1983. Fluid Dynamics. Cambridge University Press, Cambridge. Bui, H.D., 1993. Introduction aux proble`mes inverses en Me´canique des Mate´riaux. Eyrolles, Paris. Carrera, J., Neuman, S.P., 1986. Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information, 2. Uniqueness, stability and solution algorithms, 3. Application to synthetic and field data. Water Resour. Res. 22, 199 – 242. CATSIUS CLAY project, (1998). Calculation and testing of behaviour of unsaturated clay as barrier in radioactive waste repositories. Compiled by Alonso E.E. and Alcoberro J. CIMNE. European Commission. DOC XII/286/98-EN. Brussels. Edwards, A.W.F., 1972. Likelihood. Cambridge University Press, Cambridge. ENRESA, 1998. FEBEX. Bentonite: origin, properties and fabrication of blocks. Pulicacio´n Te´cnica ENRESA, 05/98, Madrid. Eykhoff, P., 1974. System Identification. Parameter and State Estimation. Wiley, Chichester. Fletcher, R., 1981. Practical Methods of Optimization: 1. Unconstrained Minimization, 2. Constrained Minimization. Wiley, Chichester. Gawin, D., Baggio, P., Schrefler, B.A., 1995. Coupled heat, water and gas flow in deformable porous media. Int. J. Numer. Methods Fluids 20, 969 – 987. Gens, A., Garcı´a-Molina, A., Olivella, S., Alonso, E.E., Huertas, F., 1998. Analysis of a full scale in situ test simulating repository conditions. Int. J. Numer. Anal. Meth. Geomech. 22, 515 – 548. Gioda, G., Sakurai, S., 1987. Back analysis procedures for the interpretation of field measurements in geomechanics. Int. J. Numer. Anal. Methods Geomech. 11, 555 – 583. Kanno, T., Kato, K., Yamagata, J., 1996. Moisture movement under a temperature gradient in highly compacted bentonite. Eng. Geol. 41, 287 – 300. Kersten, M.S., 1949. Laboratory research for the determination of the thermal properties of soils. Final Report. Corps of Engineers, U.S. Army. St Paul District. Research Laboratory Investigations. Engineering Experiment Station. University of Minnesota. Ledesma, A., Gens, A., 1997. Inverse analysis of a tunnel excavation problem from displacement and pore water pressure measurements. In: Sol, H., Oomens, C.W.J. (Eds.), Material Identification Using Mixed Numerical and Experimental Methods. Kluwer Academic Publ., Dordrecht, pp. 163 – 172. Ledesma, A., Gens, A., Alonso, E.E., 1996a. Estimation of parameters in Geotechnical Backanalysis: I. Maximum likelihood approach. Comput. Geotech. 18 (1), 1 – 27. Ledesma, A., Gens, A., Alonso, E.E., 1996b. Parameter and var-

X. Pintado et al. / Engineering Geology 64 (2002) 91–115 iance estimation in geotechnical backanalysis using prior information. Int. J. Numer. Anal. Methods Geomech. 20, 119 – 141. Maier, G., Gioda, G., 1981. Optimization methods for parametric identification of geotechnical systems. In: Martins, J.B. (Ed.), Numerical Methods in Geomechanics. Reidel, Boston, MA Nato ASI Series. Marquardt, D.W., 1963. An algorithm for least-squares estimation of non-linear parameters. J. Soc. Ind. Appl. Math. 11, 431 – 441. Menke, W., 1984. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press. Mohamed, A.M.O., Yong, R.N., Caporouscio, F., Cheung, S.C.H., Kjartanson, B.H.A., 1993. Coupled heat and water flow apparatus. Geotech. Test. J. 16 (1), 85 – 99. Mohamed, A.M.O., Yong, R.N., Onofrei, C.I., Kjartanson, B.H., 1996. Coupled heat and moisture flow in unsaturated swelling clay barriers. Geotech. Test. J. 19 (2), 155 – 163. Olivella, S., Carrera, J., Gens, A., Alonso, E.E., 1994. Non-isothermal multiphase flow of brine and gas thorough saline media. Transport Porous Media 15, 271 – 293.

115

Olivella, S., Gens, A., Carrera, J., Alonso, E.E., 1996. Numerical formulation for a simulator (CODE-BRIGHT) for the coupled analysis of saline media. Eng. Comput. 13 (7), 87 – 112. Pollock, D.W., 1986. Simulation of fluid flow and energy transport processes associated with high-level radioactive waste disposal in unsaturated alluvium. Water Resour. Res. 22 (5), 765 – 775. Selvadurai, A.P.S., 1996. Heat-induced moisture movement in a clay barrier: I. Experimental modelling of borehole emplacement. Eng. Geol. 41, 239 – 256. Tarantola, A., 1987. Inverse Problem Theory. Methods for Data Fitting and Model Parameter Estimation. Elsevier, Amsterdam. Thomas, H.R., Li, C.L.W., 1997. Modelling heat and moisture transfer in unsaturated soil using a finite difference self-implicit method on parallel computers. Int. J. Numer. Anal. Methods Geomech. 41, 409 – 419. Villar, M.V., 1994. Modelling and validation of the thermal – hydraulic – mechanical and geochemical behaviour of the clay barrier. Final Report 1991 – 1994. CIEMAT, Madrid.