3-D Rock Mass Strength Criteria—A Review of the Current Status

: The presence of complex discontinuity patterns, the inherent statistical nature of their geometrical parameters, the uncertainties involved in the estimation of the discontinuity geometrical and geo-mechanical properties and complex three dimensional (3-D) in-situ stress make the accurate prediction of rock mass strength a difﬁcult task. It has been a great challenge for the rock mechanics and rock engineering professions to develop a rock mass strength criterion in three dimensions that incorporates the effect of the minor and intermediate principal stresses and captures the scale dependent and anisotropic properties resulting from the discontinuity geometry parameters, such as the number of discontinuity sets, 3-D discontinuity intensity, and the distributions of the discontinuity orientation and size. Rock mechanics and rock engineering researchers have dealt with this topic for more than 55 years. The paper provides a critical review of the current state of the art regarding 3-D jointed rock mass strength criteria. The shortcomings of several rock mass strength criteria are discussed. The historic development of rock mass strength criteria that incorporate the effect of the minor and intermediate principal stresses and capture the scale dependent and anisotropic properties is presented. The most advanced 3-D rock mass strength criteria currently available in the literature are presented, including suggested future improvements.


Introduction
Most naturally occurring discontinuous rock masses comprise intact rock interspaced with different types of discontinuities. Fissures, fractures, joints, faults, bedding planes, folds, shear zones and dykes are different types of discontinuities that exist in rock masses. Discontinuities may be divided into major and minor using the feature size. Large features may be considered as major discontinuities, and the rest of the small features may be considered as minor discontinuities. For most of the civil and mining engineering projects at the upper most level (for the domain size considered for analysis), the rock masses contain only a few major discontinuities.
For such projects, major discontinuities can be considered as single features, and their geometry may be represented deterministically. On the other hand, due to the large number and inherent statistical nature, the geometry of minor discontinuities should be characterized statistically. Henceforth, the minor discontinuities are referred to as either "joints" or "fractures" in this paper. Due to the presence of discontinuities, the geomechanical response of discontinuous rock masses can be highly complicated under complex geology and in-situ stress systems. Therefore, mining and civil engineers face many difficulties in tackling the design and construction tasks associated with geotechnical systems that are in or on discontinuous rock masses.
Such geotechnical systems cover a wide range of surface and underground excavations made for ore extraction, tunnels for hydropower and transport, dams, deep foundations, natural and man-made rock slopes, and underground repositories for oil, gas, and hazardous waste. In dealing with these rock engineering projects, rock mass stability is the main concern for the engineers. To arrive at safe and economical designs, it is very important to have a good understanding of the rock mass strength.
Rock mass strength depends on the (a) lithology; (b) discontinuity geometry network, including the number of discontinuity sets, their intensity and the spatial distribution of the orientation, size and spacing; (c) geo-mechanical properties of the discontinuities including roughness, strength and deformation of asperities and filling material; (d) geo-mechanical properties of the intact rock; (e) in-situ stress system; (f) size of the rock mass; (g) shape of the rock mass; (h) loading/unloading stress path; (i) loading rate; (j) pore pressure in the rock mass; and (k) environmental conditions (such as the temperature, humidity etc.) of the rock mass.
The presence of complicated fracture networks, the inherent statistical nature of their geometrical parameters and the variabilities and uncertainties involved in the estimation of their geometrical and geo-mechanical properties and in-situ stress make the accurate prediction of rock mass strength a difficult, challenging task [1,2]. On the other hand, unfortunately, understanding of the mechanical behavior of rock masses is crucial to designing safe and economical structures in or on rock masses.
To obtain intact rock mechanical properties, usually 50 mm (approximately 2 inch) diameter samples according to ASTM or ISRM guidelines are tested. At this size, the distribution of micro-cracks with the location of the sample has an influence on the scatter for the considered rock block property. A further increase of sample size will include the influence of one or more fractures. When the sample size contains only a few fractures, the sample properties vary considerably from one sample to another, reflecting the statistical inhomogeneity with respect to the influence of fractures. This variability decreases as the sample size is further increased to contain more fractures and to make the sample more statistically homogeneous with respect to the fractures.
Beyond a certain minimum volume, the mass property of the rock may not change significantly (from a practical point of view) with respect to the effect of fractures. This minimum volume may be used as the element size to represent the equivalent mass property of a statistically homogeneous rock mass volume that contains a significant number of fractures. This volume may be termed as a "Representative Elementary Volume" (REV) and may be of great importance for engineering purposes. Each rock mass property may have a different REV value. For some rock masses, the REV sizes may be small compared to the size of the problem domain of interest. For such cases, the REV sizes and associated properties will be very useful for engineering applications.
For some other rock masses, REV sizes may be large compared to the size of the problem domain. In such cases, if one knows how the rock mass property of interest varies between the intact rock value and the REV value with respect to the joint geometry parameters and block size, this will be quite useful for dealing with the latter category of rock masses in engineering applications. It is important to keep in mind that the variability of the rock mass property at the REV size or greater than the REV size will be smaller compared to the rock mass property at sizes less than the REV size. If the sample size is increased further to allow major discontinuities to enter the sample, then the rock mass property will vary reflecting the influence of major discontinuities.
In summary, the REV sizes and the corresponding REV mechanical properties can be used to represent the combined equivalent continuum behavior of minor discontinuities (joints) and intact rock. To this system of the rock mass, the major discontinuities can be added as single features to complete the representation of the whole rock mass. Estimation of the REV size, the REV mechanical properties and orthotropic constitutive models for jointed rock masses at the three-dimensional (3-D) level has been reported by Kulatilake et al. [2] and Wu and Kulatilake [3]. These results clearly indicated that the strength and deformability of rock masses show very significant scale effects and anisotropic behavior at the 3-D level due to pre-existing fracture systems.
It has been a great challenge for the rock mechanics and rock engineering professions to develop a rock mass strength criterion in 3-D that incorporates the effect of important fracture geometry parameters and intermediate principal stress and to capture the scale effects and anisotropic properties of jointed rock masses. Rock mechanics and rock engineering researchers have dealt with this topic for more than 55 years. The aim of this paper is to provide the state of the art of the estimation of 3-D jointed rock mass strength criteria using physical, empirical, analytical and numerical modeling methodologies.

