An alternative coupled thermo-hydro-mechanical finite element formulation for fully saturated soils

Article history: Received 16 January 2017 Received in revised form 18 July 2017 Accepted 16 August 2017 Available online xxxx


Introduction
Significant temperature effects have been exposed in a range of geotechnical engineering problems, such as oil & gas pipelines, pavements, buried power cables, ground energy storage, and the storage of high-level radioactive waste [20]. They are also important in ground-atmosphere interaction [8]. To investigate the thermal behaviour of soils, extensive laboratory and field experimental studies have been carried out (e.g. [12,16,13,19]). It is observed that temperature changes result in the variation of stress equilibrium in the ground, as well as pore fluid flow. These problems, where mechanical, hydraulic and thermal systems in soils interact with each other, with the independent solution of any one system being impossible without simultaneous solution of the others, are defined as coupled THM problems, while the interaction between the systems is referred to as the THM coupling [34]. Further work has been carried out on THM-Chemical (i.e. THMC) coupling, with example references being those of Rutqvist et al. [29] or Seetharam et al. [30]. Despite being of interest to the modelling of some phenomena, such as chemical dissolution of dissolvable materials in soils, the present paper focuses only on the THM coupling.
To solve a coupled problem, it is necessary to first develop the coupled equations based on the governing law of each physical system [35]. A number of studies have been carried out to model the coupled THM behaviour of soils, and the most extensively used numerical tool has been the finite element (FE) method.
Aboustit et al. [1,2] adopted a general variation principle to develop the governing field equations for the problem of thermoelastic consolidation. The mechanical behaviour was modelled based on the force equilibrium of the solid-fluid mixture, while the hydraulic and thermal behaviour were taken into account using the laws of volume and energy conservation, respectively. To achieve a symmetric stiffness matrix in the FE code, Aboustit et al. [1] simply assumed the same coupling term to represent the mutual effects between the mechanical and thermal systems. It is also noted that heat advection, as well as the thermal effect on pore fluid flow, were not considered. Following this pioneering work, various approaches have been proposed for modelling the coupled THM behaviour of porous media. One commonly used approach is to first formulate the mass or volume balance equations, the momentum balance equations and the energy balance equations for each phase (e.g. solid, fluid, etc.) of a porous medium (e.g. [24]). Then the coupled THM formulation for the mixture is derived by combining the governing equations for each phase (e.g. [25,22]). Another approach takes account of the behaviour of each coupling system in the soil, thus resulting in the mechanical, hydraulic and thermal governing equations for the solid-fluid mixture (e.g. [21,10,33,18]). The difference between these two approaches is that in the first the governing equations are derived for each phase and then combined, while in the second each of the governing equations is derived for the complete mixture. It should be noted that both approaches can lead to the same governing equations. However, it is also noted that, even when the same approach is followed, the adoption of different assumptions during the derivation procedure may result in differences in the final http The new THM FE formulation for saturated soils proposed in this paper adopts the second of the above two theoretical approaches, consolidates and extends the work in the literature detailing its full derivation with clearly stated and justified assumptions. The implementation platform for the new formulation is the authors' FE software ICFEP [27,28], developed specifically for geotechnical engineering analysis. As such, it enables access to source code, which is essential for numerical developments of this type. ICFEP is already capable of performing coupled hydro-mechanical (HM) simulations (consolidation and seepage) in both fully and partially saturated soils. For THM coupling a new equation governing the heat transfer is introduced, which is also able to account for the effect of pore fluid flow and stressstrain behaviour on heat transfer. Furthermore, the existing HM coupled formulation in ICFEP is further modified to account for thermal effects on pore fluid flow and thermal deformation of soils. The adopted development approach of considering each coupling system individually ensures that any of the three systems can be disabled if they are not active in the analysis, thus reducing the full THM formulation straightforwardly to HM, TM or TH coupling. The paper further discusses the performance of the new formulation in comparison with those found in the literature. Lastly, a fully coupled THM FE analysis is performed to model the thermally induced excess pore water pressures in an undrained triaxial heating test. The results demonstrate the ability of the new formulation to predict excess pore water pressure generation due to temperature changes, without resorting to soil thermo-plastic constitutive modelling. The presented formulation adopts a tension positive sign convention, while displacement, pore fluid pressure and temperature are adopted as nodal degrees of freedom.

