Ho Wng ,Yu-Jun Cui ,Minh Ngoc Vu ,Jen Tlndier
a Ecole des Ponts ParisTech,Laboratoire Navier/CERMES,6 et 8 avenue Blaise Pascal,Marne La Vallée cedex 2,77455,France
b Andra,R&D Department,Chatenay-Malabry,92290,France
Keywords: Callovo-Oxfordian (COx) claystone Excavation damage Time-dependent behaviour Suction Viscoplastic model
ABSTRACT In order to evaluate the performance of deep geological disposal of radioactive waste,an underground research laboratory(URL)was constructed by Andra in the Callovo-Oxfordian(COx)claystone formation at the Meuse/Haute-Marne (MHM).The construction of URL induced the excavation damage of host formations,and the ventilation in the galleries desaturated the host formation close to the gallery wall.Moreover,it is expected that the mechanical behaviour of COx claystone is time-dependent.This study presents a constitutive model developed to describe the viscoplastic behaviour of unsaturated and damaged COx claystone.In this model,the unsaturation effect is considered by adopting the Bishop effective stress and the van Genuchten (VG) water retention model.In terms of the viscoplastic behaviour,the nonstationary flow surface(NSFS)theory for unsaturated soils is used with consideration of the coupled effects of strain rate and suction on the yield stress.A progressive hardening law is adopted.Meanwhile,a non-associated flow rule is used,which is similar to that in Barcelona basic model(BBM).In addition,to describe the damage effect induced by suction change and viscoplastic loading,a damage function is defined based on the crack volume proportion.This damage function contains two variables: unsaturated effective stress and viscoplastic volumetric strain,with the related parameters determined based on the mercury intrusion porosimetry (MIP) tests.For the model validation,different tests on COx claystone under different loading paths are simulated.Comparisons between experimental and simulated results indicated that the present model is able to well describe the viscoplastic behaviour of damaged COx claystone,including swelling/shrinkage,triaxial extension and compression,and triaxial creep.
In France,Andra has set up a program that includes in situ experiments at Meuse/Haute-Marne (MHM) underground research laboratory(URL),aiming at accessing the feasibility of constructing and operating a high-level and intermediate-level long-lived waste disposal facility in Callovo-Oxfordian (COx) claystone formation(Armand et al.,2013).To construct the URL,the excavation induced the damage of the host formation,and the ventilation in the galleries desaturated the host formation close to the gallery wall(Armand et al.,2016).Moreover,field monitoring showed significant time-dependent behaviour of COx claystone (Armand et al.,2014).Thereby,the benchmark ‘‘Transverse action” was proposed to develop constitutive models for characterising the claystone behaviour according to the field data (Armand et al.,2017;Seyedi et al.,2017).This benchmark aims to conduct numerical modelling at the URL scale using constitutive models adapted to COx claystone.
Up to now,several time-dependent models have been developed for unsaturated COx claystone based on different viscoplastic theories(Bui et al.,2016,2017;Bian et al.,2017;Mánica et al.,2017;Pardoen and Collin,2017).For instance,Mánica et al.(2017)adopted the rheological law proposed by Darve (1990) and the friction angle evolution method to predict the viscoplastic deformation of COx claystone and the strain softening after peak.Bui et al.(2016,2017) used the Perzyna’s viscoplasticity theory(Perzyna,1966)and damage model(Lemaitre and Chaboche,1990)to describe the relaxation and creep behaviours of COx claystone in different loading paths,for example,saturated,unsaturated,drained and undrained conditions.Bian et al.(2017) formulated a viscoplastic model for hard clay with the unified viscoplastic approach by Zhou et al.(2007) and a damage model by Mazars(1984) and then applied it to the modelling of excavation damaged zone(EDZ)in COx formation.Besides,Pardoen and Collin(2017) divided the plastic strain into two parts: an instantaneous strain and a time-dependent creep strain,and defined a viscoplastic loading surface based on the evolution of microstructure (Shao et al.,2003) to consider the creep deformation of claystone.However,these existing models have some limitations in the modelling of viscoplastic and damage behaviours.The overstress theory by Perzyna(1966)is a widely adopted method to model the creep and stress relaxation behaviours of geomaterials using different overstress functions.Nevertheless,these overstress-based models cannot meet the consistency condition as the general elastoplastic models (Liingaard et al.,2004).In this regard,Qiao et al.(2016)adopted the nonstationary flow surface (NSFS) theory in the viscoplastic modelling of saturated soils,which can meet the requirement of consistency condition and have a better convergence (Heeres et al.,2002).Few models can well predict the viscoplastic deformation under low stress in the conventional elastic zone.On the other hand,all the aforementioned models failed in explaining the damage behaviour of COx claystone,because they only account for the damage effect by using the theoretical damage model(e.g.Bian et al.,2017).Indeed,the damage mechanism of COx claystone mainly involves the generation of cracks related to the effects of suction and mechanical loading.This has been evidenced by the micro-observations on COx claystone under different hydromechanical loadings through different experimental methods(Wang et al.,2014;Menaceur et al.,2016).In this regard,as an alternative,the physics-based damage coefficient proposed by Dao et al.(2015)appears more appropriate to be used for modelling the damage of COx claystone.
In this study,a constitutive model for COx claystone based on the non-stationary flow surface theory is presented with consideration of a progressive hardening law.This model can reproduce rate-dependent and creep deformations through the change of viscoplastic strain rate with time.Besides,a damage function is proposed based on crack volume proportion identified from mercury intrusion porosimetry (MIP) results.Thus,the present model highlights the effects of viscosity and damage on the unsaturated behaviour of COx claystone.The model is further validated according to different experimental results of COx claystone under different loading paths,including swelling,rate-dependent triaxial compression/extension,and triaxial creep.
The viscoplastic model of unsaturated COx claystone is proposed and formulated in triaxial stress space.As in the viscoplastic model by Qiao et al.(2016),the strain increment contains an elastic component and a viscoplastic one,as follows:
where dεijis the total strain increment,dis the elastic one,andis the viscoplastic one.
In fact,the observations from rate-dependent loading tests(e.g.constant rate of strain (CRS) test and rate-dependent triaxial test)indicated that it is very difficult to distinguish the stress-induced strain (plastic strain) and time-induced strain (viscous strain).Thus,it is more reasonable to consider the strains induced by plasticity and viscosity as a whole.
In this model,the isotopic stressp′and the deviator stressqcan be calculated according to the effective stress tensoraccording to the Bishop effective stress concept (Bishop,1959),aiming at modelling the stresses of unsaturated soils in triaxial stress space.The specific expressions ofp′andqare simplified,as below:
where χ(Sr)is the unsaturated effective stress coefficient,depending on the degree of saturation;andsrepresents the suction.
As shown Eq.(2),the mean stress accounts for both the contributions of net stress and suction.For the latter,its contribution is quantified using the approach by Alonso et al.(2010)to model the variation of χ withSr,and the equation by van Genuchten(1980)to describe the water retention property,as presented by Eqs.(4)and(5).For the COx claystone,the model parameters related to water retention behaviour was determined based on the experimental data collected from different works,as shown in Fig.1.As for parametera,it needs to be determined according to the shear strength data and a value of 1.1 was recommended by Wang(2021).