Large Scale In-Situ Tests
Currently, the methods to estimate the strength properties of jointed rock masses can be categorized into direct and indirect methods. Direct methods include laboratory and in situ tests. Laboratory results obtained from small-sized specimens that include only micro-joints are very different from the results obtained from large-scale blocks because the laboratory samples cannot accommodate the whole spectrum of different size discontinuity network that is present in the field. A significant number of in-situ tests have been performed on various rock types to study the effect of size on the rock mass compressive strength using various sample sizes with different width to height (w/h) ratios [4][5][6][7][8][9][10][11][12][13][14][15][16][17].
Some of these investigators also performed corresponding laboratory tests on the same rock types. Bieniawski and Van Heerden [17] and Heuze [18] have reviewed the work done prior to their publications on scale effects on rock mass strength. The reported results of these investigations clearly show the reduction of rock mass strength with increasing size up to a certain size, beyond which change becomes insignificant from a practical point of view. Bieniawski [10] observed this size to be about 1.5 m for coal. Jahns [7] and Pratt et al. [15] noted this size to be about 1.0 m for diorite and iron ore, respectively. Bieniawski [10] observed the strength of the largest sample size (2 m) to be about 0.15 of the smallest laboratory sample size.
Pratt et al. [15] noted the previously mentioned ratio to be 0.1. Bieniawski and Van Heerden [17] provided valuable guidelines for experimental procedures to perform large scale in-situ tests. They also suggested various empirical equations to relate the rock mass uniaxial strength to the w/h ratio. It is important to note that the relations developed from in-situ tests in the above stated studies primarily depend on the discontinuity network of the tested rock masses. However, unfortunately, in these early investigations, no attempt was made to map the discontinuity network before subjecting the rock mass to mechanical behavior testing. Therefore, the reported relations are highly site dependent and have qualitative value only.

Jointed Block Testing with a Significant Number of Fractures in the Laboratory
To obtain realistic results for jointed rock mass mechanical properties, many large volumes of rock of different sizes with a number of different known joint configurations should be tested at significant stress levels under different stress paths. Such an experiment program is almost impossible to conduct in the laboratory. With in situ tests, such an experimental program would be very difficult, time-consuming and expensive. The results of laboratory model studies on rock-like materials [18][19][20][21][22][23][24][25][26][27] have shown that many different failure modes are possible with jointed rock and that the internal distribution of stresses and strains within a jointed rock mass can be highly complex.
Yang et al. [23], Kulatilake et al. [24][25][26] and Mehranpour et al. [27] found that the failure modes can be divided into three types, namely splitting failure through the rock blocks, sliding along joint planes and a mixed splitting-sliding mode. They also demonstrated that the failure strength exhibits pronounced anisotropy mainly due to the effects of joint orientation. Even though these jointed blocks included a significant number of joints, all the included joints were persistent joints. Only a few experimental studies have been performed using a significant number of non-persistent joints [28][29][30][31].
Prudencio and Van Sint Jan's [29] paper presents the results of biaxial tests performed on physical models of rock with non-persistent joints. The failure modes and maximum strengths developed were found to depend on, among other variables, the geometry of the joint systems, the orientation of the principal stresses and the ratio between the intermediate stress and intact material compressive strength· These test results showed three basic failure modes: failure through a planar surface, stepped failure and failure by the rotation of new blocks.
Planar failure and stepped failure were associated with high strength behavior and small failure strains, whereas rotational failure was associated with a very low strength, ductile behavior and large deformation. Chen et al. [31] investigated the combined influence of the joint inclination angle and joint continuity factor (lengths of the joints divided by the total length of the block) on the strength and deformation behavior of the jointed rock mass for gypsum specimens with a set of non-persistent open flaws under uniaxial compression.
Complete axial stress-strain curves were classified into four types, i.e., single peak, softening after multi-peak yield platform, hardening after multi-peak yield platform and multi-peak during softening. To investigate the brittleness of the specimens, the ratio of the residual strength to the maximum peak strength as well as the first and last peak strains were studied. At the same joint inclination angle, the ratios between the residual strength and the maximum peak strength and the last peak strains were observed to increase while the first peak strain decreased with the increase of joint continuity factor. At the same joint continuity factor, the curves of the three brittleness parameters versus the joint inclination angle were found to be either concave or convex single-peak or wave-shaped.
In summary, the mechanical behavior of jointed blocks with a significant number of non-persistent joints was found to be more complicated and significantly different to those with persistent joints. The mechanical behavior depended on the number of joints, joint orientation, size, density, spacing, arrangement and whether the joints were open or closed. It is important to note that, as the joint density increases, the behavior around each joint is affected by the presence of the rest of the joints in the jointed block.
He et al. [32] performed CT scanning tests to detect the pre-existing fracture systems of 36 cubical coal blocks. Systematic procedures were developed to construct the pre-existing fracture geometry of the cubical coal blocks from the CT images. Based on the constructed fracture systems, the fracture tensor components were calculated to capture the directional effects of the fracture system of the cubical coal blocks.
The effect of the pre-existing fracture system and different combinations of confining stresses on the strength of coal masses was investigated using some true triaxial test results and the computed fracture tensor components (see Kulatilake et al. [2] for details on the fracture tensor components). The strength of the cubical coal blocks in the loading direction was found to be closely related to the summation of the fracture tensor components in the two perpendicular directions to the loading direction and the level of the confining stress system.

