A gradient enhanced transversely isotropic damage plasticity model for rock ‐ formulation and comparison of different approaches

The nonlinear mechanical behavior of intact rock and rock mass is characterized by irreversible deformation, associated with strain hardening, strain softening and degradation of stiffness. In the case of stratified rock types like shales or phyllites, the mechanical behavior additionally depends on the loading direction with respect to the orientation of the principal material directions, resulting in inherent anisotropic material behavior, or in many cases more specifically in transversely isotropic material behavior. On the basis of a critical review of several approaches for modeling nonlinear anisotropic material behavior, proposed in the literature, the linear mapping of the Cauchy stress tensor into a fictitious isotropic configuration, originally proposed by Boehler [1], was identified as the best choice for the application in 3D finite element simulations. Within this framework, an isotropic damage plasticity model for rock, proposed by Unteregger et al. [2], extended to inherent transversely isotropic behavior is presented. It is implemented into a FE‐program by means of the return mapping algorithm and it is regularized by an over‐nonlocal implicit gradient enhancement for ensuring mesh‐insensitive results. For validation, numerical simulations of triaxial compression tests on Tournemire shale, documented in Niandou et al. [3], are performed at both the integration point level and the structural level. The performance of the constitutive model is assessed by its ability to predict the complex direction‐dependent mechanical behavior of transversely isotropic rock and the failure modes observed in triaxial compression tests with different confining pressures and different inclination angles of the stratification planes with respect to the direction of axial loading.


INTRODUCTION
Tunnel advance is a very demanding task in engineering practice. For the mathematical description of the material behavior of the surrounding rock mass, isotropic, linear-elastic perfectly-plastic material laws are frequently used in numerical simulations of tunneling. These material laws are based on the failure hypothesis according to Hoek et al. 4 . In the latter, discontinuities in a rock mass are considered in a simplified manner by reducing the strength parameters of intact rock by an empirical approach derived from a rock mass classification system. However, compared to the real mechanical behavior of rock mass, some essential effects of discontinuities are neglected by the respective simplifications. For example, plastic deformation in the pre-peak region of the stress-strain curves as well as strain softening in the post-peak domain, accompanied by degradation of the material stiffness, are ignored. In addition, linear-elastic behavior is incorrectly predicted for stress paths with dominating hydrostatic compression, although it is well known that in this case the material behavior depends on the hydrostatic pressure. These shortcomings significantly influence the prediction of the deformation and stress state in numerical simulations of tunnel advance. In order to overcome the mentioned shortcomings, an isotropic damage plasticity model for rock and rock mass, denoted as rock damage plasticity (RDP) model, was developed by Unteregger et al. 2 . In the RDP model the mentioned effects are considered by coupling plasticity theory with damage mechanics, the latter for describing the postpeak degradation of strength and stiffness. The flow potential and the formulations of hardening and softening, originally proposed by Grassl and Jirásek 5 , were modified for intact rock and extended for rock mass by introducing empirical reduction factors according to Hoek et al. 4 . For ensuring mesh-independent results in finite element simulations involving softening rock mass behavior, the RDP model was extended by a regularization technique based on an implicit gradient enhancement by Schreter et al. 6 .
The mentioned models for intact rock and rock mass are still restricted to rock types with isotropic material behavior. However, depending on the origin of a particular rock, primary and secondary rock structures result in inherent anisotropic material behavior of intact rock and rock mass both at the micro-and macro-level, which have a significant impact on the deformation and stress state in the respective intact rock or rock mass. Especially in a stratified rock mass, for example Innsbruck quartz phyllite or Tournemire shale, often the special cases of anisotropic material behavior, i.e., orthotropic or transversely isotropic material behavior, are encountered.
In addition to the inherent anisotropic behavior, directional dependencies can also occur due to load-induced microcracking commonly known as induced anisotropy. Since the focus of the present publication is on the representation of inherent anisotropy, an overview on different approaches for considering the latter in constitutive models is provided in the following.
The oldest approach for incorporating anisotropic material behavior in elasto-plastic models is based on the formulation of an anisotropic yield criterion. Departing from the failure criterion by von Mises, Hill 7 proposed the first anisotropic yield function for describing the influence of Lüder's bands in metals. For this purpose, the isotropic criterion was extended by six material parameters representing anisotropic behavior. Based on Hill's criterion, Pariseau 8 presented an extension for anisotropic rocks and Tsai and Wu 9 proposed a generalized anisotropic criterion for fiber-reinforced composites.
Another approach for modeling stratified materials relies on the ubiquitous joint concept, which is based on the Single Plane of Weakness Theory by Jaeger. 10 In the original ubiquitous joint (OUJ) model, introduced by Goodman, 11 joints are modeled as cracks with a given direction, distributed (or smeared) throughout the rock. In this way, a stratified material is substituted by a homogeneous anisotropic continuum. Failure may occur in the rock matrix, along the joints, or in both the rock matrix and the joints. Since in the OUJ model only strength anisotropy is considered, isotropic linear-elastic material behavior is assumed up to failure. To overcome this shortcoming, the OUJ model was further developed to the modified ubiquitous joint (MUJ) model by Ismael and Konietzky. 12 In the finite element program Plaxis 13 a jointed rock model is implemented following the ubiquitous joint concept. It is characterized by assuming transversely isotropic linear-elastic material behavior of the intact rock matrix and thus potential failure in the rock matrix is neglected. Failure is modeled by sliding of the rock along up to three stratification planes according to the Mohr-Coulomb failure criterion. However, the ubiquitous joint models presented so far, only consider linear-elastic perfectly-plastic material behavior, i.e., both hardening and softening behavior, which significantly influence the mechanical behavior, are neglected. For considering strain hardening in the pre-peak regime of the stress-strain relation and softening in the post-peak regime in ubiquitous joint models, the so-called subiquitous joint (SUJ) models, were developed. The Transubi model, presented by Ismael et al., 14 belongs to this category. In this model a bilinear Mohr-Coulomb yield surface with tension cut-off is employed for both the intact rock matrix and the stratification planes.
A further possibility of modeling inherent anisotropic behavior was published by Pietruszczak and Mroz 15 and Gao et al. 16 . It relies on the generalization of isotropic yield functions by so-called scalar anisotropy parameters. These parameters are postulated by means of a microstructure tensor, which takes into account the influence of the inclination of the stratification planes on the material parameters. Based on Pietruszczak's approach, Chen et al. 17,18 and Parisio et al. 19 proposed coupled elasto-plastic damage models.
A different continuum-based approach for the extension of isotropic failure criteria to anisotropic behavior was presented by Boehler, 1 wherein anisotropic behavior is described via a linear mapping tensor depending on the microstructure of the material. With the latter, an anisotropic failure criterion can be formulated in terms of mixed invariants of the stress tensor and a mapping tensor. In Boehler and Sawczuk, 20 this approach was used for an extension of the isotropic Mohr-Coulomb failure criterion for both, transversely isotropic and orthotropic materials and Crook et al. 21 developed an orthotropic elasto-plastic material model for weak shales by extending the Modified Cam Clay (MCC) model. Departing from the MCC model, Semnani et al. 22 also used a linear mapping tensor for deriving a transversely isotropic plasticity model to investigate the combined effects of inherent anisotropy and temperature on strain localizations.
An alternative approach for considering both, inherent and induced anisotropic material behavior, are micro-plane models. So far, they were mainly used for concrete (Bažant et al., 23 Zreid and Kaliske 24 ) and only few applications for anisotropic rock are available in the literature (Li et al. 25 ). Contrary to standard continuum models, the constitutive relations are formulated on integration planes of various orientations, representing the discretized surface of a unit hemisphere. For this purpose, the strain tensor is projected into the integration planes and the stress tensor is computed from the stress vectors acting on those integration planes by numerical integration. A similar concept is pursued by multi-laminate models (Zienkiewicz and Pande, 26 Schweiger et al. 27 ). In contrast to the micro-plane models, in the latter the stress tensor is projected into the integration planes. However, a general drawback of those models is the significantly higher computational effort compared to standard continuum based models, rendering both models computationally very expensive for large-scale simulations. For the sake of completeness, in addition to the mentioned continuum-based approaches, for modeling rock mass with single distinct stratification planes discontinuous approaches based on the discrete element method are mentioned, as discussed in Barton 28 and Wittke. 29 However, they require a high computational effort and they are usually limited by the lack of information on the structure of the discontinuities.
For the application of damage plasticity models in large-scale finite element simulations, a stable and robust formulation and implementation of the constitutive model is crucial for ensuring convergence. Hence, a damage plasticity model for transversely isotropic materials fulfilling these requirements will be presented in this paper.
For this purpose, first a review of popular approaches for considering anisotropic behavior of elasto-plastic materials is presented and main advantages and disadvantages are identified. Subsequently, the isotropic rock damage plasticity model by Unteregger et al. 2 is extended to transversely isotropic behavior. The capabilities of the new model are demonstrated by numerical simulations of triaxial compression tests on Tournemire Shale, documented by Niandou et al., 3 and a comparison of the numerical results with the respective experimental data.
The remainder of the paper is organized as follows: In Section 2, three popular approaches for describing transversely isotropic material behavior, i.e., the ubiquitous joint concept, the introduction of scalar anisotropy parameters and the linear mapping into a fictitious isotropic configuration, are critically reviewed and examined. To this end, a linear-elastic perfectly-plastic model with the simple Drucker-Prager yield criterion is extended accordingly. Based on this comparison in Section 3 the most promising approach is chosen for enhancing an existing isotropic rock damage plasticity model to account for transversely isotropic material behavior. In Section 4 the results of numerical simulations of triaxial compression tests on Tournemire shale on both the material point level and the structural level will be presented, thereby demonstrating the capabilities of the selected model. Finally, in Section 5 the paper is closed with a summary and a discussion on recommended future research activities.