Fig.1.Soil-water retention curve (SWRC) for COx claystone with experimental data from different works.
wherea(a>1) represents the material parameter governing the contribution ofsto effective stress;and αr,nandmare the model parameters related to water retention behaviour of soils.
whereKDandGDrepresent the elastic bulk and shear moduli accounting for damage effect,respectively.
The damage properties of soils are generally linked to the decrease of elastic moduli.In this regard,the adopted damage coefficientDis associated with the ratio of the current modulus to the one at intact state as shown in Eqs.(8) and (9).In this study,the determination of specificDvalue is linked to the generation of cracks inside COx claystone,as described in Section 2.3.
wherevrepresents the Poisson’s ratio,k is the elastic parameter,Eis the initial Young’s modulus,andGis the initial shear modulus.
Regarding the hardening law for unsaturated soils,the effective yield stress includes net yield stress and suction contribution(Sheng et al.,2004;Le Pense et al.,2016),as shown by Eq.(10).The former is similar to the loading collapse curve(Eq.(11))by Alonso et al.(1990),and the latter is the same as that used in the description of unsaturated mean effective stress (Eq.(2)).Eq.(11)can well describe the change of the net yield stress with suction.Meanwhile,the viscoplastic volumetric strainis adopted as a hardening variable,replacing the conventional volumetric strain,with consideration of the definition of strains in Eq.(1).
whereandare the effective and net yield stresses,respectively;v0is the initial specific volume,equal to 1+e0;and λ(0),β andrare the plastic parameters for calculating λ(s)at different suctions.These parameters can be determined based on the tests by Wang (2021) and Wang et al.(2023),as presented in Fig.2.
Moreover,to model the time-dependent deformation of COx claystone,the isotach approach for saturated soils was adopted.In the isotach approach,there exists a stress-strain-strain rate relation (?uklje,1957).The corresponding relation involving in yield stress and strain rate can be modelled using a linear function on the double logarithm plane,as shown in Fig.3 at the saturated state.The volumetric compression curve has a rate-dependent yield feature.In this study,considering the suction-dependent viscosity for unsaturated soils withCA(s),the saturated isotach approach was extended to unsaturation.The net yield stress evolves with viscoplastic volumetric strain rate and suction,as described in Fig.3.At a given suction,the relation regarding net yield stress and viscoplastic volumetric strain rate is also regarded as linear.The decrease ofCA(s)with increasing suction(Eq.(13b)) shows the reduction of viscosity with suction increase.Note that there is a reference point,at which the corresponding stress and viscoplastic volumetric strain rate can be identified based on the linear relationships at different suctions(De Gennaro and Pereira,2013).Besides,there is also a lower threshold limit of viscoplastic volumetric strain rate,below which the time-dependent deformation becomes insignificant and can be ignored (Qiao et al.,2016).Thus,Eq.(13) is proposed by combining the effects on the net yield stress:

Fig.3.Schematic illustration of changes in net yield stress with viscoplastic volumetric strain rate and suction (after De Gennaro and Pereira,2013).
wherebcontrols the evolution ofCA(s)with suction,andpatequals 0.1 MPa.
By substituting Eq.(13)into Eq.(11),the yield stress is deduced,as below:
Moreover,this study also adopts the concept of image point for defining the progressive hardening law.Fig.4 shows the identification of the image point by projecting the actual stress point on the normal yield surface.In this study,the stress ratio γ can be calculated with Eq.(15).Considering the shear effect,Eq.(16) is adopted with(Hong et al.,2016).Meanwhile,by using the viscoplastic strain,the viscosity effect on γ can also be accounted for,as below:

Fig.4.Schematic illustration of yield surfaces.

Fig.5.Pore size distributions of dried COx claystone at different suctions (experiment from Menaceur et al.,2016): (a) Cumulative void ratio,and (b) Pore size density function.

Fig.6.Pore size distributions of wetted COx claystone at different suctions (experiment from Menaceur et al.,2016):(a)Cumulative void ratio,and (b)Pore size density function.
For simplicity,the yield function adopts the modified Cam-Clay(MCC) form,the performance of which was verified by Alizadeh and Gatmiri (2017) while modelling the COx claystone behaviour.Combining the modification of the progressive hardening,the yield surfaces involved in this study are presented in Fig.4,including the inner yield surface,the normal one,the reference one and the threshold one.Different yield surfaces are also related to the definitions of yield stresses,as shown in Fig.3.This viscoplastic modelling approach is developed according to the NSFS theory,as presented by Qiao et al.(2016).The inner and normal yield functions are thus obtained,respectively,as follows:
whereMrepresents the slope of critical state line (CSL).
Besides,the adopted flow rule is non-associated,similar to that in BBM by Alonso et al.(1990).The visoplastic potential function is as follows:
where αgis used to describe the stress state under laterally confined condition.
2.3.1.Damage mechanism of COx claystone
Under the hydro-mechanical loadings,damage can be generated due to the micro-cracks in COx claystone(Wang et al.,2014,2015),which can affect the hydro-mechanical behaviour.Under hydration,the micro-cracks are created mainly along the interface of clay matrix and granular inclusions,as evidenced on environmental scanning electron microscope(ESEM)by Wang et al.(2014).Wang(2012)presented that because of the existence of inclusion,the COx claystone is locally heterogeneous.The clay matrix swells or shrinks when the relative humidity varies.However,the inclusions almost do not change and local stresses occur.Then,micro-cracks were generated and this discontinuity can further affect the elasticity of COx claystone.The phenomena were observed by Menaceur et al.(2016) according to MIP test results with the appearance of significant large pores inside COx claystone at low suction.Moreover,the loading with external stress can also generate micro-cracks,when the strength at the inclusion-matrix interfaces is lower than the external tensile stress,thus damaging the COx claystone.Furthermore,Wang (2012) and Wang et al.(2014) reported that the coarse grain inclusions can be split with increasing external stress,generating more micro-cracks.It is considered that these micro-cracks are mainly caused by plastic loading,when the external stress is larger than the yield stress.The ESEM observation by Wang et al.(2014)confirmed that local microcracks preferentially propagate along the inclusion-matrix interfaces due to plastic deformation.Similarly,Menaceur (2014)conducted numerous high stress oedometer tests on COx claystone,allowing the damage induced by compression to be investigated.The selected maximum stresses are significantly larger than the in-situ stresses.Through MIP observations on COx claystone unloaded from the maximum vertical stresses 56 MPa and 113 MPa applied,significant large pores were identified.Meanwhile,Mohajerani et al.(2011) evidenced that a higher maximum vertical stress caused a larger elastic rebounding slope k (in the ln σv′-eplane).This means that the high stress compression induced the crack generation,further increasing k,indicating the influence of crack on elasticity.The same phenomenon was also observed on a shale by Aversa et al.(1993)through the oedometer compression tests.
2.3.2.Damage modelling of COx claystone
As aforementioned,the damage mechanism involves in hydration and plastic loading.For modelling the damage behaviour,a coefficientDis adopted,regarding the ratio of φcrack(porosity of cracks) to φ (the total one) (Dao et al.,2015).It indicates thatDis physically correlated to the micro-cracks.For the determination of the micro-crack volume,a delimiting diameter 1 μm is adopted for COx claystone,which almost corresponds to the lowest positions of all the density curves (Figs.5 and 6) (Menaceur,2014).
Figs.5 and 6 show the MIP results under unconfined condition,where 34 MPa represents the initial suction of COx claystone.The detected void ratio changed significantly with the change of suction and the macro-pores(larger than 1 μm)were generated because of the swelling or shrinkage of clay minerals.Large pores appear induced by the generation of micro-cracks,indicating that the wetting/drying cycle can further damage COx claystone(Menaceur et al.,2016).Based on the micro-crack volume (Figs.5 and 6) and the total void ratio (Table 1),theDvalues at different suctions are determined.Besides,a few micro-cracks in COx claystone may close under loading,thus damage reduces (Zhang 2011;2013).The relationship betweenDand suctionsis not enough to describe this effect.For each suction,the corresponding effective stressp′can be calculated with Eq.(2),meanwhile,the damage coefficient can be determined according to the MIP data,collected from Menaceur et al.(2016).The relation betweenDandp′is thus obtained,as presented in Fig.7,and expressed by Eq.(20):