Mohr-Coulomb and Hoek-Brown Strength Criteria
One indirect way to estimate the rock mass strength is through empirical equations. Several empirical rock mass strength criteria have been suggested in the literature. The oldest one is the Mohr-Coulomb criterion. The Mohr-Coulomb criterion provides a linear prediction to the rock mass strength without consideration of the intermediate principal stress, σ 2 . On the other hand, the experimental test results indicate that the rock mass strength exhibits non-linear properties and dependence on σ 2 .
This criterion over-predicts the rock mass strength at higher confining stress levels and provides an unreasonable tensile strength value. Therefore, for practical use, this criterion is always combined with a tension cut-off. In addition to the above shortcomings, no guidelines are available to relate the two strength parameters given in the Mohr-Coulomb criterion to discontinuity geometry and discontinuity mechanical parameters. It has been used in practice mainly due to its simplicity.
Most of the rest of the rock mass strength criteria are used with a certain rock mass classification system. Edelbro [33] reviewed twenty-one different rock mass classification systems developed by different researchers in the 20th century. Each system is applied to a specific area, such as tunneling, rock support, mining, and rock slopes.
All these rock mass classification systems only provide a rough evaluation for the rock mass quality with the following shortcomings: (a) all the parameters used in these rock mass classification systems are scalar based, which is not adequate for the description of the anisotropic behavior of the rock masses; (b) the rock quality designations (RQD) used in some rock mass classification systems are highly orientation and location dependent; (c) the RQD can be considered as one dimensional joint frequency; in the rock mass classification systems, which include both RQD and joint spacing, the same factor is double counted; and (d) the most important geometric feature of discontinuities-the joint orientation-is not explicitly used in all the systems, even though some systems provide an adjustment for joint orientation.
The Hoek and Brown rock mass strength criterion [34] was introduced in 1980 as an attempt to provide input data for the analyses required for the design of underground excavations in hard rock. This is an empirical criterion developed through the trial and error curve fitting of different parabolic functions to standard triaxial test data. The choice of a parabolic function appears to have originated from one part of the equation given for Griffith's crack theory [35,36]. It also involved Hoek's experience on research results with respect to the brittle failure of intact rock [37] and Brown's experience on research results of model studies of jointed rock mass behavior [20].
The empirical parameters of the criterion were first developed for intact rock, and then those parameter values were reduced by linking them to Bieniawski's [38] Rock Mass Rating (RMR) classification system. Since 1980, the criterion has been updated several times [39][40][41][42]. The most current appears to be the Generalized Hoek-Brown criterion. The empirical parameters of this criterion are related to the Geological Strength Index (GSI) introduced by Hoek et al. [39].
The GSI is determined based on the structure of the rock blocks and the surface conditions of the joints. The number of joint sets and the fracturing level are considered in evaluating the structure of the rock blocks. The roughness/smoothness, degree of weathering, degree of alteration, and presence and type of filling are considered in evaluating the surface conditions of the joints. Cai et al. [43] modified the descriptive term "the structure of the rock blocks" to "quantitative block volume" and the descriptive term "the joint surface condition" to "quantitative joint condition factor" in estimating the GSI value for a rock mass.
The Generalized Hoek-Brown criterion is expressed by the following equation: where σ ci is the uniaxial compressive strength of the intact rock; m b is a reduced value of the material constant m i given in the Hoek-Brown equation for intact rock, and the relation between the two is given by s and a are constants for the rock mass given by the following relations: D is a factor that depends upon the degree of disturbance to which the rock mass has been subjected by blast damage and stress relaxation. It varies from 0 for undisturbed in situ rock masses to 1 for very disturbed rock masses. Guidelines for the selection of D are given in Hoek et al. [41].
Relations are also available between GSI and RMR as well as between GSI and Q system [44] as given below: (a) For better quality rock masses (GSI > 25): for RMR 76 > 18, GSI = RMR 76 and for RMR 89 > 23, GSI = RMR 89 − 5. (b) For rock masses having GSI values less than 25, Hoek et al. [40] suggested using the Q system [44] to obtain the GSI value.
Being a non-linear strength criterion, the Hoek-Brown criterion provides more reasonable predictions compared to that of the Mohr-Coulomb criterion. In addition, a lot of information is available to estimate the empirical parameters given in the Hoek-Brown criterion. Both the Mohr-Coulomb and Hoek-Brown criteria are easy to use simple empirical equations. Therefore, in practice, the Hoek-Brown criterion has been used more often than the Mohr-Coulomb criterion. However, the Hoek-Brown criterion does not incorporate the effect of σ 2 and cannot capture the anisotropy, scale effects and ductile behavior of a rock mass.
In addition, the RMR, Q and GSI systems include the RQD. For the same rock mass, infinitely many values can be obtained for the RQD by considering different directions. This creates significant uncertainty and variability in estimating a single value for RMR, Q or GSI. Due to the assumption of isotropic behavior of the rock mass, the Hoek and Brown criterion can only be applied to highly fractured rock masses with highly randomly oriented fracture systems. Some investigators extended this criterion to 3-D versions to incorporate the influence of σ 2 [45][46][47][48][49]. None of these 3-D criteria include the effect of joint orientation and joint size explicitly. Therefore, these criteria cannot capture the anisotropy and scale effect of rock mass strength.

