A plasticity model for discontinua

This article is concerned with the development and application of a simple continuum theory for rocks that may contain both randomly as well as preferentially oriented plane discontinuity surfaces. The theory stipulates that displacement discontinuities are independently activated on these surfaces as soon as an appropriate yield criterion is fulfilled; these displacement jumps account for the irreversible, ‘plastic’ part of the bulk deformation. In stress space, the critical conditions for the activation of discontinuous slip or opening displacements define an overall yield envelope that could be initially anisotropic, reflecting for example a weakness of certain orientations due to pre-existing joint sets. For the yield conditions studied in this paper, essentially a Coulomb-type friction law and a simple fracture opening condition, the inferred stress-strain response under typical triaxial loading conditions reveals the sensitivity of the two discontinuous deformation modes to the confining pressure. The incipient growth of a geological fold in such a material is modelled as a problem of plate bending. The slip- and opening-modes of deformation are found to be activated typically in the fold intrados and extrados, respectively. Under certain conditions, both modes will be activated simultaneously at the same locality and contribute to the total deformation. Field observations on a well exposed sandstone anticline are reported here, which support this conclusion. The present plasticity model for discontinua can clearly be explored in more detail for realistic distributions of faults and joints taken from field observations. It could also be improved in various ways in its description of the underlying deformation mechanisms. Apart from its interest as a mechanical constitutive model, it can also serve as a point of departure for studies of stress-sensitive, anisotropic permeability distributions in fractured formations.


Introduction
Natural fracture systems can provide a record of the history of stress, pore fluid pressure, or tectonic deformation, and an understanding of their genesis can ex plain and help to predict large-scale permeability trends in certain hydrocarbon reservoirs. Mechanical models that relate properties of natural fracture systems, such as fracture orientations or densities, to the tectonic stress and deforma tion history are therefore of fundamental interest in this context. It is hoped 1 that the simple elasto-plasticity law presented in the following may illustrate the potential of continuum theory to deal with this problem. Various classifications of natural fractures in terms of their origin have been proposed in the literature (see, e. g. , Nelson, 1985;Price & Cosgrove, 1990). A simple scheme distinguishes between tectonic fractures and regional or system atic joints. The former owe their name to the fact that they tend to accommodate tectonic deformation while being generated. Their growth results in complex 3-D geometries which often reveal an interaction of neighboring discontinuities. Lithological layering with sharp rheological contrasts and slip along bedding planes is also a key to the understanding of these tectonic fracture patterns. The origin of systematic joints, the second class of fractures, is still debated. These fractures are very regular in orientation and may cover vast areas. They are formed without significant tectonic deformation although their orientation is often considered as a marker of the paleo-stress state (Rawnsley et al., 1992;Dunne & Hancock, 1993;Petit et al., this Volume). They may also participate in tectonic deformation postdating their genesis. An ideal predictive mechani cal model should include the salient features of these discontinuities of different type and provide estimates of fracture densities and orientation based upon a knowledge of the overall deformation history.
In this paper, the link between pervasive fractures and permanent deforma tion is explored theoretically by means of a simple elasto-plasticity model in which the permanent deformation results from discontinuous (opening or slip) displacements along continuously distributed planes of weakness. The model rep resents an application to pervasively faulted rocks of ideas that may be traced back at least to the early work of Batdorf & Budiansky (1949) on the proper ties of polycrystalline materials that deform by crystal-plastic slip within single crystals. Certain of its features are also reminiscent of Reches' (1983) three dimensional model of faulted rock. It is assumed that there exists an elementary volume of the rock mass under consideration, for which a macro-scale constitutive model can be meaningfully developed. This volume element contains an infinite number of potential discontinuity surfaces with random orientation prior to de formation. Phenomenological laws are introduced that specify the resistance of a discontinuity to slip and opening: A Coulomb friction law with cohesive harden ing, to model resistance to slip, and an independent resistance law for irreversible opening displacements, which acts as a cut-off for Coulomb frictional behaviour at low confining pressures and will allow opening to occur instead of or concur rently with slip. The phenomenological laws that characterize the response of a discontinuity could in fact represent distinct physical micromechanisms, includ ing fracture nucleation and growth behaviour. Rather than giving an explicit treatment, these various effects are lumped together in the present simplified phenomenological description with the uniform elastic macro-scale response of the surroundings to localized slip or opening displacements. The micromechani cal parameters of the model could exhibit a dependence on the orientation of the discontinuities, such that the overall resistance of the rock mass to deformation would be weaker, in a sense to be defined later, in certain directions, reflecting the presence of pre-existing systematic joint sets. However, in all cases the den sity of these discontinuities in a representative vohm1e element is such that the macro-scale deformation will be continuous in space.
The idea to construct the macro-scale deformation by the superposition of micro-scale contributions from discontinuous opening or slip displacements along fractures is not new and has been pursued previously in various branches of me chanics including mechanical and civil engineering, geotechnics and geology. For this reason, the short review which follows is certainly far from complete. Oer tel (1965) suggested that the deformation observed on his soft clay models is accommodated mainly by slip along four families having orthorhombic symme try of parallel, closely spaced surfaces. The model of damage in concrete and rocks under tensile load proposed by Bazant & Oh (1985) differs from Batdorf and Budiansky's theory, apart from the choice of micro mechanism, mainly by the proposition of a uniform strain over the representative volume element. The tun nel and slope stability analysis of Zienkiewicz & Pande (1977) relies on the idea that the failure of rock masses is due to slip and opening of inherited fracture sets that act as planes of weakness. The tectonic origin of fracture sets and their ori entation with respect to the principal stress directions was addressed by Wallace (1951) & Bott (1959). A plane of fracture initiates from a heterogeneity which is then described by a micromechanism such as a penny-shape crack along which frictional sliding is accommodated (Kachanov, 1982a). Local tensile stresses can result in the propagation of branched cracks which is often the micromechanism invoked to model rock dilatancy (Kachanov, 1982b). The macroscopic stress strain relations based on such mechanisms are derived using, for example, a self-consistent scheme (Horii & Nemat-Nasser, 1983). Sliding along randomly oriented microcracks is often the only relevant mechanisms considered for rocks under overall compressive loading. The angular range of activated cracks depends then only partly on the history or loading path (Lehner & Kachanov, 1995), as was already suggested by Sanders (1954), who employed the plasticity formalism of multiple loading functions of Koiter (1953). Phenomenological plasticity mod els with the features suggesting the presence of sliding microcracks have been proposed in the past (Rudnicki & Rice, 1975), emphasizing the importance of the micromechanism of deformation on the conditions for overall rock failure in a shear-band mode. A similar influence on the critical tectonic stresses required for the initiation of folds is discussed in detail by Leroy & Triantafyllidis (this Volume).
The contents of this paper are as follows. The next section contains a small strain formulation of the plasticity model, applicable to the modeling of tectonic deformation in three dimensions. This model could easily be extended so as to allow for large deformations. The third section is devoted to a detailed analysis of a triaxial test in both extension and compression. In the latter case, a plane of weakness is also accounted for in a manner which is similar to the proposition of Jaeger (1960). � While only slip can be generated in the compressive triaxial test, both opening and slip modes are activated during extension under triaxial stress conditions. This distinction is important for a proper understanding of the initial plastic flow in a developing fold, which is studied in the fourth section as a problem of plate bending; extension generated in the extrados of the bent plate is found to result in the opening of discontinuities. In the concluding discussion, the scope of these theoretical results as an aid in the prediction and interpretation and natural fracture systems is critically evaluated for a well-exposed, densely fractured anticline of Devonian sandstone in South Morocco (Gaulier et al., 1996).

