Strength of shear bands in fluid-saturated rocks: a nonlinear effect of competition between dilation and fluid flow

This study shows the significant and nonlinear effect of the competition between dilation and fluid flow on the shear strength of constrained shear bands in fluid-saturated rocks. This effect is conditioned by the contribution of the pore pressure to the yield stress and strength. The pore pressure is controlled by the dilation of the pore space in the solid skeleton of the shear band during plastic deformation and by squeezing of pores in surrounding blocks by the dilating shear band due to the high stiffness of the host massif. A generalized equation has been derived to describe the dependence of the shear band strength on the ratio of strain rate to fluid flow rate.

A range of laboratory and full-scale geological and geophysical research suggests that irreversible deformation in rock samples and rock massifs is strongly localized in shear bands at different scales, the largest of which are tectonic faults [1][2][3] . These narrow zones not only determine the compliance of rocks in the form of localized relative shear displacements of structural blocks, but control seismic activity of rock massifs. The latter explains the current interest in the mechanical properties of fault zones and the rapid increase in the number of published works in this area.
One of the key mechanical properties is the maximum (or peak) strength of the fault zone under given stress and confinement conditions [4][5][6][7] . Reaching the maximum strength corresponds to a change in the response of the shear band from pervasive strain (strain hardening stage) to strain localization (strain softening stage) 8,9 . The behavior of the shear band at the strain softening stage is largely determined by the degree of consolidation of the gouge. In particular, for shear bands with a mature zone of unconsolidated/granular gouge the hardening-to-softening transition is typically smooth and does not lead to onset of dynamic slip as granular gouge usually exhibits velocity-strengthening frictional behavior [10][11][12] . At the same time, consolidated (lithified and indurated) shear bands behave in a fundamentally different way with rapid stress drop (corresponding to brittle rupture) and transition to dynamic slip accompanied by radiation of elastic waves (velocity-weakening behavior) 9,10,13 . The radiated acoustic/seismic power is determined by the magnitude of the dynamic stress drop, which in turn strongly depends on the peak strength of the shear band. This problem is especially important for faults or fault segments with healed (consolidated) core zones, which typically have high shear strength. Failure of such faults or fault segments has a dynamic character and is accompanied by generation of strong seismic waves [14][15][16] . Therefore, the estimation of maximum shear strength is both a fundamental and practical problem widely discussed in fault and rock mechanics [4][5][6][7] .
The conditions for onset of pervasive inelastic strain and subsequent reaching of the maximum strength of shear bands (including faults) are mainly affected by the pore structure and pore fluid pressure. The pore pressure dynamics is controlled by two interrelated processes [17][18][19][20][21][22][23][24][25][26][27] : (1) fluid flow and (2) pore volume change. The pervasive inelastic deformation of rock is often accompanied by its dilatancy [28][29][30] . The volume of the connected crack-pore space in the rock increases during pervasive shear deformation of the shear band, which leads to decrease of local pore pressure. The reduction of fluid pressure reduces the intensity of the relaxation processes associated with the formation of new discontinuities and coalescence of existing ones. This effect is called dilatancy hardening 18 . In turn, fluid inflow is able to compensate for the pore pressure drop and reduce the effect of strain hardening 27,31 . The ratio of fluid flow rate to strain rate (which governs the dilation rate) determines the specific value of shear strength of shear bands. Note that the influence of the competition between dilatancy and fluid flow on shear strength is strongly pronounced for shear bands that are surrounded by material blocks with a similar permeability to that of the shear band gouge. This particularly corresponds to healed (consolidated) faults where the difference in the porosity and permeability of the principal slip zone (of width 1-10 cm) and surrounding periphery zone (up to several meters wide) is much less pronounced than in faults with a mature zone of unconsolidated gouge 32 .
Conventionally, the effect of pore fluid on the maximum shear strength (hereinafter referred to as strength) of shear bands, including fault zones, has been studied for limiting modes of deformation (very slow and very fast) that correspond to drained and undrained hydrological conditions 33 . This is principally due to the fact that the limiting modes most closely correspond to the long-term creep and short-term dynamic modes of deformation of fault zone regions. Numerous experimental and theoretical works, starting with the classical paper of Brace and Martin 33 , show that the strength of permeable rock samples can increase significantly (up to 30-50%) in transition from the drained to undrained condition 31,[34][35][36][37][38][39][40] . This is explained by the limitation of fluid inflow to the increasing pore space of the incipient shear band and corresponding inhibition of pore pressure drop recovery.
Despite the progress achieved in studying the mechanical behaviour of shear bands (including fault zones) in permeable rocks, the influence of some essential factors is not fully understood. Among them are permeability of the shear band and the blocks and their "fluid holding capacity" (which depends on porosity and size). Another important factor is the elastic stiffness of the bulk of rock surrounding the shear band 12,41 . This factor is traditionally associated with the magnitude of effective normal stress defined as the difference between applied normal stress and local pore pressure 9,27,42 . The magnitude of elastic stiffness reaches its maximum value in the bulk of the rock massif, where the mechanical constraint by surrounding rock strongly limits movement of blocks adjacent to the shear band in the normal direction. In such conditions, dilation of the shear band or core zone of the fault causes progressive compression of adjacent block regions and is accompanied by increase in the value of compressive volumetric stress (and compressive effective normal stress) in the fault zone.
There is still no unambiguous understanding of how the strength of a shear band in the depth of constrained permeable rock massif changes in the transition region between the undrained and drained conditions. The topicality of this question is particularly related to the important and widely discussed problem of using watering to control the dynamics of shear stress relaxation in fault zones 27,43,44 . Keeping this question in mind, we studied the nature and functional form of the dependence of the shear band strength on the ratio of shear strain rate to fluid flow rate under constrained conditions corresponding to faults in rock massifs. The study was performed by numerical modelling of the shear deformation of a fluid-saturated permeable shear band using the discrete element method.

