## Abstract

A laboratory programme of uniaxial, triaxial, cyclic and Brazilian tests was conducted to investigate the anisotropic mechanical behaviour of the Tournemire argillite, with different axial loading orientations with respect to the bedding planes (i.e. loading orientation angle, θ=0°, 30°, 45°, 60° and 90°). The experimental results show that both strength and deformation of the argillite are direction-dependent. Failure occurs in a brittle manner with a sudden collapse of the material strength. The failure mode exhibits localization along distinct failure planes and also depends on the loading orientation. This paper summarizes the experimental results and describes constitutive relationships that were developed in order to simulate the stress–strain behaviour of the Tournemire argillite. A microstructure tensor approach is adopted in order to take into account the anisotropic behaviour of the argillite. The identification procedure for material function and parameters is outlined, and the model is applied to simulate the set of triaxial tests performed at different levels of confining pressure and orientation of the bedding planes. It is demonstrated that the model adequately reproduces the anisotropy, the pre-peak stress–strain response and the onset of material collapse in those tests.

Argillaceous rock formations that are being considered as potential host or cap rocks for the geological disposal of nuclear wastes are usually characterized by the presence of bedding planes, which result in anisotropy of their strength and deformation properties. The anisotropic mechanical behaviour of the argillaceous rocks could affect the extent and characteristics of the damage induced by excavation of a deep geological repository (DGR) in these rock formations, which may have implications on their containment capability for nuclear wastes (Blümling *et al.* 2007; Millard *et al.* 2009; Labiouse & Vietor 2014). Therefore, an understanding of the anisotropic mechanical behaviour of these argillaceous rocks is one of the important aspects in the development of the safety case for a deep geological repository. In Canada, Ontario Power Generation is currently proposing a DGR for the management of its low- and intermediate-level nuclear wastes in an argillaceous limestone, known as the Cobourg limestone, which is overlain by more than 200 m of shale cap rocks. Nuclear Waste Management Organization is currently under a process to select a site for the management of used nuclear fuels in Canada in which the sedimentary rock is considered as one of the potential host rocks. The Canadian Nuclear Safety Commission (CNSC), Canada's nuclear regulator, is responsible for the licensing of future deep geological repositories in Canada. In order to better understand the ability of these argillaceous rock formations to contain and isolate the wastes, and to develop a tool to independently assess the extent and characteristics of damage induced by excavation of a deep geological repository, the CNSC has collaborated with different national and international research organizations and has access to underground research laboratories (URLs) worldwide. One of these URLs is situated in Tournemire, France, in upper Toarcian argillaceous rock (Boisson *et al.* 2001). The Tournemire argillite is characterized by the presence of closely spaced bedding planes and exhibits a strong anisotropic behaviour, so that its strength, as well as its deformation properties, are directionally dependent.

Over the last few decades, an extensive research effort has been devoted to studying the mechanical behaviour of anisotropic rocks. Comprehensive references on this topic can be found in a number of review papers (e.g. Amadei 1983; Kwasniewski 1993; Ramamurthy 1993). The transverse isotropy in geomaterials has been examined mainly through triaxial tests, and has been found to be significantly important in the analysis and design of a variety of geotechnical structures, such as foundations, retaining walls and slopes (Casagrande & Carillo 1944; Arthur & Menzies 1972; Oda *et al.* 1978). Other experimentally observed features, like the onset of shear banding and the influence of the intermediate principal stress, which cannot be properly described by isotropic criteria, have also been extensively examined (Gao *et al.* 2010). The strength of sedimentary rocks was found to be strongly affected by the loading direction. A large number of triaxial compression tests were conducted on orientated samples (Donath 1961; McLamore & Gray 1967; Hoek 1968, 1983; Atwell & Sandford 1974; Lerau *et al.* 1981) and the results generally indicate that the maximum strength is associated with specimens in which the direction of major principal stress is either parallel to or perpendicular to the bedding planes, while the minimum strength has been observed for orientations of beddings between 30° and 60°.