2
The plasticity model The constitutive relations developed in this section pertain to the mechanical response of an elementary volume of a rock mass. The tractions acting on that volume are taken to be those generated by a uniform 'macrostress' field. The inelastic or 'plastic' deformation behaviour of the element at sufficiently high stress levels is assumed to exhibit a strong dependence on the mean stress, or pressure. It is this pressure sensitivity, which is best known from dry friction behaviour of sliding surfaces, that has suggested to us the idea of a plasticity model in which the actual physical micromechanisms of deformation (i.e., grain boundary sliding, growth of transgranular cracks followed by grain crushing, crack growth and interaction on the scale of many grains, etc.) appear lumped together on a continuum scale, where they give rise to a behaviour that can be represented mathematically by the concept of multiple yield-surfaces, as in the work of Koiter (1953). In such a model material, each yield surface represents the constraint imposed on the state of stress by the strength behaviour of real or imagined plane surfaces of discontinuity in the displacement 1 . Despite the somewhat ficticious nature of these pervasive discontinuities, it will be helpful to imagine a volume of rock that contains a pervasive system of potential discon tinuity surfaces of all orientations, potential, because in order to contribute its share to the bulk strain of the elementary volume, any such discontinuity must first be activated. Moreover, when taken as real, these pervasive discontinuities will suggest the right kind of mathematical continuum model for a rock mass that contains discrete sets of regional fractures or systematic joints. The range of orientations of the potential planes of discontinuity, as defined by their unit normal vectors n, extends over the hemisphere H (of unit radius) shown in Figure 1. Accordingly, every unit normal vector n defines a plane in an oriented material sample, which can accommodate jumps in tangential or normal displacement. It is assumed that the magnitudes of these jumps are not affected by the presence of neighbouring discontinuities. It is further stipulated that any increase in magnitude of these irreversible discontinuous displacements will necessitate an increase in the macrostress, so that the material element exhibits an overall work hardening behaviour. This assumption ensures a stable response on the macroscale. Its absence would limit the load to a maximum, attained at the first activation of a discontinuity, and this would immediately enforce localized faulting on the macroscale. (For further discussion of these issues, see also the article by Leroy & Triantafyllidis in this Volume). The macro-scale strain rate E is composed of an elastic and a plastic part, the former being related to the stress rate by Hooke's law. The plastic part of the strain rate is denoted by (oP. Bold quantities identify vectors and second-order tensors. A superposed dot denotes differentiation with respect to time. The rate of macro-scale plastic straining (oP results from the additive contributions (oP(n) of all active discontinuities of orientation n on the hemisphere H: ;ov = L ;ov ( n) dH . (1) In order to obtain the permanent or 'plastic' strain rate associated with the discontinuity of normal n, it is necessary to identify the force acting on that plane. It is assumed that the local stress vector, t(n), acting on the plane may be obtained directly from the macro-scale stress tensor u by the relation 2 2 A single dot is used to denote the scalar or 'dot product' of two vectors, the product of a tensor and a vector as in (2) with components ti = O"ijnj, or the product of t(n) = a· n. (2) The stress vector may be viewed as composed of a normal and a tangential component according to t(n) = tN(n)+tr(n), as shown in Figure 2. The normal stress uN(n) on the plane of orientation n, i.e., the magnitude the component tN(n) is given by 3 '4 uN(n) = n · t(n) = n ·a· n = (n ® n): a. (3) Writing tN(n) = uN(n)n = (n · t(n))n = (n ® n) · t(n), for the normal com ponent, also yields 5 for the tangential component of the stress vector. This component defines a unit tangent vector rn in the plane of orientation n by ( ) _ tr(n) rn n -T(n) ' (5) in which the resolved shear stress T(n) is defined as the magnitude of tr(n). This shear stress is conveniently expressed in terms of the macrostress and the unit vectors rn and n by 1 T(n) = rn ·a· n = 2 [rn ® n + n ® rn] : a, where the symmetry of the stress tensor has been exploited to introduce the symmetric Schmid tensor as a factor forming a scalar product with the stress. Note also that the rate + is related to the time derivative of the stress tensor ir in the same way.
two second-order tensors C=A · B with components Cij = AikBkj in a Cartesian coordinate system. Here repeated indices imply summation over the range of the index, i.e., from 1 to 3. 3 Note that compressive stresses are taken negative throughout this paper.
Having identified the forces acting on the discontinuity of normal n, we now propose a phenomenological plasticity law to describe the two mechanisms of interest: the generation of discontinuous slip and opening displacements on a plane surface of discontinuity. For slip to occur on a given discontinuity surface, we require that the equality in the Mohr-Coulomb yield criterion (7) be satisfied. Here J.L is a constant coefficient of friction and c0 is the cohesion, which is taken to be a function of the internal variable 1(n) . Note that these strength parameters could both be functions of the orientation of the plane of weakness. The internal variable is conceived here as dimensionless equivalent plastic shear strain and is defined by (8) where E� ' (n) is the deviator6 of the plastic strain E�(n) that is generated by slip on discontinuity planes of orientation n. The equivalent plastic strain rate in shear, ')t(n), for a discontinuity of orientation7 n may be determined from the so-called 'consistency condition' which calls for the continuing satisfaction of the yield criterion (7) as a condition for continuing plastic straining. The condition therefore demands that ;ps ( u, n, 1( n)) = 0. After substitution in (7) of the expressions (3) and (6) for uN (n) and T(n) in terms of the macrostress, and differentiation, this produces the following result for the plastic rate of shearing: with (9) Here hs stands for the derivative of the function c0(1). Note that 8cps j 8u defines the normal to the yield surface cps in stress space.
We must now identify the plastic deformation rate E� ( n) which contributes to the microstrain E�(n) whenever the loading conditions cps = 0 and �� : & > 0 are fulfilled. It is constructed by assuming that slip on a plane with surface normal n takes place in the direction rn of the tangential component of the stress vector and that the deviator of the resulting micro-strain rate is proportional to the Schmid tensor: At least three interpretations can be given to this expression, interpreted as the micro-scale contribution to the macro-scale strain rate introduced by (1). Thus, it is first observed that the right-hand side of (1) has a structure that is familiar from theories of crystal plastic deformation (see, e. g. , the review by Asaro, 1983) in which the integral is typically replaced by a sum over a finite number of slip systems. The tangential vector rn may however be interpreted differently in these theories as being determined by a crystallographic direction rather than the resolved macrostress.
In a second interpretation of (1) and (10), the integration over all orientations is again replaced by a finite sum, this time over preferentially oriented sets of faults in a fractured formation, as in the work of Gauthier & Angelier (1985). The present model differs from theirs mainly in that a constitutive relation is invoked for the displacement along every discontinuity.
Expression (10) also has an interpretation for slip occurring under dynamic conditions. The Schmid factor on the right-hand side of (10) then equals the seismic moment tensor density, mqp, normalized by 2G, where G is the shear modulus of an isotropic elastic medium (Molnar, 1983). The difference between the present static and the dynamic interpretation of (10) is thus that in the former ')' depends on the rheology of the sliding discontinuities rather than on the elastic properties of the surrounding rock mass. Dilatancy during slip results from the presence of stiff heterogeneities and asperities on the sliding surfaces. Following Rudnicki & Rice (1975), the vol umetric strain rate is therefore defined in terms of the equivalent strain rate in shear and a dilatancy coefficient f3 that may depend on I· In the following, this coefficient is assumed to be constant and independent of the orientation of a discontinuity. The microstrain resulting from the dilation of a discontinuity of orientation n is thus taken to be proportional to i'( n )n ® n, with /3 as the factor of proportionality. The total strain rate associated with slip on a dilating discontinuity thus becomes and this contains the essence of our plasticity model for slip along a discontinuity of orientation n. Note that for the sake of simplicity the dependence of "y(n) on the macro-scale stress state has not been made explicit in the above expressions.
Let us now consider, as a second micro-scale deformation mechanism, the opening of discontinuities independent of any slip displacement. Such a mecha nism could already be activated at small compressive normal stresses (()N(n) < 0) and its operation can be allowed for by the introduction of a second yield criterion of the form which involves only the normal stress ()N and a material function k0(d) of the accumulated opening strain, d, an internal variable that may depend on the ori entation of the discontinuity. This criterion acts as a cut-off of the Mohr-Coulomb criterion in stress space, as is illustrated Figure 3. Note that k0(d(n)) is taken to be initially negative, which excludes the state of zero stress from the initial elas tic domain. This choice is consistent with numerous observations of the damage induced in core samples retrieved from great depth. It makes allowance for the operation of one or more micromechanisms, not to be discussed here, that effect a transformation of compressive macrostresses into locally tensile states of stress. ko is taken to be an increasing function of d, reflecting an increasing resistance of the discontinuity to opening. The structure of the plastic strain produced by this second deformation mechanism resembles that associated with slip-induced dilation in (11 ), if "y/3 is replaced by d. The rate of plastic deformation associated with this opening mechanism for a discontinuity of orientation n is thus written Here d is determined from the consistency condition 4P yield criterion (12), giving in which the scalar h0 denotes the derivative of the function k0(d) . Note that the opening mechanism is activated, if the conditions qP = 0 and ( oqP I 00') : (r > 0 are simultaneously satisfied. 3. The yield locus in stress space for a potential discontinuity: A Mohr-Coulomb envelope truncated by the opening criterion a-N = ko(d).
The second mechanism, which potentially excludes the zero stress point from the admissible domain, may raise a question about the stability of the postulated model material in the sense of Drucker (1951). To dispel such concerns, it suffices to consider a single discontinuity of orientation n and a stress state a * in the elastic domain or on the loading surface. Identify, for example, a * with point P in Figure 3. The material is stable in the sense of Drucker, if the work done by an external loading process in a cycle starting at a * and resulting in plastic deformation is strictly non-negative. Drucker showed that this condition requires that (aa * ) : E� + 6-: E� > 0, in which a lies on the loading surface. The second term on the left-hand-side of this inequality is strictly non-negative if d is non-zero, given the normality of the plastic flow (13) and as long as there is positive work hardening (h0 > 0). The same property and the convexity of the yield condition (12) entail the positiveness of the first term in the above inequality. The assumed deformation mechanism of the opening of discontinuities is thus stable in the sense of Drucker. It should also be recalled here that the relevance of Drucker's postulate as a sufficient, but not necessary condition for the stability of the first mechanism of frictional slip has been repeatedly discussed in the literature (see, for example, the note by Mandel (1964) and related discussions in the same volume) and will not be addressed here. A further point of concern could be the lack of uniqueness of the normal direction to the loading surface at the intersection, in stress space, of the yield surfaces corresponding to the two distinct criteria (7) and (12) (cf. Fig. 3). However, no such problem arises with the present model in view of the independence of the two postulated deformation mechanisms (Koiter, 1953). An indeterminacy could have resulted here, if such an interaction had been allowed for, e.g. by the introduction of latent hardening as it is appropriate for some crystalline materials (Pan & Rice, 1983).
In summary, the total plastic rate of deformation, due to the activation of a planar discontinuity of orientation n by the two independent micromechanisms, is obtained as the sum of the strain rates (11) and (13): Here the scalar factors a3 ( n) and a0 ( n) take on the values 1 or 0, depending on whether or not the corresponding micromechanism is activated.
Having characterized the plastic deformation accumulated on individual dis continuities by the two independent micro-scale deformation mechanisms, we must now determine the resulting macro-scale plastic deformation E_P. This last operation involves the substitution of (15) for the integrand in (1), followed by an appropriate integration over the hemisphere H. Before carrying out such a calculation for the triaxial test and the bending of a plate, we wish however to discuss the incorporation in our plasticity model of weak discontinuities of preferred orientation. These discontinuities are taken to represent inherited sys tematic joints that were formed during earlier episodes of stressing and deforma tion and will contribute to either slip or opening. In that sense, they are weaker than any potential discontinuity of the rock mass. Systematic joints are also called regional joints because of their large extent. It is assumed, however, that their contribution to the tectonic deformation occurs by slip or opening along patches that are sufficiently well distributed to justify the assumption that the macro-scale stress field prevails on every discontinuity. The present model is built on the assumption that there exists a finite number N of joint sets of a given precise orientation, which cut through the elementary volume. The rele vant components of stress acting on a discontinuity plane are assumed to satisfy a plastic yield criterion, similar to the one proposed for the above continuously distributed potential discontinuities, but with different material parameters. The total macroscopic plastic strain rate is composed of the integral term in (1) and a sum Ef €�(na) of terms that represent the contributions to the macrostrain from N individual sets of discontinuities. It is thus given by where the €�(na) have been assumed to satisfy an expression of the form (15) .
A similar superposition of strains resulting from slip on weak discontinu ities and microstrains generated pervasively throughout a rock mass has been considered earlier by Jaeger (1960). An obvious problem arising in this context is that of a scale-dependence of certain results, such as for example the calcu lated stiffness of the material volume element under consideration. In general, one must expect quantities of this kind to depend on the size of the elementary volume by the mere fact that different volumes will sample different ensembles of systematic joints8. In the best of all situations, this size dependence will fade as the volume exceeds a certain minimum size that is still small in comparison with a characteristic overall length of the rock mass of interest. One can then speak meaningfully of a representative elementary volume (REV), a qualifica tion we have so far avoided on purpose. Indeed, a challenging task for studies of the present kind remains precisely the determination of typical length scales of REV's in the field.

