Draft Development of a viscoelastoplastic model for a bedded argillaceous rock from laboratory triaxial tests

Argillaceous rocks are candidate host and/or cap formations for the geological disposal of nuclear wastes in many countries, including Canada, France, and Switzerland. The understanding of the long-term mechanical behaviour of such rocks is an essential requirement for the assessment of their performance as a barrier against radionuclide migration. The French Institute for Radiological Protection and Nuclear Safety (IRSN) operates an Underground Research Laboratory (URL) in Tournemire, France, in a rock formation known as the Tournemire shale. Many types of experiments are conducted at the Tournemire URL in order to better understand the physical and chemical behaviour of this shale and its interaction with seal materials intended to be used in the geological disposal of radioactive wastes. The Canadian Nuclear Safety Commission (CNSC) collaborates with the IRSN and CanmetMINING laboratories to perform experimental and theoretical research on the mechanical behaviour of the Tournemire shale. Using the data from creep tests, and monotonic and cyclic triaxial tests performed at CanmetMINING Laboratories, we developed constitutive relationships for the mechanical behaviour of the Tournemire shale. The model is based on the theory of plasticity, and takes into consideration the inherent anisotropy due to the existence of bedding planes, hardening behaviour before the peak strength, and viscosity.


Introduction
Deep geological disposal is being proposed for the long-term management of nuclear wastes in many countries, including Canada, France and Switzerland.Sedimentary rock formations are one of the host media being considered.In Tournemire, France, IRSN operates an underground research laboratory in Tournemire located in an argillaceous sedimentary rock, the Tournemire shale that is found in a Mesozoic marine basin (Rejeb and Stephansson 2007).The Tournemire shale possesses strong anisotropy in stiffness, deformation and permeability (Niandou et al. 1997;Zhang et al. 2002;Zhang et al. 2004).Anisotropy of geomaterials always involves the compositional layered structures, i.e. bedding, layering, and foliation.Loading history, mineral constituents, and deposition intervals are all reponsible for the formation of anisotropy.Considering the potential risk of preferential pathways associated with the excavation induced damaged zone for radioactive contaminants, it is thus required to thoroughly understand the anisotropic behaviour of the sedimentary rocks of interest.
Experiments on sedimentary rocks have indicated that the general trend of compressive strength varies with the loading direction (Rejeb 1999).The weakest plane exists between 30-60 o of loading angle against the bedding plane, while the maximum strength is found parrallel or perpendicular to the bedding plane.Besides the strength properties, elasticity of these materials also indicates significant dependence on the loading orientation.The anisotropy in mechanical properties of sedimentary rocks leads to challenges in the numerical modelling of various underground geotechnical engineering problems within the framework of classical isotropic elastoplastic theories.Numerous efforts have been devoted to the development of appropriate constitutive models for anisotropic geomaterials.Amongst them, the microstructure tensor approach developed by Pietruszczak (2001) proved to be robust and relatively easy to use for the characterization and modelling of the constitutive behaviour of materials with transverse isotropy.The microstructure tensor approach was later applied by Nguyen andLe (2015a,2015b) and Le and Nguyen (2015) to model the behaviour of Opalinus clay during the excavation of a micro-tunnel, followed by water and gas injection in the tunnel.
Since the 1960s, the deformation and strength characteristics of brittle rocks have been studied by many researchers (e.g., Brace 1964;Bieniawski 1967;Wawersik and Fairhurst 1970;Martin and Chandler 1994).These studies come to a general consensus that with increased loading, the mechanical behaviour of brittle rocks is governed by the initiation, propagation and coalescence of cracks.The above processes can be characterized and evidenced by acoustic occurrence (Martin 1993;Popp et al. 2001;Vajdova et al. 2004;Ghazvinian et al. 2013).Argillaceous rocks formed from sedimentation are, however, considered to be semi-brittle and experimental data from laboratory testing (e.g.Abdi et al. 2015) indicate several distinctive features, as compared to brittle rocks; mainly the inherent anisotropy due to sedimentation, but also strong timedependency (Wang et al. 2012;Sun et al. 2014;Chauveau and Kaminski 2008;Zhang and Rothfuchs 2004;Fabre and Pellet 2006;Gunzburger and Cornet 2007), among other factors.Many researchers took into account the time-dependency by developing constitutive relationships within the theoretical framework of viscoelasticity, or elasto-viscoplasticity D r a f t e-Doc 4879168 3 (e.g., Simo 1987;Fafard et al. 2001;Rouabhi et al. 2007;Pellet et al. 2009;Xiong et al. 2014).The present study shows a systematic development of a constitutive model using data from a comprehensive triaxial testing program performed collaboratively by the authors.Tournemire shale samples from the IRSN's Tournemire Underground Research Laboratory (URL) in France were tested.The main physical phenomena found from the tests were inherent anisotropy due to bedding, stiffness and strength degradation, strain hardening, and energy dissipation found in cyclic tests.Most existing models would only consider one or two of the above characteristics; the present model integrates all of them together.