Mechanical equation
If only the mechanical behaviour of soils is considered, the FE formulation can be developed based on the soil constitutive behaviour, expressed as: where {Dr} and {De} are the incremental total stress and strain vectors, respectively, and [D] is the constitutive matrix. With the FE formulation derived from Eq. (1), only the extreme soil conditions, i.e. fully drained or fully undrained, can be simulated. If soil behaviour is somewhere between these two extremes, the principle of effective stress is adopted to formulate the mechanical governing equations for fully saturated soils in an isothermal HM coupled analysis [27]. Therefore, Eq. (1) becomes: where ½D 0 is the effective constitutive matrix, {Dr f } T = {Dp f Dp f Dp f 0 0 0} and Dp f is the change in pore fluid pressure. In a coupled THM problem, the effect of temperature change on the mechanical behaviour of soils is generally characterised by a thermally induced volumetric change, namely thermal expansion/contraction. To formulate the mechanical governing equation that accounts for thermal behaviour, two assumptions are introduced: (1) The temperature of the soil particles, T s , is assumed to be equal to the temperature of the pore fluid, T f . Therefore, only the overall temperature of the soil, T, is adopted here, which implies an instantaneous temperature equilibrium between soil particles and pore fluid; (2) If soil grains are in mineral-to-mineral contact and the temperature is changed, the thermal elastic strain of the soil particles, De T,p , is assumed equal to the thermal elastic strain of the soil skeleton, De T,s . This assumption agrees with the observations of Campanella and Mitchell [12]. Given this equality, in the remainder of the paper the symbol De T is used instead of De T,p or De T,s to represent the effect of temperature on the solid part of the porous medium. Conversely, the incremental thermal strain of the pore fluid is defined independently as De T,f .
Under non-isothermal conditions, the incremental total strain {De} can be expressed as the sum of the incremental strain due to stress change (mechanical strain), {De r }, and the incremental strain due to temperature change (thermal strain), {De T }: fDrg ¼ ½D 0 fDeg þ fDr f g À ½D 0 fm T gðDTÞ ð 4Þ where {m T } T = {a T a T a T 0 0 0}. The last two terms on the right-hand side of Eq. (4) illustrate the change in total stresses induced by the variation in pore fluid pressure and temperature, which represent the coupled effects of the hydraulic and thermal systems on the mechanical system. Applying the principle of minimum potential energy, and minimising the potential energy with respect to the incremental nodal displacements {Dd} nG , results in a finite element equation associated with force equilibrium for a coupled THM problem of fully saturated soils: The matrices in Eq. (5) are defined in the Appendix, where [K G ] represents soil mechanical behaviour, [L G ] is the HM coupling term for pore fluid effect on mechanical response and [M G ] is the TM coupling term for temperature effect on the mechanical response.