3
The triaxial test This section is devoted to an analysis of the permanent deformation of an elasto plastic material of the type just defined, when a representative sample of this material is subjected to standard triaxial loading conditions. Linear hardening laws are adopted for the cohesion and the opening hardening functions in order to obtain analytical expressions. The triaxial compressive test, in which only sliding discontinuities can be activated, is considered first. Triaxial extension, which permits sliding and opening displacements to occur, is discussed in the second part of this section. Assume now that the material sample is initially subjected to the isotropic compressive stress -PI, where P is a positive scalar pressure. The specimen is subsequently loaded by a monotonically increasing or decreasing axial compo nent of stress, q, in direction 2 ( cf. Fig. 4), so that the total stress state becomes u = -PI + q ez 0 ez , where q varyies monotonically with time and e2 is the base vector associated with the x2 coordinate axis.
The components of the unit normal n of a potential discontinuity plane in the Cartesian coordinate system of Figure 1, when expressed in terms of the two Euler angles e and rp, are The normal stress uN (n) , the resolved shear stress T (n) , and the vector rn (n) are given by uN (n) = -P + qcos 2 e, T (n) = l q l cose sine, rn (n) = Sign(q) (-cose cos cpe1 + sin Be 2 -cos Bsin cpe 3 ) .
These expressions are found by using (5) and (6) in (4) and (5), while observing that T = l tr · tr l 1 1 2 . Note that the normal stress and the resolved shear stress are independent of the Euler angle tp. A direct consequence of this axial symmetry is that the activation of potential discontinuities always occurs for all values of this angle. However, this symmetry will be carried over to the deformation of the sample only in the absence of any discrete sets of joints of preferred orientation. The first part of this section is concerned with the compressive triaxial test for which the scalar q in (5) is negative and monotonically decreasing. However, the equations to be presented are kept general, allowing q to be of any sign in anticipation of the subsequent analysis of the extension test. Inspection of the evolution of UN and T according to (7) and of the structure of the yield criterion (12) leads to the conclusion that it is impossible to activate the opening mode if q is negative. This conclusion is also evident from the illustration in Figure 3 of the stress path followed during compression. We shall thus focus on the sliding mode in this first part of this section. Inserting the expressions (7) for UN and T in the yield criterion for sliding mode (7) provides the following condition valid throughout the test: 13 sin(2( B + Sign( q )q)) :S [ 2c0 (r( n )) cos q) + (2P -q) sin q)]/lql , The inequality is satisfied in the elastic range of deformation and we now inspect the conditions for which the equality holds for the first time, defining first yield. Since c0 increases with /, we must examine solutions to equation (8) for the initial value c0(0). The two solutions found are is the classical Coulomb angle, that is the angle of bisection of the first activated pair of discontinuity surfaces by the largest compressive principal stress. Since 0 :S q) :S 7f/2, this falls into the ranges 7f/4, (0) :S Be :S 1rj2, (7r/4) for Sign(q) = -1, (+1). The angle B1 = B2 =Be thus represents the solution of (8) obtained for a critical load such that the argument of the arcsin function in (11) drops to the value 1 for the first time, that is for the load S . ( ) 2[c0 (0) cos q) + P sin q)] qc = 1gn q .
· 1+ Sign(q) sinq) Beyond this critical value of q, two solutions, B1 and B2 will exist for equation (8). These furnish the limits of the orientation range of activated discontinuities [B1, B2] = B2 -B1. The existence and growth of this active orientation range has been indicated already in Figure 4c. It is made more precise by Figure 5, where the limits B1, B2 are displayed in a Mohr diagram, making use of the stress 'pole' at the point (-P, 0) (see, e.g., Mandl, 1988). The fact that the orientation range [B1, B2] grows symmetrically about the Coulomb orientation Be is essentially a consequence of the assumed hardening of the slip systems. This permits a continuing increase in the load, which in turn results in the activation of less favourably oriented slip systems and a corresponding growth of the angular range of activated discontinuities. At any given instant, the most critical orientation, Be, will have hardened more than any other orientation. The state of stress on that plane is represented by the tangent point of a Coulomb envelope, which has been shifted in accord with the amount of cohesive hardening corresponding to the current value of !(Be, q), and the stress circle passing through the points (-P, 0) and (-P + q, 0). (In constructing Fig. 5, we have made use of (18) and the linear hardening law ( 17).) Since the stresses on any newly activated planes of orienation B1 and B2 remain constrained by a Coulomb condition for first slip, with c0 = c0(0) , they must be represented by the points of intersection of the initial Mohr envelope with the current stress circle.
Let us now quantify the shear deformation associated with slip on active discontinuities and from this determine the relation between the applied load and the permanent macroscopic deformation of the sample. If slip-induced dilatancy is disregarded (/3 = 0) in addition to opening mode (13) and a correspondingly simplified form of expression (15) is substituted in (16), together with (5),(6), and (7c), there results the following expression for the macro-scale plastic strain rate: where R( e, cp) stands for the Schmid tensor � ( m 0 n + n ® m ) , the component matrix of which becomes [ -sin2Bcos 2 cp (1 -2cos 2 B) cos cp -� sin2Bsin2cp l (l -2 cos 2 B) cos <p sin 2B (1 -2 cos 2 B) sin cp . -� sin 2B sin 2<p (1 -2 cos 2 B) sin cp -sin 2B sin 2 cp (25) The fact that the strain rate 'Y does not depend on the angle <p simplifies the integration of (14). Moreover, on account of (9a) and (5), !(B) depends on time only through the time-like factor q; (14) is therefore readily integrated along the loading path to provide the total macroscopic plastic strain: To render explicit the relation between this plastic macrostrain and the ap plied load, the functional dependence of the cohesion c 0 on the internal variable 1 ( n) must be specified for all orientations n of potential discontinuity surfaces. For simplicity, we consider the linear relation (27) Here c1 and c 2 are constants that are independent of the orientation n. A re lation of this type can be motivated for crack-like discontinuities which, in the absence of crack growth, exhibit a linear elastic compliance, i.e., a form of 'linear hardening' (see, e.g., Lehner & Kachanov, 1995). A shortcoming of (17) remains however the absence of an expected saturation level for work-hardening, that is a value of 1 beyond which no further appreciable increase in cohesive strength will occur.
Using (17) in (8) and rearranging, the accumulated plastic shear strain may be expressed as a function of the loading parameters P and q as follows: The terms within vertical bars are now evaluated for the limits of the active orientation range, 81 and 82, the former being subtracted from the latter. By use of ( 11) the plastic strain may then be expressed entirely in terms of the load parameters P and q.
To complete this discussion of triaxial deformation, let us now determine the contribution from a set of systematic joints or 'weak planes', initially present in the sample, to the overall plastic deformation. For simplicity, only a single joint set of orientation ( B I, c.p I) is considered. If activated, its contribution to the macro-scale deformation will be the inelastic strain (30) [ -sin2Bicos 2 cpi (1 -2cos 2 BI) cost.pi -sin2Bisin2cpJ/2 ] 1 I ( B I) ( 1 -2 cos 2 B I) cos c.p 1 sin 2B I ( 1 -2 cos 2 B I) sin c.p I , -sin(2BI) sin(2cp I) /2 (1 -2 cos 2 BI) sin c.p I -sin 2BI sin 2 cp I where !I( BI) is obtained from expression (18) for r( B) upon replacing the argu ment B by the fixed angle BI and using appropriate parameter values. The total permanent macro-scale deformation in the presence of the joint set is then given by the sum of (29 ) and (30).
In order gain some perspective, it will be of interest to evaluate the above general results for the following specific data set: A confining pressure P for the compressive test of 10 MPa, an initial cohesion c1 on all potential discontinuities of 10 MPa, and a reduced initial cohesion one half the value of c1 on a single preexiting joint set of orientation ( B I, cp I) = ( 7r j 4, 0). All discontinuity surfaces have the same coefficient of sliding friction f.1. = 0.6, and the work-hardening constant c 2 is given the uniform value of 50 MPa.
Consider first the stress-strain response in triaxial compression, as given by (29 ) for a sample without weak joints of preferred orientation. First yield occurs in this test at a critical axial load, whose normalized magnitude I 0"221 /c1 = 1 -P + qc l/c1 is found to equal 6.46 from (13) (cf. also Fig. 5). This marks the intercept with the zero-strain axis of the inelastic stress-strain response shown in Figure 6 (dashed line), in agreement with classical Mohr-Coulomb theory. The subsequent computed stress-strain response exhibits nonlinear hardening behaviour and this is clearly a consequence of the widening orientation range [ B1, B 2 ] of active slip systems with increasing axial load, since individual systems have been assumed to respond linearly; the larger the range of activated discontinuities, the softer the material response.
The evolution of the angular range of activated slip systems with increasing plastic deformation is shown in Figure 7, in which the limiting angles B1 and B 2 have been plotted versus the axial strain. The contribution of each activated plane to the macroscopic strain may be appreciated from the graphs shown in Figure 8 of the accumulated plastic shear on individual planes as a function of their orientation B. Note the increase in the range of activated systems and in the magnitude of 1 as loading proceeds.
The same triaxial loading history is now imposed on a specimen containing a set of weak planes with the above-assumed orientation and material proper ties. The computed macroscopic stress-strain response is given by the dotted line in Figure 6. Slip commences on the weak planes, but at the same normal ized critical load of 6.46, found previously in the absence of weak planes, slip on distributed discontinuities sets in with a gradually expanding orientation range of activated planes. Note that the slope of the stress-total strain curve remains constant as long as only the single set of weak joints is active (linear hardening), but decreases as soon as the growing range of distributed discontinuities be comes active and softens the response. When plotted versus total plastic strain, the orientation ranges of actively slipping discontinuities differ markedly in the absence (solid lines) and presence (dotted lines) of the weak joint set, the range being substantially diminished in the latter case (Fig. 7).
The second part of this section is concerned with a triaxial test in which the sample is allowed to undergo extension in the axial direction, so that both both sliding and opening displacements can be operative on potential surfaces of discontinuity. The representative volume element is assumed to be free of any weak joints. The macro-scale plastic deformation resulting from the opening mode is now discussed briefly and only the main results will be given.
The yield condition (12) gives the critical normal stress at which the opening mode of potential discontinuities is activated. On substituting the stress aN from (7) in this criterion, the resulting condition for the initiation of the opening mode is It predicts that the first discontinuities will open at a critical value of q equal to p + k1, at an orientation angle e = 0, i.e., in planes perpendicular to the Xz-axis ( cf. Fig. 4b ).
Since work-hardening has been assumed for the opening mode as well, further loading results in the opening of discontinuities sub-parallel to the x1, x3-plane in a manner similar to the activation of the sliding cracks. The angle range of open discontinuities is given by the following solution to (31):

20
() (J ko(O) + P) The macroscopic permanent strain associated with the opening of distributed discontinuities is found by integrating the appropriate term in expression (15) over the active orientation range and over time, giving As with the sliding mechanism, work-hardening is characterized by a linear re lation (34) in which the two constants k1 and k 2 are the same for all discontinuity orienta tions. This choice yields the expression d( ()) = q cos 2 ( ()) -P-k1 for the opening displacement. Substitution in (33) Since the two micro-scale deformation mechanisms operate independently, the total plastic deformation is obtained by adding up the individual contributions from slip (Eq. 29) and opening displacements (Eq. 35). Note again that the equa tions for slip have been written in a form that remains valid both for compression and extension. The feature of interest in the extension test is the possibility of a simultaneous activation of the two micromechanisms of slip and opening and its consequence for the stress-strain response. This is illustrated by the examples shown in Fig  ure 9. The parameters k1 and k 2 of the hardening relation (34) are set to -5 and 10 MPa, respectively. Three extension tests are considered for confining pres sures of 30, 40, and 60 MPa, corresponding to 3, 4, and 6 times the value of c 1 . At the largest confining pressure P, yield occurs first by slip. From (12), the orientation of the discontinuity that slips first is found to be ()c = 7r /4 -¢/2. Further loading increases the range of activated discontinuities until the normal stress becomes small enough to satisfy the yield criterion for the opening mode. Opening first occurs in planes normal to the x 2 -axis when the axial stress reaches the value k1, i.e., at 0' 22 /c1 = -0.5 in the present example. Continued straining is accompanied by the simultaneous activation of opening and sliding modes. If the same test is conducted at the intermediate value for P of 4c1 (40 MPa), then the order of activation of the micromechanisms is reversed: The opening mode preceeds the sliding mode. According to (13), the latter is activated as soon as the axial stress attains the value [ 2c1 cos¢ -(1 -sin¢ )P]/(1 + sin¢), i.e., when 0' 22 /c1 = -0.1785 for the above parameter values. Finally, for the smallest con fining pressure P, equal to 3c1 (30 MPa), only the opening mode of deformation is activated. Note again the nonlinearity of stress-strain response, which results from the progressive opening of discontinuities sub-parallel to the x1, x 3 -plane.
This sequential activation of the two micro-scale deformation mechanisms and its sensitivity to the prevailing hydrostatic stress state is also an essential feature of plate bending. As will be discussed next, only slip will be activated in the plate intrados (region of compression), whereas slip and opening displace ments can occur in the extrados (region of extension).

4
Bending of a plate The initial stage in the formation of a fold is analyzed in this section as a problem of plate bending for a plate of thickness 2D and infinite extent ( cf. Bending about the x 3 -axis generates principal stresses in the x2 and X 3 directions that vary linearly across the plate thickness with a gradient proportional to the load parameter q, as illustrated in (b) for 0'22· With reference to Figure , it is assumed that the deformation in every plane x 3 =const. is one of plane strain. The unit normals n of all activated discontinu ities must lie in the (.:z:1, .:z: 2 ) plane, as a consequence of which the second Euler angle rp only takes on the values 0 and 11. Moreover, the plane strain assumption and the absence of permanent deformation in the x 3 -direction obviously implies that the elastic strains in that direction must also vanish and this constrains the component of stress cr 33 . The state of stress in the plate is therefore given by: where Pis the initial isotropic pressure and v denotes Poisson's ratio. The stress is seen to vary linearly with x1 across the plate and the parameter q is assumed to increase monotonically from q = 0 during loading. The stress distribution (36) satisfies equilibrium pointwise throughout the plate as long as gravity effects are disregarded. For q 2': 0, the bending bending generates compression in the region x1 > 0 (the intrados region) and extension in the region x1 < 0 (the extrados region). These stresses in the two regions thus differ in a similar way as in the triaxial compression and extension tests.
The components of the unit normal to potential discontinuities direction are again given by (6), with the second Euler angle rp set to either 0 or Jr. From (6) and (36), one proceeds as in the previous section to derive the following expressions for the normal stress and the resolved shear stress (for q > 0): CTN (n) = -P -(qxl/ D) cos 2 e, T (n) = (q l xl l / D) sine cos e.
Substitution in the yield conditions (7) and (12)  and this clearly happens first in the extrados at x1 = -D, while the initiation of slip at x1 = + D in the intrados requires a larger value of q (exactly by a factor 3 for a friction angle of 30 degrees).
It follows from (39) that the opening mode of potential discontinuities can be activated only in the extrados. This will happen first at x1 = -D for planes running parallel to the x1, x 3 -plane as soon as q = k0(0) + P. Further loading results in the opening of sub-parallel discontinuities with orientation angles that lie within the range [0, 8 3 ]. The value of the initial pressure P determines whether opening or slip will occur or whether both mechanisms are activated, as in the triaxial extension test. This effect is now illustrated by an example computation for three different values of P.
In this example, the plate is taken to be free of any weak joint or fracture systems of preferred orientation. All potential discontinuities have a friction co efficient p, of 0.3, a zero initial cohesion c1 and a hardening coefficient c 2 of 50 MPa. The value k0(0) in (30), i.e. the constant k1, is set to -10 MPa and the hardening parameter associated with opening, k 2 , to 10 MPa. The plate half thickness D is 500 m. Results are presented for a loading parameter q equal to 30 MPa and the three values of 30, 35 and 40 MPa for the confining pressure P.
In Figure 11a two plots are shown for the smallest confining pressure of 30 MPa, one for the extrados ( -1 :::; x1 :::; -0.5) and one for the intrados (0.8 :::; x1 :::; 1.). The boundaries between the elastic core of the plate and the adjacent plastic zone are marked by a dashed line and only a portion of the core is shown. In the extrados, only the opening mechanism is found activated. Opening is initiated in a direction perpendicular to the x1, x 3 -plane at the limit of the elastic and plastic regions (x1 = -0. 67 D) and the range of activated discontinuities increases away from the center of the plate. The limit 8 3 of the activated range is largest at the top surface of the plate. In the intrados, the range [ el ' e 2 ] of slipping discontinuities exhibits a similar trend, beginning with a single critical orientation Be at elastic/plastic boundary and growing with increasing X1 .
The results obtained for the larger values of confining pressure of 35 and 40 MPa (Figs. llb and llc) remain qualitatively the same in the intrados, while differing markedly in the extrados. In the intrados, the main effect of a higher confining pressure is an enlargement of the elastic core of the plate. At 40 MPa the intrados remain entirely within the elastic range (Fig. llc). In the extrados, on the other hand, increasing confining pressures are found to bring about a gradual shift from opening (dotted lines) to the shear (solid lines) as the pre ferred mode of plastic deformation, allowing both to spread simultaneously over different sections of the plate within a certain range of pressures (Fig. 11 b). At the largest value of P, only the slip mechanism is found activated (Fig. llc). 5 Concluding discussion To obtain a first impression of the potential and the limitations of the simple model described in the foregoing, we now wish to to compare its predictions for the bending of a plate with the fracture pattern observed in a sandstone anticline outcropping in the lVIoroccan Anti-Atlas as shown in Figure 12a, which has kindly been provided by Gaulier et al. (1996). The sandstone formation seen on Figure 12 is of Devonian age and is part of a fold formed during the Hercynian compression. The exposed sequence is at most 20 meters thick and is composed of several massive sandstone layers with thicknesses ranging from 2 to 5 meters. The study of Gaulier et al. (1996) provides a classification of the joint families in relation to the structural evolution of the folding. Fracture data, i.e., fracture orientation with respect to bedding, fracture density and spacing as well as fracture length and penetration depth, have been collected at more than 30 measurement sites. The classification of the joint families yields three fracture sets, two of them being apparent on Figure 12a. A first set, considered to predate the folding event, is made up of systematic joints perpendicular to bedding with a marked orientation. They can be detected on Figure 12a on the left-hand side of the photograph thanks to the sunlight's orientation: the exposed surfaces between two shadows belong to the systematic regional set. A second family is composed of extrados fractures perpendicular to the bedding and parallel to the fold hinge that is sub-perpendicular to the plane of the cliff. A precise description of this family is now given starting from the top of the anticline and moving down the exposed strata.
Looking at the top few meters of the anticline in cross-section, we note a dense population of meter-scale extrado joints. Their orientation, as plotted in stereographic projection in in Figure 12b (pole to fractures and fault plane trace) have been measured from the top horizontal surface along a distance of 10 me ters for each diagram. The fracture dip and strike show a small scatter around the vertical and fold hinge direction, respectively. The density of meter-scale extrados joints is larger in some intermediate layers than in the top layer. This observation is often explained by a difference in mechanical properties of the various layers. Note that some of the joints are intersecting several sandstone strata. These discontinuities could result from the coalescence of fracture sur faces initiated in adjacent beds. They are often curved, dipping towards the center of the photograph. Furthermore, detailed field observations clearly indi cate that some discontinuities have been activated successively as normal and reverse faults during folding. The reverse faulting event is disregarded from the subsequent discussion, since it belongs to the late stage of the compression. A typical normal fault, shown in Figure 13, is arrested as it crosses a layer with a different lithology. The secondary fracture pattern with a horse-tail geometry in the bottom layer permits to locate the fault boundary and to appreciate the extent of the extensional area.
We now turn our attention to a comparison between the predictions made by the constitutive model proposed in this paper and the field observations that have just been summarized. To begin with, it is necessary to fix the size of the elementary volume on the length scale of the field structure. We shall assume the elementary volume to cover the 20 metres of exposed thickness of the anticline. Motivated by the occurrence of joints that have been activated in are also provided at three locations on the top of the anticline, separated by a distance of 10 metres. (Courtesy Gaulier et al. , 1996). Fig. 13. Photograph of an extrado fracture that has been reactivated as a normal fault. It terminates in a layer composed of thin beds as a secondary fracture array with a 'horse-tail' geometry.
an opening mode, we interprete the elementary volume exposed in the field as representing a location along the vertical a.xis of Figure 11 in the extrados. The predicted orientations for discontinuities activated in the opening mode are then in qualitative agreement with those recorded on the stereonets of Figure 12b. The sliding mode is also found to occur simultaneously with the opening mode in the field, suggesting a situation as found for the extrados in Figure 11 b at an intermediate pressure level. However, the Coulomb angle predicted by (12) for a friction coefficient of 0.3 is ec = n/4+¢/2 � 53°, which is very different from the subvertical orientation of the sheared discontinuities in the field. It would appear therefore, that in contrast with multiple faults of the classical Coulomb-type ( cf . Mandl, 1988, p. 120 and Fig. I.3-23), 'normal faulting' in the present exposure has resulted from the initiation and coalescence of extrado joints. Their relatively narrow spacing, which appears to be controlled by bed thickness, has allowed the activation of a 'bookshelf mechanism' (Mandl, 1987). This is still recognizable in the left portion of Figure 12b, even though the faults have subsequently been reactivated in reverse and rotated clockwise. Rotation must have been anti clockwise during 'normal faulting', however. The deformation seen in the field also comprises a ductile component in some of the lithologies present, an aspect that is clearly beyond the scope of our plasticity model in its present form.
In conclusion it may perhaps be said that both the strength and weakness of the proposed continuum model come from its simplicity. A major difficulty remains the justification of our use of simple yield criteria and hardening laws as a way of lumping together the description of complex processes of activation and growth of discontinuous or localized deformation on a finer scale. A second principal difficulty has to do with geometrical simplifications and the choice of scale. As illustrated by the example of the complex horse-tail geometry shown in Figure 13, there are always important features that can only be studied in isolation and on their own proper scale, features that fail to be representative even as an element of a larger ensemble of like elements. On the other hand, given suitable elements of this kind, such as the weak joint systems discussed in this paper, there remains the important task of addressing the scale effects that have been mentioned here only in passing, i.e., when does one have a representative elementary volume and how to deal with situations in which no statistically stationary distributions of the envisaged structural elements exist? In view of this and the above fundamental question about the lumped mechanical behaviour of discontinuities it seems clear that simple models of the present kind should possess a significant potential for further development.