Besides experimental studies, an extensive research effort has also been devoted to modelling the mechanical behaviour of anisotropic rocks. Comprehensive reviews on this topic, examining different approaches, are provided by Kwasniewski (1993) and Duveau *et al.* (1998). One of the first attempts to describe the conditions at failure in anisotropic rocks was the work reported by Pariseau (1968), which was an extension of Hill's criterion (Hill 1950). This was followed by more complex tensorial representations (Boehler & Sawczuk 1977; Nova 1980; Amadei 1983). The application of the latter criteria to practical problems is generally difficult due to a large number of independent material functions and/or parameters that appear in the formulation. A simple and pragmatic approach, which incorporates a scalar anisotropy parameter that is a function of a mixed invariant of the stress and the structure orientation tensor, has been developed by Pietruszczak & Mroz (2000, 2001). This approach was later applied to the modelling of sedimentary rocks (Pietruszczak *et al.* 2002; Lydzba *et al.* 2003; Lade 2007).

In this study, a laboratory experimental programme of uniaxial, triaxial, cyclic and Brazilian tests with monitoring of acoustic emissions (AE) was performed on the Tournemire argillite at the CANMET laboratories in Ottawa, Canada, to investigate the mechanical behaviour of the argillite (Abdi & Evgin 2013; Abdi *et al.* 2015). The results of the triaxial tests were modelled using a mathematical framework that employs the microstructure tensor approach. This paper summarizes the experimental results and presents the numerical description of the anisotropic characteristics of the Tournemire argillite.

## Description of samples and testing equipment

The Tournemire argillite contains closely spaced bedding planes dipping 5° towards the north and randomly distributed fractures. The rock samples were retrieved in two sets from the Tournemire URL (Fig. 1) in seven boreholes at different angles (0°, 30°, 45°, 60° and 90°) to the sub-horizontal bedding planes (Fig. 2, five boreholes shown only). The first set of samples was drilled between 24 May and 1 June 2011, and the second set was drilled between 5 and 22 December 2011. Boreholes were drilled using dry air (compressed air) as the drilling fluid to minimize interactions of the rock with the drilling fluid, which allows better preservation of samples for further analyses. The drilling was coupled with an industrial vacuum device with particle filtering to minimize dust. The seasonal changes in relative humidity (RH) and temperatures in 2011 were measured by Hedan *et al.* (2014). The temperature at the Tournemire URL in late May and early June 2011 was about 12°C, while the RH was approximately 70–80%. The temperature at the URL in December 2011 was between 11 and 12°C, while the RH was approximately 45–75%.

A series of 1 m-long cores were obtained using a single-tube drill bit (76 mm outer diameter) coupled with a core lifter (61.7 mm inner diameter). Core logging (lithology, palaeontology and structures) was carried out right after core extraction. Each core was then cut into 20–25 cm-long sections, which were wrapped in aluminium foil, vacuumed then sealed. Finally, each sample was bubble-wrapped and all samples were carefully emplaced within a rigid box, which was sent to the CANMET laboratory in Ottawa, Canada, by plane freight.

Rock samples were received at the laboratory as sent (i.e. sealed to preserve original moisture) less than 1 week after their departure from France. The samples were immediately transferred and stored inside the CANMET Rock Mechanics Laboratory environmental chamber at ambient temperature (22°C) with 85% constant RH. Test specimens were prepared according to the American Society for Testing and Materials (ASTM) standards D3967 (ASTM 2008*a*) and D4543 (ASTM 2008*b*), and then sealed back under vacuum (one test specimen per bag) and returned to the environmental chamber. Testing was carried out shortly afterwards according to ASTM standards D3967 (ASTM 2008*a*) and D7012 (ASTM 2007) at the ambient laboratory temperature and humidity. The moisture content and porosity of the Tournemire argillite were measured for four samples before the testing.

The testing specimens were cylindrical samples about 60 mm in diameter and 130 mm in height for uniaxial, triaxial and cyclic tests. Shorter specimens with the same diameter and a nominal length of 40 mm were used for Brazilian tests. The depths of the tested specimens from the gallery boundary are provided in Table 1. The testing equipment consisted of a Mechanical Testing Systems (MTS) rock mechanics system (model 815), with linear transducers to measure the axial displacement, a chain extensometer to measure the circumferential displacement at the centre of the test specimen, transducers to measure porewater pressures and a data acquisition module. A detailed description of the MTS system was reported by Labrie & Conlon (2008). Acoustic emissions were recorded through a measuring system developed by Physical Acoustics Corporation (Princeton Junction, NJ, USA), which was operated simultaneously but independently of the main MTS.

## Summary of experimental results