Hydraulic equation
Under isothermal conditions, the governing equations for pore fluid flow through the soil skeleton can be established by combining the continuity equation with Darcy's law (e.g. [27]). For a soil saturated with a compressible pore fluid, the continuity equation can be formulated based on the volume conservation of the pore fluid, which implies that the net volume of the pore fluid flowing into and out of a compressible element of fully saturated soil is equivalent to the total volumetric change of the soil skeleton, such that: where {v f } represents the vector of the seepage velocity, rÁ is the symbol of divergence defined as r Á H ¼ @H=@x þ @H=@y þ @H=@z, n is porosity, K f is the bulk modulus of the pore fluid, Q f represents any pore fluid sources and/or sinks, e v is the volumetric strain of the soil skeleton, and t is time. The seepage velocity {v f } in Eq. (6) is considered to be governed by the generalised Darcy's law, expressed as: where [k f ] is the permeability matrix, or hydraulic conductivity, of the soil, p f is the pore fluid pressure, frp f g represents the gradient of pore fluid pressure, the vector {i G } T = {i Gx i Gy i Gz } is the unit vector parallel, but in the opposite direction, to gravity, c f is the specific weight of the pore fluid. Under non-isothermal conditions, however, the changes in volume of the pore fluid due to the temperature change need to be taken into account, as often its coefficient of thermal expansion is different from that of the soil particles (or the soil skeleton, as per earlier assumption). The difference in the two thermal expansion coefficients can generate a volume of pore fluid flowing into or out of the soil element when there is a temperature change. If no flow can enter or leave the soil element, a variation in the pore fluid pressure is induced. Adopting the principle of volume conservation, the equation of continuity for a compressible pore fluid can be extended from Eq. (6) as: where a T,f is the thermal expansion coefficient of the pore fluid. It should be noted that the effect of buoyancy is not taken into account here. Clearly, the first term on the left-hand side of Eq. (8) is related to the flow of pore fluid into and out of the soil element, the second term denotes the changes in the volume of the pore fluid due to its compressibility, while the third term expresses the thermal expansion/contraction of the pore fluid. On the right-hand side of Eq. (8), the first term is the total volumetric strain, while the second term reflects the changes in pore space available for the pore fluid within the porous medium due to the thermal expansion/contraction of the soil skeleton. Eq. (8) can be further rearranged as: where De vT = 3a T DT. It should be noted that when the undrained heating behaviour of soils is modelled by employing Eq. (9) with a very low soil permeability and a relatively small time-step, the difference in thermal expansion coefficients between the soil skeleton and the pore fluid leads to an increase in pore fluid pressure. It should also be noted that applying the principle of mass conservation (the expression for which is shown in the next section) instead of volume conservation, leads to the same formulation as Eq. (6) for the isothermal case, and Eq. (9) for the non-isothermal case. The same hydraulic formulation was also obtained by Lewis and Schrefler [22], who followed the approach which combines the mass balance equation for the solid phase with the mass balance equation for the fluid phase. Transfer of pore liquid as well as pore moisture under temperature gradients was observed in experiments on porous materials [26,11]. For partially saturated soils the effect of temperature gradients on pore liquid flow can be significant; however, for fully saturated soils, the fluid flow due to a temperature gradient is negligible [26]. Therefore, Eq. (7) is adopted here to model the seepage velocity under non-isothermal conditions. Substituting Eq. (7) into Eq. (9) and applying the principle of virtual work as well as the divergence theorem, the finite element equation governing the pore fluid flow under non-isothermal conditions can be assembled in terms of global matrices as [35,27]: Details of the matrices in Eq. (10) are given in the Appendix, where [L G ] and [Z G ] represent the influence of stress-strain and temperature changes on the hydraulic behaviour respectively, and [U G ] and [S G ] are the permeability and compressibility matrices of the pore fluid. To solve this time-dependent equation, a time marching scheme is adopted here, which assumes that for the time interval [t 1 , t], where t 1 is the initial time and t is the current time, where b 1 is the time marching parameter for pore fluid flow, and ({p f } nG ) 1 is the initial pore fluid pressure at the given time interval.
To ensure the stability of the marching process, the value of b 1 should be between 0.5 and 1.0 [9], while the time-step should be chosen carefully to guarantee accuracy (i.e. not excessively large) and to avoid spatial oscillations (i.e. not excessively small), as highlighted in Cui et al. [15]. Substituting Eq. (11) into Eq. (10) yields:

Thermal equation
Heat transfer in soils generally involves three modes [17]: diffusion, advection and radiation. To obtain the equations governing heat transfer for THM coupled problems, the following assumptions have been made: (1) The effect of radiation is assumed negligible [23] and is therefore not taken into account; (2) Following the assumption adopted for the mechanical equation, that the volumetric strain of particles and solid skeleton due to changes in temperature are identical (i.e. De T,p = De T,s = De T ), thermal volumetric changes of the solid part in a free draining, mechanically unrestrained soil element do not change the void ratio. Indeed, conceptually, it is argued that when the particles thermally expand, volumes of both the soil skeleton and the voids change proportionally, resulting in a constant void ratio, such that: and V s0 are the initial volumes of the voids and the soil skeleton, respectively, and DV v,T and DV s,T are the corresponding incremental volumes due to the temperature change. Therefore, in this formulation, only the mechanical volumetric change, which is induced by the change in effective stresses, can affect the void ratio of the porous medium. It should be noted that under undrained heating conditions the change in the void ratio could indeed be caused by changes in pore fluid pressure due to the differential thermal expansion of the pore fluid and the soil skeleton, as discussed in the previous section.
The equation governing ground heat transfer is formulated based on the law of energy conservation [31,32]: where U T is the heat content of the soil per unit volume, Q T is the heat flux per unit volume, which includes heat diffusion and heat advection, Q T represents any heat source and/or sink, and dV is an infinitesimal volume of the soil. The heat content term, U T , is defined proportionally to soil phases, hence for fully saturated soils: where C pf is the specific heat capacity of the pore fluid and therefore associated with the pore space (i.e. porosity, n) and C ps is the specific heat capacity of soil particles, hence associated with the remaining volume (i.e. 1 À n), q f and q s are the densities of the pore fluid and soil particles respectively, and T r is a reference temperature. The heat flux Q T is written as the sum of heat diffusion q d and heat advection q c , where and [k T ] is the thermal conductivity matrix. Substituting Eqs. (14)- (16) into Eq. (13) gives: The first term on the left-hand side of Eq. (17) can be expanded to give: @f½nq f C pf þ ð1 À nÞq s C ps ðT À T r ÞdVg @t where terms (2) and (3) in Eq. (18) represent the influence of deformation of the soil on heat transfer. In effect, these two terms denote the influence of changes in the heat capacity of the porous medium, resulting from variations in porosity which arise from mechanical deformation. Naturally, the overall impact of this part of the formulation is essentially controlled by the difference in heat capacity of the fluid and solid phases of the medium. The second term on the left-hand side of Eq. (17), which represents the advective heat transfer, can be derived as: Applying the principle of mass conservation for the pore fluid (omitting the source term) and soil particles, gives [22,32]: where {v s } is the velocity of the soil particles. Clearly, using Eq. (20), it is possible to establish that term (3) in Eq. (18) and term (2) in Eq. (19) are equal in magnitude but have opposite signs. Therefore, eliminating those terms from Eq. (17) and neglecting the advective heat flux in the solid phase (the second term in Eq. (21)), Eq. (17), which governs heat transfer, can be expressed as: ½nq f C pf þ ð1 À nÞq s C ps @T @t þ q f C pf fv f g T frTg À r Á ð½k T frTgÞ À Q T ¼ 0 It is noted that Eq. (22) reproduces that proposed by Lewis and Schrefler [22] and in a fully coupled THM FE analysis it may adequately illustrate the coupled behaviour of soils. However, when a coupled thermo-mechanical (TM) FE analysis (e.g. an undrained analysis, typical of low permeability clays, where the two-phase nature of the soil is ignored and only the behaviour of the mixture is considered) is carried out, cancelling the two terms which are equal in magnitude from Eq. (17) leads to the absence of TM coupling in the heat transfer equation. Therefore, the formulation proposed here retains all the terms in Eq. (17). Introducing the void ratio, e, as a soil property that is measured in experiments, n = e/(1 + e), and noting that dV = (1 + e)dV s , where dV s is the infinitesimal volume of the soil particles, into Eq. (17) yields: @f½eq f C pf þ q s C ps ðT À T r ÞdV s g @t þ fr Á ½q f C pf fv f gðT À T r Þ À r Á ð½k T frTgÞ À Q T gð1 þ eÞdV s ¼ 0 Under isothermal conditions, dV s is generally assumed to be constant in a simulation, regardless of the change in effective stresses. Under non-isothermal conditions, however, based on the assumptions presented earlier, dV s is temperature dependent and can be written as: where dV s0 is the initial volume of the soil particles. It is noted that dV s0 is assumed to be constant here, which is different from Thomas et al. [32] who assume that dV s is constant in a coupled THM problem. Substituting Eq. (24) into Eq. (23) and eliminating dV s0 leads to: Noting that variables in Eq. (25) are e, T, e vT , q f and q s , and fv f g, the first term on the left-hand side of Eq. (25) can be further derived as: @f½eq f C pf þ q s C ps ðT À T r Þð1 þ e vT Þg @t 1 ð1 þ eÞð1 þ e vT Þ ¼ ½nq f C pf þ ð1 À nÞq s C ps @T @t ð1Þ þ q f C pf ðT À T r Þ 1 1 þ e @e @t ð2Þ þ ½nq f C pf þ ð1 À nÞq s C ps ðT À T r Þ 1 þ e vT @e vT @t ð3Þ þ nC pf ðT À T r Þ @q f @t þ ð1 À nÞC ps ðT À T r Þ @q s @t ð4Þ ð26Þ In this equation, term (2) denotes the influence of the mechanical deformation, through the variation of void ratio, on the simulated heat transfer. Clearly, the contribution of this term could potentially be significant in applications where considerable stress changes are expected, such as in the design of thermo-active structures. If the material density is assumed to be temperature dependent, substituting q = q 0 /(1 + e vT ) for the pore fluid and solid phases, into Eq. (26) leads to the elimination of terms (3) and (4) from Eq. (26), resulting in @f½eq f C pf þ q s C ps ðT À T r Þð1 þ e vT Þg @t 1 ð1 þ eÞð1 þ e vT Þ ¼ ½nq f C pf þ ð1 À nÞq s C ps @T @t ð1Þ þ q f C pf ðT À T r Þ However, if this is done, then a variable material density should also be applied to the development of the mechanical and hydraulic equations and be updated during the iteration of each increment to ensure the consistency of this approach. As far as the authors are aware, this complication has not been fully taken into account in any of the existing published THM formulations using the finite element method. Therefore, constant densities are assumed in the formulation proposed here so that comparisons can be made with results from the literature. As a result, Eq. (26) becomes: @f½eq f C pf þ q s C ps ðT À T r Þð1 þ e vT Þg @t 1 ð1 þ eÞð1 þ e vT Þ ¼ ½nq f C pf þ ð1 À nÞq s C ps @T @t Since the changes in void ratio are assumed to be only a result of the mechanical strain, it is possible to establish the following relationship in a small displacement analysis: where e 0 is the initial void ratio of the soil. If the fluid density is assumed not to change (r Á q f ¼ 0), the second term on the lefthand side of Eq. (25) can be derived as: Substituting Eqs. (7), (28)-(30) into Eq. (25) and following a procedure similar to that for the hydraulic FE formulation, gives the following finite element equations associated with heat transfer: where b 2 and a 1 are time marching parameters for heat transfer.
Note that the same concerns on the adequate selection of the time-step previously raised for the solution of the hydraulic equation apply to the solution of the thermal equation (Eq. (32)). Indeed, Cui et al. [15] demonstrates that for problems involving heat flux, the minimum time-step required to prevent the occurrence of spatial oscillations is highly affected by the chosen discretisation, the applied boundary conditions and by the relative weight of advection on the overall rate of heat transfer. Assembling Eqs. (5), (12) and (32) gives the fully coupled THM formulation in the incremental matrix form as:

Boundary conditions
To adequately model a coupled THM problem with the FE method, appropriate mechanical, hydraulic and thermal boundary conditions (BC) are required. Details of the mechanical and hydrau-lic BC for geotechnical engineering can be found in Potts and Zdravković [27], while the thermal BC, which are commonly used and have also been implemented into ICFEP, involve: (i) Prescribed temperature; (ii) Prescribed heat flux; (iii) Heat source and sink; (iv) Natural heat loss BC which enables heat flux exchange between a body and its surrounding fluid due to temperature difference between the two (for example between the ground and the atmosphere).
In particular, a new coupled thermo-hydraulic BC has been developed and implemented into ICFEP to balance the change in energy associated with the flow of pore fluid through the boundary in a coupled THM or TH analysis [14].
It should be noted that both the natural heat loss BC and the coupled thermo-hydraulic BC are non-linear, as the variables employed by these two BC can vary over an increment of the analysis.

Validation and application
To demonstrate the performance of the newly developed and implemented THM formulation, a series of benchmark analyses have been performed and some are discussed here. In the first instance ICFEP predictions are compared against analytical solutions and numerical studies found in the literature. Subsequently, a comparison is made against laboratory experimental results.

Thermo-elastic consolidation of a soil column
The 1D thermo-elastic consolidation analysis of a fully saturated column of soil was first presented by Aboustit et al. [1,2], and was subsequently used as a benchmark by Noorishad et al. [24], Lewis et al. [21], Lewis and Schrefler [22], and Gatmiri and Delage [18] for checking their FE formulations. Consequently, this example is also analysed using ICFEP.
To ensure the consistency between material properties adopted in the literature and those used in ICFEP, a plane strain analysis of 1D isothermal elastic consolidation was first performed, for which an analytical solution [7] is available. The adopted FE mesh is shown in Fig. 1, employing 8-noded quadrilateral elements, with both displacement and pore pressure degrees of freedom at all element nodes. The adopted material properties are shown in Table 1.
A vertical surface pressure of 1 Pa was applied at the top boundary of the soil column at the start of the analysis and displacement boundary conditions as shown in the figure were prescribed to ensure 1D vertical deformation of the column. The prescribed hydraulic boundary conditions allowed for drainage only through the top boundary and the initial pore water pressure is hydrostatic from this boundary. Two sets of time-steps (see Table 2) were investigated and b 1 = 0.875 was adopted for both sets as the time marching parameter in Eq. (11). Set I was the same as in Aboustit et al. [1,2] and Lewis et al. [21] and it resulted in spatial oscillations in the pore pressure at the beginning of the analysis (both in the above references and in the ICFEP analysis). The authors' earlier studies have shown that these oscillations are induced by an excessively small time-step. In this respect, Cui et al. [15] derived the minimum time-step that ensures non-oscillatory results: where h is a measure of element size. Substituting the values given for the parameters in Table 1 and with h = 0.2 m, gives a minimum time step of 0.056 s. Clearly, the initial time step of 0.01 s in the first analysis violates this constraint and explains the observed oscillations. Fig. 2 shows the evolution of the column surface settlement with time, due to the applied pressure, for both analyses. Analytical solutions from Biot [7] are also plotted. It is noted that both analyses produce very close transient and steady state results, with the second one matching the analytical solutions exactly. The latter result is then compared in Fig. 3 to those shown in the literature. A very good match is evident with the results of Lewis et al. [21] and Noorishad et al. [24], while the formulation of Gatmiri and Delage [18] differs significantly from all others.
The subsequent exercise is a non-isothermal consolidation analysis, which uses the same 2D mesh in Fig. 1, but with an additional temperature degree of freedom at all element nodes. The material parameters are listed in Table 1 and the initial temperature of the soil column is 0°C. The mechanical and hydraulic boundary conditions stayed unchanged, while a constant temperature of 50°C and a vertical surface pressure of 1 Pa are prescribed at the top surface. No heat flow is prescribed at any of the other boundaries (i.e. insulated boundary). Both time-step sets in Table 2 were attempted and Fig. 4 shows that the biggest difference between the two results is in the predicted maximum settlement. Further inspection showed that the Set I analysis missed the peak settlements due to a large time-step Dt = 1000 s between increments 41 (t = 1000 s) and 42 (t = 2000 s). Additionally, oscillations in the nodal temperature distribution along the soil column were observed in this analysis, as the initial time step is lower than the minimum one established by Cui et al. [15]. Therefore, the results from the new THM formulation (as implemented in ICFEP) with the Set II time steps are compared to those from the same references used in the HM analysis. It should be noted that Lewis et al. [21] adopted a Petrov-Galerkin FE method and upwinding techniques to avoid oscillations as the effect of advection is taken into account here. The authors have recently implemented the Petrov-Galerkin FE method into ICFEP for highly advective flows. However, their previous work [14] has shown that even the conventional Galerkin FE method can accurately simulate this problem, depending on the magnitude of the maximum Péclet number. This parameter is used to represent the ratio between the advective and conductive heat fluxes, which in this analysis is 0.002 and therefore below the critical value of 2.0 established by Cui et al. [14]. Consequently, the Galerkin FE method is adopted here in the THM analysis of the soil column. Fig. 5 compares the prediction of the variation in both temperature and excess pore fluid pressure at a depth of 1 m (point 5 in the mesh, see Fig. 1) from the new THM formulation against the    available data from Lewis et al. [21]. Good matches can be found for both temperature and pore fluid pressure. The slight difference in the temperature variation at the initial stage of the analysis is due to the different time marching schemes adopted for heat transfer in the two analyses (i.e. Set I by Lewis et al. [21] and Set II by the authors). Fig. 6 compares the prediction of the surface settlement with time from the new THM formulation against those shown in the literature. The solid line result is from the analysis with the new THM formulation (as implemented in ICFEP) in which a T = a T,f , as adopted in the analyses from the literature. As noted by Lewis et al. [21], there was probably a printing error in the results of Aboustit et al. [2] which has been corrected in Fig. 6 and shows very good agreement with the ICFEP results, although there are limited discrete results for the former publication. Also, the same steady-state was achieved by Lewis et al. [21], Noorishad et al. [24] and the current work. The difference in the peak settlement between the prediction of the new formulation and that given by Lewis et al. [21] is mainly due to the fact that the latter formulation assumes no coupling effect of temperature on the hydromechanical response. In terms of the THM governing Eq. (33) this is equivalent to the coupling matrix [Z G ] being equal to zero. In the new formulation, however, this is not the case and the [Z G ] matrix depends on the magnitudes of a T and a T,f ,. A further example of the level of difference in predicted maximum settlements from the new THM formulation is presented by the dashed line, from an analysis with different thermal expansion coefficients (a T = 3 Â 10 À7 as in Table 1, and a T,f = 2.1 Â 10 À6 ). Although differences in excess pore fluid pressure (see Fig. 5(b)) and maximum settlement (see Fig. 6) between the two analyses with the new THM formulation may not be seen as significant (10% on average for both) for this validation exercise, they could be substantial in a complex boundary value problem. Unfortunately, it is difficult to assess the reasons for differences between the results from the new THM formulation and those from Noorishad et al. [24] or Gatmiri and Delage [18], due to limited information on their numerical formulations given in the literature.