Problem Statement
We considered a 2D model sample consisting of two blocks separated by an interfacial layer (shear band) in the plane strain approximation (Fig. 1). The blocks imitate the regions of the medium adjacent to the shear band, which are less damaged than the shear band, have higher cohesion, and therefore deform elastically under the considered loading conditions. The shear band of width 2 h is a layer of an elastic-plastic dilatant material, which simulates the layer of consolidated gouge in principal slip zones of faults 9,45 . The width of the model blocks was H = 20 h. We used the following reference values of widths of the shear band and blocks: 2 h ≈ 1.5 cm, H = 15 cm. The shear band and blocks were assumed to be permeable and fluid saturated.
The model shear band with surrounding fragments of blocks was numerically simulated by the discrete element method 46-49 using a fully coupled macroscopic model of fluid-saturated porous brittle materials 50,51 . Within this model the discrete elements simulating parts of the shear band and surrounding blocks are treated as porous and permeable. The effect of the fluid contained in the crack-pore space of a discrete element on its stress state is described on the basis of the Biot linear model of poroelasticity [52][53][54] . The inelastic behavior of the permeable brittle material of the discrete element is described using a plastic flow model of rocks with a non-associated flow law and the Mises-Schleicher yield criterion (Nikolaevsky's model) 55,56 :  where β and Y are the internal friction coefficient and cohesion (yield stress under pure shear deformation) of dry material of the discrete element, P pore is the pore pressure, σ mean and σ eq are the mean stress and the equivalent stress within the discrete element volume 51 . Nikolaevsky's model assumes a linear relation between the rates of shear and bulk plastic strains with the dilatancy factor Λ. The dilation of the material of the discrete element leads to an increase in its porosity and permeability 51 . Within the formalism of discrete elements, local failure is modeled by changing the state of a pair of interacting elements from linked (or bonded) to unlinked 49 . As the failure criterion we apply the Drucker-Prager criterion in the following formulation: where λ = σ c /σ t is the strength ratio of "dry" material under uniaxial compression (σ c ) and tension (σ t ). The pressure-dependent Mises-Schleicher and Drucker-Prager criteria are generalizations (smoothing) of the Mohr-Coulomb yield and failure criteria, which are widely used in rock and fault mechanics 4,7,24,38,57 . The influence of the pore fluid pressure P pore on the conditions of onset of inelasticity and fracture is taken into account within these criteria by means of their formulation in terms of effective mean stress ). The initial values of porosity (φ 0 = 0.1) and permeability k 0 of the shear band and the blocks were assumed to be equal. This approximation is consistent with the low gradients of porosity and permeability in central zones of healed (consolidated) faults.
The model sample was simulated by a close-packed ensemble of equally sized discrete elements. Initially, all interacting elements were assumed to be linked to imitate a consolidated shear band. The element size was 5·10 −4 m. Separate tests have shown that further reduction of the element size does not have any significant effect on the simulation results (see Supplementary Materials).
We modelled constrained shear of the sample in the horizontal plane along the X axis ( Fig. 1). Periodic boundary conditions were specified on the lateral faces in the horizontal direction to simulate an infinitely long shear band. The sample was loaded in two stages. At the first stage, a normal load σ N was applied to the upper and lower sample faces. The initial fluid concentration in the pore space of the sample was chosen so as to create the specified pore pressure P pore 0 in the sample. There was no plastic deformation in the sample by the end of the first loading stage. The stress and pore pressure distributions were homogeneous. At the second stage, the sample was subject to simple shear by applying constant tangential velocity V x and zero normal velocity (along the Y axis) to the upper and lower faces to fulfil the constrained shear condition. The sample deformation proceeded until crack initiation in the shear band. The described 2D system models a horizontal cross section of a healed fault between structural blocks of a rock massif at a certain depth. Note that the initial "horizontal" (in the XY plane) stresses in the given formulation of the problem exceed the "vertical" ones. This is consistent with experimental data indicating that horizontal stresses are considerably higher than vertical ones in regions with high deformation activity 58 .
We used the following types of hydrological boundary conditions on the external surfaces of the sample: (1) impermeable boundaries (hydraulically isolated sample) or (2) perfectly permeable (pore pressure on the surfaces was permanently equal to the initial value P pore 0 ). These boundary conditions correspond to the hydrological conditions in the central regions of fault zones in the bulk of low and high permeability host rocks, respectively. Figure 2a shows examples of shear loading curves of the reference sample (2 h ≈ 1.5 cm, H = 15 cm) at different shear strain rates. Here shear stress is calculated as the ratio of shear resistance force, measured at the upper sample face, to the square of the face. It can be seen that the initially elastic deformation transitions into inelastic deformation distributed in the volume of shear band (pervasive shear strain accompanied by dilation). Upon reaching the maximum shear stress (which is a peak shear strength) the shear crack is nucleated and developed in the shear band along the direction of band strike.