A CRITICAL REVIEW ON MODELING TRANSVERSELY ISOTROPIC MATERIAL BEHAVIOR WITHIN THE FRAMEWORK OF PLASTICITY THEORY
This section is devoted to the review of three popular approaches for extending an existing isotropic elastoplastic model to a transversely isotropic elasto-plastic model, i.e., (i) the ubiquitous joint concept, (ii) the introduction of scalar anisotropy parameters and (iii) the projection of the Cauchy stress tensor into a mapped isotropic F I G U R E 1 Sketch of the global and local orthogonal bases, ( ) and ′( ) , = 1, 2, 3, respectively, for a transversely isotropic material, with ′(1) defining the axis of material symmetry and ′(2) and ′(3) the plane of isotropic material behavior configuration. To this end, the approaches (i)-(iii) will be briefly reviewed, and the key features will be subsequently assessed based on the extension of a linear-elastic perfectly-plastic model with the Drucker-Prager yield criterion to transversely isotropic behavior, serving as simple yet illustrative representatives of transversely isotropic elasto-plastic constitutive models.

Preliminaries
The focus is on transversely isotropic constitutive behavior based on coupling linear elasticity and plasticity theory. Therein, it is convenient to express the stress-strain relation in the principal material directions where the superscript prime (•) ′ designates variables formulated in the local coordinate system aligned with the principal material directions (cf. Figure 1). ′ describes the nominal Cauchy stress tensor, ′ the total strain tensor, ′ ( ) the elastic strain tensor and ′ ( ) the plastic strain tensor, considering the additive split of the total strain into an elastic and a plastic part. ℂ ′ ( ) denotes the fourth-order transversely isotropic elastic stiffness tensor, reading in Voigt-notation as with the Young's modulus 1 in ′ 1 -direction, the Poisson's ratio 12 and the shear modulus 12 in planes containing the axis of material symmetry ′(1) and the Young's modulus 2 and the Poisson's ratio 23 both in the ′ 2 ′ 3 -plane of isotropic material behavior.
The stress and the strain tensor, defined in the local coordinate system, are transformed into the global coordinate system according to = ′ (3a) and = ′ (3b) by means of the orthogonal coordinate transformation matrix = cos ∠( ′( ) , ( ) ), with the entries representing the direction cosines between the base vectors of the local and the global coordinate system. (3a) and (3b) are applied in an analogous manner for the transformation of the elastic and plastic strain tensors. Rearrangement of (3a) and insertion of (1) and (3a) into (3b) yields the stress-strain relation formulated with respect to the global coordinate system ℂ ( ) represents the fourth-order transversely isotropic elastic stiffness tensor in the global coordinate system. The stress tensor (4), referred to the global coordinate system, enters the equilibrium equations.

