Numerical analysis of thermal impact on hydro-mechanical properties of clay

Xuerui Wng,Hu Sho,Jürgen Hesser,Chunling Zhng,Wenqing Wng,Olf Kolitzb,

aFederal Institute for Geosciences and Natural Resources-BGR,Hanover,Germany

bTechnical University of Dresden,Dresden,Germany

cGesellschaft für Anlagen-und Reaktorsicherheit-GRS,Braunschweig,Germany

dHelmholtz Centre for Environmental Research-UFZ,Leipzig,Germany

1.Introduction

Highly consolidated clay formations are investigated for deep geological disposal of heat-emitting high-level radioactive waste(HLW)due to their favourable properties such as extremely low hydraulic permeability,predominant diffusive mass transport,good isolation capability,homogeneous structure and especially high sorption capacity for most radionuclides,as well as the ability to seal cracks and fissures by swelling.However,heat output from HLW will dissipate through the surrounding rocks and induce complex thermo-hydro-mechanical-chemical(THMC)processes.Because of the very low permeability and the water-saturated state in clayey rocks,the thermal responses of clay formations are significant such as thermally induced pore pressure changes and therewith swelling and shrinkage behaviours,expansion and contraction due to temperature change,thermally induced deformation and change of strength.In order to investigate the thermal impact on clay host rocks,a large number of in situ heating experiments have been performed in underground research laboratories(URLs),i.e.the HE-D experiment on the Opalinus Clay in the Mont Terri Rock Laboratory in Switzerland(Wileveau,2005;Zhang et al.,2009),the TER experiment on the Callovo-Oxfordian(COX)Clay at the Bure site in France(Wileveau and Su,2007),and the ATLAS experiment on the Boom Clay in the rock laboratory Mol in Belgium(De Bruyn and Labat,2002).Numerous in situ experiments have been carried out under defined conditions to estimate the material-specific parameters.However,there are still challenges to answer the following questions:What is the long-term thermohydro-mechanical(THM)behaviour of clay host rock in a realistic setting?Which possible changes of the material properties,such as permeability,porosity,stiffness and strength,could occur under the long-term heating and cooling processes?Is it necessary to consider the effects of gas thermal expansion,gas flow,and the phase change between gas and liquid in the pore space?What could be the consequences if micro-cracks or fractures formed?How could the integrity of the host rocks be affected by the longterm thermal impact?

The investigation of the above-mentioned issues requires a realistic numerical model which is calibrated by experimental data.This paper makes a contribution to describing the development of a coupled THM model to simulate the thermal influences on the hydro-mechanical processes.Generally,it is difficult to achieve a fully saturated state under natural or experimental conditions.Therefore,it is assumed that gas and waterexist in the porespace at the same time.Both the gas and liquid phases could expand during heating;moreover,the liquid can be vaporised inside a gas-tight system.Thus,a multiphase flow model has been applied to describe the hydraulic processes for gas and liquid as well as the phase change.This model is based on associated theoretical formulations and the multiphase flow model in the finite element programme OpenGeoSys(OGS)(Kolditz et al.,2012a).According to the special properties of clay formations and the experimental conditions,some modifications have already been made in the OGS code.The developed numerical code has been applied to simulate a laboratory heating experiment on the Callovo-Oxfordian(COX)claystone which has been performed by GRS(Zhang et al.,2010).This modelling can reasonably reflect the observed responses of the claystone to the thermal loading.

Fig.1.Interactions between THM processes in porous media(Zhang et al.,2010).

2.Physical processes and numerical approaches

Generally,thermal effects can induce a complex interaction between THM processes.For instance(see Fig.1),thermal loading can induce expansion of pore fluids and solid skeleton,and then change the pore space,leading to increase in pore pressure.Furthermore,the increase of temperature can cause a decrease in gas and liquid viscosity(T?H).It can also induce a deformation with stress variations(T?M).On the other hand,the change of pore pressure has an effect on the effective stress(H?M)and the degree of water saturation.The water/gas flow can influence the heat conductivity(H?T)because of the saturation-dependent thermal conductivity(Gens et al.,2007;Ghabezloo,2010;Zhang et al.,2010).All the above-mentioned couplings have been considered in our numerical modelling.

To analyse the THM coupling processes,the porous medium is usually assumed as a homogenous continuum(Fig.2).Generally,such porous media are composed of three species:mineral,water and gas,distributed as three phases:solid,liquid and gas(Kolditz et al.,2012b).The fluid phase contains water and dissolved air,while the gas phase is a mixture of dry air and water vapour.To present the behaviour of all three phases,a multiphase flow model is required for numerical analysis.

Fig.2.Porous system with two-phase flow(Kolditz et al.,2012b).

2.1.Balance equations

Mathematical descriptions of the physical coupled THM processes for saturated porous media have been proposed by several researchers(e.g.Booker and Savvidou,1985;Olivella et al.,1994).Generally,for the solution of a coupled THM problem,a set of balance equations for internal energy,solid(s)mass,water(w)mass,gas(g)mass,and stress equilibrium have to be solved.

(1)Solid mass balance:

(2)Water mass balance:

(3)Gas mass balance:

(4)Internal energy balance:

(5)Stress equilibrium:

In Eqs.(1)-(5),ρis the density;φis the porosity;J→is the total mass flux;q is the external mass supply per unit volume of medium;Swis the degree of water saturation;E is the specific internal energy;icis the conductive heat flux;JEis the energy flux due to mass motion;σis the total stress tensor;b is the body forces vector;and subscripts “s”,“w”,and “g”stand for the solid,liquid and gas phases,respectively.