Results
The simulation results showed that at high strain rates the magnitude of shear strength tends toward the upper limit, and at low strain rates -toward the minimal value. In the example shown in Fig. 2a the values of shear strength for curves 1 and 4 are close to the upper and lower limits, respectively. Such regularity was first observed by Brace and Martin 33 , and then reported in numerous experimental and theoretical studies of confined compression of rock samples and shear loading of model fault zones 31,[34][35][36][37][38][39] . At the same time, we found that, in general, the reduction of shear strength from the upper to the lower limit with reduction of strain rate is not monotonous. At a certain intermediate strain rate the shear strength reaches a local minimum (curve 2 in Fig. 2a). Further reduction of strain rate leads to an increase of shear strength up to a local maximum (curve 3 in Fig. 2a). At even smaller strain rates the shear strength again decreases down to the lower limit. This result was unexpected, but is confirmed by recent experimental studies 31 In the present study, we varied the shear strain rate V h h /( ) xy x  ε = + , the initial permeability of the blocks and the shear band k 0 (at fixed φ 0 = 0.1), the dynamic fluid viscosity η, and system size (2H + 2h) within several orders of magnitude: ε xy  from 5·10 −4 s −1 to 1 s −1 , k 0 from 10 −18 m 2 to 10 −13 m 2 , η from 2·10 −4 Pa·s to 2·10 −2 Pa·s (dynamic  viscosity of water at room temperature is about 10 −3 Pa·s), (2H + 2h) from 15 cm to 150 cm. We found that the parameter combination xy xy 2 0 unambiguously determines the value of shear strength of the shear band for a given initial mean stress σ mean 0 , pore pressure P pore 0 , and ratio h/H. In other words, shear band zones have the same shear strength if they are characterized by the same value of A xy , (even if the specific values of the parameters k 0 , η, xy ε  , and h differ by orders of magnitude). The parameter A xy has the meaning of strain rate normalized by fluid flow rate, because the combination k 0 /η(H + h) 2 detеrmines the fluid flow rate in conventional equations of fluid mass transfer. Figure 2b shows a typical dependence of the shear strength τ c of the modelled shear band on the parameter A xy for a hydraulically isolated system. Each point of the curve corresponds to a separate calculation at given values ). The region A xy → ∞ (region I in Fig. 2b) corresponds to combinations of k, η, xy  ε and h at which the fluid flow rate is extremely low compared to the rate of pore pressure change caused by pore volume variation. This corresponds to the hydrological conditions close to the undrained condition of the shear band. The region of low A xy values (region III in Fig. 2b) corresponds to low shear rate, low fluid viscosity or high permeability of the blocks. In this region the fluid flow rate is relatively high, and the hydrological conditions for the shear band approach a fully drained (the pore pressure distribution in the sample is close to homogeneous during the entire course of deformation).
The strength τ c of the shear band approaches the absolute maximum at A xy → ∞ (dynamic loading, highly viscous fluid or impermeable blocks) and tends to the absolute minimum at A xy → 0 where permeability of the sample is large enough to provide homogeneous distribution of pore pressure at all stages of deformation. At the same time, the dependence τ c (A xy ) in the interval of intermediate A xy values (region II in Fig. 2b) has unexpected non-monotonic character. The non-monotonicity of the τ c (A xy ) dependence can be explained by analysing the dynamics of pore pressure and effective mean stress distributions across the vertical section of the samples (Fig. 3).
At high A xy , the pore fluid pressure in the shear band drops to zero at the early stage of inelastic deformation accompanied by dilation (curves 1-4 for P pore in Fig. 3a). Due to the mechanical confinement of the system, shear band expansion leads to progressive compression of the blocks, which in turn causes increased compressive mean stress and pore pressure (at the same time, the magnitude of expansion is smaller than in an unconstrained system). Consequently, the absolute value of effective mean stress within the shear band increases during the course of deformation (effective mean stress is compressive, i.e. negative). The magnitude of the effective mean stress in the shear band is up to several times higher than in bordering regions of blocks, and this difference gradually increases with shear strain (curves 1-4 for σ mean eff in Fig. 3a). The maximum absolute value of effective mean stress in the shear band is achieved for the case of negligible fluid flow rate, which results in the maximum shear strength at A xy → ∞ (Fig. 2b).
A reduction of the parameter A xy , for example as a result of a decrease in strain rate or viscosity or an increase in the permeability, leads to increase in the rate of fluid flow from the blocks to the shear band. Moreover, at sufficiently low values of A xy , the pore pressure in the blocks decreases (instead of increasing) under shear deformation (curves 1-4 for P pore in Fig. 3b), despite the dilation of the interface and block compression. Fluid outflow is accompanied by a decrease in the linear dimensions of the blocks. The latter determines the increase in the rate of shear band expansion (curves a and b in Fig. 4) and the corresponding decrease in the absolute values of effective mean stress in the shear band and blocks as compared to an impermeable system (curves 1-4 for mean eff σ in Fig. 3b). Besides this, due to increase of fluid inflow, a substantial gradient of effective mean stress forms in the cross section of the shear band. The minimum of σ mean eff is localized at the center. The described changes in the stress state cause a decrease in the strength τ c of the shear band with decreasing A xy (region I in Fig. 2b). At a certain "threshold" combination of the filtration and deformation parameters A xy = A crit1 , the shear strength τ c reaches its local minimum (point "c" in Fig. 2b).
In the region A xy < A crit1 (region II in Fig. 2), the ratio of fluid flow rate to strain rate becomes large enough to maintain a nonzero fluid pressure in the shear band up to large values of applied shear strain ε xy (curves 1-4 for P pore in Fig. 3c), including the point of failure. In this case, the smaller is the value of A xy , the higher is the pore pressure in the shear band at the same value of applied shear strain ε xy . As the pore pressure rises, the inelastic contribution to the total strain rate increases. Consequently, the dilation rate increases with decrease of A xy . This effect is most pronounced at the early stage of inelastic strain (curves b and c in Fig. 4). In the interval of A xy corresponding to region II in Fig. 2b the increase in the rate of block compression by the dilating shear band overcompensates for the contraction of blocks caused by fluid outflow. As a result, in region II the effective mean stress in the shear band again increases with decrease of A xy (curves 1-4 for σ mean eff in Fig. 3c). Moreover, the gradient of effective mean stress in the cross section of the shear band decreases at the same time. The above considerations explain the observed increase in the strength τ c of the shear band. At some "threshold" combination of the filtration and deformation parameters A xy = A crit2 , the shear strength τ c reaches its local maximum (point "c" in Fig. 2b).
In the region A xy < A crit2 (region III in Fig. 2), which corresponds to an extremely low rate of applied shear strain or high permeability, the fluid flow rate provides nearly homogeneous distributions of pore pressure across the sample and of effective mean stress across the blocks and the shear band (Fig. 3d). A key feature of region III is that the pore pressure in the shear band remains nonzero up to the moment of crack formation and its value increases with decrease in A xy (P pore ≈ 5 MPa or about 30% of the P pore 0 value in the example shown in Fig. 3d). The threshold A crit2 corresponds to a combination of the parameters at which the pore pressure in the center of the shear band reaches zero at the moment of crack initiation. The increase of dilation rate with increasing P pore tends to saturation in region III (curves c and d in Fig. 4). Therefore, increase of the pore pressure in the shear band in region III leads to slower growth of effective mean stress during the course of deformation and hence to decrease in the shear strength. The lower limit τ c is reached at A xy → 0 (fully drained condition for the shear band).
The described qualitative change of the distribution of effective means stress also determines the position of the shear crack in the shear band (Fig. 5). At large and small values of A xy (A xy ≫ A crit1 and A xy < A crit2 , respectively) a pair of symmetrically positioned cracks is formed (formation of a pair of cracks is a consequence of the symmetry of the system). The distance between the cracks and the center of the shear band is approximately 60% of the half-width of the shear band (the upper scheme in the Fig. 5a). In the region I the crack position gradually shifts to the center of the shear band following the decrease in strength of the shear band with decrease of A xy (Figs 2b and 5b). The interval of A xy , where crack shifts from "limiting" position (60% of the half-width of the shear band) to the center of the shear band, approximately amounts to two orders of magnitude. In the region II (within the range A crit2 < A xy < A crit1 ) there is only one crack propagating along the central line of the shear band (the lower scheme in the Fig. 5a). Transition to the region III is accompanied by backward shift of crack position to the "limiting" position. This behaviour is caused by changes in the degree of inhomogeneity of effective mean stress across the shear band with changing A xy . For A xy several orders of magnitude larger than A crit1 the distribution of mean stress in the central part of the shear band is almost homogeneous, while the zone of maximum  stress gradient is located relatively close to blocks (the cracks nucleate in this zone). For values of A xy belonging to the region II (lower than A crit1 ), a narrow zone of lower effective mean stress forms in the centre of the shear band (Fig. 3b). This zone determines the position of the crack in the centre of the shear band. At A xy < A crit2 the new factor appears, which leads to backward shift of the crack position. This factor is nonzero pore pressure in the shear band during the entire course of deformation. In the samples with A xy lying in the region III this factor determines the profile of effective mean stress distribution across the shear band with maximum absolute value in the centre and the zone of maximum stress gradient located relatively close to the boundaries with the blocks.