Yudhbir et al., Shoerey et al., and Ramamurthy Strength Criteria
Yudhbir et al. [50], Sheorey et al. [51] and Ramamurthy [52] proposed different forms of σ 1 = f (σ 3 ) equations to estimate rock mass strengths. Out of the three criteria, only the Ramamurthy criterion includes factors considering the effect of the orientation of a sliding joint or joint set; however, multiple joint sets are not considered in this criterion. The joint orientation is not explicitly considered in all three criteria. It is important to note that all three criteria cannot capture the influence of σ 2 , the scale effect and anisotropic rock mass behaviors. For details on these three criteria, the reader is referred to Kulatilake [53].

Kulatilake et al. Strength Criterion
Kulatilake et al. [26] performed uniaxial and biaxial compression tests on glastone model material intact prismatic samples of size 35.6 × 17.8 × 2.5 cm to determine the strength of the intact model material. For intact prismatic samples, tensile fractures almost perpendicular to the σ 2 direction were observed for samples subjected to low levels of σ 2 /σ u ≤ 0.03. Note that σ u is the uniaxial compressive strength of the intact model material. Fracturing parallel, or slightly inclined, to the σ 1 (major principal stress)-σ 2 plane (i.e., normal to σ 3 (minor principal stress) = 0) was observed for most of the intact samples that were subjected to high levels of σ 2 (σ 2 /σ u > 0.03). A new intact rock failure criterion in 3-D given by Equation (3) was proposed by Kulatilake et al. [26]. In the equation, J 1 and J 2D represent the mean stress (one third of the first invariant of stress tensor) and the second invariant of the deviatoric stress tensor, respectively.
and A, B and C are empirical material constants. In Equation (3b,c), σ 1,I , σ 2 , and σ 3 represent, respectively, the major principal stress applied on intact rock samples at failure, intermediate principal stress, and minor principal stress. In the case of a biaxial state of stress, J 1 and J 2D reduce to 1/2 × (σ 1,I + σ 2 ) and 1/3 × σ 1,I 2 − σ 1,I σ 2 + σ 2 2 , respectively. Therefore, the new intact rock failure criterion in the case of a biaxial state of stress can be expressed as: The intact material strength data of glastone showed a good fit to the new criterion (regression coefficient of determination of 0.97) along with an almost equally strong fit for modified Wiebols and Cook criterion [54]. However, the Drucker and Prager [55] and Mogi [56] criteria showed poor fits for the glastone data.
Uniaxial and biaxial compression tests on prismatic jointed blocks of size 35.6 × 17.8 × 2.5 cm were performed using the same biaxial load frame used to test intact prismatic samples. Both persistent and impersistent smooth joint configurations were included in producing jointed model material blocks. The results indicated three failure modes for the jointed blocks: (a) tensile failure through intact material; (b) combined shear and tensile failure or only shear failure on joints; and (c) mixed failure of the above two modes depending on the joint geometry. The fracture tensor component [2] was used to quantify the directional effect of the joint geometry, including the number of fracture sets, fracture density, and probability distributions for the size and orientation of the fracture sets.
The strength results obtained for the jointed glastone blocks were used along with the strength results obtained for the intact glastone blocks in developing the following new rock mass failure criterion (see Equation (5)) for biaxial loading conditions. In Equation (5a), σ 1,b and F 22 are, respectively, the jointed rock block strength under biaxial loading and the fracture tensor component in the σ 2 direction; σ 1,I is the intact rock strength under biaxial loading conditions and is estimated using Equation (4). In Equation (5b), σ u,I is the uniaxial compressive strength of the intact model material; and ω 0 , a and b are parameters of the rock mass strength criterion.
The parameter ω 0 is obtained from the relation given in Equation (6) between the uniaxial compressive strength of the jointed block, σ u,b , F 22 and σ u,I . Figure 1 shows the comparison obtained between the predictions and observed values for the relation between ω and σ 2 /σ u,I for the glastone material. Figure 2 shows the rock mass failure criterion on the σ 1 /σ u,I versus σ 2 /σ u,I space for different F 22 values. This figure clearly shows the importance of σ 2 on the jointed rock mass strength.
( ) (6) Figure 1 shows the comparison obtained between the predictions and observed values for the relation between ω and σ2/σu,I for the glastone material. Figure 2 shows the rock mass failure criterion on the σ1/σu,I versus σ2/σu,I space for different F22 values. This figure clearly shows the importance of σ2 on the jointed rock mass strength.