2.2.Thermal model

Heat transport in the porous media is governed by conduction and advection flow of water,gas and vapour(Garitte et al.,2014).A linear heat transport process is assumed to describe the thermal process.Using Fourier’s law to describe the thermal conduction,the heat balance equation can be expressed as the simplification of the energy balance equation(Eq.(4)):

where?T is the temperature gradient,c is the heat specific capacity,v is the advection velocity,andλis the thermal conductivity which depends on water saturation and is determined by the conductivity of solid associated with the conductivities of water and gas:

An important thermal effect on hydraulic and mechanical fields is the thermal expansion of each phase(Shui et al.,2010).Generally,it can be described by a thermal expansion coefficient which is defined under undrained condition as(Ghabezloo and Sulem,2010):where pfis the pore pressure andσcis the confining stress.

Different materials have different thermal expansion coefficients.Fig.3 illustrates the thermal expansion coefficient for clay minerals,water and air as a function of temperature.

In case of an unsaturated condition,the thermal expansion of a porous medium is governed by the expansion of the solid minerals,pore water,vapour and dry air(Ghabezloo and Sulem,2010).The expansion coefficient of pore water is almost two orders of magnitude higher than the solid grains.Although air shows the largest thermal expansion,it is much more compressible compared with water.Therefore,the thermal expansion of unsaturated porous medium is dominated by the water expansion.To represent the dependency of thermal expansion of water on temperature,a curve(Fig.3)is defined in the OGS code,which is interpolated from the observed data.

Moreover,the thermal expansion of each phase has already been implemented in the extended mass balance equations(Eqs.(1)-(3)).

Fig.3.Thermal expansion coefficient of water(w),air(g)and clay mineral(s)as a function of temperature(Zhang et al.,2010;Garitte,2013).

With combination and extension of Eqs.(1)and(2),the porous medium mass balance for the solid skeleton,water and vapour is given by where pcis the capillary pressure.

Moreover,combining Eq.(1)with Eq.(3),the extended mass balance equation for solid skeleton and gas can be written as where pgis the gas pressure.

In Eqs.(9)and(10),αis the thermal expansion coefficient for each phase;“ga”and “gw”stand for dry gas and water vapour,respectively;M is the molar mass of each phase;u is the displacement vector of solid matrix;Dgwgis the effective diffusivity tensor of water vapour in dry air;Dgagis the effective diffusivity tensor of dry air in water vapour;and a is the Biot coefficient.

2.3.Hydraulic model

As is already mentioned above,with consideration of the existing gas phase in the pore space and the phase change between liquid and vapour, the multiphase flow model (Kolditz and De Jonge, 2004; Wang et al., 2010) is modified to describe the hydraulic process. The total hydraulic flux Jγkγcontains two parts: the advective part JγAKand the diffusion part JγDK:

where the indexγrefers to the different phases,i.e.liquid(l)and gas(g);the subscript“k”denotes the components of fluid:dry gas(k=a)and water(k=w).

The advective part is based on Darcy’s law and includes water flux JAlkas well as gas flux JAgk:

where kintis the intrinsic permeability,kris the relative permeability,g is the gravity acceleration,andμis the viscosity.Moreover,the water viscosity is a function of temperature and can be given by(Garitte,2013):

whereAandBare the shape parameters,A=2.1×10-6Pa s,andB=1808.5 K.

In the numerical calculation,gas pressurepgand capillary pressurepcare considered as primary variables instead of thewater pressurepwin the single phase flow.Water pressure can be then calculated as the difference between the gas pressure and capillary pressure:

The water retention is introduced by the van Genuchten function:

wherep0is the constant air entry pressure,Seffis the effective degree of water saturation,mis the pore size distribution index,andn=1/(1-m).

Moreover,the dependency of the relative permeability for gas and water on the effective degree of water saturation can be written as

where Dγkis the diffusion tensor.

To describe the phase change,only the vaporisation will be considered,however,the condensation processes are neglected(Rutqvist et al.,2002).Based on these assumptions,diffusion is only observed in the binary gas phase.So the corresponding diffusion part can be obtained as follows(Lewis and Schre fl er,1998):

wheredenotes the partial pressure induced by phase change;subscripts “a”, “w”,and “g”stand for dry gas,water,and gas mixture,respectively.The molar mass of gas mixtureMgcan then be written as

where ρgis the density of gas mixture,refers to the mass concentration of dry gas component,is the mass fromwater vapour component in gas phase.

using the physically appropriate assumption of

Furthermore,according to the Clapeyron equation for an ideal gas and Dalton’s law,the following equations can be obtained:

whereRis the universal gas constant with avalue of 8.314 J/(mol K).

In the unsaturated zone,the equilibrium water vapour pressure is given by the Kelvin-Laplace equation:

wherepgwsis the saturated vapour pressure,as given in Philip and De Vries(1957):

In order to enhance the numerical stability,the gas phase density related derivative terms in the mass balance equations(Eqs.(9)and(10))are expanded by the chain rule:

Eq.(26)represents the contribution of the change of the different variations to the mass transition from liquid phase to gas phase.

2.4.Mechanical model