Laboratory investigation of the mechanical behaviour of the Tournemire argillite was performed at the CANMET laboratories in Ottawa, Canada. The experimental programme and the detailed results were presented by Abdi & Evgin (2013) and Abdi *et al.* (2015). Some main facts and experimental results are summarized here for the purpose of numerical characterization of the mechanical behaviour of the Tournemire argillite.

The test loading pressure was applied at different angles (θ) with respect to the bedding planes (Fig. 3): that is, parallel to the bedding planes (θ=0°), perpendicular to the bedding planes (θ=90°) and at three different inclinations (θ=30°, 45° and 60°). Confining pressures of 0, 4 and 10 MPa were used for uniaxial and triaxial tests based on expected *in situ* stresses of σ_{v}≈σ_{H}≈4 MPa. The choice of loading rate was based on the nominal capacity of the test system. A 0.03 mm min^{−1} loading displacement rate was used for all tests.

### Basic physical properties of Tournemire argillite

The following properties were determined from laboratory experiments. The moisture content of the rock specimens was 3.86% with a standard deviation (SD) equal to 0.078, the bulk density was 2550 kg m^{−3} (SD=0.0244), porosity was 9.52% (SD=0.211) and the anisotropy coefficient of primary (P)-wave velocity was 2.16. The mineralogical composition of the Tournemire argillite was reported by Niandou *et al.* (1997).

### Acoustic emission data

The AE technique was used to detect the initiation and propagation of microcracks in rock specimens under mechanical loading. Based on the experimental AE data, new microcracks start to develop at a stress level where the crack initiation (σ_{ci}) is approximately 70–75% of the peak strength. At a stress level above σ_{ci}, the development and propagation of microcracks accelerate with the increase in stress level. In addition, according to the observations made in the experimental data, the crack damage stress level (σ_{cd}) is reached at a stress level where σ_{cd} is approximately 85–90% of the peak strength. The recorded data indicate that the stress levels of damage initiation and growth for the argillite are different to those of hard rocks (Martin 1993).

### Results of Brazilian indirect tensile tests

Brazilian tests were carried out on specimens obtained from different boreholes at depths ranging from 3.55 to 7.4 m (Table 1). The experimental results show that the tensile strength is influenced by the loading orientation angle, θ, and varies between 4.8 MPa (θ=0°) and 5.9 MPa (θ=90°). However, the variation of tensile strength with a loading orientation of between 30° and 60° is small (i.e. between 5.2 and 4.9 MPa).

### Results of uniaxial and triaxial tests

Uniaxial and triaxial tests were carried out on specimens obtained from depths ranging from 3.55 to 6.85 m (Table 1). Typical stress–strain curves are shown in Figures 4 and 5 for all tested loading orientations and confining pressures. Table 2 presents the experimental results of peak strength, and the axial and volumetric strains at peak strength for different values of loading orientation angles.

The experimental results show that both peak strength and axial strain depend on θ and the confining pressure. The maximum strength was obtained in tests performed at θ=0° (i.e. for loading applied parallel to bedding). The strength obtained in tests performed at 30°≤θ≤60° exhibits the lowest values. This indicates that the compressive strength is highly sensitive to the load inclination. The ratio between the principal stress difference at peak strength when θ=0° and the principal stress difference at peak strength when θ=90° increases slightly with the increase in confining pressure, but its value is nearly equal to unity. This indicates that the strength anisotropy between the two principal axes, θ=0° and 90°, is small. The ratio between the maximum and minimum peak strengths for the same confining pressure is nearly equal to 2, and decreases with increase in confining pressure, which indicates that the strength anisotropy becomes smaller for higher confining pressures.

The axial strain at peak strength increases with increase in confining pressure. Small values of θ (<30°) have no significant effect on the values of the axial strain at peak strength. However, the axial strain at peak strength increases significantly with larger values of θ (>30°). These experimental results indicate that the compressibility of the rock in the direction perpendicular to bedding is higher than the compressibility of the rock along the bedding. The experimental data also indicate that Young's modulus, *E*, decreases with increase in θ. The maximum Young's modulus value was obtained at tests performed at θ=0°; the Young's modulus value was lowest at θ=90°.

The failure of rock specimens occurs in a brittle manner with a sudden collapse of the material strength after peaking, as shown in Figures 4 and 5. Fractures observed on failed samples indicate that failure occurs in three principal modes: extension (splitting), shearing, and a combination of extension and shearing. The failure mode is a function of the loading orientation, θ, as shown in Figure 6. The range of confining pressures used in this investigation has little effect on the failure mode.