Undrained heating test
The next validation exercise involved application of the fully coupled THM formulation proposed in this paper in the analysis of a triaxial undrained heating test reported by Abuel-Naga et al. [3]. The particular interest here is the development of thermallyinduced excess pore water pressures. A fully saturated Soft Bangkok clay was used in the test with an initial temperature of 25°C. The specimen was isotropically consolidated under isothermal conditions to 200 kPa, followed by an unloading which resulted in an overconsolidation ratio (OCR) of 4.0. Subsequently, the sample was heated under undrained conditions to 90°C with increments of 10°C, and the thermally induced pore water pressure was measured after each temperature increment. Laboratory tests of the thermo-mechanical response of Soft Bangkok clay have shown that overconsolidated samples tend to behave elastically [4,5], and consequently such an assumption is adopted here.
A single element axi-symmetric fully coupled THM analysis was carried out with the present THM formulation. The relevant properties of Soft Bangkok clay (K, n o and a T ) were obtained from data provided by Abuel-Naga et al. [6,4], whereas the thermal expansion coefficient of water (a T,f ) was adopted from Abuel-Naga et al. [3] 41 42   which represents an average value within the temperature range.
All key material properties are summarised in Table 3. All boundaries of the finite element were modelled as impermeable so that an undrained condition was enforced. A total temperature increase of 65°C was applied to all element nodes with an incremental temperature of 1°C. A value of 0.8 was adopted for all time-marching parameters and the time-step was chosen arbitrarily as 1.0 s. The changes in pore fluid pressure were monitored in the analysis and plotted in Fig. 7 against the experimental results. Considering the inevitable experimental scatter, the agreement between the two is excellent. It is the coupling matrix [Z G ] in the proposed formulation that enables the prediction of excess pore water pressures in an undrained heating test, without resorting to a constitutive model specially developed to simulate such an effect.