Although thermal loading may cause plastic deformation even damage to clay sample,taking account of relative low temperature of 90°C,relative low observed pore pressure of 12 MPa in comparison with the con fining stress of 15 MPa,and relative low observed strain,only elastic deformation is speculated in our case.Specific attention is,however,paid to the thermal impact on the hydraulic field.For simplification,a linear elastic mechanical model based on Hooke’s law is used to account forde formation induced by swelling/shrinkage and temperature change:

whereσijis the effective stress;is the elastic stiffness matrix,which can be determined by the Young’s modulusEeand the Poisson’s ratio νin the linear elasticity;d3is the total strain increment;d3sis the strain increment caused by swelling or shrinkage;and d3Tis the thermally induced strain increment.The effective stress can be expressed with the Biot’s effective stress as

whereBis the Biot coefficient,I is the identity tensor,andplis the pore pressure.

It is well known that the change of temperature causes an additional strain in mechanical process.This thermal strain has been considered as an additional strain tensor in the mechanical model,related to the linear thermal expansion coefficientαand the temperature change:

Corresponding to the change of water content,swelling or shrinkage occurs with the change of the distance between solid particles of the pore space(Xie et al.,2004).To calculate the deformation caused by swelling or shrinkage,the numerical approach by Xu et al.(2013a)has been extended.It can be expressed as(Garitte,2013):

where ds=d（pg-pl）is the suction change.The increase of water content(suction decrease)leads to swelling and in contrary the decrease of water content(suction increase)causes shrinkage.Ksis the bulk modulus against suction change.Considering transversal isotropic swelling/shrinkage,a bedding direction dependent bulk modulus tensor with K||sand K⊥shas been implemented.

3.Laboratory heating experiment

To characterise the thermal properties,especially the thermal expansion of COX and therewith induced pore pressure increase,a heating experiment on a COX sample was carried out by GRS(Zhang et al.,2010).Fig.4 illustrates schematically the testing apparatus and conditions.The sample was mounted on a triaxial apparatus.During heating,axial deformation was recorded by a deformation meter mounted inside the cell between the upper and lower loading platens.A circumferential extensometer chain was mounted around the sample outside jacket at its mid-height to determine the radial deformation.The sample was loaded with an axial stress of 15.5 MPa and a confining stress of 15 MPa.Before heating synthetic COX,pore water was injected into the sample through the upper and lower porous discs until the sample was saturated.After that,the inlet and outlet valves were closed.Under the undrained conditions,heating started using an electrical heater positioned near the bottom of the sample.The water back-pressure at the bottom and the top was recorded during the experiment.

Fig.5 shows the evolution of confining stress,temperature and water back-pressure at the bottom and the top of sample as a function of test duration.The first water-saturation phase with injecting water to both end faces of the sample lasted for more than two months.The following heating phase was conducted in two steps.In the first heating phase,temperature was increased from 30°C to 60°C with a rate of 3°C/h.Then it remained constant(60°C)for about 3 d.After that,the second phase started with temperature elevation from 60°C to 90°C with a rate of 1.2°C/h.Until the end of the test,the temperature was kept constant at 90°C.

In the first heating step(from 30°C to 60°C),a rapid increase of the pore pressure from the initial value of 1.2 MPa to 4.3 MPa at the top and to 7 MPa at the bottom was measured.It was expected that the pore pressure would keep constant given a constant heat injection.However,the pressure decreased at the elevated temperature level gradually to 2.2 MPa at the top and 1.7 MPa at the bottom.The pore pressure increased rapidly to 6.5 MPa at the top and 11.7 MPa at the bottom in the second heating phase(from 60°C to 90°C).Similar to the first heating phase,the pore pressure decreased continuously when the temperature remained constant at 90°C.It is evident that the sudden pore pressure rise is probably caused by the thermal expansion of the three phases(solid,water and gas).The unexpected decrease of the pore pressure at each elevated temperature indicates that the sample sealed with rubber jacket might not be gas-tight or/and a full water-saturation was not reached,so that the testing condition was not absolutely“undrained”.The subsequent testing(indicated by numbers 3,4,and 5 in Fig.5)shows that the pressure increase at the top did not yield any pressure response at the bottom,suggesting that the sample was actually unsaturated.To prove these assumptions,several numerical calculations have been done to simulate the described laboratory test and to analyse the possible THM processes in the clay sample(Sections 4 and 5).

Fig.4.Triaxial testing apparatus(left)and applied conditions(right)of the heating experiment(Zhang et al.,2010).

4.Numerical simulation and analysis

4.1.Model setup

In order to analyse the laboratory results and to get a better understanding of the induced processes,a coupled non-isothermal multiphase flow mechanical(TH2M)model has been applied.A two-dimensional axisymmetric model with a height of 100 mm in the direction of the rotational axis(y-direction)and a width of 25 mm(x-direction)has been constructed.Fig.6 illustrates the model geometry,the boundary conditions,and the finite element meshes consistingof 212 nodes and 359elements.The bedding is in the direction perpendicular to the rotational axis and represents the transversal isotropy of the clay sample.

As for the initial conditions,the temperature was set to 30°C and the pore pressure was set to 1.2 MPa based on the laboratory experiment(Fig.5).The initial degree of water saturation was assumed to be 98%.This assumption of the unsaturated condition of the sample is based on the observation in the water injection test that no pressure response at the bottom was generated by water injection at the top.Temperature changes were set on the bottom of the model.The lateral boundary was assumed that no heat flux and no flow fluxexist because of good thermal isolation and tightness of the rubber jacket.Displacements in vertical direction are disabled at the bottom and a compressive load of 15.5 MPa was set at the top of the model.A compressive confining stress of 15 MPa was applied on the lateral boundary.