## Numerical simulations of triaxial test results

The experimental results of triaxial tests were used for numerical investigation of the anisotropic mechanical behaviour of the Tournemire argillite. The constitutive relationship was developed to describe the stress–strain behaviour of the argillite. In the simulation, we used an elastoplastic framework. In the elastic range, Hooke's law for transversely isotropic materials was used with five parameters obtained from the experimental results to account for the bedding. In the plastic range, a microstructure tensor approach (Pietruszczak & Mroz 2001) was adopted in order to take the anisotropic behaviour of the argillite into account. The material functions/parameters and results of the numerical simulation are presented.

### Constitutive relationship based on the microstructure tensor

In order to define the anisotropy parameter(s), the formulation employs a generalized loading vector that is defined as:
(1)
where *, j*=1, 2, 3, are the base vectors that specify the preferred material axes. Thus, the components of *L _{i}* represent the magnitudes of traction vectors on the planes normal to the principal material axes. Note that:
(2)
so that the traction moduli can be expressed as mixed invariants of the stress and microstructure-orientation tensors. The unit vector along

*L*is given by: (3) A microstructure tensor

_{i}*a*is now introduced that is a measure of material fabric. While different descriptors may be employed to quantify the fabric, the eigenvectors of this operator are said to be collinear with . The projection of the microstructure tensor on

_{ij}*l*becomes: (4)

_{i}Here, the scalar variable η, referred to as anisotropy parameter, specifies the effect of load orientation relative to material axes, and is defined as the ratio of the joint invariant of stress and microstructure tensor *a _{ij}*σ

*σ*

_{ik}*to the stress invariant σ*

_{jk}*σ*

_{ij}*. It is a homogeneous function of stress of the degree zero, so that the stress magnitude does not affect its value. Note that equation (4) can be expressed as: (5) where*

_{ij}*A*)

