INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS, VOL.

22, 49—71 (1998)

AN EXAMINATION OF THE VALIDATION OF A MODEL

OF THE HYDRO/THERMO/MECHANICAL BEHAVIOUR
OF ENGINEERED CLAY BARRIERS

H. R. THOMAS1*, Y. HE1 AND C. ONOFREI2

1 Cardiff School of Engineering, University of Wales Cardiff, P.O. Box 925, Cardiff, CF2 1YF, U.K.
2 Whiteshell Laboratories, Pinawa, Manitoba R0E 1L0, Canada

SUMMARY
This paper focuses attention on the development of a numerical model of the hydro/thermo/mechanical
behaviour of unsaturated clay and its consequent verification and validation. The work presented describes
on-going collaboration between the Cardiff School of Engineering and Atomic Energy of Canada. The
model development, which was carried out at Cardiff, can be described as being based on a mechanistic
approach to coupled heat, moisture and air flow. This is then linked to a deformation analysis of the material
within a ‘consolidation’ type of model. The whole is solved via the finite element method to yield a computer
software code named COMPASS (COde for Modelling PArtly Saturated Soil). Some aspects of verification
and validation of the model have been addressed in-house. However, the purpose of current AECL work is
to provide an independent, rigorous, structured programme of validation and the paper will also explore the
further validation of COMPASS within this context. ( 1998 by John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
(No. of Figures: 11 No. of Tables: 0 No. of Refs: 25)
Key words: unsaturated soil; heat transfer; moisture transfer and stress—strain behaviour; model and
validation

1. INTRODUCTION
A comprehensive analysis of coupled heat and moisture transfer in a deformable unsaturated soil
remains a major research problem of interest in nuclear waste disposal. In particular, in order to
obtain an improved understanding of the behaviour of the engineered clay barrier, heat transfer,
moisture migration, air transfer, stress equilibrium and stress—strain behaviour need to be
considered. Since the interrelated effects of these various phenomena are too complicated for an
analytical solution to be achieved, an effective numerical approach has to be employed for the
modelling work.
Coupled heat and moisture transfer in a rigid porous media under gradients of temperature
and moisture content has been extensively studied based on a model originally proposed by
Philip and de Vries.1 Adopting a mechanistic approach based on the microscopic phase interac-
tion of the liquid, vapour and porous structure, de Vries2 developed the theory further and

*Correspondence to H. R. Thomas, Cardiff Echool of Engineering, University of Wales Cardiff, P.O. Box 925, Cardiff,
CF2 1YF, U.K.

CCC 0363—9061/98/010049—23$17.50 Received 30 March 1997

( 1998 by John Wiley & Sons, Ltd. Revised 4 April 1997
50 H. R. THOMAS E¹ A¸.

presented a detailed formulation for coupled heat and moisture transfer in a rigid porous medium
under the influence of gradients of temperature and moisture. The approach has established itself
in the literature and a significant number of papers on heat and mass transfer in porous media
have subsequently appeared based on solutions of the model and its further development and
extension.3~7
The stress—strain behaviour of unsaturated soil has been the subject of numerous experimental
and theoretical investigations. A number of constitutive relationships have been proposed to
describe soil behaviour,8~10 incorporating the so-called state surface approach and more recently
elasto-plastic type models.
The concept of a state surface was first suggested by Bishop and Blight.11 Fredlund12 proposed
constitutive equations to describe a state surface using the logarithmic form. Further develop-
ment of the state surface approach were reported by Lloret and Alonso.13
The advantage of the state surface approach is that both wetting collapse and swelling
characteristics due to the effect of suction and stress interaction can be accommodated.
The uniqueness of the state surface for the void ratio and degree of saturation has been
experimentally verified by Matyas and Radhakrishna14 for the case of monotonic changes in net
stress and suction. However, experimental observations indicate that this unique relationship is
lost if soil is subjected to loading/unloading or wetting/drying cycles, which give rise to hysteresis
effect.
An elasto-plastic model was proposed by Alonso et al.9 formulated within the framework
of strain hardening plasticity using net mean stress and suction as its primary stress variables.
The model is able to represent many features of unsaturated soil in a consistent and unified
manner. Combined with a soil liquid flow balance equation and a stress equilibrium equation, it
can be applied to simulate the wetting-loading collapse and drying shrinkage behaviour of
unsaturated soil. When saturation is reached, the model becomes a conventional critical state
model.
Liquid, vapour, air and heat transfer occur simultaneously in soil. The approach presented
here treats each flow independently and relates the velocities of flow to the gradients of
relevant potentials using the generalized laws of Darcy and Fourier. Both a non-linear elastic
state surface approach and an elasto-plastic constitutive model have been incorporated in the
work performed.
Having developed a model of the hydro/thermo/mechanical behaviour of unsaturated soil, its
accuracy, both in terms of software integrity and its ability to describe the relevant physical
phenomena, is matter of interest. Within the context of the safe disposal of high-level nuclear
waste, such issues are clearly of paramount importance. This paper therefore also addresses same
aspects of verification and validation of the model, both within conventional academic considera-
tion and also the Canadian Nuclear Fuel Waste Management Program.

2. VALIDATION OF NUMERICAL MODELS FOR USE IN THE EVALUATION

OF AN ENGINEERED CLAY BARRIER SYSTEM—GENERAL CONCEPTS
The Canadian Nuclear Fuel Waste Management Program (CNFWMP) is evaluating the concept
of disposal of nuclear fuel waste in an engineered vault at a depth of 500 to 1000 m in the plutonic
rock of the Canadian Shield. In the engineered barrier system designs being developed, the waste
would be contained within durable containers that, in turn, would be isolated from the host rock
by clay-based materials.

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 51

Both the design and the performance assessment rely on experiments performed on physical
models of vault elements over relatively short times and on information inferred from calculations
(mathematical models) that simulate the probable behaviour of the system in the space—time
domain of interest. Within CNFWMP, one of the more important goals is to determine whether
the simulation models used are adequate tools that represent the probable behaviour of the real
engineered barriers.
The general procedure used in the analysis of an engineering barrier system is based on system
analysis, where the starting-point of the analysis is a conceptual system, which, by definition, is an
abstraction. A simulation model that simplifies the conceptual system is developed, in which
numerical values are given to the general statements of size, magnitude and influences. After the
simulation model is judiciously tested (validated) it is used as a direction to build the new product,
the concrete object. The end product is in a sense a model of the conceptual system, and the
simulation model sits between the concept and physical reality. Obviously, a comparison of the
simulation model outputs with measurements of the concrete object before the object is built is
not possible. Difficulties encountered in the validation of a simulation model are generally
resolved by developing physical models of the conceptual system or of components of the system,
and testing the model by comparing the outputs from a simulation with the outputs of a physical
model.
The focal goal of the validation activity is to establish in a transparent fashion the process by
which the repository developers will demonstrate a level of confidence in models used to estimate
the probable behaviour of engineered clay barriers. The validation term, with its pragmatic
meaning, is used in the programme to define the activity of testing these models that will lead to
a reasonable assurance that the simulation results are acceptable. In this context, validation is
concerned with the aspects of appropriateness and plausibility of a model to be used in solving
practical problems. Results from a series of tests performed within this activity will serve as
tangible evidence regarding the success of the model in representing the system of interest. If
a simulation model is to be successful, that is, used within a real-world situation, it must give
information or predictions that are clearly better, in some way, than the mental image or other
abstracted model that would be used instead. The validation phase is the interface between model
development and model application. Results from a series of tests performed within this activity
will serve as tangible evidence regarding the success or the degree of success of a model in
representing the system of interest.
COMPASS and several other models are included in the validation programme. Several
facets of model validity are considered in this programme. Two types of validity, the so
called event validity and face validity, are the most convincing tests regarding the usefulness
of simulation models. Event validity, i.e. ‘validation through model testing with experiments’
is essential since a model derived from theory, no matter how elegant in itself, is unlikely to
be of use solving a practical problem unless there is a practical way of assessing how well the
theory agrees with the observed data. An event in this interpretation is an occurrence related to
a concrete object; to each event belongs its place co-ordinates and time values. The tests included
in this state are based on comparisons of model outputs with the outputs (measurements) of
physical models of engineered barriers. In order to secure the integrity of these tests, a completely
independent data set, other than data sets used in calibration and sensitivity analysis,
are employed. The experiments within this activity are essentially of the same type as those
used in the model development phase: (a) bench laboratory models, (b) large-scale laboratory
models, and (c) in situ models. COMPASS preliminary tests results (event validity) using in situ

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
52 H. R. THOMAS E¹ A¸.

test of the isothermal experiment installed at AECL Underground Research Laboratory (URL)
are presented in this paper.
Face validity that is concerned with the impression of realism that the simulation makes on the
participants in the validation and application stage is also important since a simulation model,
and in general any mathematical model, is an adjunct, a supporter to professional judgement, and
never a substitute. Tests performed at this stage are primarily based on the intuition of engineers.
Therefore they are, to some extent, subjective. However, it is important to emphasize that the
results of all tests performed within the validation phase are checked continuously against the
engineers’ intuition. In this manner, face validity consolidates all other facets of the validation
phase. The validation phase is designed to start and to end with face validity. Agreement between
mathematics (simulation results) physical confirmation and intuition is regarded as a strong
indication that the simulation is successful. Alternatively, disagreement is regarded as a signal for
a need for the improvement of the existing models.

3. HYDRO/THERMO/MECHANICAL BEHAVIOUR
OF UNSATURATED CLAY—BASIC THEORY

3.1. Stress—strain constitutive relationship and governing equation

When employing a state surface approach, total strain is assumed to consist of components due
to net stress, suction and temperature changes. This can be given in an incremental form, without
loss of generality, as
de"de #de #de (1)
p s T
where the subscripts p, s and ¹ refer to the net stress, suction and temperature, respectively.
The net mean stress and suction are chosen as the primary variables of the state surface for void
ratio. Lloret and Alonso13 examined a number of forms of state surfaces of void ratio and
concluded that the following formulation gives an accurate description of soil behaviour:
e"e #a@ ln(p)#b@ ln(s )#c@ ln(p)ln(s ) (2)
0 3 3
where e is void ratio, s is the suction at the reference temperature ¹ , a@, b@ and c@ are constants.
3 3
The stress state variables p and s can be defined as10
p"1 (p #p #p )!u (3)
3 x y z !
s"u !u (4)
! -
where u is the pore liquid pressure and u is the pore air pressure.
- !
In the model developed here, temperature influence on the volumetric strain includes the
contributions of (i) thermal expansion and (ii) suction change due to temperature effect via surface
energy.
The elastic deviatoric strain induced by the deviatoric stress change q will be evluated through
a shear modulus G.
de "dq/3G (5)
q
Combining the above equations, the stress—strain relationship may be expressed as
dp@@"D(de!de !de ) (6)
6 6 64 6T
Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 53

where p@@ is the net stress which is defined as

dp@@"d(p!u ) (7)
6 6 !
and D is the elastic matrix.
The governing stress equilibrium equation for non-isothermal deformation in an isotropic
elastic unsaturated porous medium subjected to small deformation is given as follows:
L(dp )
ij #db "0 (8)
Lx i
j
where b is the body force. Substituting equations (6) and (7) into equation (8) yields
i
L Ldu
(D(de!A ds!A d¹ ))# !#db "0 (9)
Lx 6 4 T Lx i
j j
where A is a suction matrix derived from the state surface of void ratio, A is a thermal matrix
4 T
that depends on the thermal expansion coefficient of soil and the surface energy of the soil—water
system.
Experimental work indicates that significant irreversible deformation may be induceed when
wetting/collapse takes place in unsaturated soil. Furthermore large increases in suction may also
induce plastic strain. An elasto-plastic model is therefore required to accommodate these features.
In an elasto-plastic approach, total strain is assumed to include plastic components due to net
stress, suction and temperature changes, that is
de"de% #de%#de% #de1 (10)
p 4 T
where the superscripts e and p refer to the elastic and plastic strains respectively. The stress—strain
relationship can be expressed as
dp@@"D(de!de %!de % !de1) (11)
6 6 64 6T 6
Appropriate yield equations and plastic potentials need to be defined for the evaluation of the
plastic strain. In the work presented here, the elasto-plastic constitutive relationship proposed by
Alonso et al.9 is employed to determine the volumetric plastic strain.
Substituting equation (11) into equation (8) yields
L Ldu
(D(de!A ds!A d¹!de1))# !#db "0 (12)
Lx 6 4 T 6 Lx i
j j
Equation (12) is the governing equation for the elasto-plastic model. Unlike the previous elastic
approach, the method requires an iterative procedure to obtain the correct irreversible deforma-
tion.

3.2. Moisture mass transfer

The volumetric moisture content in unsaturated soil is defined as the sum of the volumetric
liquid content and the volumetric vapour content.
h"h #h (13)
- 7
( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
54 H. R. THOMAS E¹ A¸.

where h and h are the volumetric content of liquid and vapour respectively. Moisture transfer in
- 7
unsaturated soil can therefore be considered in two parts, liquid transfer and vapour transfer.
Considering the combined effects of the movement of liquid, the movement of vapour due to
vapour diffusion together with the movement of vapour in the pore air, the law of conservation of
mass for the moisture dictates that
L(h o ) L(h o )
- - # ! 7 #+ · (l o )#+ · (l o )#+ · (l o )"0 (14)
Lt Lt 6- - 67 - 6! 7
where o is the density, t is the time and l is the velocity. The subscripts l, a and v refer to liquid, air
and water vapour, respectively. h and 6 h are defined by
- !
h "nS (15)
- -
h "n!h (16)
! -
where n is the porosity and S is the degree of saturation with respect to the pore liquid.
-
In this model the soil is deformable, therefore the variation of porosity must be included. It is
necessary therefore to rewrite equation (14) in terms of porosity and degree of saturation as
follows
L(nS o ) L(nS o )
- -# ! 7 #+ · (l o )#+ · (l o )#+ · (l o )"0 (17)
Lt Lt 6- - 67 - 6! 7
The motion of pore liquid can be defined according to a generalised Darcy’s law.14
The degree of saturation of pore water may be assumed to be dependent exclusively on the
suction and net mean stress in unsaturated soil.8,10
S "S (p, s, ¹ )
- -
The state surface approach is again employed to relate the degree of saturation of pore liquid to
net mean stress, suction and temperature. Noting that suction is a function of temperature via
surface energy and differentiating the degree of saturation with respect to time yields
LS LS Lp LS Ls LS L¹
-" - # - # - (18)
Lt Lp Lt Ls Lt L¹ Lt
Considering vapour diffusion, flow is assumed to occur under a vapour density gradient. An
extended vapour velocity equation proposed by Thomas and King6 is employed, i.e.

G C DH
nD l Lo (+¹ ) Lo
l "! !5. 7 7 +t# ! 7 +¹ (19)
67 o
-
Lt +¹ L¹
where l is the velocity of vapour flow, n is porosity, D is the molecular diffusivity for vapour
7 !5.
through6 air, l is the mass flow factor, o is the density of the water vapour and t, the capillary
7 7
potential, is defined as
t"(u !u )/c (20)
- ! -
The velocity of vapour can be written in a more convenient form as follows:
l "!MK +u #K +¹#K +u N (21)
67 7- - 7T 7! !
where K , K and K are functions of the capillary potential and temperature.
7- 7T 7!
Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 55

By the application of the generalized Darcy’s law for multiphase flow in unsaturated soil,15 the
velocity of pore air is assumed to be governed by
l "!K +(u ) (22)
6! ! !
Substituting (18), (21) and (22) into equation (17), the equation of moisture transfer can be
finally written in the form of primary variables
Lu L¹ Lu Lu
-#C !#C "+ [K +u ]#+(K +¹ )#+[K +u ]#J
-6 LtN
C #C (23)
-- Lt -T Lt -! Lt -- - -T -! ! -

where u is the vector of nodal displacements and C , K and J are coefficients of the equation
lj -j
N a, u; with u in this case representing displacements -
( j"l, ¹, in general).

3.3. Dry air mass transfer

Using Henry’s law to take account of dissolved air in the pore water, conservation of dry air
flow dictates that
L[o (h #H h )]
$! ! # - #+ · [o (l #H l )]"0 (24)
Lt $! 6 ! #6-

where H is Henry’s coefficient of solubility, o , the density of the dry air, is assumed to be given
# $!
by
u R
o " ! ! 7 o (25)
$! R ¹ R 7
$! $!
where R is the specific gas constant for dry air.
$!
As stated previously, the velocity of pore air in the continuous air phase form is assumed to be
described by equation (22)16.
Substituting equations (15) and (16) into equation (24) yields
LS LS Lo Ln
no !#H no -#n (S #H S ) $!#o (S #H S ) "+ · [o (l #H l )] (26)
$! Lt # $! Lt ! # - Lt $! ! # - Lt $! 6 ! #6-

The governing equation of dry air transfer in primary variable form therefore becomes
Lu L¹ Lu Lu
-#C !#C "+[K +u ]#+[K +u ]#J
!6 LtN
C #C (27)
!- Lt !T Lt !! Lt !- - !! ! !

where C , K and J are coefficients of the equation ( j"l, T, a, u).

!j !j !

3.4. Heat transfer

From considerations of conservation of energy, the governing differential equation for heat
transfer can be written as
L'
#+ · Q"0 (28)
Lt

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
56 H. R. THOMAS E¹ A¸.

where Q is the global heat flux per unit volume and ', the heat capacity of the soil per unit
volume is defined as

'"n(o S C #o S C #o S C )(¹!¹ )#(1!n)o C (¹!¹ )#no S ¸ (29)

- - 1- 7 ! 17 ! ! 1! 3 4 14 3 7 !
C , C , C and C are the specific heat capacities of pore liquid, pore water vapour, pore air
1- 17 1! 14
and the solid particles respectively, ¹ and ¹ are the temperature and reference temperature
3
respectively, o is the density of the solid particles and ¸ is the latent heat of vaporization.
4
The heat flux per unit volume, including conduction, convenction and transfer of latent heat in
vapour form, is defined according to

Q"!j +¹#(o l #o l )¸#(C o l #C o l #C o l #C o l ) (¹!¹ ) (30)

T - 67 76! 1- - 6 - 17 - 6 7 17 7 6 ! 1! ! 6 ! 3
where j is the intrinsic thermal conductivity of the soil. Expanding equation (28) yields
T
L¹ LH Ln LS Lo
H # (¹!¹ )#¸S o #¸no !#¸nS 7!+ · (j +¹ )#¸o +(l )#¸+ · (l o )
Lt Lt 3 ! 7 Lt 7 Lt ! Lt T - N7 6! 7

#+ · [(C l o #C l o #C l o #C l o )(¹!¹ )] (31)

1- 6 - - 17 6 7 - 17 6 ! 7 1$! 6 ! $! 3
where H, the specific heat capacity of the unsaturated soil, is defined as

H"(1!n)C o #n(C S o #C S o #C S o ) (32)

14 4 1- - - 17 ! 7 1$! ! $!
Equation (31) can be finally written in primary variable form as

Lu L¹ Lu Lu
-#C !#C "+(K +u )#+(K +¹ )#(+K +u )
T6 Lt6
C #C
T- Lt TT Lt T! Lt T- N - TT T! !

#» · +(¹ )#(¹!¹ )+ · (K +(u )#K +(¹ )#K +(u ))#J (33)

TT1 3 T-2 - TT2 T!2 ! T
where C , K , » and J are coefficients of the equation, ( j"l, T, a, u).
Tj Tj TT1 T
It is necessary to set up the initial and boundary conditions to solve the set of coupled
equations. Boundary conditions take various forms as follows.

1. Prescribed primary variable /"/* on boundary S , where / represents whole unknown.

$
This is a Dirichlet condition.
2. Prescribed flux q"q* on boundary S . This is a Neumann condition.
/
3. Prescribed convection condition on boundary S , for example, q"h(¹!¹ ) for temper-
# 0
ature boundary. This is a Cauchy condition.
4. For the stress equilibrium equation, there can be loads from surface traction applied to
boundary S .
&
With appropriate initial and boundary conditions, a set of coupled, non-linear governing
equations is defined for the problem of heat and moisture transfer and deformation in
unsaturated soil in terms of pore water pressure, pore air pressure, temperature and displace-
ments.

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 57

4. HYDRO/THERMO/MECHANICAL BEHAVIOUR OF UNSATURATED

CLAY—NUMERICAL MODELLING
There are several methods of formulating the finite element discretization of a set of governing
differential equations, e.g. the variational method and the weighted residual method. The
Galerkin weighted residual method is employed here. Quadratic quadrilateral elements are coded
for pore water pressure, pore air pressure, temperature and displacement.
In the case of the mass or the energy balance equations, it is assumed that the approximation
polynomial is given in terms of a shape function N . For simplicity, the transport equation is
.
expressed in a brief form,

L(M )
j #+ · (q )"0 (34)
Lt 6j
where M refers to mass or enthalpy, q refers to the flux term, j"l, T, a.
6
Through the Galerkin method, equation (34) can be discretized as follows:

P) CN5. D
L(M )
Lt
j !+N5 · (q ) d)#
. 6j P! N5. q6 *j d!"0
2
(35)

where q* is the flux prescribed at boundary ! (Von Neumann boundary condition).

2
The stress equilibrium equation (9) for the non-linear elastic approach is discretized via the use
of shape function N
&

P) [B5D(B du#A4 du-!A4 du!!AT d¹ )] d)!P) N5& [+du!#dbi ] d)!P! N5&q6 d!"0 2

(36)

where B is the strain matrix and q is the surface traction.

Similarly, discretization of the equilibrium equation (12) for the elasto-plastic approach yields

P) [B5D(B du#A4 du-!A4 du!!AT d¹ )] d)!P) N5& [+du!#dbi ] d)!P! N5& q6 d! 2

P) 6
! [B5D de1] d)"0 (37)

The coupled governing equations, namely, equations (23), (27), (33) and (9) or (12) can thus be
discretized spatially to produce, in matrix form

C DG H C DG H G H
K K K 0 u C C C C uR f
-- -T -! - -- -T -! T6 - -
K K K 0 ¹ C C C C ¹Q f
T- TT T! # T- TT T! T6 " T (38)
K 0 K 0 u C C C C uR f
!- !! ! !- !T !! T6 ! !
0 0 0 0 u C C C C uR f
6 6- 6T 6! 66 6 6
where K and C represent the corresponding matrices of the governing equation, (i, j"
ij ij
l, T, a, u).

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
58 H. R. THOMAS E¹ A¸.

For simplicity, it is convenient to re-write equation (38) as

G H
Lu
K(uN#C "M f N (39)
Lt

where u refers to the global unknown, that is Mu ¹ u u NT.

- !
To solve equation (39), one of the useful methods for 6 transient problems is direct temporal
integration. A general form of two level difference method is employed to discretize the governing
equation temporally.17 Therefore equation (39) can be written as

Kn`1@2[sMun`1N#(1!s)MunN]#Cn`1@2[un`1!un]/*t"(1!s)MRnN#sMRn`1 N (40)

where s is the integration factor, which defines the required time interval (s3[0, 1]).
To solve the highly non-linear problem considered here, the integration factor s is recommen-
ded as 1.0, K, C and R are evaluated at the mid-interval value of the primary variables. The
scheme thus becomes the implicit mid-interval backward difference algorithm.
Re-writing equation (40) in alternate notation yields

AMun`1N"MFN (41)

where

A"Kn`1@2#Cn`1@2/*t (42)

and

MFN"MRn`1N#Cn`1@2MunN/*t (43)

Therefore a solution of un`1 is achievable, providing the matrix A and the vector F can be
obtained. This is actually achieved by the use of an iterative solution procedure.
At the beginning of each time interval, the first estimate of un`1 is assumed to be chosen as
the value at last time interval, i.e. un. Thus the value of un`1 for the first iteration can be obtained
by

Mun`1 N"A~1 MF N (44)

- - 1
where A and F are calculated at the value of un.
1 1
After the first iteration, it is possible to correct the mid-interval value of un`1 and evaluate the
matrix A and the vector F at time level n#1 . The iteration will continue until the values of un`1
2
converge.
Convergence is monitored between every successive solution and is deemed to have been
achieved when the following criterion is satisfied.

K K
un`1!un`1
i`1 i )m
un`1
i
where i is the iteration level and m is the relative tolerance.
The time step increment is controlled by two factors, maximum iterations and minimum
iterations. Should the actual number of iterations for convergence exceed the maximum specified,

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 59

the time step size is reduced. Likewise, should the iteration number be less than the minimum,
the time step size will be increased. This procedure enables a variable time step size to be
employed, which will benefit the analysis of heat and moisture transfer taking place over a
long period of time but with more rapid variations taking place during the initial stages of the
problem.

5. VERIFICATION EXERCISES
A series of exercises were designed to verify both the complete code COMPASS and various
sub-sets of the overall theory. Details of the work performed are presented below.

5.1. One-dimensional transient liquid flow with Dirichlet boundary

The first exercise was aimed at confirming the ability of the program to solve a transient liquid
flow problem. An example of pore water flow in a linear system of semi-infinite length porous
medium column was simulated.
The initial liquid pressure of the porous medium is u , at time t"0. The pressure on the face of
0
the porous medium at x"0 is raised to u , and maintained at u for all times greater than zero.
1 1
The initial and boundary conditions can be expressed as follows:

u (x, 0)"u (45)

- 0
u (0, t)"u (46)
- 1
u (R, t)"u (47)
1 0
The solution of the liquid flow equation with initial condition (45) and boundary conditions
(46) and (47) is

A A BB
x
u "u #(u !u ) 1!erf (48)
- 0 1 0 J(4at)

where erf is error function, a is the coefficient of the equation.

To achieve the finite element solution of the problem, a 20 cm long column of 100 eight noded
isoparametric elements is employed. The time step size used for time integration is 1 s. Since the
finite element mesh has a finite length, it cannot satisfy boundary condition (47). However, the
correct solution is still obtained before the transient response reaches the boundary.
The results of u /u at time"1000 s is given in Figure 1. Comparison of results indicates a very
- 1
good match.
A similar exercise to confirm the ability of the program to solve a transient air flow problem
was performed and the same conclusion reached.

5.2. Axisymmetric elastic loading-deformation test

The second exercise was designed to verify loading-settlement on an axisymmetric soil column.
A commercial FE program, PAFEC, was used to run the same test example. The final results of
both programs were the same as shown in their output files.

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
60 H. R. THOMAS E¹ A¸.

Figure 1. Comparison of FE and analytical solution—transient liquid flow

5.3. Axisymmetric heat transfer test

COMPASS was used to run a heat transfer test, the third exercise and compared with another
in-house FE program (HEAT), which is a FE program for non-linear heat flow analysis of plane
or axisymmetric problems. Both programs gave the same results for the test example.

5.4. Axisymmetric moisture transfer test with Neumann boundary

The fourth exercise was designed to check the ability of the program to solve moisture flow in
an axisymmetric soil sample with Neumann boundary condition. Another in-house FE program
was employed to run the same test example and both programs gave the same results.

5.5. Axisymmetric wetting-loading test—non-linear elastic state surface model

The fifth exercise was designed to verify the non-linear elastic state surface model when
simulating wetting-loading in an axisymmetric soil column. The soil column is soaked from
suction at 200 kPa to 100 Pa, and then the vertical stress is increased from 10 to 110 kPa keeping
the suction constant. The results of volumetric strain and void ratio are plotted against vertical
stress and suction respectively in Figures 2 and 3. During the wetting procedure, the soil sample
shows swelling behaviour as illustrated in Figure 3. In the latter loading stage, the soil sample is
compressed under the vertical load in Figure 2. The finite element results have been compared
with those obtained directly from the analytical formulation. A very good match can be seen to
have achieved.

5.6. Axisymmetric drying loading test—elasto-plastic model

COMPASS has been used to simulate a drying-loading case in unsaturated soil. The elasto-
plastic model is employed. Two stress paths involved drying steps are considered here. The first

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 61

Figure 2. Comparison of FE and analytical solution—volumetric strain versus vertical stress

Figure 3. Comparison of FE and analytical solution—void ratio versus suction

one is loaded to 600 kPa under saturated conditions and then dried to a suction of 200 kPa.
The second one is dried to suction of 200 kPa first, then loaded to 600 kPa. Due to the stiffness
induced by suction increase, the second sample produces less volumetric deformation. The results
of specific volume versus net mean stress is given in Figure 4. A very good match is obtained when
the result from COMPASS is compared with that from the analytical formulation.

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
62 H. R. THOMAS E¹ A¸.

Figure 4. Comparison of FE and analytical solution—specific volume versus net mean stress

Apart from the above exercises, a number of numerical tests have been performed using plane
stress and plane strain models. The results have been compared with either theoretical solutions
or those obtained from a commercial FE package. All tests verified that the model is able to
provide correct solutions. The model has also been used to analyse the consolidation of
unsaturated soil, which includes the movement of moisture and deformation.18 The results
obtained were compared with another in-house FE program. Again a very good match was
achieved.

6. APPLICATIONS

6.1. Isothermal moisture transfer at a constant air pressure

This section presents the application of the model to simulate an isothermal hydration experi-
ment carried out by SCK/CEN, which involves the use of computerized X-ray tomography to
measure hydration profiles in compacted plugs of powdered Boom clay.
The material employed in the research, namely Boom clay, is being considered as a suitable
backfill and buffer material for the long-term safe disposal of high-level nuclear waste. Extensive
experimental work has been carried out to determine the soil’s thermo/hydraulic physical
properties.
The Boom clay samples were produced by uniaxial compaction of the clay powder in the
permeameter cell. After compaction the cell was left for 24 h to allow the sample to reach
equilibrium conditions. To hydrate the sample the lower face of the clay plug was connected to
a 1 m vertical column of water and the upper face was in contact with the atmosphere, allowing
air to escape freely. During the hydration process X-ray tomographs were taken at 3 mm intervals
along the clay plug and the mean water content each cross-section determined.

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 63

Figure 5. Comparison of numerical and experimental results—SCK/CEN hydration test at a constant air
pressure

To apply the numerical model to the hydration experiment, thermo/hydaulic physical rela-
tionships, which were determined experimentally by SCK/CEN, are incorporated into the
code.18
The domain was discretised into a one-dimensional mesh comprising 35]2 mm isoparametric
elements and a constant time step size of 900 s was employed. The experimental and numerical
results of the hydration experiment are given in Figure 5. The results of wetting profiles and
penetration rate predicted by the numerical analysis are in good agreement with the experimental
observations.

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
64 H. R. THOMAS E¹ A¸.

Figure 6. Comparison of numerical and experimental results—a two-phase infiltration test at varying air pressure

6.2. Coupled heat and moisture transfer at constant air pressure

Coupled heat and moisture transfer at constant air pressure has been investigated by compari-
son of numerical results with experimental results from a series of laboratory controlled heating
experiments on medium sand. Excellent correlations were achieved.5

6.3. Isothermal moisture and air transfer

The isothermal air and moisture transfer subset of the model has been explored by comparison
of numerical and experimental results of ponded infiltration into coarse sand20 in Figure 6. Good
agreement between the experimental and numerical results was achieved.21

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 65

The model has also been employed to simulate the coupled transport of energy, water and
dry air in unsaturated alluvium by comparison with an alternative approach.4 Good agree-
ment between the two sets of results was observed for both the temperature and moisture
distributions.

6.4. Coupled heat, moisture and air transfer in a deformable soil

The experimental work for this exercise was performed on an unsaturated Spanish montmoril-
lonite clay by CIEMAT. The montmorillonite clay sample was uniaxially compacted at an initial
water content of 12·4 per cent to a dry density of 1·62 g/cm3 inside a stainless steel cell which is
14·6 cm high and 15cm in inner diameter.22 In the upper part of the cell, a heating device of 1·5 cm
in height and 1·0 cm in diameter is placed along the axis of the cylinder and heated up to 100°C.
The temperature distribution within the sample was measured by 9 thermocouples at different
levels.
For the experiment reported here, steady state of temperature was reached around one
hour after the heating started. At the end of the experiment, the clay sample was taken out
and systematically cut to examine the final water content and dry density through the block.

Figure 7. Comparison of numerical and experimental results—CIEMAT thermal test, plot of temperature

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
66 H. R. THOMAS E¹ A¸.

An axisymmetric finite element mesh, 146 mm high and 75 mm wide, was chosen to represent
the test sample. The temperature is fixed at 100°C at the nodes of the heater. The temperature
along the outside boundaries is fixed at 28°C due to the thermo-shower. Initial conditions for the
simulation were estimated from measured experimental data.23
Examining the experimental and numerical results given in Figure 7, it can be seen that
generally reasonable correlation is obtained. The overall pattern of temperature distribu-
tions are the same, with values at eight out of the nine experimental points quite well
matched. The temperature response throughout all regions of the sample, both those near
the heater and those near the boundary, is affected by a combination of the heat source itself
and the applied boundary conditions. However regions near the boundary are more
directly influenced by the temperature of the thermo-shower than the temperature of the heat
source.
The results presented in Figures 8 and 9 are also claimed to be encouraging. A reduction in the
degree of saturation and an increase in the void ratio near the heater occurred in the experiment.
This pattern has been matched in the numerical simulation.24,25 It is clear, however, that work is
now required to further develop the numerical simulation so as to obtain improved correlation
between numerical and experimental results.

Figure 8. Comparison of numerical and experimental results—CIEMAT thermal test, plot of degree of saturation

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 67

Figure 9. Comparison of numerical and experimental results—CIEMAT thermal test, plot of void ratio

6.5. Preliminary test results of the simulation of an in situ test at

AECL’s Underground Research Laboratory (URL)
In this exercise the model simulated a physical model, the in situ test of the isothermal
experiment installed at URL. This experiment is being carried out at the 240 m level, i.e. 240 m
below the surface, in crystalline rocks of the Canadian Shield. A 5 m deep, 1·24 m diameter
borehole forms the basis for the experiment and is used to examine the water uptake by the buffer
under in situ boundary conditions. The buffer’s hydraulic and mechanical interactions with the
surrounding rock and a concrete restraining plug is being monitored. This experiment is in
progress. Typical results, comparing observed and simulated values of the pore water pressure
head, using preliminary data, are presented in Figures 10 and 11 for 100 and 200 days, respective-
ly, from the start of the experiment.
A reasonable agreement between predicted and measured data has been achieved at each time
considered. The gradual resaturation of the buffer observed was also indicated by the simulated
values. Differences exist, however, between observed and simulated values, particularly regarding
the duration of the suction in the rock. The numerical analysis indicates that significant negative
pressure heads can be generated in the granite. The magnitude of such pressure heads and
the extent of their penetration into the granite then influences the rate of saturation of the

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
 Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)

68
H. R. THOMAS E¹ A¸.
( 1998 by John Wiley & Sons, Ltd.

Figure 10. Comparison of numerical and experimental results—AECL isothermal experiment at URL after 100 days from the start of the experiment
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 69
Figure 11. Comparison of numerical and experimental results—AECL isothermal experiment at URL after 200 days from the start of the experiment
( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)
70 H. R. THOMAS E¹ A¸.

bentonite-sand barrier. The simulation results obtained are strongly dependent on the hydrau-
lic properties of the granite. High suction in the granite results in very low hydraulic con-
ductivity in the near-field zone creating a nearly impermeable region. The modelling results
serve to highlight the importance of buffer/rock hydraulic interaction and indicate the need
for further investigation of the flow characteristics of the rock in the disturbed region near
excavation.
Although it is still early in the validation programme, it is AECL’s judgement that COMPASS
seems to be a very useful tool that can successfully estimate the behaviour of clay engineered
barriers. Their view is that the model is theoretically sound, includes the major coupled processes
that take place in the near-field, and the results of preliminary tests suggest that the simulated
values compare favourably with observed values from physical models.

7. CONCLUSIONS
The work presented in this paper describes first a model for the analysis of the coupled transport
of heat, moisture and air transfer, which is applicable to a deformable unsaturated soil. The
theoretical formulation and numerical implementation have been presented, which can accom-
modate either a thermoelastic constitutive relationship based on the state surface approach or an
elasto-plastic model appropriate to describe unsaturated soil. A computer programme COM-
PASS for the solution of coupled thermo/hydro/mechanical boundary problem of unsaturated
clay has been developed at Cardiff.
A number of tests have been carried out at the verification stage to ensure the accuracy
of the computer code. Confidence in the numerical accuracy of the software has thus been
obtained. A number of applications of the model have been performed to assess its ability to
describe the physical phenomena involved. Further confidence in the model has therefore been
achieved.
It is however recognized that further work is now required to investigate further the ability of
the model to describe accurately the complex fully coupled phenomena under consideration.
Good quality experimental data is limited and further work is necessary to provide such
information, within a structured programme of validation tests.
The model described above has been developed from a mechanistic approach, combining
together the various phenomena in an interrelated coupled manner. As such, the importance of
validation of the new model is highlighted as a feature of importance. The development of the
theoretical model and its numerical solution is seen as being fundamentally linked to a compre-
hensive on-going validation programme, in order to ensure the accuracy and practical usefulness
of the work.

ACKNOWLEDGEMENTS

The work presented here carried out at Cardiff has been part funded from a research programme
supported by the Commission of the European Union under Contract F12W-CT90-0033. This
support is gratefully acknowledged, together with the collaboration of CIEMAT and SCK/CEN
who performed some of the experiments. The Canadian Nuclear Fuel Waste Management
Program is jointly funded by AECL and Ontario Hydro under the auspices of the CANDU
Owners Group.

Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998) ( 1998 by John Wiley & Sons, Ltd.
 HYDRO/THERMO/MECHANICAL BEHAVIOUR OF ENGINEERED CLAY BARRIERS 71

REFERENCES
1. J. R. Philip and D. A. de Vries, ‘Moisture movement in porous materials under temperature gradients’, ¹rans. Am.
Geophys. ºnion, 38, 222—232 (1957).
2. D. A. de Vries, ‘Simultaneous transfer of heat and moisture in porous media’, ¹rans. Am. Geophys. ºnion, 39(5),
909—916 (1958).
3. M. Geraminegad and S. K. Saxena, ‘A coupled thermoleastic model for saturated-unsaturated porous media’,
Geotechnique, 36, 539—550 (1986).
4. D. W. Pollock, ‘Simulation of fluid flow and energy transport processes associated with high-level radioactive waste
disposal in unsaturated alluvium’, ¼ater Resour. Res., 22, 765—775 (1986).
5. J. Ewen and H. R. Thomas, ‘Heating unsaturated medium sand’, Geotechnique, 39, 455—470 (1989).
6. H. R. Thomas and S. D. King, ‘Coupled temperature/capillary potential variations in unsaturated soil’, J. Engng.
Mech., ASCE, 117 (11), 2475—2491 (1991).
7. H. R. Thomas and C. L. W. Li, ‘Modelling transient heat and moisture transfer in unsaturated soil using a parallel
computing approach’, Int. J. Numer. Anal. Methods Geomechanics, 19, 345—366 (1995).
8. E. E. Alonso, A. Gens and D. W. Hight, ‘Special problem soils-general report’ (seession 5), Proc. 9th European Conf.
Soil Mech. and Found. Engng., Dublin, 1987, pp. 1087—1146.
9. E. E. Alonso, A. Gens and A. Josa, ‘A constitutive model for partially saturated soils’, Geotechnique, 40, 405—430
(1990).
10. D. G. Fredlung and H. Rahardjo, Soil Mechanics of ºnsaturated Soil, Wiley, New York, 1993.
11. A. W. Bishop and G. E. Blight, ‘Some aspects of effective stresses in saturated and partly saturated soils’, Geotechnique,
13(3), 177—197 (1963).
12. D. G. Fredlung, ‘Appropriate concepts and technology for unsaturated soils’, Canad. Geotech. J., 16, 121—139 (1979).
13. A. Lloret and E. E. Alonso, ‘State surfaces for partially saturated soils’, 11th I.C.S.M.F.E. Vol. 2, San Francisco, 1985,
pp. 557—562.
14. E. L. Matyas, and H. S. Radhakrishna, ‘Volume change characteristics of partially saturated soils,’ Geotechnique, 18,
432—448 (1968).
15. J. Bear and A. Verruijit, Modeling Groundwater Flow and Pollution, Reidel, Dordrecht, 1987.
16. E. E. Alonso, F. Batlle, A. Gens and A. Lloret, ‘Consolidation analysis of partially saturated soils-application to
earthdam construction’, Proc. 6th Int. Conf. on Num. Methods in Geomech., Innsbruck, 1988, pp. 1303—1308.
17. R. D. Cook, Concepts and Applications of Finite Element Analysis, Wiley, New York, 1981.
18. H. R. Thomas, Z. M. Zhou and Y. He, ‘Analysis of consolidation of unsaturated soils’, Proc. Second Czechoslovak
Conf. on Num. Meth. in Geomech. vol. 1, Prague, 1992, pp. 242—247.
19. H. R. Thomas, M. R. Sansom, G. Volckaert, P. Jacobs and M. Kumnan, ‘An experimental and numerical investigation
of the hydration of compacted powered boom clay’, in: I. Smith (ed.), Numerical methods in Geotechnical Engineering,
A. A. Balkema, Rotterdam, 1994, pp. 135—142.
20. J. Touma and M. Vauclin, ‘Experimental and numerical analysis of two-phase infiltration in a partly saturated soil’,
¹ransport in porous media, 1, 27—55 (1986).
21. H. R. Thomas and M. R. Sansom, ‘Fully coupled analysis of heat, moisture and air transfer in unsaturated soil’,
J. Engng. Mech., ASCE, 121(3), 392—405 (1995).
22. M. V. Villar, J. Cuevas, A. M. Fernández and P. L. Martin ‘Effects of the interaction of heat and water flow in
compacted bentonite’, Int. ¼orkshop on ¹hermomechanics of Clays and Clay Barriers, Bergamao, 1993.
23. M. V. Villar and P. L. Martin, ‘Suction controlled oedometric tests in montmorillonite clay’, Proc. 29th Annual Conf.
of the Engineering Group of the Geological Society of ¸ondon, 1993, pp. 337—342.
24. H. R. Thomas, Y. He, A. Ramesh, Z. Zhou, M. V. Villar and J. Cuevas, ‘Heating unsaturated clay — an experimental
and numerical investigation’, in: I. Smith (ed.), Numerical methods in Geotechnical Engineering, 1994, A. A. Balkema,
Rotterdam, pp. 181—186.
25. H. R. Thomas and Y. He, ‘An analysis of coupled heat, moisture and air transfer in a deformable unsaturated soil’,
Geotechnique, 45, 677—689 (1995).

( 1998 by John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., Vol. 22, 49—71 (1998)