Ubiquitous joint models
In ubiquitous joint models, two different failure criteria are combined, one representing failure of the rock matrix ("nonsliding mode"), denoted as ( ) , and the other one failure along the stratification planes ("sliding mode"), denoted as ( ) . The latter is defined in terms of the traction vector ( ) = acting on an infinitesimal area of the stratification planes with the normal vector , leading to the incorporation of direction-dependent behavior. For introducing multiple failure criteria in ubiquitous joint models, the constitutive framework is extended to multi-surface plasticity as described by Simo. 30 According to the latter, the set of admissible yield functions is defined as with the subset of active yield functions According to the ubiqituous joint concept, an arbitrary number of stratification planes can be considered, whereas for transversely-isotropic materials one joint set is sufficient. For the evolution of the plastic strain, the non-associated flow rulė is employed witḣ( ) denoting the plastic strain rate, the subset of active yield criteria, ( ) ( ) the plastic potential function anḋ( ) the corresponding consistency parameter. In the following, the isotropic Drucker-Prager yield

Introduction of scalar anisotropy parameters
As a commonly employed representative for capturing inherent anisotropy in plasticity models by the direction-dependent formulation of strength-related material parameters, the framework by Pietruszczak and Mroz 15 is selected for this review. The layered fabric of a material is taken into account by a second-order microstructure tensor formulated with respect to the principal material axes, and 1 , 2 and 3 denoting the respective principal values. To take into account the influence of the loading direction with respect to the principal material directions, a generalized unit loading vector (cf. Figure 2) is proposed as Based on (12) and (13), a scalar anisotropy parameter is formulated as projection of the microstructure tensor in the direction of the generalized unit loading vector F I G U R E 2 Illustration of the generalized unit loading vector defined in the principal material directions It is convenient to split the microstructure tensor into a spherical part 0 = ∕3 and a deviatoric part , i.e., Inserting (15) into (14) yields the final representation of the scalar anisotropy parameter with̄= ∕ 0 . Thus, anisotropic behavior is incorporated into an existing isotropic model by replacing the directionindependent material parameter 0 by . Since in case of isotropic material behavior the principal values of the microstructure tensor are equal, the deviatoric part of the microstructure tensor vanishes and (16) degenerates to = 0 . For the present case of transversely isotropic material behavior with the stratification planes corresponding to the ′ 2 ′ 3 -plane according to Figure 1, the principal values 2 and 3 of the microstructure tensor are equal, i.e., 2 = 3 = . Thus, the deviatoric part of the microstructure tensor can be expressed in terms of one independent parameter , fulfilling the condition = 0: According to (16), only a monotonic increase or decrease of a material parameter in terms of the loading direction can be modeled. This contradicts, for example, the direction-dependent behavior of the uniaxial compressive strength of layered rocks. Hence, Pietruszczak and Mroz 31 proposed an extension of (16) by taking into account also higher-order terms up to degree + 1 with , = 1, … , , representing additional material parameters. For extending the Drucker-Prager yield criterion to transversely isotropic behavior, a direction-dependent friction parameter is derived from (18). Replacing and 0 in (18) by and 0 , and taking = 1, yields the transversely isotropic friction parameter It is characterized by two additional parameters, that is, in̄and , which replaces 1 in (18). Inserting (19) into (8) leads to the transversely isotropic Drucker-Prager yield function with 0 as an additional material parameter. In case of isotropic rock, the last two terms on the right hand side of (19) are zero and, thus, degenerates to 0 representing the isotropic friction parameter. Furthermore, it is assumed that the cohesion remains direction-independent. The transversely isotropic plastic potential function ( ) is obtained by omitting the last term on the right hand side in (20) and by replacing 0 by the dilatancy parameter 0 :