Experimental study of the mechanical behaviour of Tournemire shale
At CanmetMINING Laboratories in Ottawa, an experimental program consisting of monotonic triaxial, cyclic triaxial, and creep tests was performed using cylindrical shale samples of about 60 mm in diameter and 130 mm in height.CanmetMining Laboratories follow the ASTM test standards, which recommend a Height/Diameter ratio to be comprised between 2.0 and 2.5 (ASTM D4543).H/D = 2.17 is the ratio usually targeted by the laboratory.The shale samples were obtained from boreholes oriented at different orientations with respect to the bedding planes.The experimental results were previously reported by Abdi et al. (2015).The main features are discussed here in order to provide a conceptual basis for the ensuing development of a constitutive model.The geometry of the samples being tested is represented in Figure 1.The bedding angle β is the angle between bedding planes and the x-direction, while z-direction is for the axial loading.

Monotonic triaxial tests on Tournemire shale
Triaxial tests on brittle rocks usually show five phases in the stress-strain response (see e.g.Martin, 1993;Chandler and Martin, 1994).In phase I, the applied stress results in a closure of pre-existing cracks in the sample.This is followed by phase II of essentially linear elastic behaviour.The onset of microfracturing, call the crack initiation threshold, starts at the beginning of phase III.In phase III, crack growth is stable, meaning that crack propagation would stop if the load is removed.The microcracks are mostly intragranular, and oriented in the direction of the major principal stress.The onset of phase IV is usually called the crack damage threshold.In phase IV, unstable crack growth would start.This phase is marked by the onset of dilatancy.In Figure 2, the volumetric strain curve starts to show negative slopes at the onset of phase IV, indicating dilatant incremental volumetric strain.In addition to intragranular cracks parallel to the direction of the applied load, transgranular inclined cracks are formed.In the post-peak phase (phase V), macrocracks develop and eventually lead to the collapse of the sample By contrast to brittle rocks, a typical plot of the deviatoric stress and the volumetric strain versus the axial strain for Tournemire shale is shown in -The crack closure phase is absent for this particular test suggesting that, for the Tournemire shale, any initial cracks that exist in the sample are either inexistent or healed.Evidence of a crack closure phase is only found for unconfined tests with bedding angle β=0°.This suggests that for those tests, the crack closure phase would actually be due to closure of the bedding planes.-The crack damage phase seems to start at more than 95% of the peak.This is corroborated by the acoustic emission record, as reported by Abdi et al. (2015).
Apart from the above differences, similarly to brittle rocks, there is a linear portion of both the deviatoric stress and volumetric deformation curves up to approximately 30% of the peak.From that point, crack initiation and propagation are inferred up to the peak, resulting in nonlinearity of the curves.After the peak, collapse of the material occurs, with the localization of strain along narrow failure bands.
Figure 4 shows the strong influence of bedding plane orientations on the stress-strain and dilatancy behaviour of Tournemire shale.This influence is reflected in the stiffness and peak strength of the material.The Tournemire shale is found to be stiffer and stronger when loading is parallel to the bedding planes (β=90 o ).

Cyclic triaxial test
Figure 5 shows the stress-strain and volumetric variations of the Tournemire shale for a typical cyclic triaxial test.Hysteretic loops are observed from unloading-reloading sequences.The slope of the loops decreases with the loading level, indicating a degradation of the stiffness of the rock.Similar type of behaviour is also found for brittle rocks.However, argillaceous rocks such as the Tournemire shale exhibit larger energy dissipation as evidenced by the larger areas of the hysteretic loops.Furthermore, the areas of the loops increase with the stress level.This phenomenon is characteristic of a material where viscosity might play an important role in its stress-strain behaviour.

