Transient and Steady‐State Friction in Non‐Isobaric Conditions

The frictional properties of faults control the initiation and propagation of earthquakes and the associated hazards. Although the ambient temperature and instantaneous slip velocity controls on friction in isobaric conditions are increasingly well understood, the role of normal stress on steady‐state and transient frictional behaviors remains elusive. The friction coefficient of rocks exhibits a strong dependence on normal stress at typical crustal depths. Furthermore, rapid changes in normal stress cause a direct effect on friction followed by an evolutionary response. Here, we derive a constitutive friction law that consistently explains the yield strength of rocks from atmospheric pressure to gigapascals while capturing the transient behavior following perturbations in normal stress. The model explains the frictional strength of a variety of sedimentary, metamorphic, and igneous rocks and the slip‐dependent response upon normal stress steps of Westerly granite bare contact and synthetic gouges made of quartz and a mixture of quartz and smectite. The nonlinear normal stress dependence of the frictional resistance may originate from the distribution of asperities that control the real area of contact. The direct and transient effects may be important for induced seismicity by hydraulic fracturing or for naturally occurring normal stress perturbations within fault zones in the brittle crust.

In its simplest form, the frictional strength of rocks is often described by a linear dependence on effective normal stress (Amontons, 1699;Coulomb, 1785) where τ y is the yield strength, μ0 is the static coefficient of friction, σ is the effective normal stress accounting for the Terzaghi effect of pore fluid pressure, and c0 is the cohesion.Equation 1 is often used as a static failure criterion in engineering or geophysical applications, oftentimes with c0 = 0 (Baus & Hansen, 1980;Beeler et al., 2000;Freed, 2005;Harris & Simpson, 1998;Im et al., 2021;Stein et al., 1997;Toda et al., 1998).However, the yield strength for a wide range of rocks from atmospheric pressure to gigapascals indicates a nonlinear dependence on normal stress (Barton, 1973(Barton, , 1976;;Byerlee, 1978;Handin, 1969;Hoskins et al., 1968;Jaeger et al., 2009;Murrell, 1965), such that the linear form of Equation 1 is only approximately valid for a limited range of normal stress (Figure 1).Most rocks, with the notable exception of clays, exhibit an effective friction coefficient of 0.6 ± 0.05 for normal stress above 500 MPa.However, the ratio of shear to normal stress shows increasing scatter at lower confining pressures, with values ranging from 0.25 to 8 (Byerlee, 1978).As faults are thought to operate at low stress due to the presence of over-pressurized fluids (e.g., Faulkner & Rutter, 2001;Sibson, 1994;Suppe, 2014), understanding the nonlinear dependence on normal stress is paramount.Another puzzling observation is the transient behavior of friction associated with perturbations in normal stress (Hong & Marone, 2005;Kilgore et al., 2012Kilgore et al., , 2017;;Linker & Dieterich, 1992;Prakash & Clifton, 1993;Richardson & Marone, 1999;Shreedharan et al., 2019;W. Wang & Scholz, 1994).The friction coefficient directly decreases upon a sudden increase in normal stress, followed by a transient recovery toward steady-state.The micromechanics underpinning these effects are currently unknown, hindering accurate representations of dynamic changes in normal stress in dynamic models of induced or naturally occurring seismicity.
In this study, we derive a constitutive model of the slip-rate and state dependence of fault friction that incorporates the nonlinear effects of normal stress.We show that the nonlinear strength versus normal stress relationship revealed by Byerlee (1978) and the direct and transient responses to normal stress perturbations first highlighted by Linker and Dieterich (1992) both emanate from the sensitivity of the real area of contact to normal stress imaged by Dieterich and Kilgore (1996).In Section 2, we describe the physical assumptions of the model and the implications on frictional strength.In Section 3, we compare the model with laboratory observations, explaining first the range of apparent friction coefficients at low confining pressure for a large collection of rocks (Barton, 1973;Byerlee, 1978;Handin, 1969;Hoskins et al., 1968), and, next, the direct and evolutionary effects associated with sudden changes of normal stress for Westerly granite in bare contact (Linker & Dieterich, 1992) and synthetic gouge of quartz and smectite (Hong & Marone, 2005).In Section 4, we compare the proposed physical model with previous friction laws.The constitutive model augments and complements previous formulations limited to isobaric conditions (Barbot, 2019a(Barbot, , 2022(Barbot, , 2023)), enabling an increasingly comprehensive representation of fault mechanics during seismic cycles.

Constitutive Model
The constitutive behavior of frictional interfaces is governed by the microphysics of solid contacts.For bare contacts, the rough topography of natural surfaces creates contact junctions that concentrate the shear and normal loads (Bowden & Tabor, 1950, 1964;Engelder & Scholz, 1976;Scholz & Engelder, 1976).As a result, a small proportion of the nominal area of contact is in real contact (Ben-David et al., 2010;Ben-David & Fineberg, 2011;Dieterich & Kilgore, 1994, 1996).These results can be generalized for gouge friction, where grain-to-grain contacts occur within a thin volume instead of a single surface (Barbot, 2019a).The macroscopic frictional strength is controlled by the real area of contact (Baumberger et al., 1999;Rubinstein et al., 2004;Selvadurai & Glaser, 2015, 2017), which depends directly on the size and shape of micro-asperities and the ambient shear and normal stress conditions (Maegawa et al., 2015;Sahli et al., 2018;Sleep, 2006).
Theoretical estimates of the real area of contact depend on the assumed texture of the fault interface, which dictates the spatial arrangement and size of micro-asperities.The real area of contact for an elastically deforming isolated perfect sphere of radius d under normal stress σ is given by A r ∼ σ 2 3 d (Hertz, 1881;Timoshenko & Goodier, 1951), indicating a linear relationship with the size of micro-asperities and a nonlinear relationship with normal stress.For a rough sphere with a fractal distribution of asperities of decreasing size, the real area of contact becomes A r ∼ σ 26 27 d 1 9 (Archard, 1957), indicating an almost linear dependence on normal stress, with a deviation of the power exponent from linearity by a factor 1/27. Incorporating a large-scale roughness of the surface, so that  (Byerlee, 1978).Coefficients of friction of 0.7 and 0.6 are shown for reference.(b) Yield strength of rock for effective normal stress up to 400 MPa.There is no cohesion at vanishing normal stress.Coefficients of friction of 1.0, 0.7, and 0.5 are shown for reference.(c) Effective friction coefficient μ = τ/ σ for effective normal stress up to 2 GPa.A wide range of values between 0.25 and 1.5 can be found at low confining pressure.(d) Effective coefficient of friction for effective normal stress up to 400 MPa.Values below 0.5 and above 1.0 can be found below 200 MPa.In all panels, the superlinear (long dashed red line), linear (solid red line), and sublinear (short dashed red line) yield strength profiles correspond to steady-state friction with pβ > qα, pβ = qα, and pβ < qα in Equation 10, respectively.identical spheres come in contact starting from random heights, leads to A r ∼ σd 1 2 (Greenwood & Williamson, 1966), indicating a linear dependence of the real area of contact on normal stress on average.
In general, if the micro-asperities at contact junctions have a rough surface, vary in size and shape, and enter into grain-to-grain contact from different starting positions due to long-wavelength roughness of the fault zone or the randomness of force chains within a gouge layer (e.g., Daniels & Hayman, 2008), the area of contact becomes asymptotically linear with normal stress.We define the real area of contact density A = A r / A 0 as the ratio of real to nominal area of contact.To capture the second-order effects of effective normal stress, we can assume the following general expression where μ 0 = χ/χ n is a reference friction coefficient defined in adhesion theory as the ratio of the plowing hardness to indentation hardness (Bowden & Tabor, 1950, 1964), d and d 0 = 1 μm are the size of micro-asperities at contact junctions and a reference size, respectively, with the power exponent α ≪ 1.The possible nonlinear relationship with normal stress is captured by introducing the deviation of normal stress from a reference value σ 0 appropriate for the experimental conditions.The power exponent 0 ≤ β ≤ 1 is an adjusting variable with β = 0 corresponding to the Greenwood and Williamson end-member, β = 1/27 corresponding to the Archard model, and β = 1/3 corresponding to a Hertzian contact.For the area of contact to increase with normal stress places the hard bound β < 1.In principle, the area of contact density saturates at sufficiently large normal stress, which is not captured by Equation 2. However, saturation occurs at confining pressures found well below the brittle-ductile transition, so this is not a concern for geophysical applications.More complex expressions that capture the saturation at large confining pressure are used in engineering (e.g., Barton, 1976;Persson, 2016;Persson et al., 2004;Yastrebov et al., 2015).Equation 2 also ignores the effect of shear stress on the area of contact (Maegawa et al., 2015;Sahli et al., 2018).
The connection between the micromechanics of contacts and the macroscopic shear strength of the frictional interface is obtained by assuming the yield strength τ y = χA, following adhesion theory (Bowden & Tabor, 1950).The rate-dependent failure criterion in isothermal conditions can be written where V and V 0 = 1 μm/s are the instantaneous slip-rate and a reference value, respectively, and n ≫ 1 is the stress power exponent.Equation 3 predicts rapid sliding when the shear stress reaches or exceeds the yield strength, vanishing slip-rate at shear stress much lower than yield strength, and stationary contact in the absence of shear stress.Combining Equations 2 and 3, we obtain the slip-rate and state-dependent friction law where the state variable is the apparent size of micro-asperities at contact junctions.The parameters V 0 , d 0 , and σ 0 in Equations 2-4 are merely introduced to make the power-law terms dimensionless.These parameters define the choice of unit or relevant order of magnitude for the respective dynamic parameters.Accordingly, these constants do not embody constitutive parameters per se and cannot be constrained individually.The size of micro-asperities increases at stationary contact or sufficiently low sliding velocity, leading to frictional healing, but decreases with shear due to contact rejuvenation.The competition between healing and weakening enables seismic cycles.
Healing of the interface occurs by compaction creep under normal loading.The presence of pore space surrounding contact junctions induces a large differential stress that drives local viscoelastic deformation dominantly in the fault perpendicular direction, even in stationary contact.Hence, the rate of healing depends on normal stress.The evolution law that captures these effects in isothermal, non-isobaric conditions for a single healing mechanism is given by the additive formulation where G = (1 μm) p /s is a reference rate of growth, p is the power exponent for the size sensitivity of healing, and q is the power law exponent for the sensitivity of healing to normal stress.Equation 5does not capture the minimal particle size associated with the grinding limit (e.g., Sammis & Ben-Zion, 2008).The terms G and σ 0 are introduced to satisfy the physical units of the various terms and cannot be determined independently of each other.
The comminution rate is proportional to plastic strain rate, which is controlled by the shear zone width h, relative to a reference strain 1/λ.The constitutive framework simplifies to the isothermal and isobaric case for β = 0 and q = 0. Equation 5 simplifies to the aging law of Ruina (1983) in isobaric conditions, as shown in Appendix A, so it is referred to as the aging-law end-member evolution law.An alternative expression uses the difference of the logarithm of the healing and weakening terms, leading to the multiplicative form Equation 6 corresponds to the slip-law end-member evolution law as it simplifies to the slip law of Ruina (1983) in isobaric conditions, as explained in more details in Appendix A. Equations 5 and 6 feature different evolutionary effects (Beeler et al., 1994;Bhattacharya et al., 2017;Nakatani, 2001) and different stability conditions upon perturbations (Perfettini et al., 2001), but produce the same steady-state solution for ḋ = 0. Considering either the evolution law of Equation 5or Equation 6, the steady-state size of micro-asperities obeys the relationship showcasing a velocity and normal stress dependence, and where we have defined the reference micro-asperity size The size of micro-asperities initially does not depend on normal stress.However, the steady-state size of microasperities, after undergoing plastic deformation, depends on the applied load.Hence, the direct and steady-state frictional dependence on normal stress may differ.
The stability of frictional sliding can be inferred by examining the steady-state frictional strength.The instantaneous frictional strength can be obtained by inverting Equation 4, to get the multiplicative form The steady-state friction coefficient is obtained by combining Equations 7 and 9 and dividing by the normal stress, to obtain Steady-state velocity-weakening behavior is obtained for p < αn, velocity neutral for p = αn, and velocity strengthening for p > αn.Even though the normal stress power exponent can be positive or negative, enforcing an overall increase in frictional strength with normal stress requires qα/p β < 1.The steady-state friction coefficient is independent of effective normal stress for pβ = qα.Otherwise, there is a weak sensitivity of the friction coefficient to normal stress, with superlinear and sublinear strength envelopes for pβ > qα and pβ < qα, respectively (Figure 1).

Comparison With Laboratory Observations
The constitutive framework described in Section 2 describes the evolution of fault friction in isothermal conditions for a limited range of slip-rates in which the same healing and deformation mechanisms dominate.In this section, we describe the evolution of the area of contact and steady-state friction as a function of normal stress and then investigate the transient effects following normal stress change.We first consider the area of contact as a function of normal stress and surface roughness for transparent materials (Dieterich & Kilgore, 1996).The area of contact for acrylic, calcite, glass, and quartz is positively correlated with normal stress, but also affected by the grit number used during surface preparation (Figure 2).We analyze these data using Equation 2 assuming an inverse relationship between grit number and the representative size of micro-asperities.As the relationship is assumed up to an unknown scaling factor, the data can only be used to determine the power-law exponents α and β, which are obtained by a linear least squares procedure.The best-fitting exponents are α = 0.12 and β = 0.10 for acrylic, α = 0.25 and β = 0.14 for calcite, α = 0.16 and β = 0.41 for glass, and α = 0.12 and β = 0.22 for quartz, with uncertainties indicated by the number of significant digits.The measured area of contact generally follows the predictions of Equation 2 for these materials, confirming the nonlinear effect of normal stress on the contact area.These results generalize a previous formulation based on a linear approximation (Barbot, 2019a).
Next, we consider the frictional strength of igneous, metamorphic, and sedimentary rocks curated by Byerlee (1978).The available data on maximum shear strength correspond to finely ground surfaces, initially totally interlocked surfaces, or irregular faults produced in initially intact rocks from atmospheric pressure to 2 GPa (Figure 1).The mechanical data can be approximated by piecewise linear forms such as Equation 1 with different coefficients (Byerlee, 1978).Inspection of the effective friction coefficients, defined as μ = τ/ σ, reveals a diverging pattern across the range of normal stress relevant to crustal conditions, below 400 MPa.Some of the rocks exhibit an increasing effective friction coefficient with diminishing normal stress, with coefficients as high as μ = 1.4 at 50 MPa normal stress.In contrast, other rocks show a decrease in the effective friction coefficient with diminishing normal stress, with coefficients as low as μ = 0.35 at 20 MPa normal stress.The effective friction coefficient seems to converge to μ = 0.6 at normal stress above 500 MPa.
We explain the range of observations using the frictional resistance at steady state given by Equation 10 at a nominal velocity of V = V 0 corresponding to the onset of failure.Accordingly, the friction coefficient simplifies to which involves two degrees of freedom, as μ 0 and σ0 trade-off with each other and cannot be determined individually.We obtain three end-member behaviors based on the choice of constitutive parameters.For pβ > qα, the friction coefficient decays with increasing normal stress, with unbounded values at vanishing normal stress, converging to μ ss = μ 0 at σ = σ 0 .For pβ < qα, the friction coefficient increases with normal stress.The steadystate friction coefficient is invariant for qα = pβ.The envelope of the frictional data can be explained by the superlinear and sublinear strength profiles with σ 0 = 200 MPa.The upper bound corresponds to αq/p β = 0.16.The lower bound is given by αq/p β = 0.078.The mechanical data can be explained by intermediate values of the constitutive parameters, which likely represent the natural variability in nature.
We now examine the yield strength for individual rock types.We begin with the mechanical data for Bilsthorpe mudstone and silty mudstone, Hucknall shale, and Ormonde siltstone obtained under triaxial compression (Barton, 1973).All samples exhibit a nonlinear relationship between the maximum strength versus effective normal stress between atmospheric pressure and 60 MPa revealing a decrease in the apparent friction coefficient with increasing normal stress.These measurements are compatible with the predictions of Equation 11with βp > αq (Figure 3).Using a reference normal stress σ 0 = 10 MPa, we use least squares to obtain the remaining best-fitting parameters assuming a 5% uncertainty on the yield stress measurements.We explain the individual datasets with μ 0 = 0.70 ± 0.08 and αq/p β = 0.19 ± 0.04 for Bilsthorpe mudstone, μ 0 = 0.86 ± 0.08 and αq/ p β = 0.30 ± 0.03 for Bilsthorpe silty mudstone, μ 0 = 0.91 ± 0.08 and αq/p β = 0.25 ± 0.03 for Hucknall shale, and μ 0 = 0.93 ± 0.08 and αq/p β = 0.31 ± 0.03 for Ormonde siltstone.The inferred reference coefficients of friction and power-law exponents are negatively correlated.
We continue with the mechanical data for Wombeyan marble, trachyte, and Gosford sandstone at low confining pressure, up to 7 MPa, in a double-direct shear setup (Hoskins et al., 1968).All samples exhibit a diminishing effective coefficient of friction (Figure 4).We consider a reference normal stress of σ 0 = 10 MPa for these experiments at low pressure.The strength of Wombeyan marble, with an effective coefficient of friction varying from 1.65 to 1.0 can be explained with μ 0 = 0.76 ± 0.23 and the power exponent αq/p β = 0.36 ± 0.20.The data for trachyte is explained by μ 0 between 0.51 and 0.59 and αq/p β between 0.29 and 0.25 depending on surface roughness.The data for Gosford sandstone corresponds to μ 0 = 0.49 ± 0.24 and the power exponent αq/ p β = 0.18 ± 0.23.The inferred reference coefficients of friction and power-law exponents exhibit a negative correlation coefficient between 0.83 and 0.91.
Next, we consider the high-pressure mechanical data for Tennessee limestone and Know dolomite based on triaxial compression tests from 50 to 620 MPa (Handin, 1969).In this pressure range, both rocks showcase a mild decrease in the effective coefficient of friction from 0.75 to 0.68 for Tennessee limestone and from 0.55 to 0.40 for Know dolomite (Figure 4).We explain these data using the same procedure using σ 0 = 200 MPa.For the Tennessee limestone, we obtain μ 0 = 0.72 ± 0.7 and αq/p β = 0.056 ± 0.14.For Know dolomite, we get μ 0 = 0.52 ± 0.7 and αq/p β = 0.15 ± 0.14.The inferred parameters are highly negatively correlated due to the lack of data at low pressure.At high pressure, the constitutive parameters fall within the bounds that apply to the Byerlee (1978) dataset.
We now discuss the transient effects upon perturbations of normal stress (Figures 5 and 6).We first focus on the double-direct shear experiments conducted on Westerly granite at 5 MPa in nominally dry conditions, except for atmospheric humidity (Linker & Dieterich, 1992).The set of experiments illuminates the constitutive behavior of granite in initial bare contact during velocity steps and normal stress steps, allowing us to directly constrain many constitutive parameters independently.Considering that V 0 , d 0 , σ 0 , and G are merely scaling parameters, the remaining unknowns are the 8 constitutive parameters μ 0 , n, α, p, q, β, λ, and the loading stiffness k, which offer various sensitivity to mechanical data.The gouge thickness h is directly measured in the laboratory setting.First, velocity steps are conducted by alternating a loading rate of 1 and 0.1 μm/s at a constant normal stress of 5 MPa (Figure 5).The mechanical response includes a direct effect only mildly smoothed by the machine stiffness followed by a recovery to steady state compatible with near-neutral weakening.These experiments allow us to constrain the overall strength, the parameters controlling the direct and steady-state responses, as well as the weakening distance, and the loading stiffness, namely μ 0 , n, the ratio αn/p, and k.To explain these observations, we conduct numerical velocity step experiments with a spring-slider assembly governed by the constitutive framework devised in Section 2. We load the spring at the prescribed rate and compare the responses of the slider using the evolution law of Equation 5 and of Equation 6.The details of the calculation are described in Appendix D. To accommodate the nonlinear nature of the constitutive equations, we explore the best-fitting model parameters and associated uncertainties using a grid search.Laboratory measurements directly provide μ 0 = 0.67 and we use σ 0 = 5 MPa.Using the aging-law end-member of Equation 5, the model cannot reproduce the symmetric response to up and down velocity steps using constant constitutive parameters, as previously demonstrated (Ampuero & Rubin, 2008;Bhattacharya et al., 2015;Rathbun & Marone, 2013).In contrast, using the slip-law end-member of Equation 6 provides a satisfactory fit to the observations for up and down velocity steps at constant constitutive parameters.Specifically, the velocity step data require the stress power exponent n = 42 ± 5 and the micro-asperity size exponents α = 0.03 ± 0.005 and p = (0.95 ± 0.02) αn.The loading stiffness is estimated at 1.2 MPa/μm, close to the value of 0.5 MPa/μm estimated by Linker and Dieterich (1992).Assuming a narrow shear zone thickness of h = 0.5 mm, the weakening distance during the evolutionary phase is explained with the characteristic strain 1/λ = 0.4%.The constitutive framework captures the frictional behavior of rocks in isobaric conditions, as previously established (Barbot, 2019a).
Next, we consider the normal stress step experiments conducted at a uniform loading rate of 1 μm/s.Starting at a normal stress of 5 MPa, the interface is subject to normal stress increases of 1%, 2%, 5%, 10%, 20%, and 40%.Each up-step is conducted from the base level of 5 MPa, producing a sequence of increasing and decreasing normal stress.We model the steps in normal stress using the slip-law end-member (Figure 6) using the numerical procedure described in Appendix D. The mechanical data is reproduced with the normal stress exponent q = 10 and β = qα/p for the healing and direct effects, respectively.Ignoring the normal stress dependence of healing with q = 0 grossly misfits the data.Using the intermediate value q = 5 improves the quality of fit, but is suboptimal.The mechanical data upon steps in velocity or normal stress are explained for a single set of constitutive parameters shown in Table 1.The long-wavelength misfit is caused by the small hardening that accrues with cumulative slip.The shear stress increases slightly during the experiment due to slip hardening, but is assumed constant in the model, explaining the slight misfit in amplitude for the 40% normal stress step.However, there is a tradeoff between model parameters.The mechanical data can be explained equally well using β = 0.8 qα/p and μ 0 = 0.63 or with β = 1.2 qα/p and μ 0 = 0.71, all other parameters being the same.The constitutive model with the evolution law of Equation 6 captures the dynamic effects associated with steps in normal stress and loading velocity consistently at constant parameters for Westerly granite in bare contact.
Finally, we consider the normal stress step experiments conducted on synthetic gouge in double direct shear (Hong & Marone, 2005).The first experiments are conducted on 3 mm-thick synthetic quartz gouge based on Ottawa sand (>99% SiO 2 ) at 45% humidity (Hong & Marone, 2005, experiment number p287).A series of velocity step experiments alternating loading rates of 10 and 100 μm/s allow us to constrain some key physical parameters independently.The same sample is then subject to two sequences of varying normal stress of 20.5, 21.5, 23, 25, 28, 32, 37, 40, 37, 32, 28, 25, 23, and 20 MPa at a constant loading rate V L = 60 μm/s (Figure 7a).We first attempt to model the mechanical data with the aging-law end-member of Equation 5.The velocity step data show a symmetric response to up and down velocity steps that cannot be captured by the aging-law end-member.Furthermore, we find a strong trade-off between the constitutive parameters needed to fit the loading and unloading phases.When the weakening distance matches the increasing velocity steps, the model over-predicts the shear stress during the loading sequence and under-predicts the shear stress during the unloading sequence.
We then model the velocity steps and the ensuing normal stress steps using the evolution law of Equation 6.The velocity-step data can be reproduced assuming n = 90 ± 5, α = 0.03 ± 0.002, and p = (0.95 ± 0.05) αn, capturing the symmetry between up and down steps.The weakening distance during the evolutionary phase constrains the  (Linker & Dieterich, 1992) is shown in black and the numerical simulation based on the slip-law end-member of Equation 6 is shown in the red dashed line.The simulation corresponds to the constitutive parameters listed in Table 1.
characteristic strain 1/λ = 2.5%, taking the shear zone width h = 3 mm.Using the slip-law end-member allows us to reproduce the increasing and decreasing velocity steps using the same constitutive parameters.The model can explain the sequence of normal stress steps using μ 0 = 0.62, q = 10, and β = qα/p, including the smooth transition of shear stress during normal stress increase and the dynamic overshoot during decreasing normal stress steps.
Ignoring the normal stress dependence of healing with q = 0 produces a slight overshoot during the loading phase.
The constitutive parameters are summarized in Table 1.The model misfit during the two final stages of the unloading sequence is attributed to evolving frictional strength, as the shear stress is not identical under similar loads before and after the sequence.However, there are once again tradeoffs between model parameters and the same data can be explained with β = 0.8 qα/p using μ 0 = 0.583 or with β = 1.2 qα/p and μ 0 = 0.66, all other constitutive parameters being the same.
Lastly, we consider the experiments for a mixture of 50% quartz and 50% smectite gouge based on 25% humidity (Hong & Marone, 2005, experiment number p373).Using the aging-law end-member, we find a trade-off between the constitutive parameters needed to fit the loading and unloading sequences.Using the slip-law endmember, we can reproduce the observations consistently.The mechanical data can be explained by β = qα/p using μ 0 = 0.41, including the symmetric gradual rise and fall of shear stress upon increasing and decreasing normal stress steps (Figure 7b).The remaining parameters are listed in Table 1.Although most of the sequence is reproduced quantitatively well by the numerical simulations, the model does not perform well for the largest normal stress of 40 MPa corresponding to +100% of the reference level, presumably due to simplifying assumptions in the constitutive framework.A similar fit to the data can be obtained for β = 0.8 qα/p using μ 0 = 0.37 or with β = 1.2 qα/p using μ 0 = 0.45, all other parameters otherwise identical.
The constitutive model explains the variability of effective friction coefficients and apparent cohesion from atmospheric pressure to several gigapascals.The model also reproduces the transient effects for bare contact and synthetic gouge of various compositions.The normal stress step experiments are not sufficient on their own to indicate whether the steady-state shear strength is linear or nonlinear with normal stress.All of the sublinear, linear, and superlinear strength profiles implied by different values of β are compatible with the normal stress step data.Therefore, the nonlinear dependence of shear strength on normal stress highlighted by Byerlee (1978) is broadly compatible with the transient effects documented in various subsequent experiments (Hong & Marone, 2005;Linker & Dieterich, 1992, and references therein).Note.The constants V 0 = 1 μm/s, d 0 = 1 μm, σ 0 = 5 MPa, and G=(1 μm) p /s are scaling factors indicative of the relevant choice of physical units for the corresponding dynamic variables.The reference friction coefficient defined as μ 0 = χ s /χ n is a material property.

Comparison With Alternative Friction Laws
The constitutive model provides a physical basis for commonly used empirical friction laws in isobaric and non-isobaric conditions (Dieterich, 1979a;Linker & Dieterich, 1992;Ruina, 1983).A dimensional analysis of the evolution Equations 5 and 6 indicates that the size and age of microasperities relate to each other following where θ is the age of contact and d is the size of contact (Barbot, 2019a).Considering the change of variable of Equation 12 and a Taylor series expansion of the powers in terms of logarithms provide the friction law in additive form where a = μ/n, b = μα/p, and c = μβ.A detailed derivation, including the definition of the characteristic weakening distance L = 2h/(λp), is provided in Appendix A. Equation 13 simplifies to the classical formulation of Ruina (1983) in isobaric condition assuming σ = σ 0 .Using the change of variable in Equation 12, the evolution law of Equation 5 becomes θ = ( which simplifies to the aging law defined by Ruina (1983) in isobaric condition.Similarly, Equation 6 can be expressed in terms of age of contact, leading to the expression which simplifies to the slip law of Ruina (1983) in isobaric conditions.Equations 13-15, previously proposed by Sleep (2005) in the special case of q = c/b, highlight the direct effect of normal stress on the friction coefficient in Equation 13 and the effect of normal stress during the evolutionary phase in Equations 14 and 15.In contrast, the formulation of Linker and Dieterich (1992) uses the friction law coupled to either the modified aging law or the modified slip law The formulations of Equations 13-15 and 16-18 are strictly equivalent for q = c/b based on a change of variable, as derived in Appendix B. However, the direct effect of normal stress change is not obvious in Equations 16-18, where the direct effect is represented by an instantaneous change in the state variable.Furthermore, Equations 17 and 18 are ill-posed for step changes in normal stress, for which the rate of change is unbounded.When using the evolution laws of Equation 17 or 18, another equation is needed to describe the direct effect (e.g., Equation 11 of Linker and Dieterich (1992)).It is more consistent to capture all the direct effects in the friction law and the timedependent effects separately in the evolution law.As Equations 13-15 follow this strategy, they involve fewer dynamic variables and readily accommodate instantaneous changes in normal stress.
Regardless of notational differences, the formulations of Equations 13 and 16 are ill-defined, as they give rise to unbounded negative friction coefficients for vanishing velocity or vanishing normal stress.A form of regularization of these equations with a hyperbolic sine function has been proposed (Baumberger et al., 1999;Briscoe & Evans, 1982;Rice & Ben-Zion, 1996;Rice et al., 2001).In contrast, the multiplicative form of Equations 4-6 is valid for the entire range of slip-rate and normal stress, without the need for further regularization.In addition, the constitutive framework of Equations 4-6 generalizes the previous formulations of Linker and Dieterich (1992) and Sleep (2005) by allowing a normal stress dependence of the apparent friction coefficient, consistent with laboratory observations (Barton, 1973;Byerlee, 1978;Dieterich & Kilgore, 1996;Handin, 1969;Hoskins et al., 1968).
The proposed constitutive model provides a tentative explanation for the nonlinearity of shear strength with normal stress and the transient effects on rock friction observed upon perturbations of normal stress that may be relevant to the above-mentioned phenomena.The model seems broadly applicable to rocks in bare contact and to gouge layers, which provide more realistic proxies for natural fault zones.The nonlinear effects of normal stress originate fundamentally from the real area of contact of asperities that populate the fault interface.The distribution of contact junctions and fault roughness dictate the sensitivity of fault strength to normal stress.The variability of textural layouts at the micro-scale may explain the wide range of effective friction coefficients within the bounds of normal stress relevant to the brittle crust.The rate-independent failure criterion of Equation 11 captures the evolution of fault strength for a wide range of normal stress at constant constitutive parameters (Figure 8a).The nonlinear static yield criterion can be approximated by the linear terms of its Taylor series expansion, leading to an expression such as Equation 1, as shown in Appendix C.However, the resulting rock failure criterion is only valid for a limited range of normal stress in the neighborhood of the point of expansion.The effective friction coefficient is only apparent, varying greatly depending on the location of the expansion.The apparent friction coefficient is unbounded at vanishing normal stress for pβ > αq, allowing values as high as μ0 = 8 for small, yet finite values, as observed experimentally (Byerlee, 1978, Figure 3).The cohesion term c0 is merely an artifact of the Taylor series expansion.In reality, there is no residual strength at diminishing normal stress.Hence, the cohesion term in Equation 1 does not represent physical bonding or interlocking of asperities, as would apply for the failure of intact rocks in tension.
The nonlinear strength versus normal stress relationship applies at steady-state and during transient effects, particularly upon perturbations of normal stress.To illustrate the constitutive behavior during normal stress steps, we compute the frictional response of a spring-slider in the limit of infinite stiffness (Figure 8b).As the normal stress is suddenly increased, the apparent friction coefficient decreases as a direct effect, followed by a transient recovery to a new steady-state, which may be higher, identical, or lower than the previous value for pβ < αq, pβ = αq, and pβ > αq, respectively.The proposed constitutive framework for variable normal stress is compatible with previous formulations for the evolution of friction in non-isothermal conditions (Barbot, 2019a(Barbot, , 2022(Barbot, , 2023)).Hence, the model can be readily augmented to represent the slip-rate, state, temperature, and normal-stress dependence of friction, capturing the frictional response under an increasingly wide range of physical conditions relevant to faulting in the brittle crust.

Conclusions
We propose a constitutive framework that explains the hitherto unknown connection between the nonlinear effects of normal stress on transient and steady-state friction in isothermal conditions and nominally dry conditions, except for atmospheric humidity.Fault strength depends nonlinearly on normal stress, leading to effective friction coefficients that decrease, remain constant, or increase with normal stress, depending on constitutive properties.This effect results from the competing effects of the enhanced healing with increasing normal stress that leads to strengthening of the interface, and the direct effect of normal stress that decreases the effective friction coefficient.The relative strengths of these effects lead to superlinear, linear, and sublinear strength profiles depending on interfacial properties.In addition, the normal-stress dependence of healing causes an evolutionary response upon normal stress perturbations.
The constitutive framework generalizes previous formulations (Linker & Dieterich, 1992;Sleep, 2005), connecting nonlinear transient and steady-state effects.The transient phase is best captured by an evolution law that augments the slip law end-member of evolution laws.The multiplicative form is well-posed for any slip-rate and normal stress and does not require further regularization.The friction law is based on the area of contact and the evolving size of contact junctions, compatible with previous formulations that capture the competition between healing and deformation mechanisms (Barbot, 2019a(Barbot, , 2022(Barbot, , 2023)).As a result, the constitutive framework captures the slip-rate, state, temperature, and normal-stress dependence of fault friction for a wide range of conditions amenable to seismic cycles in the brittle crust.The physical model applies to many interconnected problems of natural and induced seismicity, including aftershock production, static and dynamic triggering of earthquakes, and rupture propagation on non-planar faults.A remaining unknown is the role of reactive fluids that affect fault strength beyond the Terzaghi effect.
where we used and we defined the characteristic weakening distance Using a Taylor series expansion of Equation A2, as in x y = exp( y ln x) = 1 + y ln x + O ln 2 x) , we can write the additive form of the friction law which is Equation 13 in the main text.Given Equation 2, for the real area of contact to increase with normal stress requires β < 1, implying c < μ 0 .When β = 1 or, equivalently, c = μ 0 , there is no direct strength increase upon step increase in normal stress, only a transient effect associated with the normal-stress dependence of healing.
We now discuss different formulations of the evolution laws.The correspondence between size and age of contact in Equation A1 implies Hence, the evolution law of Equation 5 can be written as a function of the age of contact which is Equation 14 in the main text.Equation A9 simplifies to the aging law defined by Ruina (1983) in isobaric condition for σ = σ 0 .Using the same change of variable, Equation 6 can be written which is Equation 15 in the main text.Equation A10 simplifies to the slip law defined by Ruina (1983) in isobaric condition for σ = σ 0 .The combination of the friction law of Equation A7 and the evolution laws in the form of Equations 14 and A10 were derived by Sleep (2005) by modifying the constitutive framework derived by where V L is the imposed loading rate.The constitutive framework assumes the friction law of Equation 4 and the evolution law of Equation 5 or Equation 6. Together with Equation D1, we obtain a nonlinear system of coupled ordinary differential equations.We solve for the evolution of the dynamic variables numerically using a Runge-Kutta method (Press et al., 1992).We initiate the system at steady-state whereby the state variable satisfies Equation 7. We simulate velocity steps by imposing a new loading velocity.For a sudden change of normal stress from σ1 to σ2 at constant stress, the direct effect on sliding velocity is where V 1 and V 2 are the slider velocity before and immediately after the normal-stress step.Simulations of velocity and normal-stress steps are shown in Figures 5-7.
We also conduct numerical simulations of loading and unloading sequences at infinite stiffness to illuminate the rheological behavior (Figure 8).In this case, the calculation is simplified because the sliding velocity is prescribed, that is, V = V L .At infinite stiffness, for a sudden change of normal stress from σ1 to σ2 at constant velocity, the direct effect imposes the instantaneous change in shear stress where τ 1 and τ 2 are the shear stress before and immediately after the step, respectively.

Figure 1 .
Figure 1.Normal stress dependence of the yield stress and effective friction coefficient of rocks.(a) Yield strength of rocks (circles) for effective normal stress up 2 GPa(Byerlee, 1978).Coefficients of friction of 0.7 and 0.6 are shown for reference.(b) Yield strength of rock for effective normal stress up to 400 MPa.There is no cohesion at vanishing normal stress.Coefficients of friction of 1.0, 0.7, and 0.5 are shown for reference.(c) Effective friction coefficient μ = τ/ σ for effective normal stress up to 2 GPa.A wide range of values between 0.25 and 1.5 can be found at low confining pressure.(d) Effective coefficient of friction for effective normal stress up to 400 MPa.Values below 0.5 and above 1.0 can be found below 200 MPa.In all panels, the superlinear (long dashed red line), linear (solid red line), and sublinear (short dashed red line) yield strength profiles correspond to steady-state friction with pβ > qα, pβ = qα, and pβ < qα in Equation10, respectively.

Figure 2 .
Figure 2. Area of contact density for acrylic (white triangles), calcite (white squares), glass (black circles), and quartz (black triangles) for varying normal stress and surface roughness and best-fitting model following Equation 2. (a) Linear plot.(b) Log-log plot emphasizing the power-law relationship.The model assumes an inverse relationship between grit number and representative micro-asperity size.

Figure 3 .
Figure 3. Normal stress dependence of the yield strength and effective friction coefficient of single rocks types from atmospheric pressure to 60 MPa.(a) Yield strength of Bilsthorpe mudstone (black circles) and Hucknall shale (red squares)(Barton, 1973) and corresponding model using Equation11(solid black and dashed red lines, respectively).(b) Yield strength of Bilsthorpe silty mudstone (black circles) and Ormonde siltstone (red squares)(Barton, 1973) and corresponding model.(c) Effective coefficient of friction of Bilsthorpe mudstone (black circles) and Hucknall shale (red squares)(Barton, 1973) and corresponding model using Equation11(solid black and dashed red lines, respectively).(d) Same for Bilsthorpe silty mudstone and Ormonde siltstone.

Figure 4 .
Figure 4. Normal stress dependence of the yield strength and effective friction coefficient of single rock types.(a) Yield strength of Wombeyan marble, trachyte, and Gosford sandstone up to 7 MPa (Hoskins et al., 1968) and corresponding model using Equation 11.(b) Effective coefficient of friction of Wombeyan marble, trachyte, and sandstone and corresponding model.(c) Yield strength of Tennessee limestone (red circles and error bars) and Know dolomite (black squares and error bars) and corresponding models up to 700 MPa (Handin, 1969).(d) Corresponding effective friction coefficient of Tennessee limestone and Know dolomite.

Figure 5 .
Figure5.Velocity-step experiments on Westerly granite with double-direct shear at 5 MPa normal stress(Linker &  Dieterich, 1992).(a) Comparison between the measured (black line) and simulated (red dashed line) shear stress using the aging law end-member of Equation5.(b) Comparison between measured and simulated shear stress using the slip law endmember of Equation6.The modified slip law better reproduces the symmetric response upon increasing and decreasing velocity steps.

Figure 6 .
Figure 6.Normal stress step experiments on Westerly granite with double-direct shear at a constant loading rate V L = 1 μm/s.(a) Increasing and decreasing steps of normal stress of 1%, 2%, 5%, and 10%.(b) Normal stress step of 20%.(c) Normal stress step of 40%.The mechanical data(Linker & Dieterich, 1992)  is shown in black and the numerical simulation based on the slip-law end-member of Equation 6 is shown in the red dashed line.The simulation corresponds to the constitutive parameters listed in Table1.

Figure 7 .
Figure 7. Normal stress step experiments on synthetic gouge with double-direct shear.(a) Normal stress loading and unloading sequence for a 3 mm-thick synthetic quartz gouge based on Ottawa sand (>99% SiO 2 ) at 45% humidity (Hong & Marone, 2005, experiment p287).Results from velocity step experiments alternating loading rates of 10 and 100 μm/s conducted previously in the same sequence are shown with the corresponding simulation.(b) Normal stress loading and unloading sequence for a 3 mm-thick mixture of 50% quartz and 50% smectite gouge based on 25% humidity (Hong & Marone, 2005, experiment p373) at a constant loading rate V L = 60 μm/s.The normal stress sequence is shown on top.The resulting measured (black line) and simulated (red dashed line) shear stress is shown in the bottom plot.