Faulting of rocks as a critical energy process

4 Faulting of rocks is a dominant earth process that governs small-scale fracturing, formation of 5 tectonic plate boundaries, and earthquakes occurrence1–4. Since the 18th century, the mechanical 6 settings for rock faulting were commonly analyzed with the Coulomb criterion5 that offers 7 empirical, useful tools for scientific and engineering applications1,6–12. Here we revisit the 8 processes of rock faulting by an alternative approach that incorporates elastic energy, strain-state, 9 and three-dimensional deformation; these mechanical fundamentals are missing in Coulomb 10 criterion. We propose that a stressed rock-body fails as two conditions are met: (1) The elastic 11 energy generated by the loading system equals or exceeds a critical energy intensity that is required 12 for the faulting process; (2) The internal strain of the stressed rock-body due to slip and dilation 13 along the developing faults equals the strain-state created by the loading system to maintain 14 physical continuity13,14. Our simulations reveal that meeting these energy and strain conditions 15 requires an orthorhombic, polymodal fault geometry that is similar to natural and experimental 16 fault systems15–20. The application of our formulation to hundreds of rock-mechanics 17 experiments11,21–28 provides a new, comprehensive benchmark for rock-faulting. 18


INTRODUCTION
Since Charles-Augustin de Coulomb addressed the issue of soil stability in 1773, his failure criterion has become the prevailing analytical tool of soil and rock mechanics 1,5,6,8 .The rationale of Coulomb criterion is straightforward: a body of soil or rock fails by shear fracturing when the stresses acting on a given fracture reach its strength: where  and n are the acting shear and normal stresses, C is the fracture cohesion, and , indicates the fracture strength dependency on the normal stress, which is known as the coefficient of 'internal friction' 29 .This simple and elegant expression was applied to countless rock failure experiments and adopted for the interpretation of fault systems in the field 4,6 .But the Coulomb criterion cannot explain the observed orthorhombic fault patterns in experiments and field (Fig. 1) [15][16][17][18][19][20] , apparently due to omitting strain and energy considerations that are paramount components of mechanical processes.
In 1913, von Mises proposed that ductile materials yield once the deviatoric stresses reach a critical level 30,31 .This criterion, which is related to the elastic strain energy 32 , was successfully applied to 3-D experimental yielding of ductile materials, mostly metals.However, this criterion ignores the effects of the mean-stress that has strong influence on yielding and shear of brittle solids like rocks 1,9,33 .
Later, Alan Arnold Griffith developed a theory of tensile fracturing by analyzing the stress distribution around flaws in a brittle solid 34 .He showed that the propagation of a critically stressed flaw is controlled by the balance between the elastic energy release and the energy dissipation by the new surface-area of the growing fracture.With this energy-controlled theory, Griffith established the field of fracture mechanics 35,36 .A few attempts to adopt this theory to shear faulting of rocks under compression were partially successful 25,37,38 .Following Griffith's energy-based mechanism, we envision that rock faulting is constrained by critical energy level, and develop here a new approach that incorporates strain-state and elastic energy limit.The applications of this concept to a range of rock-failure experiments provide an energy-based rock failure criterion, and naturally reveal the orthorhombic structure of in-situ observed faults.

HYPOTHESES
We consider an intact, homogeneous, linear elastic, brittle rock, at a quasi-static compressional loading system.The bulk elastic strain energy of the tested rock sample, U, is a product of the applied three-dimensional stresses and the elastic moduli of the rock sample.For simplicity, the analysis refers to the principal stress coordinate system, in which i, and i (i = 1-3), are the principal stresses and the associated strains, respectively.
It is suggested that the rock failure by shear along one or more faults (Fig. 1) is controlled by the conditions formulated by these hypotheses: A. Critical energy: Faulting (i.e., shear fracturing) occurs when U ≥ UC, where U is the elastic energy density generated by the loading system, and UC is the critical energy density necessary for dissipation by the faulting processes.We employ the total elastic energy that incorporates both the mean-stress and the distortion-energy because the failure of brittle rocks strongly depends on the mean-stress and normal stress 1,4,9 .B. Strain compatibility: The internal strain-state is compatible with the external, elastic strainstate, i that is generated by the external loading stresses 13,39 , i.The internal strain is accommodated by slip along the developing faults and by dilation normal to them.