Table 1Testing conditions of MIP samples under unconfined hydration.

Fig.7.Evolution of D(p’)with mean effective stress p’.
where αDand βDare the model parameters governing the changing rate ofDwithp′.
Besides,Wang et al.(2014) showed that the micro-cracking in COx claystone is more pronounced after a wetting-drying cycle,which is evinced by Menaceur et al.(2016).However,the related effect on damage is difficult to be quantified due to the limited experimental data.Thus,it is ignored in this study.
Fig.8 presents the pore size distributions of saturated oedometer samples,loaded to different maximum stresses and then unloaded.Table 2 shows the testing conditions in detail.Based on the MIP results,the damage coefficients corresponding to different maximum stresses were determined.In terms of the mechanical loading effect,the plastic volumetric strain is used to model the change ofDat a given suction.Correspondingly,the plastic volumetric strains are calculated based on the oedometer compression curves by Menaceur(2014).In terms of the elastic parameter k,the value 0.005 is taken for calculating plastic volumetric strains(Alizadeh and Gatmiri,2017).The determinedDincreases with increasing,which indicates the damage induced by plastic loading.The change ofDwithis presented in Fig.9,and expressed by Eq.(21):

Table 2Testing conditions of MIP samples under saturated oedometer compression.

Fig.8.Pore size distributions of COx claystone unloaded from different stresses (experiment from Menaceur et al.,2016): (a) Cumulative void ratio,and (b) Pore size density function.