Fig.5.Temperature,con fining stress and measured pore pressure at the top and the bottom of the sample(Zhang et al.,2010).

4.2.Model case studies

Only the two heating phases were numerically studied.The initial saturation phase of the test sample before laboratory heating test was neglected.As is already mentioned in Section 3,the injection tests during the second heating phase may indicate that the sample is still not fully saturated.The sample was therefore assumed to be quasi-saturated with a water saturation of 98%.Furthermore,the water injection tests are not considered in the numerical simulation.The general numerical exercise thereafter consists of three parts:coupled hydro-mechanical simulation of sample consolidation before the first temperature increase in order to obtain initial hydro-mechanical equilibrium state.After that,the simulations of the two heating phases were done,using coupled thermal multiphase flow mechanical model.The temperature conditions for the simulations are listed in Table 1.To analyse the effect of gas phase on the total pressure development,an additional simulation based on Richards’formulation was carried out as the reference of the two-phase flow model.In the Richards’formulation,gas pressure is considered constant and equal to atmospheric pressure.

Fig.6.Sketch of the numerical model:geometry,mesh and the boundary conditions.

4.3.Model parameters

The parameters used in the model are given in Table 2.The parameters for solid density and porosity come from experimental database.The mechanical parameters of the sample(Young’s modulus and Poisson’s ratio)were determined by corresponding tests on other samples(Zhang et al.,2010)which were taken from the same core as the one used for the heating experiment.Even so there are still some uncertainties related to a few parameters because of the insufficient laboratory evidence.Thermal,hydraulic and mechanical transversal isotropies of the bedded clay formation were considered and represented by transversal isotropic values(Table 2).The water and gas densities were assigned with the values at the initial temperature of 30°C and their temperature dependency was considered according to the thermal expansion coefficient. Moreover,the temperature-dependent thermal expansion of water was calculated using a curve,while the thermal expansion of solid and gas was assumed to be constant.

Table 1Phases of numerical analyses.

5.Results and analyses

5.1.Thermal results

In the experiment,only the temperature at the bottom near the electrical heater was monitored.This measured temperature was used as boundary conditions on the bottom of the model.The temperature distribution in the sample is unknown for the lack of measurements.Due to low permeability of clay,temperature transport by advection is negligible,and therefore the only significant heat transport mechanism is heat conduction(Gens et al.,2007).Thus,the temperature distribution in the clay sample is dominated by the thermal conductivity of solid and water which has been proven by using both multiphase flow model and single phase model.Fig.7 illustrates the distribution of the calculated temperature at the end of the second heating phase.In Fig.7,a temperature gradient from the bottom to the top can be seen,due to the smaller thermal conductivity in the perpendicular to bedding direction than that in the parallel direction(1 W m-1K-1).The temperature compensation in the lateral direction is more than two times faster than that in the axial direction because of the higher thermal conductivity in this direction.The significant temperature gradient has a certain impact on the pore pressure distribution in the sample.

5.2.Hydraulic results

Focus of the present study is to understand the pore pressure evolution under thermal loadings in the experiment.The observed behaviour of pore pressure in the experiment will be numerically analysed with consideration of possible influential factors.

The pore pressure increase caused by the thermal loading can be well simulated if the thermal expansion of fluids and solids is considered.However,the unexpected decrease of the pore pressure during the phase of constant temperature cannot be simulated under no- flux boundary conditions.Theoretically in an undrained system,a state of equilibrium for pore pressure will appear in the whole system when there is no temperature change.This viewpoint has also been improved through the first simulation using the basic model(undrained condition),as shown in Fig.8.In the first calculation,only the first heating phase was considered with a simulated heating duration of 60 h.As shown in the numerical results after the maximal temperature was reached,the pressure at the top of the sample increased continuously till the same pressure at the bottom.While the temperature was kept constant,the pore pressure remained unchanged.However this was in contrast to the experimental findings,where a decrease of pore pressure was measured when the temperature remained constant.Therefore,additional factors impact on the pore pressure behaviour should be applied to the basic model.

Table 2Model parameters.

5.2.1.Leakage effects

In Fig.5,a significant pore pressure decrease in the phases of constant temperature can be seen.There are two possible reasons for this laboratory finding.One could be the thermally induced extension of the pore space.An analytical calculation was firstly carried out in order to verify if the pore pressure decrease is significantly induced by pore space thermal expansion.Because the thermally induced expansion of pore space cannot be measured in the experiment,the pore space enlargement is equal to the measured expansion of sample volume.Fig.9 shows that the maximal volume expansion is found at the end of the first heating phase with an increment about 0.05%.Moreover,the pore space is regarded as fully saturated.Therefore,the maximal possible pore pressure decrease due to the pore space enlargement can be supposed according to the Tait equation:where Kwis the water bulk modulus,and Bhis the heat capacity ratio.For water under natural condition,Bh≈0.3214 GPa and χ≈7.With these assumptions,a possible maximal pore pressure decrease of 1.12 MPa due to the pore space expansion can be overestimated.However,from the measurement a pore pressure drop about 6 MPa at the bottom and more than 2 MPa at the top after the first heating phase can be observed.Thence it is postulated that the leakage effects should be the main reason of the pore pressure drop.

Fig.7.Calculated temperature distribution at the end of the second heating phase for t=96 h(multiphase flow).