Discussion
The derived nonlinear and non-monotonic character of the dependence of the shear strength of mechanically constrained shear band on the parameter A xy as well as the explanation of the reasons for this dependence are well supported by the results of recent experimental studies on triaxial deformation of different sandstone 31 and edifice rock 40 samples as well as on cone penetration tests in consolidated silts [59][60][61][62] . In particular, Duda and Renner 31 showed in their detail study that non-monotonic dependence of the strength of permeable fluid-saturated samples on axial strain rate under constant value of nominal effective confining pressure is determined by a change in the internal drainage state. They concluded that effective drainage during deformation is determined by the interaction between the mechanical (deformation of the skeleton) and hydraulic (fluid redistribution) evolution. An important issue supported by issues of other cited experimental papers is that "bulk hydraulic properties rather than local inhomogeneities control the condition of drainage" 31  where Δε inelast denotes an interval of axial strains associated with significant volume changes, s is specific storage capacity, l is a characteristic length of fluid redistribution. Considering ε  crit as an analogue of A crit1 , we can estimate the values of the parameter A crit1 for sandstones studied by Duda and Renner 31 . Substituting the expression (4) for  ε crit into the expression (3) for A xy and applying the values of Δε inelast and s for three experimentally studied sandstones we obtain an estimation of ~Ã s / 10 10 GPa This estimation corresponds within the order of magnitude to the numerically derived value A crit1 for the considered shear band. The all above said confirms our idea that non-monotonic profile of the dependence of strength on the parameter A xy characterizing the ratio of strain rate to fluid flow rate is general for fluid-saturated rock samples and fragments of rock massifs where fracture has a localized pattern in the form of shear band formation or activation.
We have to emphasize that the key factor determining non-monotonic variation of the shear strength of a dilatant shear band at varying ratios of strain rate to fluid flow rate is mechanical constraint by the surrounding bulk of massif. Mechanical constraint in combination with dilation of the shear band leads to increasing compression of adjacent regions of rock and consequently to increase of effective stress. The rate of increase of effective stress determines the normal elastic stiffness of the surrounding massif. As shown in Fig. 2b, the rate of increase of effective stress, and consequently the elastic stiffness, nonmonotonously change with decreasing ratio of strain rate to fluid flow rate, and this change determines the shear strength variation.
The curve shown in Fig. 2b has three characteristic regions. In each region the change of strength is monotonous, which implies the presence of a dominant mechanism, which determines the direction of the change. When passing into a different region, both the dominant mechanism and the direction of the trend are changed.
In region I the dominating mechanism lies in the decrease of the linear dimensions of the blocks due to fluid outflow to the shear band (poroelastic contraction). Decrease in the value of A xy is accompanied by the inflow of a large amount of fluid into the shear band and hence by reduction of the constraint imposed on the shear band by the compressed blocks (effective normal stiffness of the blocks decreases). This mechanism determines the decrease of the shear strength in region I as A xy decreases.
In region II the trend-determining mechanism is tied to the increase of the dilation rate of the shear band with decreasing value of A xy due to slowing of pore pressure reduction and maintaining nonzero pore pressure during most of the shear process. This mechanism provides an increase in the absolute value of effective mean stress in the sample and hence increase in the strength of the shear band in region II as A xy decreases.
In region III the trend-determining mechanism is linked to the fact that pore pressure in the shear band remains non-zero during the entire shear process. In this region pore pressure in the shear band is higher in the samples characterized by lower values of A xy . Because of this fact an absolute value of effective mean stress in a shear band is also lower in the samples with lower A xy . Decrease of an absolute value of effective mean stress leads to gradual decrease in the shear strength down to the absolute minimum at A xy → 0.
The described three parts of the curve τ c (A xy ) have sigmoid profiles, which is the result of the competition between shear band dilatancy and poroelastic contraction of blocks due to fluid outflow. Analysis of the obtained result allowed to formulate the following general dependence of the shear strength of constrained shear band zones on the parameter A xy . The dependence is expressed as a sum of a constant and three sigmoid contributions (curve 1 in Fig. 6a): The third contribution is determined by the influence of non-zero pore pressure in the shear band on effective mean stress. A value of τ 3 can be obtained from a constrained shear test at A xy → 0, measuring pore pressure in the shear band. Special numerical analysis at different P pore 0 and the fixed value of 35 mean 0 σ = − MPa has shown that the amplitude τ 3 is proportional to the value of pore pressure P pore max in "fully drained" state of the shear band (A xy → 0) at the moment of failure. This linear dependence is in good quantitative agreement with the analytical estimate derived from failure criterion (2): λ τ = .
− P 0 5 3 ( 1) pore 3 max (see Supplementary Materials). The second sigmoid contribution is determined by the dependence of the dilation rate of a shear band material on the magnitude of pore pressure. Having experimentally estimated a value of τ 3 , it is possible to calculate τ 2 in the following way: 3 . Another six parameters can be estimated on the basis of fitting the curve to experimental data on the shear strength in intermediate range of A xy . More promising way is a development of an analytical approach allowing estimation of these parameters based on material properties. This is a subject of further research. However, we should mention that analytical estimation ~ε Δ A s / crit inelast 1 makes determining of parameters of the first and the second sigmoidal contributions into Eq. (5) much easier.
The specific parameter values of the Eq. (5) depend on mechanical characteristics of the skeleton of shear band and blocks, bulk modulus of fluid, the h/H ratio, initial mean stress mean 0 σ , fluid content in the blocks and hydrological conditions. An example of constant and sigmoid contributions to shear band strength for the particular case shown in Fig. 2b are drawn in Fig. 6b. Figure 6a shows examples of the τ c (A xy ) dependences for different pore pressure values and hydrological boundary conditions (curves 1-3) approximated by Eq. (5). Comparison of the curves 1 and 2 shows the determining role of fluid content for the non-monotonic pattern of variation of τ c . Decrease in initial pore pressure leads to "degeneration" of the curve to the simple sigmoid function (τ 3 approaches zero, τ 1 and τ 2 decrease) and then to trivial τ c = const at → P 0 pore 0 . The strength of a shear band also strongly depends on the permeability of the host massif enclosing the shear band and adjacent regions of blocks. At a relatively high host massif permeability (curve 3, Fig. 6a) modelled by means of "perfectly permeable" boundary conditions, the amplitudes τ 2 and τ 3 increase, compared with the same parameters for a hydraulically isolated system (curve 1, Fig. 6a), while the parameters of the constant and first contributions remain unchanged. This leads to a sharp increase in the nonmonotonicity of the curve including a pronounced increase in the value of the local maximum and considerably lower value of the minimum at A xy → 0. Finally, we stress that in this study we assumed the permeability of the blocks to be same as shear band permeability, which is reminiscent of the situation in the central zones of lithified fault segments. At the same time, it is a widespread situation that the permeability of the shear band is higher than permeability of surrounding blocks. In particular, permeability of unconsolidated shear band may be several orders of magnitude higher than block permeability. Separate analysis of such systems has shown that the form of the relationship (5) and the values of its parameters do not change if the combination A xy is defined as follows: is block permeability (Fig. 7). So, the dependence τ c (A xy ) is general for shear band zones, where block permeability does not exceed shear band permeability (at the same φ 0 ). This allows to state that it is the permeability of the blocks that controls the shear band strength in rock massifs, all other factors being equal. The shear band strength is directly connected with the length of the stable stage of pervasive shear (Fig. 2a). The results of this study show that the length of this stage can change (and potentially can be managed) as a result of changes in shear rate, fluid viscosity and permeability of the shear band zone and surrounding massif. Depending on the magnitude and sign of change of the parameters that determine the instantaneous value of A xy , this change can be positive (elongation) or negative (shortening). The latter may be especially important for control of the point of transition of the regime of shear of consolidated fault segments from stick to dynamic slip.