Conclusions
This paper presents a complete finite element formulation accounting for fully coupled thermo-hydro-mechanical behaviour of saturated soils. While the formulation is independent and based on a set of clearly stated assumptions, the final equations are similar to, but differ in some important aspects from those in the literature. The following are the key characteristics of the proposed formulation: (1) The formulation is fully coupled and all terms, including the cross-coupling terms, are independently defined. (2) Full derivation of the coupled THM formulation identifies cross-coupling terms which can be cancelled out in the heat transfer equation. This has been exploited by others in the literature. However, the proposed formulation shows that this results in an error in a thermo-mechanical analysis.
Consequently, it retains all the terms. (3) The difference in the thermal expansion coefficient between the pore water and the soil skeleton is explicitly accounted for. This enables simulation of undrained heating which results in the generation of excess pore water pressures without having to use a specially developed thermo-plastic constitutive model to simulate this effect. Comparison of predictions using this facility against experimental results shows excellent agreement.
Z Vol ½nq f C pf þ ð1 À nÞq s C ps 1 þ 3a T ðT À T r Þ 1 þ e vT ½N T T ½N T À q f C pf ðT À T r Þ Z Vol q f C pf ðT À T r Þ½E T T ½k f fi G gdVol i fDF G g ¼ fÀ½n G þ Q f þ ½U G ðfp f g nG Þ 1 gDt fDH G g ¼ fÀ½n T G þ Q T þ ½X G ðfp f g nG Þ 1 À ½C G ðfTg nG Þ 1 gDt ½E p ¼ @N P @x ; @N P @y ; @N P @z T ½E T ¼ @N T @x ; @N T @y ; @N T @z T