A constitutive model for the Tournemire shale
The comprehensive test data obtained from the laboratory test program summarized above allow us to formulate a constitutive model that would capture the main physical features found during those tests: -The existence of three regions of the stress-strain curves as described above.
-The directional dependence of the mechanical behaviour with respect to the bedding orientation.-The time-dependent behaviour that is evidenced from creep tests and the existence of important energy dissipation observed in cyclic triaxial tests.
It is also assumed that the stress-strain response from the triaxial tests represents the drained behaviour of the shale.The samples were obtained at natural water contents that correspond to D r a f t e-Doc 4879168 degrees of saturation higher than 90%.Therefore, due to the increased compressibility of the equivalent pore fluid (mixture of air and water) pore pressure generated during loading is assumed to be negligible (Nguyen and Selvadurai, 1995).Tests performed on oven-dried samples (Abdi et al. 2015) shows the significant influence of suction on the mechanical behaviour of Tournemire shale.However at higher saturation level (larger than 90%), the influence of suction seems to be small, and therefore was not considered in this study.Future experimental program should consider undrained tests on fully saturated samples with the measurement of pore pressure in order to verify the above assumption.
Similarly to Haghighat and Pietruszczak (2015) and Nguyen and Le (2015a), the model is developed within the framework of elasto-plasticity, using the Mohr-Coulomb yield criterion that is formulated to include directional dependence and strain-hardening of the yield parameters.Before yielding starts, the material is assumed to be viscoelastic.
The constitutive relationship is an equation that relates the stress increment to the strain increment: where dσ is the increment of the stress tensor (written as a vector); dε is the increment of total strain tensor (written as a vector); dε p is the increment of the plastic strain tensor (written as a vector); dγ is the increment of the viscous strain tensor (written as a vector); and D is the elastic stiffness tensor (written as a matrix).
In the following sections, we describe how D, dε p and dγ are derived in order to fully define the proposed constitutive relationship.

Determination of the elastic parameters of the stiffness tensor D
The secant modulus E 50 , determined from a secant line between the beginning of loading up to 50% of the peak load in a triaxial test, is used in this study to estimate the initial Young's modulus.That modulus varies with the direction of the bedding plane β.For a transversely isotropic material, the following relationship could be derived from rotation of the reference axes (Niandou et al. 1977): where E β is the modulus at bedding angle β; E L and E T are the in-plane and transverse moduli, respectively; ν LT is the transverse Poisson's ratio; and G LT is the transverse shear modulus.
Provided that the stiffness is known for a series of bedding angles, as is the case in this study, the transverse Poisson's ratio and shear modulus can be estimated by curve fitting using the above equation.Figure 6 shows the best-fit curves using equation (2).Due to the small difference of values for different confining stresses, the average value was taken as the required parameter for modelling purposes.
The cyclic triaxial test results show that the stiffness degrades with increasing loads.Similarly to Martin (1993), we found that this decrease could be expressed as a function of the cumulative irreversible deformation.Using the secant slopes of the stress reversal loops, a normalized curve (Figure 7) of stiffness degradation versus the effective plastic strain γ p is derived.The effective plastic strain is used as a measure of cumulative irrecoverable deformation, as detailed in section 3.3.The experimental data in Figure 7 shows very little disparity between different bedding angles, and thus a unique exponential function can be fitted to the experimental points.
Table 1 summarizes the elastic parameters for the Tournemire shale.The in-plane Poisson's ratio is found to be independent of the confining pressure.In order to estimate the in-plane Poisson's ratio, the bulk modulus K was first calculated by plotting the mean stress against the volumetric strain (ε v ).The mean stress is related to the volumetric strain by the stress strain equation: p=Kε v.The experimental values of p were plotted against the experimental values of ε v in the elastic range.A straight line passing through the origin was then plotted to best fit those experimental points.The slope of the line gives the value of K.The in-plane Poisson's ratio can be computed (