Fig.8.Temperature and pore pressure evolutions during the first heating phase(base case model).

Fig.9.The measured volume strain compared with the measured pore pressure.

Fig.10.Distribution of pore pressure(a)and water saturation(b)at the end of the first heating phase(t=71 h).

According to the description of the experiment equipment,a possible leakage is assumed at the edge of both upper and lower surfaces of the cylindrical sample.To simulate the possible leakage effect on the rapid pore pressure decrease,atmospheric pressure was imposed on this edge as an additional boundary condition.That means the simulated system is not gas-tight,especially the gas and vapour can flow out of the system through the edge boundary.Fig.10 illustrates the calculated distribution of pore pressure and water saturation in the sample at the end of the first heating phase.These results were calculated with the base case model considering a possible leakage.The area around the leakage position shows a significant lower pore pressure and lower water saturation compared with other areas inside the model.

With this modification,reasonable pore pressure evolutions at the measurement points have been obtained in the model.In the first heating phase,the simulated pore pressure fitted very well with the measured data in trend and in magnitude(Fig.11).However,in the second heating phase,the simulated maximum pore pressure was almost three times higher than the measured one.It can be concluded that with consideration of the leakage points,the trend of pressure evolution can be well captured in the model but not the magnitude in the second heating phase.The difference between the measured and simulated pore pressure development during the second heating phase may indicate that some hydraulic properties of the sample have been changed under the relatively high temperature,e.g.permeability,which has a significant impact on the pore pressure evolution.Under the assumption that the higher temperature can induce a larger pore pressure increase which may cause a pore space enlargement and accordingly the permeability increases,a pore pressure dependent permeability model should be developed to capture the pore pressure behaviour under higher thermal loading.

Fig.11.Temperature and pore pressure evolution during the first heating phase(model with consideration of leakage effect).

5.2.2.Pore pressure dependent permeability model

As mentioned in Section 5.2.1,the permeability may have been changed under thermal loading.For those possible mechanisms,for instance,pore space increase induced by pore pressure increase,fracturing process due to plastic failure should be considered.In consideration of relatively low temperature and high confining stress,the increase of permeability may merely cause the enlargement of pore space induced by the increase of pore pressure.Some researchers(Kozeny,1927;Ghabezloo,2010)have proved that increase of pore pressure may lead to an increase in porosity and therefore to an increase in permeability.The phenomenon has been observed in the laboratory as well as in situ and has been investigated theoretically and numerically(e.g.Walls,1982;Li et al.,2009;Xu et al.,2011,2013b).Xu et al.(2011)had proposed a linear relationship between permeability change and gas pressure in the pore space,with the following two assumptions:(1)Relatively low pore pressure can only cause extension or compaction of pore space that has limited effects on permeability;(2)If gas pressure oversteps a threshold value(e.g.the minimal principal stress),micro-or macro-fractures can be generated which cause a significant permeability increase.The permeability is assumed to increase slightly if gas pressure is lower than this threshold.Otherwise the permeability increases rapidly after the gas pressure exceeds the threshold.Based on the permeability model proposed by Xu et al.(2011),a pore pressure dependent permeability model(Eq.(32))was modified and implemented in this study with the assumption that the permeability begins to increase if the pore pressure is higher than a critical value.The increase of permeability is linearly dependent on the actual pore pressure with a factorβ(Fig.12).The factorβ(Eq.(32))and the critical pore pressure pcriare empirical parameters assumed in this study.Using the trial-and-error method,the critical pore pressure pcriis set to 7.2 MPa and the permeability enhancement factor β is given in Fig.12 in this study.

The actual permeability in the model can be calculated using the empirical pore pressure dependent permeability model:where β is the permeability enhancement factor and determined by a linear function as shown in Fig.12.

Fig.12.Permeability enhancement factor as a function of pore pressure.

The validity of the proposed model is still not proven for a general application because there is still less of evidence on the micro-fracture in case of dilatancy-controlled flow regime from the laboratory experiment.Some further investigations should be carried out to estimate the behaviour of pressure-controlled permeability in such situation.

5.2.3.Effects of gas phase

The amount of gas in the above-mentioned experiment is very low.Concerning temperature load,there is also only a small amount of fluid phase that could contribute to a vaporisation process.Therefore,it is justified to state the question for the necessary of consideration of the gas effects in the hydraulic model:How does the presence of gas affect the THM interactions?

To demonstrate the effect of the gas phase,a reference simulation with Richards’ flow model has also been carried out.The difference between the Richards’ flow model and multiphase flow model is the handling of gas pressure.In the Richards’model,gas pressure is assumed to be constant as atmospheric pressure with a value of 0.1 MPa.However,in the multiphase flow model,both the liquid and gas phases have a defined volume fraction and velocity field.Moreover,phase change between gas and liquid can take place under certain conditions.Fig.13 illustrates the comparison of the calculated results using the Richards’ flow model with the measured results.In this model the leakage effects as well as the pore pressure dependent permeability model have been considered.The numerical results show that if the gas effects(e.g.thermal expansion of gas,gas flow,change of gas pressure)are neglected,the computed increase of pore pressure by temperature increase and the decrease of pore pressure at the constant temperature are more rapid than the measured results.From the results we can conclude that gas has a so-called “buffer”effect on pore pressure evolution.There are two possible reasons for the better matching of the observed pore pressure evolution using the multiphase flow model(Fig.14)than using Richards’ flow model(Fig.13).On one hand,the thermal expansion of gas in comparison with water is small(Eq.(10)),but it can affect the increase of pore pressure in both magnitude and trend.On the other hand,the phase change from liquid phase to gas phase can equilibrate the change of liquid field(Eqs.(9)and(26)).Therefore,the computed change of pore pressure using the multiphase flow model taking the gas pressure into consideration is more reasonable than a single phase Richards’flow model.It is verified that the gas effects should also be considered in the model.