CRITICAL-ENERGY-OF-FAULTING-FAILURE (CEFF)
General elasticity theory indicates that loading of a linear, elastic solid by stress-state, i, can be translated to elastic strain-state, i: where  and E are the Poisson's ratio and Young modulus, and the other two principal strains (2 and 3) are determined by the 2 nd and 3 rd permutations 8 .The elastic energy per unit volume, namely the elastic-energy-density, is defined as one half the product of the stress and the associated strain 8 , U = ½  .Under three-dimensional deformation, U is effectively presented by the stress invariants 32 : where I1 = (1+ 2 + 3) and I2 = (1 2 + 2 3 + 3 1) are the first and second invariants in the principal stress system (Supplementary).The elastic energy, U, depends on the three principal stresses, and we found that it is both simple and unique to relate UC, the critical elastic energy density at failure, to the maximum stress,1 (Methods, Extended Data Fig. 1).
The critical energy hypothesis implies that the elastic energy density at yielding stage, Uc, defines a failure envelope that delineates the transition from intact to faulted states.This relationship is demonstrated for the Shirahama sandstone with 37 failure experiments 11 (Extended Data Table 1, Fig. 2A).The elastic energy density, which is calculated from the experimental stress data and the elastic moduli, displays power relationships, UC = a (1) b , where a and b are constants of 7.3610 -5 and 1.87, respectively.The power law relationship is further demonstrated for 13 published datasets of experimental rock faults (Methods, Extended Data Table 1, Fig. 2B).Most of the experiments were conducted under the general, polyaxial loading , 1> 2 > 3, also called 'true-triaxial ' 23,40 ; the rest of the experiments are of triaxial loading, 1> 2 = 3, and uniaxial loading, 1> 2 = 3 =0 (Extended Data Table 1).The systematic, unique relations of the critical energy to 1 (Methods, Figs.2A, B), demonstrates that UC can be considered as a robust, energy- based faulting criterion.We also noted a clear separation between the eight sets of brittle, strong rocks and the five sets of ductile, relatively weak faults.This separation reflects the effect of the relatively smaller Young moduli of the ductile, weak rock resulting in larger elastic energy for a given 1.We abbreviate the Critical-Energy-of-Faulting-Failure as CEFF, and later discuss its physical source.