Determination of the viscous component
The stress tensor and strain tensor can both be decomposed into volumetric and deviatoric components.
where the mean stress p and the volumetric strain can be written as: For an elastic material, can be expressed as a function of by the mean and deviatoric components.For example, the isotropic elastic material obeys the following relationship 4879168  7 where K and G are respectively the bulk and shear moduli, which are scalars for an isotropic material, but tensors for an anisotropic material.
For a viscoelastic material, the deviatoric stress is not only dependent on the deviatoric strain but also on its time variations.The deviatoric stress is usually written in an integral form as follows (Simo and Hughes, 1998 ݐ‬ is a relaxation function, which can be approximated by a Prony series where N is the total number of viscous branches, and τ is relaxation time ߬ = ఎ ீ . A simple rheological model, i.e. the simplified general Maxwell model (also known as the standard linear solid model) that consists of a spring and dashpot in sequence, was implemented in the study of the rheological behaviour of the Tournemire shale.It is assumed that viscosity only influences the deviatoric components of stress and strain, as is the common assumption adapted in geomechanics.Viscosity under hydrostatic loading is rarely measured in the laboratory and is often poorly defined in the literature (Hasanov et al. 2016).A limited number of creep tests on shales under hydrostatic loading were recently reported.Villamor Lora and Ghazanfari (2015) tested unsaturated shale specimens and found that creep strain under hydrostatic loading is relatively small (approximately 20% ) compared to creep strain for similar deviatoric loading levels.Sone (2012) attempted to perform creep tests under hydrostatic loading, but found that the data were difficult to interpret due to experimental uncertainty caused by erratic instrument responses to the change in confining pressure.For the above reasons, the authors adopted the simplifying assumption that viscosity is solely dependent on the deviatoric component of stress.
As only one viscous branch is considered, the integral deviatoric stress can be rewritten as The above equation reduces to another form in the absence of viscosity, i.e.
which suggests that the total stress is imposed solely onto the bulk elastic material when the viscosity diminishes.
As shown in Figure 8, the total deviatoric stress σ d equals the sum of both the bulk stress and the viscous stress (σ q ) in the form of D r a f t e-Doc 4879168 The stress on the spring-dashpot branch takes the following form, where G 1 and q 1 are respectively the shear modulus and shear strain tensor of the viscous branch: η 1 and γ 1 are respectively the viscosity and the viscous strain tensor of the dashpot.
The above equation gives rise to the following formula where τ 1 is the relaxation time.
The total deviatoric elastic strain thus consists of two parts 1) the viscous strain ߛ ଵ that is induced by the dashpot deformation, and 2) the elastic strain q 1 that is attributed to the spring.Then we get Then the viscous strains are defined by the following equations:

Creep test and the parameterization of viscoelastic model
The time-dependent behaviour of the Tournemire shale is evidenced from creep tests.Figure 9(a) shows the variation of axial strain with elapsed time for an unconfined creep test on Tournemire shale.The experiment was carried out under unconfined conditions, e.g. the uniaxial compression loading mode, in a stepwise incremental manner as shown in Figure 9(b).Stress was increased by increments of 1 MPa and each increment was maintained for 1 hour.Samples were retrieved from the same boreholes as those for the triaxial tests.Since the samples were unsaturated but with degree of saturation higher than 90%, we excluded the effects of porewater pressure on their time-dependent behaviour, and assumed that the latter was due to viscosity.deemed to be sufficient, taking into account the time constraints, the extent and the objectives of the test program As shown in Figure 9(a), the loading orientation appears to influence the creep strain significantly.The initial strain of the sample with β=0 o is about 3 times of the one observed at β=45 o .When loading is parallel to the bedding plane (β=90 o ), the instantaneous strain and longterm creep strain are both at their minimum.The behaviour appears to be anisotropic in creep, as suggested in previous studies (Rejeb 1999;Pietruszczak, Lydzba and Shao 2004).
The differential equations of the Maxwell models could be solved either numerically or semianalytically.The authors used a finite difference method implemented into a spreadsheet to solve the Maxwell's equations under creep test conditions (e.g.constant loads).The viscoelastic constants used as input of the Maxwell model are the shear modulus and the relaxation time.The output is the time evolution of strain under constant load conditions, from which the creep rate could be determined.The authors used a trial and error approach, and vary the values of the shear modulus and the relaxation time until a visually good match is achieved between the numerical output of the spreadsheet results and the experimental data.We found that the calibrated viscoelastic constants could be considered as isotropic, as shown in Table 2, while the orientation dependent strain is attributable to the anisotropy in plasticity.