Fig.13.Comparison of the computed pore pressure using the Richards’ flow model(single phase)with the measured data.

Fig.14.Pore pressure evolution at the bottom and top of the sample with the applied temperature.

5.2.4.Results of the multiphase model

Based on the above-mentioned observations from the preliminary calculations,a simulation using multiphase flow model has been carried out,in which gas effects,leakage effects and pore pressure dependent permeability evolution have been considered.The comparison between the calculated pore pressure at the bottom and the top with the measured data is depicted in Fig.14.The numerical results agree fairly well with the experimental data,especially in the first heating phase.Both the simulated and the measured results show that in the first heating phase,the effect of thermal expansion on solid and water is dominant.Therefore,the pore pressure increases to a maximal value of 7.2 MPaat the bottom and 4.3 MPa at the top in the first heating phase.Afterwards,under constant temperature conditions,a pore pressure decrease to 2 MPa is observed due to the leakage effect in the system,even when the pressure is still higher than the atmospheric pressure.In the second heating phase,a much steeper pressure increase to a maximal value of 11 MPa at the bottom and 6.7 MPa at the top was measured and simulated.Unlike the behaviour in the first heating phase,the pressure dropped down more rapidly in the case of constant temperature in the second heating phase.Keeping the leakage properties unchanged,the rapid pressure decrease cannot be simulated by the same model(Fig.11).Due to high pore pressure,possible pore space expansion may happen.This additional expansion of pore space can be simulated using the abovementioned pressure dependent permeability model.The simulated results indicate that the drop of pressure is not only due to the leakage but also due to the increase of permeability.The simulated decrease of pressure is considered as the results of the combined effects of leakage and permeability increase.At the time of 100 h,an additional small increase of temperature from 81°C to 89°C took place.However,this increase of temperature induced merely a small pore pressure rise at the bottom both simulated and measured.It can be seen that the main feature of the observed pore pressure evolution has been captured in the modelling.During the second heating phase,the calculated pore pressure is generally higher than the measured value as the temperature was kept unchanged.The reason for this may be that the pore pressure induced permeability enlargement could immediately get back in the model(Eq.(32))when the pore pressure drops down.In reality,however,the permeability response to the decrease in pore pressure may have a retardation effect.

5.3.Mechanical results

The changes of strains and stresses as well as material properties result from the combining effects of temperature change,pore pressure generation and dissipation.The evolution of the calculated radial and axial strains associated with the measured strain and temperature is shown in Fig.15(positive value means extension,negative refers to compression).Only relatively small strain due to thermal expansion was observed for both measured and calculated results because of the relative high confining stress of 15 MPa.Moreover,the calculated strain evolution shows a similar tendency compared with the observed strain behaviour.As the temperature increases,the sample extends in the radial direction mainly due to the thermal expansion.However,when the temperature remains constant,a compaction effect of the sample occurs in all directions.The compaction results from the leakage induced decrease of pore pressure and the shrinkage effects.It can be also seen that while radial direction extends,compaction prevails in axial direction.The difference in strain behaviour between radial and axial directions reflects the strong transversal isotropic properties of clay considered by the orientation dependent stiffness in the direction and thermal properties.To be specific,the sample has a smaller thermal conductivity,a weaker shrinkage/swelling behaviour and a smaller stiffness in the direction perpendicular to bedding(axial direction)than that in the direction parallel to bedding(radial direction).

In Fig.15,the measured data from the laboratory experiment and the calculated results are depicted.In the first heating phase,the calculated results match the measured data quite well,especially the axial strain.However,a major difference between numerical and experimental results can be observed at the beginning of the second heating phase.This difference might be caused by different pore pressure evolutions between the numerical and observed data(Fig.14).

With the temperature increase,a thermal extension and a pore pressure increase occur.Obviously,both effects lead to an extension in radial direction.Moreover,in the period of constant temperature,leakage effects cause a decrease in pore pressure and also shrinkage behaviour of the sample.Therefore,after the first heating phase,while the temperature is kept constant,the radial extension decreases and the compaction in axial direction increases.In the second heating phase,the radial extension increases again and there is a greater increase in the axial compression.When the temperature remains constant at 90°C during the second heating phase,the measured pore pressure decreases more rapidly.Within this phase of constant temperature,the measured pore pressure reaches much smaller values compared with the calculated results(Fig.14).This smaller pore pressure could be the reason for the more compressive deformation behaviour in the laboratory experiment.Another reason accounting for the difference between the measured and observed axial strains is perhaps due to the implemented mechanical constitutive model. As already mentioned,a pure elastic mechanical constitutive model with additional swelling/shrinkage strain was implemented in the numerical analyses.Indeed,plastic deformation can be observed under such heating experiment.During the second heating phase when the temperature remained constant,the observed compression behaviour of the sample was more rapidly than in the first heating phase.That may indicate a plastic material behaviour during the second heating phase which has caused an irreversible permeability enhancement.