Fig.9.Evolution of D()with plastic volumetric strain .
whereD0is the initial value of damage coefficient with=0 at a given suction,and γDis the model parameter governing the changing rate ofD.
Thus,effective stressp′and plastic volumetric strainare used to model the change ofD.Thereby,both the effects due to hydration and mechanical loading can be considered.Combining Eqs.(20)and (21),Eq.(22) is thus deduced.According to the proposed damage function,the damage variable depends on both the volumetric plastic strain and the mean effective stress (including the influences of suction and confining stress).Thus,even in the elastic phase,the damage state can actually change,if the loading induces the changes of suction and confining stress.
In addition,to consider the time effect,Eq.(20) is revised by using the viscoplastic volumetric strain,as shown by Eq.(23).
Then,by substituting Eq.(23) into Eqs.(8) and (9),the damage effect is accounted for in the viscoplastic modelling.
The viscoplastic strains can be calculated with Eq.(24) and the corresponding dilatancy equation is presented as Eq.(25).
The consistency condition can be deduced as follows:
From Eq.(26),it is found thatcan change withwhen stress and suction keep constant.Besides,can also evolve with time without the changes of volume and suction,resulting in the typical stress relaxation phenomenon.Moreover,the introduction of progressively hardening law allows the viscoplastic behaviour under low stress to be reproduced.
The proposed model includes 22 parameters and their determinations include four aspects: (1) k,ν,Mand λ represents the are the elastic and plastic parameters,which can be determined based on unsaturated triaxial and oedometer test results and also allow β,randpcto be determined;(2)The water retention related parameters αv,mandncan be fitted according to VG model,by using the experimental data.In this model,the hysteresis behaviour of SWRC is not considered;(3) The viscous parametersCA,b,can be obtained according to unsaturated CRS tests;(4) The damage-related parameters αD,βDand γDwere determined according to the MIP results.The sensitivity analysis method can be used to obtainkanda.The determined parameters are listed in Table 3.

Table 3Model parameters for COx claystone.
In this part,the simulated and experimental results are compared.The proposed model is evaluated in terms of the swelling,shear and creep behaviours.
Fig.10 presents the comparisons between the experimental and the simulated volumetric strains with/without consideration of damage.In Fig.10a,the experimental data of COx claystone were collected from Vales(2008)under unconfined condition(wetted or dried from the same initial state).The initial state corresponds to the suction 35.4 MPa in a relative humidity of 76%.The swelling simulations with and without damage are found to be different.The latter underestimates the volume change,no matter at low or high suctions,similar to the modelling of Boom clay by Le Pense et al.(2016) and bentonite by Mokni et al.(2020).By contrast,the simulation with damage is better comparable with the experimental data under low suction range,while underestimating the volume change under high suction.The contribution ofsto χ under high suction is more difficult to be estimated (Alonso et al.,2010).This may induce the more significant differences in high suction range.

Fig.10.Comparison between simulated and experimental volumetric strains at different suctions (experimental data from Vales,2008;Pham et al.,2007): (a) Wetting or drying from the initial state,and (b) Drying-wetting cycle from a low suction.
Moreover,Fig.10b shows the experimental and simulated results accounting for damage.The experimental data by Pham et al.(2007)present the changes of volumetric strain of COx claystone in a drying-wetting cycle from an initial suction of 2.35 MPa under unconfined condition.The increase of swelling strain with increasing suction is well described by the model in the drying process.However,during the wetting,this model overestimates the volumetric strain.This difference is due to the determination of damage-related parameter,which is based on limited experimental data without considering the effect of drying-wetting cycle.Besides,the hysteretic effect on water retention behaviour was not considered in the present model,which may affect the modelling of volume change due to the drying-wetting cycle.
Fig.11 presents the simulated and experimental volumetric strains with the degree of saturation under unconfined condition.The tested samples have an initial degree of saturation of 98.3%(Cariou,2010).Compared with the simulation without damage,the simulation with damage can better reproduce the change of volumetric strain with the degree of saturation,although the difference is not significant.This is consistent with the evolution of volumetric strain with suction.Nevertheless,the increase of deformability with damage is well reproduced by the present model.It is also similar to that by Mohajerani et al.(2011) on COx claystone.