Determination of the plastic strain
The Mohr-Coulomb yield criterion can be written in the following form where I 1 is the first invariant of the stress tensor; J 2 is the second invariant of the deviatoric stress tensor; θ is the Lode angle; c is the cohesion; and φ is the friction angle.
or in a more compact form, where θ is the Lode angle in the range of (0, π/3).
When the yield criterion is reached, plastic strain occurs.The plastic strain is derived from the plastic potential equation.Here we use a nonassociative flow rule, with the plastic potential taking the form of a Drucker-Prager relation: In this regard, the plastic strain rate is given by where λ is the consistency parameter that can be determined by classical elasto-plastic theory with appropriate rearrangement of the differentiated form of the yield function (Simo and Hughes, 1998).
The effective plastic strain can be defined as (Hashiguchi, 2014) where ߜ is the Kroenecker delta.
In triaxial conditions, the above formula can be expressed in another form, where ߝ ଵ is cumulative axial plastic strain; ߝ is the cumulative radial plastic strain; and ߝ ௩ is the cumulative volumetric plastic strain.

Initial yielding
Here we hypothesize that initial yielding starts at the crack initiation threshold, at a certain percentage of the peak strength (σ 1p ), ߸(%), and the elastic modulus decays gradually with accumulated damage.Then we obtain the initial yield stress where ߪ ଵ is the initial yield stress and ߪ ଵ is the peak value of the major principal stress.
In case of stress higher than the crack initiation threshold value, the plastic strain in the principal stress direction starts to grow with respect to the crack initiation threshold value where ߝ is the plastic strain in principal direction; ߝ is the strain in principal direction and ߝ ధ is the strain corresponding to crack initiation.

Hardening law -a normalized hardening model
A normalized hardening yield surface based on the Mohr-Coulomb criterion is proposed here.A normalized damage parameter is first defined in terms of ω as an important variable contributing to the strain hardening process.It relates the effective plastic strain to the plastic strain at peak strength and is transformed into a unified range between (0, 1).Similar processing of the effective plastic strain has been carried out by Martin (1993) in the analysis of strength of the Lac du Bonnet granite.
where ω is the normalized damage factor, ߝ is the plastic strain corresponding to the peak strength that marks the end of strain hardening and the start of softening.For the Tournemire shale, a linear relationship is found to exist between the plastic strain at the peak axial stress and the confining pressure, as shown in Figure 10.
As discussed above, the crack initiation threshold of the Tournemire shale is determined to be 30 % of the peak strength.When that threshold is exceeded, plastic yielding is assumed to occur with hardening until the peak strength.As shown in Figure 11 The parameters of the above equations, i.e. k, m, h and l, are directionally dependent as shown in table 3. Figure 12 shows the evolution of friction and cohesion as calculated from equations ( 31) and ( 33), with β=90°.It is interesting to compare the evolution of cohesion and friction of the Tournemire shale to the one for the same parameters of brittle rocks.The friction angle of the shale continuously grows to the peak value with increasing damage, and is in good agreement with Martin's (1993) conceptualization.However, the cohesion of the Tournemire shale does not decrease as compared to brittle rocks, but instead increases up to the peak strength.Loss of cohesion in brittle rocks is attributed to crack initiation and propagation.For the Tournemire shale, that phenomenon is offset by an increase in cohesion due to an increase in the mean stress when the load is increased.The effects of the mean stress on shear strength parameters are typically more pronounced for softer rocks like shale as compared to harder rocks like granite.

Microstructure tensor approach for yield parameters
The existence of bedding in Tournemire shale results in inherent anisotropy of its yield parameters. Considering only the bedding angle β=0°, the anisotropy reduces to a case of transverse isotropy, with two principal directions parallel to bedding, and the third being perpendicular to it.As in Nguyen and Le (2015), we used the microstructure tensor approach of Pietruszczak and Mroz (2001) in order to take into account the transverse isotropy of the yield parameters.In that approach, a generalized loading vector is first defined from the stress tensor as: where ݁ (ఈ) , α=1,2,3, are the base vectors in the principal directions of the transverse isotropy.
The unit vector along l i is given by A microstructure tensor a ij , which is a measure of the material fabric, is introduced.The projection of the microstructure tensor on l i becomes The scalar variable η, referred to as the anisotropy parameter, specifies the effect of load orientation relative to material axes.The above equation can also be expressed as where A ij is a symmetric traceless operator.This relation may be generalized by considering higher order tensors, i.e.,