Fig.15.Comparison of the calculated axial and radial strains with the measured values.

6.Conclusions

Thermal effects on the hydro-mechanical properties and therefore the integrity of clay host rocks are very important issues for repositories of heat-emitting radioactive waste in deepgeological formations. Experimental studies have provided a basicdatabase for the characterisation of thermal behaviour and the determination of material parameters. To investigate the long-termbehaviour of clay for a repository of heat-emitting waste, the development of a numerical model is necessary. Therefore, a coupled thermo-mechanical multiphase flow model has been developed and implemented in the finite element programme OGS to analyse the thermal effects on clay behaviour. A two dimensional axis symmetrical model has been constructed for the interpretation and explanation of the observations in the laboratory heating experiment. By stepwise inclusion of influential factors in the numerical model, an inside view of the experimental results can be obtained. Based on the experimental observation, the developed numerical model is able to analyse the coupled THM behaviour. The basic achievements can be summarised as follows:

(1)Generally,the increase of temperature in the porous system can

result in thermal expansion of solid,liquid and gas.Therewith pore pressure will be increased due to thermal expansion.To represent this phenomenon,the specific thermal expansion coefficients of solid,water and gas(Eqs.(9)and(10))have been implemented in the model.Moreover,the thermally induced stress change has also been considered(Eq.(30)).The comparison of the numerical results with the measured data(Fig.8)indicates that the pore pressure increase can be well interpreted when taking the thermal expansion of solid and liquid into account in the model.However,the decrease in pore pressure cannot be captured by applying the no- flux hydraulic boundary(undrained conditions).

(2)During the experiment,it can be observed that at elevated temperatures,the pore pressure decreases gradually.This is also a common phenomenon by the in situ heating experiment in the underground laboratory resulting from the dissipation of the thermally mobilised pore water seepage outward in drained conditions.However,in our study,pore pressure decreases even in undrained experimental conditions.That may suggest that vapour or gas flows out of the experimental system through some leakage area.Using a model with leakage points at the no- flux boundary,the decrease in pore pressure at elevated temperatures can be well represented(Fig.11).

(3)Fig.11 shows that the calculated pore pressure increase in the second heating phase is almost three times the observed pore pressure increase.Therefore,it is assumed that the permeability has been changed under the relatively high temperatures.Special approaches for the interpretation of the experimental observations especially in the second heating phase have been developed.By considering the special experiment conditions and based on researches from other authors,a linear relationship between the pore pressure and permeability(Fig.12)has been developed.In this pore pressure dependent permeability model,when the pore pressure exceeds a critical value,the permeability increases linearly with an enhancement as the pore pressure increases.The critical value and the enhancement factor are determined through back calculation by matching the curve to the measured pore pressure.Only by using the modified permeability model can the observed pore pressure evolution in the second heating stage be numerically interpreted well.

(4)Normally under natural or experimental conditions,gas and water exist in the pore space at the same time.Similar to water,gas phase also expands during thermal loading,which has a certain contribution to pore pressure increase.Simultaneously,the temperature and pressure change can induce a phase change between water and gas.To investigate the gas effects,a reference modelling using single phase flow/Richards’ flow model has also been carried out.The comparison of the numerical results using Richards’ flow with the measured results(Fig.14)indicates that gas expansion,gas flow,and phase change between liquid and gas have a certain impact on the pore pressure evolution in both magnitude and trend.Moreover,it also verifies that the modelling using multiphase flow model in consideration of the gas effects is closer to reality.

(5)The triggers of displacement and strains under thermal loading contain three parts:change of pore pressure,swelling/shrinkage behaviour,and thermal expansion.The swelling/shrinkage behaviour is simulated by adding an additional increment of strain which is caused by the change in water saturation.Because of the relatively high confining stress,only a very small thermally induced strain with a maximum value of 0.05%was measured and computed.Using the current framework of the scientific code,OGS,the laboratory heating experiment on the COX clay stone can be numerically interpreted in most hydraulic and mechanical aspects.However,there are still spaces left in the presented model for improvement.A three dimensional model will be developed for better interpretation of the transversal isotropic effects on a clay formation and to simulate large-scale in situ experiment or real repository.Moreover,a more complex elasto-plastic mechanical constitutive model related to thermal effects,swelling/shrinkage behaviour,and transversal isotropy should be developed to better analyse mechanical behaviour.Future work will also focus on the thermal effects on material properties,e.g.the change in porosity and permeability,anisotropic thermal expansion of solid,the thermally induced change of material stiffness and strength,etc.

Conflict of interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Acknowledgement

This work was supported by BMWi(Bundesministerium für Wirtschaft und Energie,Berlin).We are grateful to Ms Jakobs(UFZ,Germany)for her reviewing of the paper concerning English language.

Booker JR,Savvidou C.Consolidation around a point heat source.International Journal for Numerical and Analytical Methods in Geomechanics 1985;9(2):173-84.

De Bruyn D,Labat S.The second phase of ATLAS:the continuation of a running THM test in the HADES underground research facility at Mol.Engineering Geology 2002;64(2-3):309-16.

Garitte B,Gens A,Vaunat J,Armand G.Thermal conductivity of argillaceous rocks:determination methodology using in situ heating tests.Rock Mechanics and Rock Engineering 2014;47(1):111-29.

Garitte B.FMT FE-experiment:1D benchmark definition for the modelling groups;2013[unpublished].