Analytical Modeling Applications to Estimate Rock Mass Strength
The second indirect approach to study rock mass strength behavior is through the analytical decomposition technique. Jennings [57] expressed the combined strength of joint and rock bridges through the following equation by a simple linear weighing of the strength contributed by each fraction of material: (6) Figure 1 shows the comparison obtained between the predictions and observed values for the relation between ω and σ2/σu,I for the glastone material. Figure 2 shows the rock mass failure criterion on the σ1/σu,I versus σ2/σu,I space for different F22 values. This figure clearly shows the importance of σ2 on the jointed rock mass strength.

Analytical Modeling Applications to Estimate Rock Mass Strength
The second indirect approach to study rock mass strength behavior is through the analytical decomposition technique. Jennings [57] expressed the combined strength of joint and rock bridges through the following equation by a simple linear weighing of the strength contributed by each fraction of material:

Analytical Modeling Applications to Estimate Rock Mass Strength
The second indirect approach to study rock mass strength behavior is through the analytical decomposition technique. Jennings [57] expressed the combined strength of joint and rock bridges through the following equation by a simple linear weighing of the strength contributed by each fraction of material: In Equation (7a), (c j , φ j ) and (c r , φ r ) represent the cohesion and friction angles of the joint and intact rock, respectively, and k is the joint continuity factor given by where L j and L r are the length of the joint and rock bridge, respectively. Note that Equation (7) ignores the interaction between the joints and intact rock and assumes si-multaneous failure of the intact rock and joints. This means it disregards the possibility of progressive failure. Amadei [58] also used the analytical superposition theorem to describe the strength of a regularly jointed rock mass under the full range of states of stress that are likely to be encountered in situ with one or several stress components being tensile. The intact rock was modeled using the Hoek and Brown criterion. The joints were modeled with no tensile strength and the shear strength described by the Coulomb criterion. Again, the interactions between the intact rock and joints were ignored in producing the derivation. In these derivations, simplified fracture systems are used with constant deterministic spacing and deterministic size and orientations for joint sets. In reality, joints are of finite size, and all the joint geometry parameters are inherently statistical.
Therefore, even though this early research has given some analytical thinking with respect to developing rock mass strength criteria, the assumptions made in these models are quite difficult to satisfy for most in-situ rock masses. Bekaert and Maghous [59] and Pouya and Ghoreychi [60] also have used analytical approaches, selecting suitable intact rock and rock discontinuity strength criteria and applying simplified methods to combine them to develop rock mass strength criteria. These methods are rarely applicable in dealing with field rock masses, which are typically more complicated than the assumed simplified models.

Numerical Modeling Applications to Estimate Rock Mass Strength
The third indirect approach available to estimate rock mass strength is the numerical decomposition technique. This technique calculates the strength of rock masses using a numerical model, incorporating the strength and deformability properties of the intact rock and joints. This technique allows the incorporation of any joint network to the jointed rock mass and also includes interactions between the joints and intact rock. Mainly, the following three types of numerical modeling have been used to model rock mass strength: (a) Finite Element Modeling (FEM), (b) Distinct Element Modeling through 3DEC [61] and (c) Particle Flow Modeling through PFC3D [62,63]. The following sections provide a summary of the literature available on the aforementioned numerical modeling categories.

FEM-Based Numerical Modeling Applications to Estimate Rock Mass Strength
A pioneering scale effect research performed using this approach at the two-dimensional level using the finite element method [1] with Goodman's joint element [64] showed anisotropic, scale dependent mechanical behavior for jointed rock masses. The REV and the equivalent continuum behavior for mechanical properties have been investigated at the 2-D level using the finite element method [1,60,65,66]. However, since the finite element method is based on continuum mechanics, the simulation of large displacements and large rotations is difficult with the method, even though they may occur in jointed rock masses.