D r a f t
e-Doc 4879168 The above representation is rather complex in terms of implementation and identification.Therefore, it is convenient to use a simplified functional form such as: It could be shown that (Pietruszczak and Mroz, 2001): where A 1 is the principal component of A ij in the bedding direction; l 3 is the component of the loading unit vector in the direction perpendicular to bedding (i.e.x 3 axis in Figure 1).
Nguyen and Le (2015a) expressed the Mohr-Coulomb strength parameters, e.g.c and φ, as a function of the microstructure tensor and loading orientation using the functional form of equation ( 39).In this study, c and φ are described by the hardening functions in equations ( 31) and ( 33).The parameters of those equations are then expressed in the functional form of equation ( 39), with best fit parameters as shown in Figure 13.The resulting parameters for the crack initiation threshold and the peak strength are shown in Figure 14.

FEM model for triaxial test
A 3-D FEM model for the triaxial test was developed with the commercial finite element software COMSOL Multiphysics (Version 5.1) in order to solve the equations of static equilibrium.The FEM model has a height of 133 mm and a diameter of 61 mm in order to exactly represent the dimensions of the test samples (Abdi et al. 2015).The FEM mesh of the model is shown in Figure 15.Tetrahedral elements with quadratic shape functions were used.
The bottom plane was assigned with fixed displacement in all directions while the side surfaces were subjected to a confining pressure Pc.The upper boundary was prescribed with timedependent displacement in the vertical direction.The evolution of vertical stress was monitored throughout the simulation for comparison with the test data.Modelling with different confining pressures, e.g.0, 4 and 10 MPa, and different bedding orientations, e.g.0°, 45°, and 90 o , was performed.In this study we focussed on the behaviour of the Tournemire shale up to the peak stress.In the pre-peak state, when the crack initiation threshold is exceeded, damage of the sample is due to the initiation and propagation of microcracks.The authors adopted Martin's (1993) hypothesis by using the accumulated effective plastic strain as a macroscopic quantitative measure of damage in order to simulate the degradation of the elastic parameters (Figure 7) and the evolution of the parameters of the yield function (in this case the cohesion and friction angle).

Viscoelastoplastic modelling of static triaxial test
After the peak, cracks and microcracks coalesce along narrow zones and the sample usually fails in a brittle manner along those zones.The modelling of the post-peak behaviour is currently being investigated by the authors and is outside of the scope of the present paper.

Stress-strain behaviour
Figure 16(a-c) shows the modelling results of stress-strain curves in comparison with test data.Good agreement between modelling and test data can be observed.In particular, the model can correctly simulate the influence of confining stress and the direction of loading with respect to the bedding orientation.The linear hardening law with which the model parameters were determined seems sufficient to reflect the gradually curved growth of axial stress with increasing axial strain.The peak strength fits well with the estimated Mohr-Coulomb strength envelope.

Volumetric strain versus axial strain
Figure 17(a-c) shows the calculated volumetric strain in comparison with the test data.The trends for the volumetric strain evolution are well reproduced by our model.In particular, consistently with the experimental data, the modelling results show that: -The elastic linear portion of the volumetric deformation is restricted to stress levels below the crack initiation threshold, which corresponds to approximately 30% of the peak axial load.-The point of reversal of the volumetric versus axial strain curve identifies the crack damage threshold.Consistently with the experimental data, the model indicates that the crack damage threshold for the Tournemire shale is very close to the peak axial stress.
The calculated values for volumetric compression are consistently higher than the experimental values for all angles of bedding orientation.Therefore, dilation as estimated by the model is lower than indicated by the test data.After a thorough re-examination of the test data and procedure, the authors believe that the way volumetric strain was estimated is probably inaccurate.During the tests, circumferential strain was mechanically measured with a ring wire.The volumetric strain was then estimated from that circumferential strain, assuming it was uniform along the height of the sample.Using that method, the volumetric deformation as shown in Figure 19 suggests that the pre-peak dilatancy angle would be approximately 60 o .This value is approximately three times higher than the peak friction angle as shown in Figure 12, and is deemed to be unrealistic.Niandou et al. (1997) performed triaxial tests on smaller Tournemire samples (37 mm in diameter, 75 mmm in height, versus 60 mm and 130 mm in this study).Tests performed by Medhurst (1996) showed that smaller samples should result in higher dilation; yet, D r a f t e-Doc 4879168 the test data reported by Niandou et al. (1997) showed lower dilation compared to the ones estimated from the present experimental data.Furthermore, the modelled volumetric deformation reported in this study compares well with the experimental results from Niandou et al. (1997), for all loading orientation, giving additional confidence in the adequacy of the model and the input data .