Methods
We numerically simulated constrained shear of a model shear band using one of the representatives of the group of discrete element methods (DEM), namely the explicit method of simply deformable elements. The term "explicit" means that evolution of an ensemble of interacting discrete elements is defined by a numerical solution of the system of classical equations of motion using an explicit numerical scheme (the velocity Verlet algorithm in our case) 48 . We use the well-known approximation of equivalent discs in the 2D problem statement for interpreting the element shape when numerically solving equations of motion (this kind of implementation of DEM is often called a distinct element method) 47,48 . This approximation allows to use simplified (Newton-Euler) equations of motion and to divide the interaction between two discrete elements into two independent components: central (oriented along the line connecting the centres of mass of the elements) and tangential (in the plane transverse to the mentioned line). In the simply deformable element approximation the state of a discrete element is determined by the average stress tensor σ αβ and average strain tensor ε αβ , which are calculated through interaction forces, pair overlaps and relative shear displacements in the interacting pairs of the element with its neighbours 47,49 . One important feature of our implementation of the simply deformable DEM is the use of the original many-body formulation of the element-element interaction 49,63 . It allowed us to implement coupled models of poroelasticity, poroplasticity and fluid affected failure within the DEM. Detailed description of the features of the numerical method and model implementation can be found in the papers 49,50 .
The problem of modelling the mechanical behaviour of permeable fluid-saturated material by an ensemble of interacting discrete elements is divided into two subproblems: (1) description of the mechanical behaviour of the enclosing solid skeleton with interstitial fluid; (2) description of fluid transfer in the filtration volume (system of connected channels, pores and microcracks) of the solid. In the considered model we assume the characteristic thickness of pores, channels and microcracks to be smaller than the discrete element size. Therefore, the discrete elements are treated as porous and permeable. The influence of pore fluid in the volume of the element on its mechanical properties and response is described implicitly with use of the models of poroelasticity and poroplasticity.
We used Biot's linear model of poroelasticity to describe the mutual relation between average stresses and strains and fluid pore pressure in the volume of the discrete element. The formulation of the constitutive equation for a fluid-filled porous isotropic elastic solid is used 54 : pore mean where α, β = x, y, z; G and K are the shear and bulk elastic moduli of the material of the element; δ αβ is the Kronecker delta; a is the poroelastic constant proportional to the ratio of bulk moduli of porous and nonporous material, P pore is the average value of pore pressure in the volume of the element. The behavior of the material of the discrete element beyond the yield stress is simulated using Nikolaevsky's rock plasticity model with a nonassociated flow rule and modified Mises-Schleicher yield criterion (1), which handles effective mean stress in the volume of the discrete element σ = σ + P mean eff mean pore instead of conventional mean stress in solid skeleton σ mean . A feature of Nikolaevsky's model is that it postulates a linear dependence between the rates of the bulk and shear plastic strain components with a proportionality factor Λ, termed the dilatancy factor. Nikolaevsky's model is implemented within the DEM formalism with use of the radial return algorithm of Wilkins. According to the algorithm, a numerical solution of an elastic-plastic problem at each step of integrating the equations of motion consists of: (1) incremental solution of an elastic problem and (2) subsequent reduction of the average stress tensor σ αβ in the volume of the discrete element to the yield surface Y if the inequality Φ = βσ + σ > Y / 3 mean eff eq holds true 49 : mean mean where σ′ αβ are the reduced values of the components of the average stress tensor in the volume of the discrete is the reduction factor of stress deviators, N = KΛ(Φ − Y/ (KΛβ + G)) is the correction of the mean stress value calculated in an elastic approximation. Nikolaevsky's plasticity model of brittle materials is a typical representative of nonassociated flow models. Its main limitation is the neglect of possible pore collapse. Therefore, this model is typically applied to brittle materials with relatively low initial porosity (less than 15-20%) at mean stresses below the brittle-ductile transition threshold 64,65 .
The rheological properties of the dry material of a discrete element are assigned through the dependence Y(ε ms ), where ε = ε β + ε K G /3 / 3 ms mean eq is a combination of mean strain ε mean and equivalent strain ε eq (here we imply total values of strains including both elastic and irreversible parts) 49 . The strain parameter ε ms can be conventionally called "modified Mises-Schleicher strain parameter"). The form of parameter ε ms ensures the equality Φ/ε ms = 3G within the region of elastic behavior (ε ms is an analog of equivalent strain in the conventional models of plasticity of metals with von Mises yield criterion). The dependence Y(ε ms ) for dry material (P pore = 0) is easily calculated from available or assigned stress-strain curves under shear or uniaxial compression. The strain hardening modulus Π is a derivative Π = dΦ/dε ms .
The modified failure criterion of Drucker and Prager (2), which takes into account local pore fluid pressure is used as the criterion of bond failure (linked pair → unlinked pair) in pairs of interacting fluid-saturated discrete elements 49 .
Note that the yield and failure criteria (1-2) are generalizations (smoothing) of the Mohr-Coulomb yield and failure criteria, which are widely used in rock and fault mechanics 4,7,24,38,57 . The main distinction of criteria (1-2) is that they take into account the intermediate principal stress along with the major and minor principal stresses. Typically the estimates of yield stress and strength by (1)(2) are quite close to Mohr-Coulomb criteria, and main limitations and working ranges of (1-2) as well as sensitivity of the results to the values of criterion parameters are similar for both kinds of criteria.
The filtration volume (also called the open porosity) Ω pore of the discrete element is considered to be the sum of two components: the volume of the initial system of connected pores, cracks and channels (Ω pore,el ) and the volume of new pores (Ω pore,pl ) formed due to material dilatancy during inelastic deformation (i.e., due to irreversible opening of microscopic discontinuities). Reversible changes of Ω pore,el during elastic deformation of porous material are described within the framework of Biot's linear model of poroelasticity 54 : pore el mean pore 0 0 where Ω 0 is the initial volume of the element, φ 0 is the initial porosity. The irreversible change of pore volume Ω pore,pl due to material dilatancy is expressed through the plastic component of bulk strain of the discrete element θ plast : pore pl plast e last xx yy zz mean pore , 0 0 0 where θ and θ elast are the total bulk strain of the element and its reversible (elastic) part, ε αα are diagonal components of the average strain tensor in the element volume. Permeability k of the material of a discrete element is treated as a function of porosity φ = Ω pore /Ω 0 and the so called "characteristic diameter" of filtration channels d ch In the general case the parameter d ch is not the average transverse size of discontinuities, but the size of the narrowest channels in the crack-pore space, which determine the rate of fluid flow through the volume of the discrete element. , where K fl is the bulk modulus of interstitial fluid, ρ pore = m fl /Ω pore is fluid density in the pore space of the element (m fl is mass of fluid in pore space), ρ 0 is the equilibrium value of fluid density under atmospheric conditions. We assume P pore = 0 if ρ pore ≤ ρ 0 . The fluid pressure gradient is assumed to be a driving force of filtration. We used the classical equation of fluid mass transfer in the "micropore" space 66  This equation is numerically solved using the Euler method on a grid formed by an ensemble of interacting discrete elements. Here grid nodes are the centres of mass of discrete elements. There is no mass transfer between grid nodes in which ρ pore ≤ ρ 0 .
The described model was verified in a series of numerical tests including uniaxial compression of samples of dry model materials with different values of the dilatancy factor Λ, fluid filtration through a thin layer of permeable material and fluid discharge from a sample under rapid uniaxial compression (without failure). Verification results are shown in Supplementary Materials. Data availability. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.