Gens A,Vaunat J,Garitte B,Wileveau Y.In situ behaviour of a stiff layered clay subject to thermal loading:observations and interpretation.Geotechnique 2007;57(2):207-28.

Ghabezloo S,Sulem J.Effect of the volume of the drainage system on the measurement of undrained thermo-poro-elastic parameters.International Journal of Rock Mechanics and Mining Sciences 2010;47(1):60-8.

Ghabezloo S.Effect of porosity on the thermal expansion coefficient:a discussion on the paper‘Effects of mineral admixtures on the thermal expansion properties of hardened cement paste’by Z.H.Shui,R.Zhang,W.Chen,D.Xuan,Constr.Build.Mater.24(9)(2010)1761-1767.Construction and Building Materials 2010;24(9):1796-8.

Kolditz O,Bauer S,Bilke L,B?ttcher N,Delfs JO,Fischer T,G?rke UJ,Kalbacher T,Kosakowski G,McDermott CI,Park CH,Radu F,Rink K,Shao H,Shao HB,Sun F,Sun YY,Singh AK,Taron J,Walther M,Wang W,Watanabe N,Wu N,Xie M,Xu W,Zehner B.OpenGeoSys:an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical(THM/C)processes in porous media.Environmental Earth Sciences 2012a;67(2):589-99.

Kolditz O,B?ttcher N,G?rke UJ.Theory.In:Kolditz O,G?rke UJ,Shao H,Wang W,editors.Thermo-hydro-mechanical-chemical processes in porous media.Heidelberg:Springer;2012b.p.10-54.

Kolditz O,De Jonge J.Non-isothermal two-phase flow in low-permeable porous media.Computational Mechanics 2004;33(5):345-64.

Kozeny J.Ueber kapillare Leitung des Wassers im Boden.Sitzungsgeber Akad Wissensch Wien 1927;136(2a):271-306(in Germany).

Lewis RW,Schrefler BA.The finite element method in the static and dynamic deformation and consolidation of porous media.2nd ed.Chichester:John Wiley&Sons Ltd.;1998.p.57-65.

Li M,Xiao WL,Guo X,Zhang LH,Zheng LL.Laboratory study of the effective pressure law for permeability of the low-permeability sandstones from the Tabamiao area,Inner Mongolia.Chinese Journal of Geophysics 2009;52(6):1402-13.

Olivella S,Carrera J,Gens A,Alonso EE.Non-isothermal multiphase flow of brine and gas through saline media.Transport in Porous Media 1994;15(3):271-93.Philip JR,De Vries DA.Moisture movement in porous material under temperature gradient.Transactions,American Geophysical Union 1957;38:222-32.

Rutqvist J,B?rgesson L,Chijimastu M,Kobayashi A,Jing L,Nguyen TS,Noorishad J,Tsang CF.Thermohydromechanics of partially saturated geological media-governing equations and formulation of four finite element models.International Journal of Rock Mechanics and Mining Sciences 2002;38(1):105-27.

Shui Z,Zhang R,Chen W,Xuan D.Effects of mineral admixture on the thermal expansion properties of hardened cement paste.Construction and Building Materials 2010;24(9):1761-7.

Walls JD.Effects of pore pressure,con fining pressure and partial saturation on permeability of sandstones. PhD Thesis. Stanford, USA: Stanford University; 1982.Wang WQ,Rutqvist J,G?rke UJ,Birkholzer JT,Kolditz O.Non-isothermal flow in low permeable porous media:a comparison of Richards’and two-phase flow approaches.Environmental Earth Sciences 2010;62(6):1197-207.

Wileveau Y,Su K.In situ thermal experiment carried out in Opalinus clay and Callovo-Oxfordian claystones by Andra-experiment set-up and measurement results.ANDRA report.Paris,France;2007.

Wileveau Y.THM behaviour of host rock:(HE-D)experiment:progress report.Part 1.Technical Report TR 2005-03;2005.

Xie M,Agus SS,Schanz T,Kolditz O.An upscaling method and a numerical analysis of swelling/shrinking processes in a compacted bentonite/sand mixture.International Journal for Numerical and Analytical Methods in Geomechanics 2004;28(15):1479-502.

Xu WJ, Shao H, Hesser J, Wang WQ, Kolditz O, Popp T. Simulation of dilatancy controlled gas migration process in saturated argillaceous rock. In: Proceeding of the 2nd International Symposium on Computational Geomechanics: COMGEOII. Rhodes: International Centre for Computational Engineering; 2011. p. 693-703.Xu WJ, Shao H, Hesser J, Wang WQ, Schuster K, Kolditz O. Coupled multiphase flow and elasto-plastic modelling of in-situ gas injection experiments insaturated claystone (Mont Terri Rock Laboratory). Engineering Geology2013b;157(8):55-68.

Xu WJ,Shao H,Marschall P,Hesser J,Kolditz O.Analysis of flow path around the sealing section HG-A experiment in the Mont Terri Rock Laboratory.Environmental Earth Sciences 2013a;70(7):3363-80.

Zhang CL,Czaikowski O,Rothfuchs T.Thermal-hydro-mechanical behaviour of the Callovo-Oxfordian clay rock.Final report of GRS;2010.p.121-36.

Zhang CL,Wieczorek K,Rothfuchs T,Armand G,Lebon P.Responses of the Opalinus clay to heating during the HE-D experiment at Mont Terri.In:Proceedings of the European Commission TIMODAZ-THERSA International Conference.Luxembourg;2009.