Viscoelastoplastic modelling of creep test
The calibrated viscoelastic constants shown in Table 2 were used as input to the finite element model shown in Figure 15 to verify that good agreement was achieved between the COMSOL results and the experimental data from the creep tests.Figure 18 shows the viscoelastoplastic modelling of the unconfined creep behaviour of Tournemire shale with different loading orientations.It is found that the modelling reproduces the test data very consistently.The relaxation time and the Maxwell shear modulus were determined to be 500 s and 3.5 GPa, respectively, for all three loading orientations.Our results show that a unified set of parameters for the viscoelastic model applies to Tournemire shale, indicating the isotropy in viscoelasticity.This is in agreement with other reports on rheology of anisotropic clay rocks (Zhang et al. 2002).
For instance, after extensive experimental investigation of Callovo-Oxfordian argillites, no obvious anisotropy effect on the creep behaviour of the clay rock was found to be correlated to the mineralogical composition and water content factors (Zhang et al. 2002).
It is also worth mentioning that the creep rate at the end of the test averages 8x10 -9 (s -1 ) in this study.Fabre and Pellet (2006) reported a smaller creep rate for Tournemire shale at 4x10 -10 (s -1 ), while Zhang et al. (2002) determined the steady state creep rate of Oxfordian argillite at 5x10 -11 (s -1 ).These experiments were conducted with different time duration, i.e. the longer the test lasts, the lower the creep rate becomes.Zhang et al. (2002)'s creep test was finished over a year, while Fabre and Pellet (2006)'s and ours were done within a day.The creep rate may accelerate if cyclic loading takes place and the level of shear stress overcomes a certain threshold value (Fabre and Pellet 2006).

Viscoelastoplastic modelling of cyclic triaxial test
Figure 19(a-f) shows the stress-strain behaviour of Tournemire shale under cyclic loading conditions.The modelling is able to reproduce the hysteretic cycles as observed in the test data with good agreement.It is noted that the viscoelastic constants that are used for the modelling of cyclic loading tests are identical to those calibrated from the creep test, suggesting that the rheological model and identified parameters are representative of the material properties.The simulated stress appears to be higher than the observed value before the first stress reversal.This could be attributed to the closure of cracks under increasing stresses, which is not considered in the modelling.When the axial stress is reversed, there is always a certain amount of residual strain even at very low stress levels.Abdi et al. (2015) interpreted such accumulated residual strain as plastic strain, which may not be appropriate according to our analysis.The residual strain should be attributed solely to viscosity when stress is lower than the crack initiation threshold, under which conditions no plastic strain is likely to take place.Compared to previous studies on the mechanical behaviour of Tournemire shale (Niandou et al. 1997), the compressive strength, axial and volumetric strains reported by Abdi et al. (2015) are all in the lower end, which has been thought to be a scale effect due to the difference in sample sizes.The residual strain after removal of the deviatoric stress is between two-thirds and one-half of the total strain in Abdi et al. (2015), while Niandou et al. (1997) reported a value of one-half.This could also be caused by the different loading-unloading rate.According to the viscoelastic model, a higher loading-unloading rate would induce a higher stress increment, and thus a lower strain change under the same range of stress variations.Zhang et al. (2002) reported that the creep behaviour of argillites is not influenced by scale effects even for a factor of two in size difference of the tested specimens.This is in agreement with the inference of our present rheological analysis of the Tournemire shale.