Fig.11.Comparison between simulated and experimental volumetric strains at different degrees of saturation (experimental data from Cariou,2010).
Fig.12 compares the simulated and experimental results under triaxial condition.The triaxial tests were conducted at the same axial strain loading rate and the same initial water content with different confining pressures and an initial relative humidity of 90%.More details about these tests can be found in Armand et al.(2017).Results show that with increasing confining pressure,the shear strength significantly increases.Good comparisons are obtained on the changes of deviator stress with axial and lateral strains between the simulated and experimental results prior to the maximum deviator stress.However,the softening behaviour is not reproduced.This is attributed to the changes in frictional angle and cohesion after the occurrence of strain localization,which is not accounted for in the proposed model.

Fig.12.Comparison of simulations and experiments under different confining pressures (experimental data from Armand et al.,2017).
Fig.13a displays the comparisons between simulations and experiments with and without damage at a confining pressure of 12 MPa.The triaxial tests were also conducted by Armand et al.(2017) with three unloading-reloading cycles.It is found that there are some differences for the simulated results with and without consideration of damage,although not significant.In Fig.13b,the evolution of normalizedEwith axial strain is obtained based on the unloading-reloading cycles and compared with the simulated ones with damage.The decrease ofEwith increasing axial strain is well reproduced in the simulation,which also indicates the evolution of damage under triaxial condition.

Fig.13.(a)Comparison of tests and simulations with and without damage at the confining pressure of 12 MPa,and(b)Evolution of normalized Young’s modulus with axial strain(experiment from Armand et al.,2017).
Fig.14 compares the experimental and simulated results under triaxial condition.Liu and Shao(2016)conducted the three triaxial tests and applied a confining pressure (12.4 MPa) and a constant axial strain rate (5 × 10-6s-1) under different constant water contents.The samples were loaded along the axial direction(perpendicular to the bedding plane).The three triaxial samples have different initial suctions (37.8 MPa,21 MPa and 4.2 MPa)corresponding to different relative humidity (76%,85% and 98%)prior to triaxial tests,respectively.In Fig.14,a smaller initial suction leads to a smaller stress at the same axial and lateral strains.In addition,it seems that at the same deviator stress,the difference between simulations and experiments for the axial strain ε1is larger than that for the lateral strain ε3,whatever the relative humidity.This phenomenon is also observed in Fig.12.Zhang et al.(2019) indicated that COx claystone has a rather anisotropic mechanical behaviour.In the direction parallel to bedding,the stiffness is larger,thus more difficult to deform (Zhang et al.,2015).It is inferred that the selected elastic parameter in this study is closer to the one in the direction parallel to bedding.Thus,the contribution of elasticity to the increase of deviator stress in simulations is relatively larger than that in experiments at the same axial strain.This may induce the larger differences in the axial direction than the lateral one.The difference can be reduced when considering anisotropy with different elastic parameters in two directions.However,the present model does not consider the anisotropy.

Fig.14.Comparison of simulations and experiments under different initial relative humidities (experiment from Liu and Shao,2016).
Fig.15 depicts the comparisons between experiments and simulations under triaxial condition following a constant mean stress path.For this loading path,the mean stressp=(σ1+2σ3)/3 keeps constant when increasing the axial stress σ1.The triaxial tests were performed at constant water content condition and different mean stresses (12 MPa and 14 MPa) (Liu et al.,2019;Wang et al.,2022).The comparisons show a good model performance in predicting the shear behaviour under constant mean stress.

Fig.15.Comparison of simulations and experiments under different mean stresses(experiment from Liu et al.,2019;Wang et al.,2022).
In Fig.16,the experimental results and simulations of COx claystone under triaxial compression condition at different axial strain rates are presented.The triaxial tests were also conducted under the mean stress of 12 MPa,same loading path as that in Fig.15.The initial water content of tested samples is 6.2% corresponding to the degree of saturation of 90.4%.They were loaded in the direction perpendicular to bedding.As presented in Fig.16,the deviator stresses for both simulations and experiments increase with the increase of axial strain rates,showing the rate dependent shear behaviour.In this aspect,the differences in Fig.16,no matter for experiment or simulation,are not obvious.That is because COx claystone is quite stiff and the viscosity feature and the ratedependency are thus not significant.However,considering the long-term nuclear waste disposal (up to several thousand years),COx claystone as the host formation may experience a very large time-dependent deformation.Thus,the rate-dependency,although insignificant,cannot be ignored,which is closely related to the time-dependent deformation at a large time scale.In addition,the comparison between experiments and simulations shows large difference of deviator stresses.This may be related to the yield function adopted for COx claystone.Hong et al.(2016)reported that the MCC yield function is not the best one for natural stiff clay.This may result in the differences between experiments and simulations.