Distinct Element Based Numerical Modeling Applications to Estimate Rock Mass Strength
The distinct element method introduced by Cundall [67] and further developed by Lemos et al. [68], Cundall [69] and Hart et al. [70] is a powerful technique to perform stress analyses in discontinuous blocky rock masses. In this method, the rock mass is modeled as an assemblage of rigid or deformable blocks. Discontinuities are considered as distinct boundary interactions between these blocks, and joint behavior is prescribed for these interactions. The distinct element algorithm includes not only continuum theory representation for the blocks but also force displacement laws, which specify the forces between blocks, and a motion law, which specifies the motion of each block due to unbalanced forces acting on the block.
By taking into account the interaction of intact blocks and joints, the distinct element method can effectively calculate the mechanical behavior of block systems under different stress and displacement boundary conditions. The method employs an explicit solution procedure. An advantage of the explicit method is that, because matrices are never formed, large displacements, rotations and complex constitutive behavior for both the intact material and joints are possible with no additional computing effort. For details on the distinct element theory, the readers are referred to the publications mentioned in this paragraph.
To use the distinct element method for stress analysis, first the problem domain should be discretized into polygons in 2-D and polyhedra in 3-D. To achieve this for rock masses with finite size actual joint configurations, Kulatilake et al. [71] and Wang and Kulatilake [72] introduced some fictitious joints that behave as intact rock to the domain to interact with actual joints. After performing a detailed study under different stress paths, Kulatilake et al. [71] provided recommendations to select the proper mechanical property values for these fictitious joints to reflect the intact rock behavior.
Using these techniques, pioneering detailed investigations have been performed to study the effect of finite size joint geometry networks on the deformability and strength of jointed rock blocks, REV size and equivalent continuum behavior at the 2-D [73] and 3-D levels [2,74]. The results of these studies have shown anisotropic, scale dependent mechanical behavior for jointed rock masses.
An incrementally linear elastic, orthotropic constitutive model has been suggested at the 3-D level to represent the pre-failure mechanical behavior of jointed rock blocks [2]. This constitutive model has the capability to capture the anisotropic, scale-dependent behavior of jointed rock blocks. In that model, the effect of the joint geometry network in the rock mass is incorporated in terms of fracture tensor components, which captures the effect of all the joint geometry parameters-the number of joint sets, joint density and distributions of orientation and size.
Wu and Kulatilake [3] applied similar procedures as stated above to investigate the effect of a fracture system that existed in a limestone rock mass located at the Yujian River Dam site, China on the rock block strength and deformability, REV sizes for mechanical properties and equivalent continuum behavior by performing numerical stress analyses in three perpendicular directions (x, y and z) on different block sizes ranging between 5 m and 50 m. Figure 3 shows the relation obtained between the jointed rock block strength in the i direction/intact rock block strength (S i /S I ) and block size for the investigated limestone rock mass. This plot shows the rock mass strength scale effect, anisotropy, REV behavior and equivalent continuum behavior.
Geotechnics 2021, 1, FOR PEER REVIEW 11 jointed rock mass at a refined level. In summary, these numerical studies showed the possibility of developing relations between the rock mass strength and rock mass joint geometry properties.      Section 3.3 discussed the Kulatilake et al. [26] rock mass strength criterion developed for biaxial loading conditions. Recently, He et al. [32] extended the Kulatilake et al. [26] criterion to the polyaxial stress condition based on laboratory and 3DEC numerical modeling on jointed coal blocks with non-persistent discontinuities. This was the first time any research group in the world developed a 3-D coal mass strength criterion incorporating non-persistent joints, both σ3 and σ2, and the main fracture geometry parameters (orientation and size distributions and 3-D intensity) explicitly. However, in the conducted research, because 3DEC modeling was used, the effect of fracturing was not incorporated in developing the rock mass strength criterion. He et al.'s [32] criterion is given below through Equations (8) and (9).
In Equation (8), JCMS and ICS are the jointed coal mass strength and intact coal strength under the same combination of σ2 and σ3, respectively; F22 and F33 are the fracture tensor components in the directions of σ2 and σ3, respectively. In Equation (9), ω incorporates the effect of σ3, as well as σ2; ω0 is the ω value for the uniaxial compression condition; ICSu is the intact coal strength under the uniaxial stress condition; and a, b, c and d are empirical coefficients. Figures 5 and 6 show a few typical relations obtained between the JCMS and the intermediate and minor principal stresses. For further details on distinct It is important to note that, in Figure 4, the calculated values obtained from the three perpendicular directions are given. An incrementally linear elastic, orthotropic constitutive model has also been suggested to represent the pre-failure mechanical behavior of the jointed rock mass at a refined level. In summary, these numerical studies showed the possibility of developing relations between the rock mass strength and rock mass joint geometry properties. Section 3.3 discussed the Kulatilake et al. [26] rock mass strength criterion developed for biaxial loading conditions. Recently, He et al. [32] extended the Kulatilake et al. [26] criterion to the polyaxial stress condition based on laboratory and 3DEC numerical modeling on jointed coal blocks with non-persistent discontinuities. This was the first time any research group in the world developed a 3-D coal mass strength criterion incorporating non-persistent joints, both σ 3 and σ 2 , and the main fracture geometry parameters (orientation and size distributions and 3-D intensity) explicitly. However, in the conducted research, because 3DEC modeling was used, the effect of fracturing was not incorporated in developing the rock mass strength criterion. He et al.'s [32] criterion is given below through Equations (8) and (9).
In Equation (8), JCMS and ICS are the jointed coal mass strength and intact coal strength under the same combination of σ 2 and σ 3 , respectively; F 22 and F 33 are the fracture tensor components in the directions of σ 2 and σ 3 , respectively. In Equation (9), ω incorporates the effect of σ 3 , as well as σ 2 ; ω 0 is the ω value for the uniaxial compression condition; ICS u is the intact coal strength under the uniaxial stress condition; and a, b, c and d are empirical coefficients. Figures 5 and 6 show a few typical relations obtained between the JCMS and the intermediate and minor principal stresses. For further details on distinct element method applications, the readers are referred to the publications mentioned in Section 5.2. element method applications, the readers are referred to the publications mentioned in Section 5.2.