_{ij}=dev(a_{ij}*/*η

_{0}is a symmetrical traceless operator. Equation (5) may be generalized by considering higher-order tensors, i.e.: (6)

The above representation is rather complex in terms of implementation and/or identification. Therefore, it is convenient to use a simplified functional form obtained by replacing the higher-order tensors by dyadic products of second-order tensors, i.e.:
(7)
where *b*_{s} (*s*=1, 2, …) are constants.

In view of the consideration above, the failure function can be expressed in the general form:
(8)
where *I*_{1} is the first stress invariant, while *J*_{2}, *J*_{3} are the basic invariants of the stress deviator. As mentioned earlier, the parameter η is typically identified with a relevant strength descriptor, the value of which is then assumed to depend on the orientation of the sample relative to the direction of loading. Thus, the existing criteria can be easily extended to an anisotropic regime by assuming that the strength parameters vary according to equation (7). In this work, a linear form of equation (8) has been adopted, which corresponds to the well-known Mohr–Coulomb criterion, i.e.:
(9)

Here:
and
(10)
where θ is the Lode's angle, and ϕ and *c* are the angle of friction and cohesion, respectively.

An extension to the case of inherent anisotropy can be accomplished by assuming that the strength descriptors, in this case η_{f} and *C*, are orientation-dependent and have a representation analogous to that of equation (7). Note that, in the context of the equation (9), the parameter *C* is associated with a hydrostatic stress state. The latter is, in fact, invariant with respect to the orientation of the sample. Thus, the effects of anisotropy can be primarily attributed to a variation in the strength parameter η_{f}, i.e.:
(11)
The general plasticity formulation can be derived by assuming that the yield/loading surface in the form is consistent with equation (11), i.e.:
(12)
where is the deviatoric part of the plastic strain increment, *A* and ζ are the material parameters, and κ is the plastic distortion. According to the hardening rule, for κ→∞, there is , where ζ>1. The parameter ζ is introduced here in order to define the transition to localized deformation, which is assumed to occur at . Note that the latter equality implies that *f*=*F*, so that the conditions at failure are consistent with the Mohr–Coulomb criterion (equation 9). The plastic potential can be chosen as:
(13)
where η_{c}∝η_{f}, so that η_{c}=η_{c}(*l _{i}*). Following now the standard plasticity procedure (Pietruszczak 2010) – i.e. invoking the consistency condition

*df*=0 – the constitutive relationship can be expressed as: (14) where: where is the elastic compliance, the representation of which, once again, depends on the type of material anisotropy.

### Procedure for identification of material functions/parameters

The identification process employed here is based on the experimental data of a series of triaxial compression tests presented in the previous sections. The first step involves the specification of the conditions at failure which, in turn, requires the identification of material function η_{f} and the parameter *C*, as defined in equation (11). The specification of coefficients appearing in η_{f}=η_{f}(*l _{i}*) entails the information on conditions at failure in samples tested at arbitrary orientation β (β=1−θ) with respect to direction of loading. Consider, for this purpose, the response of the sample under axial compression at a confining pressure

*p*

_{0}. Refer the geometry of the sample to the coordinate system in which are the material axes (see Fig. 3). In this case, (15) where β defines the orientation of the bedding planes. Given the representation above, the function η

_{f}=η

_{f}(

*l*) (equation 11) reduces to: (16)

_{i}Note that in the absence of confinement, *p*_{0}=0, there is *l*_{2}=cosβ, so that the anisotropy parameter is an explicit function of the deposition angle β, i.e.:
(17)

Given now the stress parameters at failure, for each specific orientation of the sample, the anisotropy parameter can be determined – i.e. η_{f}=*q*/(*p*+*C*) – together with the corresponding value of *l*_{2} (equation 15). The results can then be plotted in the affine space to generate a set of data describing equation (15). It should be noted that the above procedure requires an estimate of the value of the constant *C*, which can also be obtained based on the results of standard triaxial tests.

The specification of strength parameters for Tournemire argillite employs the results plotted in Figures 7, 8, 9, 10. Figure 7 shows a linear Mohr–Coulomb approximation of the conditions at failure in samples tested at different orientations of the bedding planes (experimental data points are shown as symbols). Based on these results, the value of *C* was estimated to be *C*=10.6 MPa, which represents an orientation-average for tests conducted at β=0°, 30°, 45°, 60° and 90°. Given the value of the constant *C*, the spatial distribution of the anisotropy parameter η_{f}=η_{f}(*l _{i}*) was established by following the procedure described above.

Figure 8 shows the experimental results projected in the affine space . Again, the data incorporate the tests conducted at different confining pressures, viz. 0, 4 and 10 MPa, and different orientations of the bedding planes. The best-fit approximation is based on equation (16). In this case, the second-order approximation appears to be adequate and corresponds to the following set of coefficients: , *A*_{1}=0.17034 and *b*_{1}=5.4957.

Figure 9 shows the plot of η_{f} against the loading angle , which incorporates the parameters given above.

The next step of the identification process pertains to specification of the material parameters and functions that govern the deformation characteristics. The elastic properties were obtained based on the results of unconfined cyclic tests as follows: *E*_{1}=12.5 GPa, *E*_{2}=21.0 GPa, ν_{21}=0.16, ν_{13}=0.08 and *G*_{21}=4.57.

Again, these values are referred to the coordinate system associated with principal material axes, as shown in Figure 3.

The evaluation of the plastic properties requires the identification of the potential function, equation (13), which incorporates the parameter η_{c}. The latter defines the transition from plastic contraction to dilatancy and its value can be assessed by examining the respective volume change characteristics. For the Tournemire argillite, the results of experimental tests as described previously indicate that the transition to dilatancy occurs at deviatoric stress intensities that approach those associated with unstable strain-softening responses. Therefore, it may be adequate to approximate this material function as: η_{c}(*l _{i}*)≈0.99η

_{f}(

*l*).

_{i}The specification of the hardening function involves the identification of the material parameter *A*, namely equation (12). In order to accomplish this, the mechanical characteristics should be replotted in {κ,η=*q*/(*p*+*C*)} space. Note that constructing such characteristics requires the value of the elastic shear modulus so that the plastic deviatoric strain can be determined via the additivity postulate. The deformation properties in the plastic range are likely to be direction-dependent. Thus, the best-fit procedure should be repeated for different orientations of the sample in order to obtain a set of corresponding values of the parameter *A*. In order to maintain the simplicity of the formulation, it may then be sufficiently accurate to assume that *A* represents a material constant, the value of which is identified with the orientation average.

It needs to be pointed out here that the experimental results provided are not exactly adequate to follow the identification procedure. This stems from the fact that the deformation measurements comprise only the vertical and volumetric strain, which for samples tested at orientations other than β=0° (i.e. horizontal bedding planes) is not sufficient to uniquely define the value of the plastic distortion κ. Note that at β=90°, the lateral strain will not be the same in the radial directions. The deformation in different radial directions is, however, not measured during the laboratory tests. At the same time, for inclined samples, the shear strain will develop under a stress-controlled regime, as the material is anisotropic.

Given the restrictions mentioned above, parameter *A* was determined by a trial-and-error procedure that involved fitting the deviatoric characteristics *q* v. axial strain for samples with horizontal bedding planes. Figure 10 shows the plots corresponding to the confining pressures of 0, 4 and 10 MPa. It is important to point out that the unstable strain-softening characteristics reported in experimental tests are associated with localized deformation: that is, the formation of macrocracks. Therefore, they do not represent a material response and should be analysed within the context of a boundary-value problem. Thus, the characteristics reported in Figure 10 are restricted to a stable range associated with strain hardening. The experimental response is characterized by an abrupt transition to a localized mode. Within the current mathematical framework, the onset of localized deformation is associated with and is defined by means of parameter ζ, which appears in the hardening function (equation 12). In general, the constraint ζ>1 ensures that the transition occurs along the ascending branch of the deviatoric characteristic. The simulations reported in Figure 10 were conducted assuming ζ=1.25, and resulted in the value of *A*=0.0012. Thus, the material parameters governing the plastic response of the Tournemire argillite were taken as: , ζ=1.25 and *A*=0.0012.

Finally, it should be pointed out that the results of the triaxial tests for inclined samples need to be taken with caution. This stems from the fact that the samples will have tendency to distort under the applied load: however, in a triaxial set-up, such distortion is kinematically constrained by the presence of the loading platens. The latter will lead to a non-uniform stress distribution within the specimen, so that the results are representative of a boundary-value problem rather than a material test. This is, in fact, the main argument why, at the present stage, a simplified procedure for the identification of material parameters in the plastic range has been adopted.

### Numerical simulation

The results of numerical simulations with the finite-element code Abaqus are shown in Figures 11, 12, 13, 14, 15. The figures present both the deviatoric characteristics *q* v. axial strain (*E*_{11}), as well as the corresponding evolution of volumetric strain (*E*_{v}). It is again noted that the simulations were terminated at the instant the criterion for the onset of localization, defined as , was satisfied. As mentioned earlier, the subsequent unstable response (indicated by an arrow pointing downwards) constitutes a boundary-value problem and cannot be examined within the context of the point-integration scheme. It appears that, in spite of simplifications embedded in the identification procedure, the results of simulations are quite reasonable in terms of depicting the basic trends. This is certainly the case given the doubts concerning the validity of experimental data for samples tested at orientations other than those involving the horizontal bedding planes.

The most significant discrepancy between the experimental data and the numerical simulations occurs at β=90° (Fig. 15). Here, the conditions at failure are predicted quite accurately. However, the stiffness characteristics are markedly different. The reason for this discrepancy is probably because the strain-hardening effect before failure for β=90° (i.e. the axial loading is parallel to the bedding plane) is stronger than that for other β angles, as shown in the experimental results. Based on the hardening function represented by equation (12), a smaller value of *A* would reflect a stronger hardening effect, which means that parameter *A* might also be orientation-dependent. However, this hypothesis would need more experimental and numerical work to support it. Therefore, the predictions can be improved by accounting for the fact that the hardening characteristics are likely to be influenced by the orientation of the sample. This is illustrated in Figure 16, which shows the results for the simulations corresponding to a value of *A*=0.0002. It is evident that a significant improvement in terms of predictive abilities is achieved.

## Conclusion

The Tournemire argillite is an anisotropic material, and its strength and deformation behaviour are direction-dependent. Stress–strain relationships are non-linear, and failure occurs in a brittle manner with a sudden collapse of the material. The failure mode exhibits localization along distinct failure planes and also depends on the loading orientation.

The identification procedure for material functions and parameters was outlined, and the model was applied to simulate the set of triaxial tests performed at different levels of confining pressure and orientation of the bedding planes. It is demonstrated that the model adequately reproduces the anisotropy, the pre-peak stress–strain response, and the onset of material collapse in those tests.

- © 2017 The Author(s). Published by The Geological Society of London

**Gold Open Access:** This article is published under the terms of the CC-BY 3.0 license.