Fig.16.Comparison of simulations and experiments under different axial strain rates(experiment from Liu et al.,2019).
Fig.17 presents the simulations of the triaxial extension tests under constant mean stress on COx claystone at the same relative humidity of 90%.The triaxial extension tests were performed under different mean stresses (12 MPa and 13 MPa),with axially unloading at a strain rate 1×10-6s-1(ε1<0 and ε3>0)(Armand et al.,2017).Note that the model does not account for the effect of Lode angle in this study.Therefore,the critical state ratioMein Table 3 was adopted in the triaxial extension simulation.Results indicate that with a higher mean stress,a smaller deviator stress(q=σ1-σ3<0)is obtained for both simulation and experiment at the same axial and lateral strains,highlighting the influence of the confining pressure on the mechanical behaviour.The simulations and experimental results are close,validating the proposed model in describing the behaviour under extension.

Fig.17.Comparison of simulations and experiments under different mean stresses in triaxial extension (experiment from Armand et al.,2017).
Figs.18 and 19 show the comparisons between the simulated and experimental results under triaxial creep condition at different deviator stress levels with the confining pressures of 2 MPa and 12 MPa.In this study,the deviator stress level represents the ratio of the applied deviator stress to the corresponding peak deviator stress.Each tested sample was first subjected to a designed confining pressure,then sheared till the target deviator stress.Afterwards,the deviator stress kept constant when creep started.
Fig.18 presents the comparisons of axial and lateral strains between simulations and experiments at different deviator stress levels (50%,75% and 90% of peak).The applied confining stress is 2 MPa and the maximum deviator stress is 22.7 MPa.Obviously,with increasing the deviator stress level,the increases of axial and lateral strains with time are more significant for both experiments and simulations.From the comparison it appears that the axial creep strains are well reproduced regardless of the deviator stress level,while for the lateral strains,the simulations overestimate the developments with time after a long term at the larger deviator stress.This can be attributed to the anisotropic behaviour,which is not accounted for.

Fig.18.Evolutions of axial and lateral strains with time in triaxial creep tests at different deviator stress levels under the confining pressure of 2 MPa (experiment from Armand et al.,2017).
In Fig.19,the experimental and simulated strains with time are compared at different deviator stress levels (50%,75% and 90% of peak).The corresponding confining stress and peak deviator stress are 12 MPa and 34.9 MPa,respectively.The continuous increases of strains with time are observed in both the axial and lateral directions.Similar to the observations from Fig.18,the comparisons of axial strains are also more satisfactory than those of lateral strains.This difference may also be attributed to the anisotropy of COx claystone.

Fig.19.Evolutions of axial and lateral strains with time in triaxial creep tests at different deviator stress levels under the confining pressure of 12 MPa(experiment from Armand et al.,2017).
In this study,an unsaturated rate-dependent constitutive model was formulated to simulate the viscoplastic behaviour of COx claystone considering damage.The model also considered the suction effect on the time-dependent deformation by adopting a suction-dependent linear relationship with respect to viscoplastic volumetric strain rate and yield stress.Besides,a damage function considering the effect of micro-cracks was developed based on MIP observations,aiming at describing the evolution of damage under hydro-mechanical loadings.In combination with the progressively hardening MCC yield functions,the model is able to well predict the shear behaviour of COx claystone in triaxial stress space.
Then,the constitutive model was validated using the experimental data from swelling,rate-dependent triaxial extension and compression,and creep tests.Comparisons under different loading paths indicated the good performance of this model.Based on that,the proposed model can be used for the assessment of excavation at the URL scale after implementation in a finite element code.However,it is also found that the anisotropic and strain softening properties of COx claystone affect the predictions of axial/lateral strains in triaxial shear and creep tests,which are not considered in the present study.Future work should be conducted to consider these features for improving the prediction of the time-dependent behaviour of COx claystone.
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Journal of Rock Mechanics and Geotechnical Engineering2024年1期