Conclusions
The understanding and prediction of the mechanical behaviour of argillaceous rocks are important in many areas of industrial activity in this rock type, including fracking for oil and gas exploitation, carbon dioxide sequestration, underground mining and geological disposal of radioactive waste.In this work, the authors took advantage of a comprehensive set of experimental data on the Tournemire shale (Aveyron, France) in order to develop a constitutive relationship for argillaceous rocks.The experimental data were generated from three sets of triaxial tests: monotonic compression, creep and cyclic loading-unloading tests, using samples with different bedding orientations with respect to the axial loading.The constitutive relationship being proposed is based on the viscoelastoplastic theoretical framework.It includes directional dependence of the elastic and plastic parameters, strain-hardening, and viscosity in order to simulate the time-dependent behaviour.A unique set of model parameters could be determined from the experimental data.That set of parameters was then used to perform the modelling of the triaxial tests, using the finite element method to solve the equation of static equilibrium.The ability of the model to reproduce the main physics that prevailed in all tests, using a unified set of experimental data, provides confidence in its robustness.
In order to validate the proposed constitutive relationship and its practical application, it is currently being implemented in a finite element model to simulate the shape and extent of the EDZ around tunnels excavated in 1996, 2003 and the century-old tunnel at the Tournemire Underground Research Laboratory.

Derivation of the hardening laws
In the hardening stage, the major principal stress can be approximated by the following equation: where ω is a normalized damage parameter, ߪ ଵ is the stress level corresponding to the initiation of plastic strain, and the parameter α is a model constant to be determined by experimental data as shown in Figure A1.
Actually ߪ ଵ determines the crack initiation threshold that we assumed represents the initial yielding.Instead of focusing on the peak strength, the elasto-plastic transition stress state has to be identified in this study in order to calculate the strength parameters.
Tests with three confining pressure levels have been conducted on the Tournemire shale.These led to the estimation of a set of Mohr-Coulomb parameters noted as ܿ and ߶ .
By analyzing the test data, an empirical linear function of the confining pressure is found to be valid for the parameter α in such a form as where the model constants k and m can be determined by fitting the experimental data for the hardening stage.As shown in Figure A2, a good linearized form of the test data can be obtained by this equation.The slope and intercept are simultaneously determined.This model takes into account both the initial yielding and the peak failure by a normalized damage parameter ω.Let us look at the two extreme cases with ω=0 and ω=1 representing the initial yielding and the peak failure, respectively.
The transition between the initial yielding and the peak failure can be modelled smoothly and seamlessly with a series of varying value of ω.The evolution of the Mohr-Coulomb envelope from initial yielding to peak failure can also be reflected with the proposed hardening model as shown in Figure 13.
This enables the derivation of the following function

Figure 3 .
Compared to Figure2for brittle rocks, there are several differences: matches the Mohr-Coulomb yield criterion at compressive meridians.The Drucker-Prager surface is smooth and its use as a yield potential avoid the singularity points at the corners of the Mohr-Coulomb surface, reducing the possibility of numerical divergence associated with those points.Even if the flow potential is based on a Drucker-Prager surface, one can define an equivalent dilatancy angle at a given Lode angle value.Equation (23) specifies that this dilatancy angle is the same as the friction angle of the Mohr-Coulomb surface for compressive meridians.
, the hardening stage encompasses a gradual change of Mohr-Coulomb strength parameters that can be related to the normalized damage factor ω. The development of analytical expressions for the hardening law is detailed in the appendix.The expressions are as follows:

Figure 1 .Figure 2 .Figure 8 .Figure 18 .
Figure CaptionsFigure1.Geometry of the samples used in triaxial tests, with indication of bedding angle (between the bedding and the x-direction)

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 7 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 14 .Figure 16 Figure 17 Figure 18 .Figure 19 .Figure
Figure 1.Geometry of the samples used in triaxial tests, with indication of bedding angle (between the bedding and the x-direction) Load duration was limited to 1 hour mainly because of time constraints imposed by the laboratory schedule and the test program.Nevertheless, in the authors' opinion, this duration is sufficient to exhibit both the full elastic and primary creep response of the rock, and reach the steady-state secondary creep stage, which was the main objective of these tests.In total there were 11 creep tests, with 7 tests on samples with β=0 o , and one test each for β=30, 45 and 60 and 90 o .The variability of the material is minimal, less than 10% -hence, this number of tests was

Table 2 .
Viscoelastic model constants best-fitted by Maxwell model for Tournemire shale

Table 3 .
Orientation dependent hardening parameters determined from curve fitting on test data