INTERNAL STRAIN
The external loading of a rock-body generates an elastic strain-state (Fig. 2C, Equation 2), and hypothesis B implies that at failure this strain is internally accommodated by slip along faults (Fig. 2D).Further, the internal strain by faults' slip should equal the externally applied strain while the energy dissipation along the faults is constrained by the available critical energy (hypothesis A).A general, 3-D strain-state (Fig. 2C) has orthorhombic symmetry (three mutually perpendicular principal strain axes).To accommodate such 3-D strain by fault slip requires four faults in orthorhombic (=polymodal) symmetry 14,41 .This requirement is analogue to the five slip systems in metal plastic deformation 13 .Notably, due to the symmetry of the principal stress axes, an orthorhombic fault system can be reproduced by mirror and rotation actions on a single preferred fault 39 .Similar rotation and mirror actions are commonly utilized to define the conjugate set of Coulomb criterion, and other failure criteria, e.g., von-Mises, or Tresca 42 .
We search for a fault system that satisfy the above two hypotheses.The search is based on simulations of a single preferred fault that is defined by its normal angles (1, 2, 3) with respect to the three-dimensional axes (X1, X2, X3) of the principal stress state (1, 2, 3, Fig. 3A) (Methods).At failure, the rock-body is internally deformed by fault slip along the DS axis and dilation DD normal to it.The faulting energy dissipation is the work done by DS and DD under the corresponding stresses applied on the preferred fault.The preferred fault geometry is derived here by a multi-parameter minimization of the differences between the external elastic deformation (strain, and energy), and the internal strain and energy dissipation by the fault (Methods).Our dedicated Matlab code, OrthoFault, and a few experimental datasets, are available in GitHub (https://github.com/nadavwetzler/OrthoFault.git).
The predicted faults orientations under polyaxial loading (1> 2 > 3) are demonstrated for the Shirahama sandstone that includes 37 failure experiments (Extended Data Table 1, Fig. 3B).The calculated preferred faults are shown by red great circles with black dots for the slip axes.Next, we detail the results of one selected experiment with strain-state of maximum vertical shortening and horizontal extensions (Fig. 3C, D).The orthorhombic faults of this experiment are reproduced by rotation and mirror operations on the preferred fault about the principal planes, and are displayed by blue circles (Fig. 3B) and 3-D view (Fig. 3C).We further calculated the preferred fault orientations in 14 sets with 344 individual experiments (Extended Data Table 1).These calculations indicate that the slip axis of a given fault is not necessarily dip-parallel, strike-parallel, or parallel to a principal stress axis (Fig. 3B, C, Extended Data Fig. 2, Fig. S4).These features contrast the classical Coulomb criterion that predicts a single conjugate set (two pink great circles in Fig. 3B) with fault planes that parallel the 2 axis 6,8 .

ENERGY DISSIPATION DURING FAULTING
We now examine the physical significance of the proposed critical energy density, CEFF, (Fig. 2A, B) by following the rationale of Griffith's theory 34 .According to this theory, tiny flaws within a brittle solid may locally amplify the stresses beyond the local strength, however, a flaw would grow unstably only when the released elastic energy exceeds the dissipated energy by the growing fracture.The dissipated energy is defined as the 'fracture energy', G, and it quantifies the energy required to form new surface-area during growth of a tensile fracture 35,36 .Similarly, the fracture energy of a shear fault quantifies the energy dissipated by the processes of fragmentation, microfracturing, and heating during faulting of an intact rock 43 .The fracture energy of faulting was measured by the integration of the post-failure stage of experimental stress-strain curve 44 (Fig. S1) revealing G = 0.3-1010 4 Jm -2 for experimental faulting 28,44,45 .
In a given experiment, the dissipated-energy-density, WD, can be evaluated by WD = GA/V ( 4 ) where A is the fault area, and V is the sample volume.The dissipated-energy-density, WD, is calculated here for a series of nine triaxial tests on Stripa granite samples 28,46 (Methods).The range of measured fracture energy during faulting in these experiments is 0.3-8.1210 4 Jm -2 (Extended Data Table 2) similar to typical experimental values 44,45 .We calculated WD for these experiments by using the observed fault area A (Insert Fig. 4A, Methods, Supplementary), the reported G, and the sample volume, V.The WD and UC energy parameters, which are independently calculated, show similar, correlated trends when plotted relatively to 1 (Fig. 4A).We argue that this correlation between UC and WD indicates the fundamental significance of hypothesis A: The critical energy generated by the loading system, CEFF, facilitates the failure by providing the energy to be dissipated by the faulting.
The last statement raises a central question: what happens if the local stresses in an experiment exceed the local strength while the total energy is below the critical energy?This question cannot be addressed by Coulomb's criterion, but Griffith's approach may offer an explanation.It is likely that an incipient fault would grow if the local stresses exceed the local strength, however, if U < UC this failure will remain local without crossing the sample.Indeed, local yielding events before rock failure are very common, and experimentally manifested by countless microfracturing, distributed damage and micro-faults as monitored by acoustic emission 47 and microscopic analyses 3 .We show here one example of local faulting during a triaxial experiment on Mount Scott granite (Extended Data Table 1).Fig. 4B displays the stress-strain relationship and the elastic energy density, U, of this experiment.Both curves show two faulting events recognized by distinct stress and energy drops.The first event, with small stress drop, indicates partial, local faulting within the sample, but during this first event, U < UC.We argue that additional damage and energy are needed to rupture the whole sample by a crosscutting fault.

APPLICATION OF FAULTING SIMULATIONS
The OrthoFault code calculations of a given experimental series reveal a linear relationship between the normal stress, n, on the preferred fault, and the shear stress, S, along its slip axis (Methods, Extended Data Table 1, Fig. 4C), Sn, where A and B are constants.This function defines a new strength criterion that is based on strain compatibility and energy constrain.We also noted distribution similarity between shear strength, S, the energy dissipation, CEFF, and the stress/strain invariants (Methods, Extended Data Fig. 4) which indicate that the new strength criterion characterizes the energy dissipation during faulting.
The new strength envelope has similar form to Coulomb criterion (equation 1) and Byerlee law 9 , allowing useful comparisons.First, the B parameter of ( 5) is functionally equivalent to , the mysterious 'internal friction' parameter in (1) 6,7,29 .We thus argue that fault dissipative energy, which is analyzed here (Fig. 4A), is the physical basis of 'internal friction'.Second, our analyses of 13 experimental sets show that Coulomb envelope is systematically higher than our strength envelope (column 5, Extended Data Table 1 and Fig. 2).For example, the Shirahama sandstone set (red line, Fig. 4C) predicts shear strength that is 26-46% lower than the Coulomb strength for these experiments (black line, Fig. 4C) (the later was calculated directly from the experimental 1and 3 8 ).Thus, the CEFF model predicts more conservative rock strength than Coulomb criterion, and this result could bear major implications to rock-engineering designs.
We finally examine the orientations of the preferred faults relatively to experimental observations.For this comparison, we use the results of the Berea sandstone set 22 that is one of the very few sets reporting fault patterns under polyaxial loading (Fig. S4).In this set, the fault systems developed in intact cubic samples that were loaded under polyaxial strain state with principal axes perpendicular to the cubes faces (Fig. S4).The eight faults, which were mapped on the exterior of three samples, trend diagonally to the principal axes.This observation is consistent with orthorhombic fault patterns in nature and other experiments (Fig. 1).The eight faults and the predicted preferred faults for the same experimental conditions (by OrthoFault) share similar orientation as displayed in a combined stereonet and rose diagram (Figs.4D, S4D).This similarity further supports the concept of orthorhombic faulting under polyaxial loading and the ability to predict their orientations by using OrthoFault.
We analyze here rock faulting processes by integrating the mechanical fundamentals of energy, stress, and strain.The analysis application to multiple experimental series revealed a systematic dependence of rock failure state on the critical intensity of the elastic energy.This energy is essential for the faulting process, and it is partly or fully dissipated by slip and dilation along the developing faults.Our analysis provides the theoretical basis for faulting, and predicts the observations of orthorhombic fault patterns in the field and experiments.Lastly, the newly derived strength envelope offers a conservative strength prediction that can be applied to a range of rock-engineering situations.

ACKNOWLEDGMENTS
We benefitted from discussions and comments by Vladimir Lyakhovsky, Einat Aharonov, Yossi H Hatzor, Brett Carpenter, Zonghu Liao, and Amotz Agnon.We thank Alexey Ryabov for providing the Matlab code for the subplot lettering AddLetters2Plots https://www.mathworks.com/matlabcentral/fileexchange/95478),MATLAB Central File Exchange.Retrieved July 9, 2021.The research was partially supported through the project "Rock Mechanics of Reservoir Rock" provided by China University of Petroleum-Beijing (ZR), and the Israeli Science Foundation #363/20 (NW).for the 37 failure experiments (black dots) of polyaxial loading (1> 2 > 3) of the Shirahama sandstone 11 as a function of the maximum normal stress, 1, with fitted power relationship, UC = a (1) b , with a and b as constants 7.3610 -5 and 1.87.B. The calculated experimental relationship of UC and the corresponding 1 for 13 experimental datasets, all of which display similar power relationships (Extended Data Table 1).Note the clear separation between the eight sets of brittle, strong rocks and the five sets of ductile, weak faults.C. Schematic display of the stress and strain states during the loading of a rock block.The externally applied polyaxial stresses (1> 2 > 3) deform the initial sample (dotted prism) into the pre-failure strained sample with strains 1 2  3, that depend on the elastic properties and stresses magnitudes (equation 2) (gray prism).D. Faulting at failure form an orthorhombic pattern with multiple planes (dashed lines), sub-parallel to the principal orientations, deforming the rock-body by the same strain magnitudes of the external loads (3-D green arrows in C) (hypothesis B).Fig. 3. Determination of the expected fault patterns according to the requirements of hypotheses A and B: The internal strain by fault-slip and fault-dilation equals the external elastic strain, and that the dissipation energy during faulting is constrained by the available elastic energy.A. The geometry of a preferred fault (gray triangular) in the principal stresses coordinate system.The fault is defined by the angles 1 2 and 3 of its normal to the principal stress axes X1, X2, X3.  2).Calculations are detailed in Methods.The inset shows a crosscutting fault in the center of the sample (yellow ellipse), and the largest possible fault (dashed red ellipse); see Fig. S2 and Supplementary.B. Stress-strain relations (blue) and contemporaneous elastic energy (red) in a triaxial experiment of Mount Scott granite 24 .Note two faulting events (text).C. The strength envelope as a shear-stress vs normal-stress line (red) that was calculated by OrthoFault code for the Shirahama sandstone experiments 11

Elastic strain energy density, U, and maximum stress,  1.
The critical elastic energy density at failure, CEFF, is a function of the three-dimensional loading stresses.We present CEFF as function of 1 alone (Figs. 2A, B, 4A) as if the other two principal stresses may be ignored.This choice provides unique relationship of critical energy and the stress-state in addition of allowing simple and intuitive presentation.For the demonstration, we consider the experimental stage of rock failure (hypothesis A) in the series of Shirahama sandstone under three-dimensional loading.These experiments display power relations between CEFF and 1 (Fig. 2A), and we examine here the uniqueness of these relationships.
Extended Data Fig. 1 displays, in addition to CEFF, the strain energy density, Ui, for a range of hypothetical triaxial experiments under constant maximum principal stress of 1= 200 MPa, and confining pressures that ranges from CPi =0 (uniaxial) to CPi = MPa.The elastic energy density, Ui, corresponding to each of the hypothetical tests is calculated by using the 1 = 200 MPa, and CPi values, increasing from right-to-left in the secondary coordinate system.The plot (Extended Data Fig. 1) displays the CPi (green dot) and the associated Ui (blue dots with horizontal line).The lowest energy of Ui = 0.93 MJm -3 is associated with CPi = 120 MPa, and it is clearly within the intact rock state.At this state, the energy density level is insufficient to promote failure.Decreasing CPi, under the constant 1, increases Ui until under CPi = 35 MPa the energy density reaches Ui = UC = 1.35 MJm -3 , and both green and blue dots are on the failure curve.It is predicted that under these conditions of 1= 200, 2 = 3 = 35 MPa, a sample of Shirahama sandstone will fail by faulting.For lower confining pressures, the elastic energy is above the critical value of the failure envelope indicating that this rock cannot support this stress state.
Obviously, there are many other paths from the intact-state to the faulted-state, but it is clear that all these paths must cross the UC curve at a single value of 1.Thus, the critical elastic energy density, CEFF, is uniquely defined by the 1 value, and this uniqueness allows using the CEFF-1 power relations in the main text (Figs.2A, B, 4A, S1).

The OrthoFault Code for CEFF analysis
The present critical energy model of rock faulting is based on two hypotheses and the assumption that these hypotheses can be realized by slip and dilation along preferred faults.The preferred faults are determined by numerical search with a dedicated Matlab program.The input file for a given dataset includes a table of the three principal stresses at experimental failure, and the Young modulus and Poisson's ratio parameters of the tested rock (Extended Data Table 1).The calculations are conducted in the following steps: I. Calculate and plot the general input data: principal stresses, the principal strains 1, 2 and 3 (equation 2) at failure, the stress invariants, I1 and I2, and the elastic energy density during failure, UC (equation 3) (Figs.2A, B, Extended Data Fig. 2).II.Guessing an initial preferred fault orientation with slip intensity, DS and dilation, DD.This fault is defined by directional cosines of its normal (N1, N2, N3), and the slip axis is defined by its direction cosines (S1, S2, S3) in the principal stress coordinates.One can show that N1 and N2 are independent, and the other four directional cosines can be determined from N1 and N2 and the strain-ratio 2/1 39 .III. Calculation for the initial preferred fault: a.The normal stress across the initial fault, n, and the shear stress, s, in the slip direction based on the known stress-state (1, 2, 3) (input), and the initial fault parameters (Ni, Si) were defined in step II.b.The internal strain-state due to the slip, DS, and dilation, DD, along the initial fault.c.The dissipated work, W, by the fault slip under s for distance DS, and work against n by dilation DD.IV.Steps II and III are repeated in a grid scanning the full range of N1, N2 (0-1), in search for the minimum DS, and DD, for the best-fit preferred fault.We used the Levenberg-Marquardt nonlinear least-squares method provided by MATLAB.The best-fit fault has these features: (1) the slip and dilation generate an internal three-dimensional strain-state that equals the elastic strain-state generated by the principal stresses; (2) the dissipated work by the slip and dilation is smaller than the critical energy, CEFF; and (3) the slip intensity, Ds, is minimized.V.The calculations of step IV produced one fault, but due to the symmetry of the stress system, four faults in orthorhombic symmetry (Figs. 1, 3C) have the same stress and strain of the preferred fault, and have equal likelihood to form.This situation is equivalent to the development of a conjugate set of faults under 2D strain loading (Reches, 1983a, b).VI.Steps I-V are for one experiment with known stress-state, and calculated strain-state (Fig. 3C) and preferred faults (Fig. 3D).The calculations are then repeated to all experiments in the dataset to produce plots of: (1) the preferred faults with slip axes (Fig. 3B); and (2) a Mohr diagram with s vs n of the preferred faults.This plot represents the strength of the analyzed rock according to the present theory.

Fracture energy and fault slip during failure
We explored the relationships between the critical elastic density, CEFF (hypothesis A, text), and the energy dissipation during experimental faulting 28,44 .In rock failure experiments, the dissipated work during faulting is evaluated by integrating the area under the stress-displacement curve after the peak stress (Fig. S1).This measured mechanical work is considered the shear fracture energy 28,44 .We examine here the results of Hakami 46 and Hakami and Stephansson 28 who presented all the data needed for dissipation calculations: sample dimensions, fault angle, shear fracture energy, G, and fault slip during failure, u (Extended Data Table 2).The triaxial experiments were conducted on Stripa granite while using intact samples and samples with well healed fracture (column 2, Extended Data Table 2).The determined fracture energy in these experiments ranges G = 0.36-8.1210 4 Jm -2 (column 5, Extended Data Table 2) which is similar to the values reported by Wong (1982) of G = 0.3-7.310 4 Jm -2 for experimental measurements (there Table 2).
Hakami and Stephansson 28 used cylindrical samples of 21.8 mm diameter and twice the height.In the experiments the confining pressure was first increased to an initial selected level (column 3A, Extended Data Table 2), and then the confining chamber was locked allowing the pressure to rise with the increasing axial load (initial and final confining pressures in column 3, Extended Data Table 2).The images of the faulted samples (Figs.5.18-5.23 in Hakami 46 ) indicate that the faults cut across the central part of the samples (Fig. S3A-E).Thus, the area of the faults is in the range of 0.5-0.75 of the largest possible fault (Fig. S3F).We calculated the fault area (column 9) by using the fault angle (column 4, Extended Data Table 2), and sample dimensions; the 0.5-0.75range was used to estimate the area error.The total dissipated energy during each experiment (column 10A, Extended Data Table 2) was calculated by using the fracture energy, G (column 5, Extended Data Table 2), and the estimated fault area.Then the dissipated energy per the volume is calculated by using equation 4 (text), WD = GA/V where V is the sample volume.The faulting dissipated energies are listed as total energy during failure (column 10A, Extended Data Table 2), and as energy density (column 10B, Extended Data Table 2).Fig. 4A (text) displays the elastic energy density, CEFF, (column 8, Extended Data Table 2), and the dissipated energy density, WD, (column 10B, Extended Data Table 2) due to the formation of the sample-crossing fault.The significance of the good agreement between the two sets of values, which were independently calculated, is discussed in the main text.
Hakami and Stephansson 28 also reported the faults' slip during failure, u (their Table 1, column 6A in Extended Data Table 2).This central measurement can be compared to our calculated slip along the preferred faults.Our calculations are based on hypothesis B (fault strain equals the elastic strain), and the preferred fault satisfies this strain requirement with the best-fit minimum slip, DST (Methods).The slip magnitude that we calculated for the nine experiments indicate DST = 1.78-6.78mm.However, as our strain and energy calculations are in SI units, the results are implicitly related to slip within a 1m 3 rock-body.Thus, we scaled these slip values to the dimensions of a given sample by using the relations, Expected fault slip (uexpected) = calculated slip (DST)  (sample height) / 1m The scaling to the 0.0436 m sample height of Stripa experiments returned values of uexpected=0.08-0.30mm (column 6B, Extended Data Table , 2).The average measured slip is crudely twice the calculated slip (column 6, Extended Data Table 2), and both the measured and calculated slip display irregular scatter with a weak increasing trend with 1 (Fig. S3).We claim that this agreement, within a factor of two, between the measured and predicted slip values supports the validity of our approach.More refined agreement may be achieved in dedicated experimental work.

CEFF strength criterion
It is shown in the main text that the OrthoFault calculations allow to define a new strength criterion which is based on strain compatibility and energy constrain.This criterion has the linear relationship (Extended Data Table 1, Fig. 4C), Sn, where n is the normal stress, S is the shear stress on the preferred fault along its slip axis, and A and B are constants.We examined here the relationships between the shear strength, s, of the strength criterion, the energy dissipation during failure, CEFF, and the stress/strain invariants by plotting these parameters for the experimental series of Dunham dolomite (Extended Data Fig. 3).The clear similar dependence of both shear strength and energy dissipation on the stress/strain invariants support the notion that the physical control of the strength function (equation 5) is the energy dissipation by faulting during failure.Further, this result offers an explanation for the elusive term of 'internal friction' of Coulomb criterion ( in equation 1).This term has the form of friction, however, common 'friction' between two solid blocks cannot exist before a fault is formed (Anderson).Handin 7 stated: "In cohesive material, however, internal friction is a fictitious quantity that cannot be measured directly."Our analysis, e.g., Extended Data Fig. 3, indicates that 'internal friction' ( in equation 1 and B in equation 5) actually represents the intensity of the energy dissipation during faulting (Fig. Extended Data 3), and it is related to 'friction' only by formula similarity.Extended Data Fig. 1.The relationship between the total elastic strain energy and the magnitude of principal stress is shown for the 37 experiments of the Shirahama sandstone (red dots).
Extended Data Fig. 2. The graphic output of OrthoFault (Methods) simulations of the 523 Extended Data Fig. 3. OrthoFault calculated relations of faulting failure in the experimental series of Dunham dolomite.The first invariants of stress and strain are linearly related to each other, and for this rock A =41.3, and B =0.66.The power fit functions are: Shear stress = 7.5510 5 I1 1.17 , r 2 = 0.75 and Dissipation energy = 3.5510 5 I1 1.98 , r 2 = 0.79 Extended Data Table 2. Experimental and calculated parameters of the triaxial failure experiments on Stripa granite samples.Experimental data from Hakami 46 and Hakami and Stephansson 28 , expected parameters were calculated in the present study.We explored the relationships between the critical elastic density, CEFF (hypothesis A, text), and the energy dissipation during experimental faulting 28,44 .In rock failure experiments, the dissipated work during faulting is evaluated by integrating the area under the stress-displacement curve after the peak stress (Fig. S1).This measured mechanical work is considered the shear fracture energy 28,44,45 .We examine here the results of Hakami 46 and Hakami and Stephansson 28 who presented all the data needed for dissipation calculations: sample dimensions, fault angle, shear fracture energy, G, and fault slip during failure, u (Table S1).The triaxial experiments were conducted on Stripa granite while using intact samples and samples with well healed fracture (column 2, Table S1).The determined fracture energy in these experiments ranges G = 0.36-8.1210 4 Jm -2 (column 5, Table S1) which is similar to the values reported by Wong 44 of G = 0.3-7.310 4 Jm -2 for experimental measurements (there Table 2).
Hakami and Stephansson 28 used cylindrical samples of 21.8 mm diameter and twice the height.In the experiments the confining pressure was first increased to an initial selected level (column 3A, Table S1), and then the confining chamber was locked allowing the pressure to rise with the increasing axial load (initial and final confining pressures in column 3, Table S1).The images of the faulted samples (Figs.5.18-5.23 in Hakami 46 ) indicate that the faults cut across the central part of the samples (Fig. S2A-E).Thus, the area of the faults is in the range of 0.5-0.75 of the largest possible fault (Fig. S2F).We calculated the fault area (column 9) by using the fault angle (column 4, Table S1), and sample dimensions; the 0.5-0.75range was used to estimate the area error.The total dissipated energy during each experiment (column 10A, Table S1) was calculated by using the fracture energy, G (column 5, Table S1), and the estimated fault area.Then the dissipated energy per the volume is calculated by using equation 4 (text), WD = GA/V where V is the sample volume.The faulting dissipated energies are listed as total energy during failure (column 10A, Table S1), and as energy density (column 10B, Table S1).Fig. 4A (text) displays the elastic energy density, CEFF, (column 8, Table S1), and the dissipated energy density, WD, (column 10B, Table S1) due to the formation of the sample-crossing fault.The significance of the good agreement between the two sets of values, which were independently calculated, is discussed in the main text.
Hakami and Stephansson 28 also reported the faults' slip during failure, u (their Table 1, column 6A in Table S1).This central measurement can be compared to our calculated slip along the preferred faults.Our calculations are based on hypothesis B (fault strain equals the elastic strain), and the preferred fault satisfies this strain requirement with the best-fit minimum slip, DST (Methods).The slip magnitude that we calculated for the nine experiments indicate DST = 1.78-6.78mm.However, as our strain and energy calculations are in SI units, the results are implicitly related to slip within a 1m 3 rock-body.Thus, we scaled these slip values to the dimensions of a given sample by using the relations, Expected fault slip (uexpected) = calculated slip (DST)  (sample height) / 1m The scaling to the 0.0436 m sample height of Stripa experiments returned values of uexpected=0.08-0.30mm (column 6B, Table, S1).The average measured slip is crudely twice the calculated slip (column 6, Table S1), and both the measured and calculated slip display irregular scatter with a weak increasing trend with 1 (Fig. S3).We suggest that this agreement, within a factor of two, between the measured and predicted slip values supports the validity of our approach.More refined agreement may be achieved in dedicated experimental work.
Table S1.Experimental and calculated parameters of the triaxial failure experiments on Stripa granite samples.Experimental data from Hakami 46 and Hakami and Stephansson 28 ; expected parameters were calculated in the present study.We compare the experimentally observed fault orientations during failure of intact cubic samples of Berea sandstone 22 (Fig. S4 A, B) with the orientations of the preferred fault calculated with OrthoFault (text).In these experiments, intact cubic samples were loaded under polyaxial strain state with principal axes perpendicular to the cubes faces (Fig. S4 A) were loaded under polyaxial stress state with principal axes perpendicular to the cubes faces.We mapped eight faults on the exterior faces of three faulted cubes (Fig. S4 B), and calculated the fault inclinations in the principal axes coordinates.The orientations of the preferred faults predicted by OrthoFault for the three experiments are displayed together with the experimental observations in a combined plot of stereonet and rose diagram (Figs.S4 C, 4D).The critical-energy concept leads to another interesting deduction.The analysis refers to the elastic-energy-density that is total elastic energy, UT, per unit volume.Consider a rock-body with characteristic dimension, L (Fig. S2F).The elastic energy density is proportional to the applied stresses, however, while he total energy depends on the volume, UT  L 3 , the dissipated energy on the fault surface is proportional to L 2 .Consider now doubling the sample characteristic size from L to 2L, and maintaining the same stresses, shape and material.The elastic energy density remains the same, the total energy increases by a factor of 8, but the area of a cutting fault increases by a factor of 4; thus, the total elastic energy available for faulting increases by a factor of two.Therefore, as relatively more energy is available for faulting in the larger sample, it is expected that failure could occur under lower stresses.This strength dependence on sample size was recognized by Sih 48 who stated: "Specimen volume or size has been known to affect the mechanical and fracture properties of metals.Tests have shown that the apparent strength of a larger size specimen is lower than a smaller size specimen of the same geometry and material."The inverse relationship between strength and size was also demonstrated in rock-mechanics experimentally (Fig. S5); yet, one should note the experimental difficulty to experiment with large rock samples, e.g., L > 20 cm 1 .
We also note that the inverse strength-size relation may be explained differently.Following Griffith 34 , it is assumed that fracture grows from the most critically stressed flaw, and that the local stress amplification is proportional to the flaw size.As already demonstrated by Griffith, a larger sample is likely to contain larger flaws, and thus it is reasonable that larger samples would fail under lower stresses.We thus conclude that the strength-size relations can be attributed to two effective mechanisms: larger flaws in larger samples and more available energy in the larger sample.Figs.S6-9 display the graphic output of OrthoFault (Methods) calculations of Westerly granite, KTB amphibolite, Maha-Sarakham salt and Mount-Scott granite, respectively.The legend of all figures is the same as for Extended Data Fig. 2.
A. Stereonet plot of the preferred faults (red great circles) with the related slip axis (black dot).The orthorhombic fault set (blue circles) is displayed with the corresponding slip direction (light blue dots) is produced by rotation and mirror of the calculated preferred one (text).The conjugate fault set of Coulomb criterion is shown by two black circles.
B. The 'critical energy', UC, as a function of the maximum normal stress, 1, with fitted power relationship, UC = a (1) b , and the calculated dissipated energy by fault slip and dilation.
C. The calculated strength envelope as a shear-stress vs normal-stress (red line) for the 37 experiments (red dots).The Coulomb strength envelope for these experiments (black line) was calculated by using 1vs 3 relations shown in I. D. Mohr plot of the loading principal stresses at failure.E. The orthorhombic fault set (four dashed blue circles) of for one selected preferred fault (blue great circle) with the corresponding slip direction (fink arrows).
F. The calculated orthorhombic fault set (dashed blue in E) is shown in 3-D block view with corresponding strain-state displayed by three pairs of arrows in the center, and the corresponding strain and stress Mohr diagrams (G).
Fig. 4. Investigation of the critical energy model of rock-faulting, CEFF. A. The elastic energy density at failure, CEFF, in triaxial experiments of Stripa granite as function of 1 (red dots), and the dissipated-energy-density, WD, due to faulting (equation4) in these experiments (blue squares) (Methods, Extended Data Table2).Calculations are detailed in Methods.The inset shows a crosscutting fault in the center of the sample (yellow ellipse), and the largest possible fault (dashed red ellipse); see Fig.S2and Supplementary.B. Stress-strain relations (blue) and contemporaneous elastic energy (red) in a triaxial experiment of Mount Scott granite24 .Note two faulting events (text).C. The strength envelope as a shear-stress vs normal-stress line (red) that was calculated by OrthoFault code for the Shirahama sandstone experiments11 (red dots).The Coulomb strength envelope for these experiments (black line) is shown for comparison (text).Note that the present predicted shear strength that is 26-46% lower than the Coulomb strength.D. Observed faults in three samples of Berea sandstone under polyaxial strain loading (Fig. S5), and calculated faults by OrthoFault code.Insert is a block-diagram of the experimental faults in sample RD-42.Fault orientations are displayed on a combined stereonet plot and rose diagram (NW quarter).The stereonet shows poles to eight observed faults in three experiments (red dots), and poles to preferred faults predicted by OrthoFault for these three experiments (purple dots).The rose diagram shows the dip-direction of the same data (experimental in blue sectors, and preferred faults in red sectors).
(red dots).The Coulomb strength envelope for these experiments (black line) is shown for comparison (text).Note that the present predicted shear strength that is 26-46% lower than the Coulomb strength.D. Observed faults in three samples of Berea sandstone under polyaxial strain loading (Fig. S5), and calculated faults by OrthoFault code.Insert is a block-diagram of the experimental faults in sample RD-42.Fault orientations are displayed on a combined stereonet plot and rose diagram (NW quarter).The stereonet shows poles to eight observed faults in three experiments (red dots), and poles to preferred faults predicted by OrthoFault for these three experiments (purple dots).The rose diagram shows the dip-direction of the same data (experimental in blue sectors, and preferred faults in red sectors).
Fig. S1.Method of measuring fracture energy in faulting experiment.Left: Schematic view of the stress-strain relations of a sample in a typical triaxial failure test.Right: The fracture energy, G = 0.5 (p -f)  u) in the integrated (grey) area; the peak (p) and final (f) shear stresses vs the fault displacement, u) (converted from the sample curve, left).

Fig. S4 .
Fig. S4.Fault systems that developed during failure of intact cubic samples of Berea sandstone that were loaded under polyaxial stress state with principal axes perpendicular to the cubes faces 22 .A. Photos of the exterior faces of three faulted cubes.All faults in the 12 rock faces trend diagonally to the sides of cubes, excluding faults at the cubes' edges due to edge effect.B. Traces of the eight main faults in these experiments that were used to calculate the fault inclinations in the principal axes coordinates.C. A combined plot of stereonet and rose diagram (NW quarter) of the observed faults (in A and B) and the preferred faults predicted by OrthoFault.The stereonet shows: Poles to observed faults in three experiments (A) (red dots), and poles to preferred faults predicted by OrthoFault for these three experiments (purple dots).The rose diagram shows the poles to dip-direction of the same data (experimental faults in blue sectors, and preferred faults in red sectors).
Fig. S5.Sample size effect on rock strength.Experimental data indicate up to 70% weakening per decade of length increase of the tested sample (after Fig. 9 in Lockner 1 ).