Stress projection into a fictitious isotropic configuration by means of a microstructure tensor
A further approach for considering transversely isotropic elasto-plastic material behavior, based on the notion of a fictitious isotropic configuration, was proposed by Boehler 1 and applied, for example, in Semnani et al., 22 Borja et al. 32 and Bryant. 33,34 Based on a fourth-order super-symmetric linear mapping tensor ℙ, the Cauchy stress tensor is projected to a mapped stress tensor * as follows: * = ℙ .
Variables defined in the fictitious isotropic configuration are designated by a superscript asterisk (•) * . ℙ is defined as incorporating transversely isotropic material behavior by means of the microstructure tensor with the components determined from the normal vector to the stratification planes (cf. Figure 1). An expression of the mapping tensor ℙ in terms of the Walpole algebra is provided in Appendix A. The constants 1 , 2 and 3 denote parameters for considering transversal isotropy. According to Semnani et al. 22 a constraint for the latter parameters is required for restricting the otherwise infinite number of parameter combinations leading to an identical mechanical response to a unique set. In the present paper this constraint is chosen as A detailed description of the constraint (25) together with an illustrative example is given in Appendix B. Following the assumption of Boehler, 1 the anisotropic yield function in terms of the stress tensor is replaced by an isotropic yield Example of applying the mapping tensor ℙ to a transversely isotropic material: (A) uniaxial stress 33 , (B) resulting multiaxial stress state * in the fictitious isotropic configuration, (C) resulting strain state function * in terms of the mapped stress tensor * : As shown, for example, in Betten, 35 the objectivity condition * ( ℙ , is met for arbitrary orthogonal transformation matriceŝ. In Figure 3, the stress tensor for a transversely isotropic continuum with stratification planes, defined by the normal vector = ⌊0 1∕ √ 2⌋ , subjected to uniaxial tension in the loading direction 3 , is projected into the fictitious isotropic configuration by the mapping tensor ℙ. Applying (22) to the respective stress tensor with the only non-zero component 33 results in the stress state shown in Figure 3B. In the mapped state, in addition to the two axial stress components, also shear stresses are present in the 2 3 -plane. However, since the normal vector lies within the 2 3plane, the components * 12 , * 13 and * 11 are zero. The resulting strain tensor, illustrated in Figure 3C, is obtained in two ways: (i) either from the original uniaxial stress state acting in the transversely isotropic material or (ii) from the mapped stress tensor * acting in an equivalent isotropic material.
According to (26) the transversely isotropic Drucker-Prager yield function is simply expressed as with the first invariant of the mapped stress tensor * 1 = ℙ and the deviatoric part of the mapped stress tensor * = ℙ , calculated by the deviatoric part of the mapping tensor ℙ = ℙ − 1 3 ℙ . Compared to the isotropic yield function, only the stress invariants are replaced by the respective invariants of * . Applying the chain rule to the non-associated plastic strain rate (7) giveṡ defining the plastic potential function of the transversely isotropic Drucker-Prager model.

Evaluation of the selected approaches for considering transversely isotropic material behavior
The performance of the reviewed approaches for considering transversely isotropic material behavior, i.e., (i) the Ubiquitous Joint Concept (cf. Section 2.2), (ii) Pietruszczak's approach (cf. Section 2.3) and (iii) Boehler's approach (cf. Section 2.4) is evaluated on the example of the presented extensions of the Drucker-Prager yield criterion. Based on this comparison the most suitable approach is selected for developing a gradient enhanced transversely isotropic rock damage plasticity model in Section 3. At first, the ability of the three approaches to reproduce the direction-dependent uniaxial compressive strength of a transversely isotropic rock is investigated on the basis of the experimental data taken from Pietruszczak and Mroz 31 for Tournemire shale, which can be viewed as a typical representative of layered rocks.
For this purpose, the material parameters for the three selected models, summarized in Table 1, are determined from the test data on Tournemire shale aiming at the best fit approximation of the direction-dependent uniaxial compressive strength using a nonlinear least-squares optimization algorithm. of the compressive strength is represented well at a slightly different corresponding inclination angle of = 45 • . In contrast, the Ubiquitous Joint Model yields an entirely different behavior: Whereas the experimentally determined minimum value of the uniaxial compressive strength is perfectly matched, for a wide range of inclination angles, i.e., from 0 • to 40 • and from 77 • to 90 • , the uniaxial compressive strength remains constant, i.e., direction-independent, corresponding to the isotropic strength of the rock matrix. The value of the constant compressive strength for a quite large variation of is determined from the material parameters of the rock matrix and y . The range for which the latter becomes active is influenced by both, the parameters of the rock matrix and the parameters and of the yield function (10) for sliding failure along the stratification planes, which only becomes active if high shear stresses in the stratification planes are present or high tensile stresses perpendicular to the latter occur. This leads to an overestimation of the uniaxial compressive strength for inclination angles between 15 • and 40 • by a factor of up to 1.4. In other words, the uniaxial compressive strength is overestimated for a wide range of inclination angles, for which the isotropic failure criterion for the rock matrix is dominant and the uniaxial compressive strength is independent of the inclination of the stratification planes. This shortcoming could be remedied by formulating an anisotropic failure criterion for the rock matrix, however at the expense of losing the simplicity of the Ubiquitous Joint Model. For the evaluation of the strength of a transversely isotropic rock, subjected to multi-axial stress states, predicted by the three selected modeling approaches, sections through the yield surfaces in the principal stress space obtained for different orientations of the stratification planes are plotted in Figure 5. In Figure 5 A-C, sections through the yield surface are plotted in the 1 2 -plane at 3 = 0 for the three selected inclination angles = 0 • , 30 • and 90 • with respect to the 1 3plane. For biaxial stress states the transversely isotropic yield surfaces for = 0 • and = 90 • are mirror-symmetric with respect to the equibiaxial-line 1 = 2 for all three approaches (cf. Figure 5A and C). In Figure 5 D-F, sections through the yield surface in the Rendulic plane are shown for the stratification planes inclined by the angle with respect to the 1 2 -plane. They represent axisymmetric stress states characterized by 1 = 2 , occuring, e.g., in conventional triaxial compression tests.
Summarizing, from Figure 5 follows that the Ubiquitous Joint Concept offers a simple approach for introducing transversely isotropic (and in general anisotropic) plastic behavior by adding failure criteria for considering sliding along stratification planes, complementing to the one of the isotropic matrix. Thus, a physical interpretation of the additionally introduced parameters and consequently the direct measurement of the latter in experiments is possible. Nevertheless, from Figure 4 to Figure 5 follows that for a wide range of stress states the failure criterion of the rock matrix is dominant, leading to an over-or underestimation of the direction-dependent strength.
One of the main advantages of Pietruszczak's approach is its simplicity: An isotropic yield criterion can be used for introducing transverse isotropy with the only difference that the material parameters are replaced by direction-dependent ones. Furthermore, for certain sets of material parameters, a direct fit of the anisotropic distribution function might be possible. However, according to Figure 5 some regions of the yield surface exhibit non-convexity, which is traced back to the unintended stress-dependence in addition to the direction-dependence of the material parameters, which is also discussed in Fuchs. 36 Thus, in view of applying a transversely isotropic constitutive model in finite element simulations, in which arbitrary stress paths may occur, this approach is questionable.
According to Boehler's approach, anisotropic behavior is considered by projecting the physical stress into a mapped stress, entering the isotropic failure criterion. This way, a direction-dependent strength for arbitrary multiaxial stress states is obtained. Two additional independent material parameters 1 and 3 are introduced for describing the influence of transversely isotropic behavior. However, they cannot be directly measured, they rather have to be determined by means of inverse parameter identification, explained in detail in Section 3.4.
Hence, for extending an existing isotropic yield criterion to transverse isotropy within the theory of plasticity, Boehler's approach, based on the projection of the Cauchy stress to a mapped stress, is considered superior to both the Ubiquitous Joint Model and Pietruszczak's approach. Thus, the gradient enhanced transversely isotropic rock damage plasticity model, to be presented in the following section, will be based on Boehler's approach.

TRANSVERSELY ISOTROPIC ROCK DAMAGE PLASTICITY MODEL
For modeling both the distinct direction-dependent and the highly nonlinear material behavior of layered rock, in the following a transversely isotropic rock damage plasticity (TI-RDP) model is proposed. It is based on a combination of (i) the isotropic rock damage plasticity model by Unteregger et al. 2 for representing the highly nonlinear material behavior of isotropic intact rock up to failure, and (ii) the approach by Boehler 1 (cf. Section 2.4) for capturing the transversely isotropic response.

Stress-strain relation
The TI-RDP model is formulated based on the theory of plasticity in the effective stress space combined with the theory of continuum damage mechanics. To this end, (4) is reformulated by replacing the Cauchy stress tensor by the effective stress tensor̄(force per undamaged area) leading to the effective stress-strain relation expressed in terms of the global coordinate system̄= The nominal Cauchy stress tensor (force per total area) is linked to the effective stress tensor by with the scalar isotropic damage variable ranging from 0 (undamaged material) to 1 (completely damaged material).

Plasticity formulation
For describing transversely isotropic elasto-plastic behavior, by analogy to (22) the effective stress tensor̄is projected to a mapped effective stress tensor̄ * = ℙ ∶̄. Hence, the Haigh-Westergaard coordinates * m =̄ * 1 ∕3 ,̄ * = √ 2̄ * 2 ,̄ * = 1 3 arccos of the mapped stress tensor̄ * , i.e., the mean stress̄ * m , the deviatoric radius̄ * and the Lode anglē * replace the respective quantities in the yield function of the original isotropic RDP model, proposed by Unteregger et al. 2 . Thus, according to (26) The equivalent isotropic yield surface exhibits a cap for predominant hydrostatic loading of the mapped stress tensor and evolves depending on a stress-like hardening variable 0 ≤ h ≤ 1. The parameters * cu and * 0 represent the equivalent isotropic uniaxial compressive strength and the friction parameter, respectively, and (̄ * , ) designates the Willam-Warnke function 37 for controlling the shape of the yield function in deviatoric planes with the eccentricity parameter ranging between 0.5 (triangular shape) and < ≤ 1 (circular shape) for ensuring a convex shape. In case of uniaxial compression perpendicular to the stratification planes, characterized in Voigt-notation by * = [− 0 0 0 0 0] , the equivalent isotropic yield surface reduces to * = ( ∕ * cu ) 2 − 2 h . It can be seen that h controls the uniaxial compressive stress normalized by the uniaxial compressive strength.
The evolution of the plastic strain is governed by the non-associated flow rule (29) with the non-associated plastic potential function of the original isotropic rock damage plasticity model, however, formulated in terms of the mapped effective stress tensor: For obtaining a realistic volumetric plastic strain response for rock, two equivalent isotropic dilatancy parameters * 1 and * 2 are introduced. Considering (25), * 1 can be determined from compression tests, loaded perpendicular to the stratification planes and triaxial compression tests, loaded in the principal material directions described in Unteregger, 38 whereas * 2 is a dependent parameter. The stress-like hardening variable with * cy denoting the equivalent isotropic compressive yield stress, is formulated in terms of the strain-like internal variable p . The initial yield surface is obtained for p = 0, with increasing values of p the yield surface is expanding until at p = 1 the ultimate strength surface is attained, corresponding to the smooth version of the Hoek-Brown failure criterion for intact rock, proposed by Menétrey and Willam. 39 In Figure 6 the evolution of the yield surface according to (34) is illustrated during strain hardening in the principal effective stress space for three different inclination angles of the stratification planes = 0 • , 45 • and 90 • , respectively. Since the stratification planes are rotated by around the 2 -axis, the uniaxial strength for loading in 2 -direction does not depend on . For = 0 • , the uniaxial strength in the 1 2 -plane is constant for arbitrarily chosen directions in this plane, as the latter represents the stratification planes. In contrast, the uniaxial strength in 3 -direction, perpendicular to the stratification planes, differs from the one in the stratification planes. Rotating the transversely isotropic rock by 90 • around the 2 -axis results in the yield surface for = 90 • . For = 45 • , the axes 1 and 3 are equally inclined with respect to the stratification planes, hence, the respective uniaxial strengths for loading in the directions of the respective axes are equal but differ from the uniaxial strength for loading in 2 -direction. Note that due to the transversely isotropic behavior the apex of the initial yield surface is not located on the hydrostatic axis anymore.
In (36) p follows from the hardening laẇp with the hardening function expressed in terms of the mapped stress tensor. It is obtained by replacing the effective stress in the isotropic RDP model by the mapped effective stress and controls the degree of strength mobilization with increasing plastic strains represented by the norm of the derivative of the plastic potential function with respect to the stress tensor. The term in the parentheses in (38) represents the Lode angle dependence factor f . For mapped stress states located at the compressive meridian (̄ * = ∕3), f = 1 and for mapped stress states located at the tensile meridian (̄ * = 0), f = 4. In case of mapped hydrostatic stress states (̄ * = 0) the dependence on the Lode angle vanishes, i.e., f = 1. The equivalent hardening ductility measure * controls the ductility in the hardening regime. h , h , h and h are independent hardening parameters and are parameters, depending on the former. The transition point which depends on the orientation of the stratification planes and the effective stress, ensures a more ductile pre-peak response in compression than in tension similarly to the original isotropic RDP model by Unteregger et al. 2 h depends on the model parameter 0 < h ≤ 1∕3 and, particularly for the TI-RDP model, on a modified mapping tensor for projecting the effective stress tensor into a mapped stress tensor ℍ ∶̄.
In (42) h is an additional anisotropy parameter. The modified projection tensor ℍ governs the direction-dependence of the ductility measure * h , which is required for replicating a realistic hardening behavior observed in experiments on transversely isotropic rocks. For the special case of isotropic material behavior h = 1, ℍ ∶̄degenerates tōyielding the original formulation for the isotropic RDP model.
In Figure 7 (left), the impact of the hardening parameter h in (42) on the equivalent isotropic hardening ductility measure * h is illustrated in terms of the inclination angle of the stratification planes for a given uniaxial stress̄= [0 0 − 25 MPa 0 0 0] , serving as a representative for an arbitrary compressive stress̄3 3 .
For h = 0.05 the inverse of the hardening ductility measure increases monotonically with increasing inclination angles , whereas for h = 1.0 a constant hardening ductility measure is obtained, i.e., * h (ℍ ∶̄) −1 is independent of the orientation of the material. For h = 2.0, * h (ℍ ∶̄) −1 decreases with increasing . In Figure 7 (right) the impact of h on the stress-strain relation is depicted for a uniaxial compression test, loaded in 3 -direction, for three different values of the inclination angle .
For inclination angles of ≠ 0 • , the influence of values h < 1 on the stress -strain curves is manifested in a more brittle hardening behavior and a more ductile hardening behavior for h > 1. It is also shown, that due to the formulation of a modified mapping tensor ℍ the stress -strain curve for loading perpendicular to the stratification planes, i.e., = 0 • , is not affected by the parameter h . Note that the uniaxial compressive strength is not influenced by h and its dependence on the loading direction with respect to the stratification planes is only controlled by the parameters 1 and 3 .

Damage formulation
Damage is initiated when the stress-like hardening variable h ( p ) attains its maximum value of 1. The damage evolution law of the TI-RDP model is assumed as an exponential relation controlled by the strain-like softening variable d and the equivalent isotropic softening modulus * f . The rate of the strainlike softening variablėd is given in terms of the volumetric plastic strain ratė( ) vol = tr(̇( ) ) and the softening ductility measure s . According to (29), the plastic strain rate is dependent on the mapping tensor ℙ. Hence, the evolution of the strain-like softening variable d and consequently the evolution of the damage parameter is influenced by the direction-dependent material behavior. The softening ductility measure, reading as A detailed description of the numerical implementation of the TI-RDP model, i.e., the stress-update algorithm, is given in Appendix C.

Proposal for determining the material parameters from experimental data
For determining the parameters entering the plasticity formulation, the global effective stress tensor̄is transformed into the local coordinate system aligned with the principal material directions as shown in Figure 1. Thus, the mapping (22) together with (23) between the local stress tensor̄′ and the mapped stress tensor is expressed in Voigt-notation as i.e., the mapping tensor ℙ boils down to a diagonal tensor. Five different loading scenarios are now considered: (1) For uniaxial loading perpendicular to the stratification planes, i.e., uniaxial loading in ′ 1 -direction, the only nonzero effective stress component is̄′ 11 . According to (48) together with the constraint (25) the only non-zero mapped effective stress component is̄ * 11 =̄′ 11 . Hence, in the limit state of̄ * 11 attaining the equivalent isotropic uniaxial compressive strength * cu , the latter corresponds to the uniaxial compressive strength (1) cu of a layered rock, loaded perpendicular to the stratification planes, with the superscript (1) indicating the loading direction, i.e., from̄ * 11 =̄′ 11 follows * cu = (1) cu . The same holds for the equivalent isotropic uniaxial compressive yield stress * cy = (1) cy as well as the equivalent isotropic uniaxial tensile strength * tu = (1) tu .
(2) For uniaxial loading parallel to the stratification planes, let us consider, e.g., loading in the direction of the ′ 2 -axis. The effective stress component̄′ 22 yields the only non-zero mapped effective stress component̄ * 22 = 1̄′ 22 . Thus, the limit state of attaining the uniaxial compressive strength parallel to the stratification planes is characterized bȳ ′ 22 = (2) cu with the superscript (2) denoting an arbitrarily chosen direction in the isotropic plane. Hence, the equivalent isotropic uniaxial compressive strength is obtained as * cu = 1 (2) cu . Combining the equivalent isotropic uniaxial compressive strength * cu = (1) cu from the previous loading scenario (1) with the present one yields the ansiotropy parameter Thus, 1 represents the ratio of the uniaxial compressive strengths in directions perpendicular and parallel to the stratification planes. (3) For simple shearing in the isotropic stratification planes (i.e., in the ′ 2 ′ 3 -plane), characterized by the only non-zero effective stress component̄′ 23 , the mapped stress component̄ * 23 = 1̄′ 23 is obtained from (48). The limit state is described by attaining the shear strength̄′ 23 = (23) sh in the respective plane, yielding the equivalent isotropic shear strength * sh = 1 (23) sh . (4) For simple shearing in planes perpendicular to the stratification planes let us consider, e.g., shearing in the ′ 1 ′ 2 -plane. The only non-zero effective stress component is̄′ 12 , which according to (48) is related to the mapped effective stress component̄ * 12 = ( 1 + 3 ∕2)̄′ 12 . The limit state is obtained by attaining the shear strength̄′ 12 = (12) sh in the respective plane, corresponding to the equivalent isotropic shear strength * sh = ( 1 + 3 ∕2) (12) sh . Combining the equivalent isotropic shear strength * sh = 1 (23) sh from the previous loading scenario (3) with the present one yields the anisotropy parameter Finally, with 1 and 3 at hand, 2 is obtained from the constraint (25) as with (23) bu denoting the biaxial compressive strength in the isotropic stratification planes.
Summarizing, the parameters * cu , * tu , * bu , * cy , 1 and 3 can be calibrated by means of experiments comprising uniaxial compression and tension tests as well as simple shear tests, characterized by loading directions perpendicular and parallel to the stratification planes and biaxial compression tests with biaxial loading in the stratification planes. According to the formulations for the eccentricity parameter and the isotropic friction parameter 0 , proposed by Grassl and Jirásek, 5 the eccentricity parameter of the equivalent isotropic yield surface can be calculated as and the equivalent isotropic friction parameter can be determined as * for the TI-RDP model. In the absence of respective experimental data, i.e., cu , tu and (23) bu , Unteregger et al. 2 suggested a default value for = 0.51. In this case, the equivalent isotropic friction parameter * 0 needs to be determined via inverse parameter identification.
For calibrating the equivalent isotropic softening modulus * f loading scenario (1) is considered with uniaxial tension perpendicular to the stratification planes. In this case, due to the formulation of the plastic potential function (35), the only non-zero plastic strain rate component controlling the evolution of d in (44) iṡ′ ( ) 11 . Thus, together with the constraint (25), the equivalent softening modulus * f corresponds to the softening modulus for uniaxial loading perpendicular to the stratification planes (1) f . The latter can be calibrated by the mode I fracture energy (1) fI for loading perpendicular to the stratification planes.

NUMERICAL STUDY
The aim of the numerical study is to evaluate the performance of the transversely isotropic rock damage plasticity model to predict the direction-dependent, highly nonlinear mechanical behavior of layered intact rock. For this purpose triaxial compression tests, documented in the paper by Niandou et al. 3 on specimens of Tournemire shale, are considered. In the respective experimental program cylindrical specimens, characterized by a height of 75 mm and a diameter of 37 mm, were mounted in a triaxial test apparatus. The specimens were tested at different confining pressures 0 and inclination angles of the stratification planes in a displacement controlled manner by monotonically increasing the axial displacement at the top surface of the specimens. Most of the experimentally determined stress -strain curves for confining pressures of 0 = 5, 10, 30, 40 and 50 MPa and inclination angles of = 0 • , 45 • and 90 • are documented in Niandou et al. 3 The capability of the transversely isotropic rock damage plasticity model to reproduce the mechanical behavior of the triaxial compression tests is evaluated by (i) single element tests and (ii) three-dimensional structural finite element simulations. The latter were carried out for capturing the complex, non-uniform deformation pattern of the specimens due to strain localization in the post-peak regime.

Parameter identification
The experimental results for two different confining pressures, i.e., 0 = 5 MPa and 0 = 40 MPa, and three different inclination angles of the stratification planes, i.e., = 0 • , 45 • and 90 • , are used for the identification of material parameters. The remaining experimental results are used for validation. An overview of the experimental tests is given in  1. Based on the ultimate strength surface, i.e., the yield surface (34) at h = 1, the parameters (1) cu and * 0 as well as the two independent anisotropy parameters 1 and 3 are identified for a best fit with measured strength values using a least-squares optimization algorithm. The distribution of the maximum deviatoric stress dev = 33 − 0 in terms of the confining pressure 0 and the inclination angle , computed from the already identified material parameters, is compared in Figure 8 with the experimental data by Niandou et al. 3 . The experimental data, plotted as orange dots, was used for parameter identification, whereas the test data, plotted as blue dots, served for validation. 2. The elastic parameters 1 , 2 , 12 , 23 and 12 and the parameters (1) cy , * 1 , h , h , h and h are calibrated for a best fit with the measured stress-strain curves from the tests conducted at the confining pressures 0 = 5 MPa and 0 = 40 MPa, again using a least-squares optimization approach. For the remaining parameters h , h , s , s and the default values suggested by Unteregger et al. 2 are employed. Figure 9 contains a comparison of the computed stress-strain curves, used for paramteter identification, with the respective experimental data.

Single element tests
The TI-RDP model is validated by computing the stress-strain curves for the tests, which were not used for parameter identification, on the basis of the material parameters, summarized in Table 3. To this end, a single 3D finite element with linear shape functions, resulting in a uniform stress and strain state throughout the specimen, is used. Figure 10 depicts a comparison of the computed stress-strain diagrams with the respective experimental data.

F I G U R E 9
Stress-strain curves (deviatoric stress versus axial strain) for the triaxial compression tests, used for calibrating the material parameters: experimental results by Niandou et al. 3 (discrete data points) and numerical results (continuous lines) From Figure 10 follows that by the TI-RDP model the triaxial compressive strength is predicted well for different confining pressures and loading directions with respect to the stratification planes and the computed stress -strain curves are in good agreement with the experimental data.
As a non-uniform strain distribution, which typically develops in the post-peak regime, cannot be captured by the employed single finite element with a uniform strain distribution, in Figures 9 and 10 the computed stress-strain curves are only shown for the pre-peak regime.

Three-dimensional finite element simulations
The application of the TI-RDP model for predicting the complex deformation states of specimens in the post-peak region of tests due to the formation of shear bands, is demonstrated by 3D finite element simulations of the triaxial compression tests, carried out by Niandou et al. 3 . The geometry of the FE-model together with the prescribed boundary and loading conditions is illustrated in Figure 11. Since the stratification planes are inclined with respect to the horizontal plane by the angle , defined as a rotation around the 2 -axis, failure modes symmetric with respect to the 1 3 -plane are expected. Hence, the respective symmetry is exploited by discretizing only one half of the specimen by a structured mesh, consisting of reduced integrated 20-node quadrilateral finite elements with an element size of 2.5 mm, and applying the respective boundary conditions in the plane of symmetry. Following the test procedure, in the first step of the numerical simulations the confining pressure is applied, followed by application of the axial loading in the subsequent steps by imposing a monotonically increasing uniform vertical displacement of the top surface of the specimen. To trigger localized failure, slightly reduced values for (1) cy and (1) cu are assigned to some elements at the center of the specimen, as indicated in Figure 11. The respective load-displacement curves are computed in a mesh-insensitive manner by means of an implicit overnonlocal gradient enhanced regularization of the TI-RDP model, following the approach described in Schreter et al. 6 for the isotropic damage plasticity rock model. For this purpose, the internal strain-like softening variable d is related to its nonlocal counterpart̄d by the screened Poisson equation with denoting an internal length scale parameter. The discretized weak form of (55) together with the discretized weak form of the equilibrium equations results in a coupled system of equations, which is solved for the nodal displacements and the nodal values of the nonlocal field̄d (cf. Schreter et al. 6 ).
In order to overcome the spurios localization of the plastic strain, an over-nonlocal softening variablêd is defined according to Vermeer and Brinkgreve 40 aŝd with denoting a weighting parameter for controlling the regularization. In the context of an over-nonlocal approach > 1 is chosen. The damage variable is finally calculated by replacing the local variable d in (43) by the over-nonlocal variablêd For the FE -simulations the additional parameters for the over-nonlocal implicit gradient enhanced formulation are chosen to = 5 mm, = 1.05 and * f = 4.75 ⋅ 10 −4 . Apart from the notion of a mapped stress tensor in the TI-RDP model, the theoretical framework for the damage formulation of the proposed model is identical to the one of the isotropic RDP model (Schreter et al. 6 ). Therein, it has been already confirmed that complex shear failure modes can be predicted in a mesh-insensitive manner using the gradient enhanced formulation. For the sake of conciseness, an additional mesh refinement study based on the TI-RDP model is omitted.
In Figure 12 the computed load-displacement curves are plotted for the specimens with the inclination angles of the stratification planes = 0 • , 45 • and 90 • and the three confining pressures 0 = 10, 30 and 50 MPa. In the respective diagrams Δ denotes the axial compression force, acting on a specimen in addition to the confining pressure. The predicted influence of the confining pressure on the structural behavior of the specimens is summarized as follows: with increasing confining pressure the peak value of the axial force and the corresponding displacement are increasing. A similar trend applies to the post-peak region, for which an increasingly ductile behavior is predicted for higher values of the confining pressure. The influence of the inclination angle of the stratification planes on the load-displacement curves is manifested by an increase of the initial stiffness with increasing values of the inclination angle, accompanied by an increasingly brittle behavior of the specimens.
For all specimens the peak values of the deviatoric stress, computed from the peak values of the axial force Δ and the initial cross -sectional area 0 = 1075 mm 2 of the specimen, are in good agreement with the respective peak values, documented in Niandou et al. 3 and plotted in Figures 9 and 10.
The axial strain documented in Niandou et al. 3 was measured by strain gauges attached to the surface of the specimens halfway between the top and bottom surface. As the strain is localizing in a shear band, a non-uniform strain distribution develops. Thus, the top displacement of a specimen cannot be calculated from a single measured strain value and compared to the computed one. The computed deformation patterns of the specimens, plotted in Figure 12 for an advanced stage of localization of the displacements, can be compared to the failure modes sketched in Figure 20 in Niandou et al. 3 . Good agreement can be stated for computed and observed failure modes.

CONCLUSION AND OUTLOOK
In this paper, the extension of an isotropic rock damage plasticity (RDP )model to transversely isotropic material behavior was presented. It is denoted as transversely isotropic rock damage plasticity (TI-RDP) model. After a review of approaches for modeling elasto-plastic anisotropic material behavior a detailed comparison of three popular approaches for considering transversely isotropic material behavior, comprising (i) the Ubiquitous Joint Concept, (ii) the introduction of scalar anisotropy parameters and (iii) the linear mapping of the Cauchy stress tensor into a fictitious isotropic configuration revealed the superiority of the latter approach which originally was proposed by Boehler. 1 In this approach shortcomings of the other two approaches are not present and, moreover, it is well suited for a stable and robust implementation into a FE-program on the basis of a return mapping algorithm. Hence, a quadratic rate of asymptotic convergence at the structural level is ensured, which is essential for applications in large-scale finite element simulations.
Two approaches for determining the material parameters of the TI-RDP model were proposed. In the first one they are determined from standard tests on rock samples with different orientations of the stratification planes with respect to the loading direction. Alternatively, inverse parameter identification procedures can be employed. The latter approach was demonstrated on the example of test results on Tournemire shale proposed by Niandou et al. 3 , as some of the required results from standard tests were not available.
The TI-RDP model was validated by numerical simulations of triaxial compression tests on Tournemire shale, which were not used for parameter identification. They comprised tests on rock samples at different confining pressures and different orientations of the stratification planes with respect to the direction of axial loading (Niandou et al. 3 ). In a first step the stress-strain curves were simulated numerically up to the peak stress at the integration point level. As in the post-peak region strain softening, resulting in the formation of shear bands, is observed, 3D FE-simulations of the rock samples were performed at the structural level. For ensuring mesh-insensitive results an implicit over-nonlocal gradient enhanced regularization of the TI-RDP model was employed according to Schreter et al. 6 .
The numerical simulations at the integration point level and at the structural level confirmed the capability of the TI-RDP model of representing the nonlinear behavior of transversely isotropic rock both in the pre-peak and post-peak domain. The predicted response of the triaxial tests up to failure corresponds well to the observed behavior in the tests, documented in Niandou et al. 3 . In particular, for advanced stages of localization of the displacements good agreement of the predicted deformation patterns of the specimens with the failure modes reported in Niandou et al. 3 can be stated. This holds for different confining pressures as well as different inclination angles of the stratification planes with respect to the direction of axial loading.
Summarizing, the transversely isotropic rock damage plasticity model is well suited for modeling the inherent directiondependent highly nonlinear behavior of stratified rock types in a realistic manner. The successful application to 3D numerical simulations of triaxial compression tests pave the way to the application of the constitutive model in 2D and large-scale 3D finite element simulations of tunnel advance in transversely isotropic rock mass and numerical simulations of anchor pull-out tests in stratified rock mass.

A C K N O W L E D G M E N T S
The authors express their gratitude to Dr. Benjamin Fuchs, who investigated the approach proposed by Pietruszczak and Mroz in his doctoral thesis.
Partial financial support of this research project by the Tyrolean Science Fund (F.18712) is greatfully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.  It can be seen, that the first term is affected by the material parameters = ( 1 + 2 + 3 ) and * cu . That is why the same material response can be obtained for different combinations of the parameters 1 , 2 , 3 , and * cu . Thus, is set to 1 by the constraint (25) and the axial behavior perpendicular to the stratification planes corresponds to the fictitious isotropic behavior, described by the parameters of the yield surface. Hence, the equivalent isotropic uniaxial compressive strength * cu corresponds to the uniaxial compressive strength (1) cu perpendicular to the stratification planes. In case of uniaxial compression parallel to the stratification planes, the yield function simplifies to * = cu according to (49) it can be seen, that in this case the uniaxial compressive strength parallel to the stratification planes (2) cu enters the yield function.