PFC Based Numerical Modeling Applications to Estimate Rock Mass Strength
Some investigators have resorted to particle flow codes (PFC 2D and PFC 3D ) [62,63,75] to model jointed rock behavior under uniaxial loading [25,[76][77][78][79][80][81][82], bi-axial loading [26] and polyaxial loading [27]. PFC 3D allows one to study the mechanical interaction behavior between intact rock and joints incorporating a significant number of joints without making unrealistic assumptions about the surrounding medium around each joint. In addition, it allows failure through both the intact rock and joints under both tensile and shear modes leading to progressive failure, which usually occurs in jointed blocks with non-persistent joints.

PFC Based Numerical Modeling Applications to Estimate Rock Mass Strength
Some investigators have resorted to particle flow codes (PFC 2D and PFC 3D ) [62,63,75] to model jointed rock behavior under uniaxial loading [25,[76][77][78][79][80][81][82], bi-axial loading [26] and polyaxial loading [27]. PFC 3D allows one to study the mechanical interaction behavior between intact rock and joints incorporating a significant number of joints without making unrealistic assumptions about the surrounding medium around each joint. In addition, it allows failure through both the intact rock and joints under both tensile and shear modes leading to progressive failure, which usually occurs in jointed blocks with non-persistent joints.

PFC Based Numerical Modeling Applications to Estimate Rock Mass Strength
Some investigators have resorted to particle flow codes (PFC 2D and PFC 3D ) [62,63,75] to model jointed rock behavior under uniaxial loading [25,[76][77][78][79][80][81][82], bi-axial loading [26] and polyaxial loading [27]. PFC 3D allows one to study the mechanical interaction behavior between intact rock and joints incorporating a significant number of joints without making unrealistic assumptions about the surrounding medium around each joint. In addition, it allows failure through both the intact rock and joints under both tensile and shear modes leading to progressive failure, which usually occurs in jointed blocks with nonpersistent joints.
Kulatilake et al. [25] performed pioneering research in providing a realistic calibration procedure for micro-mechanical parameters of PFC 3D for a contact bonded particle flow model. Using this model, they studied the jointed rock behavior of blocks with persistent joints under uniaxial loading. They included spherical particles to model both intact material and joints. In other words, they considered closed flaws. They mentioned the difficulties they had in modeling the joint behaviour of closed flaws. Their focus was on the global mechanical behavior of jointed blocks and possible failure modes.
Lee and Jeon [77] and Zhang and Wong [78] have used PFC 2D to study crack initiation, propagation and coalescence using one or two flaws. In their models, they removed particles to simulate open flaws. They reported the successes, failures and difficulties that they encountered in comparing their PFC results with experimental results.
By selecting appropriate micro-mechanical parameter values through a trial and error procedure, Fan et al. [82] used PFC 3D to study the macro-mechanical behavior of jointed blocks with multi-non-persistent joints with high joint density under uniaxial loading. The focus was to study the effects of joint orientation, the size and joint mechanical properties on jointed block strength, deformability, the stress-strain relation and failure modes at the jointed block level. The uniaxial compressive strength of the block, UCS B , was found to depend heavily on the joint dip angle, β, and joint continuity factor, k (see Figure 7).  [82]).
Based on laboratory and PFC 3D numerical modeling, Kulatilake et al. [26] proposed a rock mass strength criterion for biaxial loading conditions. The developed rock mass strength criterion is given through Equations (5) and (6) in Section 3.3.
In He et al.'s [32] criterion, stated in Section 5.2, for a set of constant values of σ2 and σ3, ω is a constant for a specified rock mass irrespective of the directions of σ2 and σ3. When the σ2 and σ3 directions rotate around the vector normal to the plane of σ2 and σ3 (i.e., around the major principal stress, σ1, direction) the fracture tensor component in the The jointed blocks resulted in four failure modes as follows: (1) splitting failure; (2) plane failure; (3) stepped path; and (4) intact material failure. These failure modes were compared with the corresponding experimental results obtained from Chen et al. [30,31].
The main features of each failure mode and possible relations between the failure modes, UCS B and β and k values are given in Section 5 of Fan et al.'s [82] paper.
Based on laboratory and PFC 3D numerical modeling, Kulatilake et al. [26] proposed a rock mass strength criterion for biaxial loading conditions. The developed rock mass strength criterion is given through Equations (5) and (6) in Section 3.3.
In He et al.'s [32] criterion, stated in Section 5.2, for a set of constant values of σ 2 and σ 3 , ω is a constant for a specified rock mass irrespective of the directions of σ 2 and σ 3 . When the σ 2 and σ 3 directions rotate around the vector normal to the plane of σ 2 and σ 3 (i.e., around the major principal stress, σ 1 , direction) the fracture tensor component in the σ 1 direction, F 11 , remains as a constant. As the first invariant of the fracture tensor (F 11 + F 22 + F 33 ) is always a constant, F 22 + F 33 also remain as constants. Therefore, under the above-mentioned conditions, Equation (8) (10) and (11).
where λ 2 , λ 3 , p 2 , p 3 , q 2 and q 3 are empirical coefficients. In Equations (10) and (11), σ J and σ I are the jointed block strength and intact block strength, respectively. The number of empirical coefficients is high in Equation (10). To reduce the number of coefficients, the second rock mass strength criterion is given with less empirical coefficients, a 2 , a 3 , b 2 and b 3 , as given in Equation (11).
For the details of the used procedures and extensive rock mass strength criteria fitting results, the reader is referred to Mehranpour et al. [27]. For further detail on the PFC3D theory and applications, the readers are referred to the publications mentioned in Section 5.3.

Needed Improvements for the Currently Existing Rock Mass Strength Criteria
Note that Mehranpour et al.'s [27] criterion was developed for a plaster-based synthetic rock mass using persistent joints. The use of only persistent joints was a major limitation of this rock mass strength criterion.
Rock masses may fail under the following modes: (a) failure through intact rock under the tensile or shear mode or a combination of the two modes; (b) separating on discontinuities under tensile mode; (c) sliding on discontinuities under shear mode; (d) a combination of the modes (b) and (c); and (e) different combinations of the modes (a), (b) and (c). It will be an extremely difficult task to develop a 3-D rock mass strength criterion that can capture failures under all the previously mentioned modes.
Even though He et al. [32] and Mehranpour et al. [27] provide the most advanced rock mass strength criteria available in the rock mechanics literature, it is likely that both of those may not capture rock mass failures resulting from all of the aforementioned modes. Therefore, in the future, attempts should be made to extend the Mehranpour et al. [27] rock mass strength criteria to develop a new deterministic 3-D rock mass strength criterion that is capable of capturing failures resulting from all the above-mentioned modes for jointed blocks with non-persistent joints.
It is recommended that future investigations should be performed to include the joint weathering and infilling conditions and ground water conditions in the developed 3-D rock mass strength criteria. The developed 3-D rock mass strength criteria should be applied to different types of rock masses of sedimentary, igneous and metamorphic origins with both orthogonal and non-orthogonal fracture systems to generalize the applicability of the included coefficients of the developed 3-D rock mass strength criteria.
Variability and uncertainty are unavoidable in estimating rock discontinuity geometry networks [2,83,84] and the geo-mechanical properties of intact rock and rock discontinuities [85][86][87]. These aspects were not addressed in developing the He et al. [32] and Mehranpour et al. [27] rock mass strength criteria because those attempts were limited to developing a deterministic 3-D rock mass strength criterion. However, those were found to be huge successes. Now, the time has come to, first, improve the deterministic He et al. [32] and Mehranpour et al. [27] criteria and, then, to extend the deterministic criterion to a probabilistic version by incorporating the aforementioned variability and uncertainty. It will be the first time in the world any research group attempting to develop a probabilistic 3-D rock mass strength criterion with the features mentioned above.

Conclusions
The intermediate principal stress plays a major role with respect to the rock mass strength. However, most of the rock mass strength criteria available in the literature, including the most widely used criterion in rock engineering, do not include the intermediate principal stress. The majority of the available rock mass strength criteria, including the most widely used criterion in rock engineering, assume isotropy.
Therefore, those have a chance of applying only to rock masses that are heavily fractured and do not have preferred discontinuity sets and orientation directions. On the other hand, most rock masses contain a few discontinuity sets, which provides distinct discontinuity orientation directions. Such rock masses exhibit anisotropic rock mass strength. The majority of the available rock mass strength criteria do not include joint size/block size as a parameter in the rock mass strength criterion. These rock mass strength criteria do not have the capability to predict changes of the rock mass strength with the block size.
This author feels that a proper rock mass strength criterion in three dimensions should incorporate the intermediate principal stress, the intensity of all the discontinuity sets and the probability distributions of the orientations and size for all the discontinuity sets explicitly in the strength criterion. Significant progress has been made in this direction (see [2,3,26,27,32,74,82]). The 3-D rock mass strength criteria developed by He et al. [32] and Mehranpour [27] already include the intermediate principal stress and the aforementioned discontinuity geometry parameters explicitly through the fracture tensor components to capture the scale dependent and anisotropic rock mass strength.
However, the He et al. [32] and Mehranpour et al. [27] rock mass strength criteria still require further modifications at the deterministic level to increase their capability in predicting failures resulting from all the different modes for jointed blocks having nonpersistent joints. In addition, the improved deterministic He et al. [32] and Mehranpour et al. [27] criteria require extension to a probabilistic version by incorporating the variability and uncertainty associated with estimation of the discontinuity geometry parameters and discontinuity mechanical properties. Data Availability Statement: Some of the results given in the paper are based on the data available in the corresponding references, which are listed in the reference section.