Reliable determination of rate- and-state friction (RSF) parameters depends on achieving steady-state (SS) friction conditions before and after experimental velocity-stepping friction tests. This operation, through nonlinear least squares fitting, is commonly preceded by the removal of any overall slip weakening/hardening after friction velocity steps (VSs) through a sufficiently large window of slip displacement at SS ( = linear detrend). However, to date, the identification of SS and thus the correct linear detrend is dependent on the user, which potentially results in differing RSF outputs from the same data set. Here, we demonstrate that large errors in the determination of the fitted RSF parameters can result if SS conditions are not reached before and after VSs. Such errors can be particularly relevant for materials characterized by long evolution of frictional resistance with slip, such as clay-rich gouge layers, in which identifying SS after VSs is not always obvious. To this end, we propose a methodology to accurately and consistently identify where SS is achieved after VSs. This methodology is coded into a new MATLAB-based routine, steadystate. We show the key features of the methodology, as well as how to use steadystate and read its output. We also illustrate the broad applicability of the approach to friction data with different noise levels and sampling frequencies referenced to slip velocity, by reviewing observations from synthetic data sets and specific examples of experiments from different laboratories involving various sheared materials.

The rate-and-state friction (RSF) formalism pioneered by J.H. Dieterich, A.L. Ruina, and J.R. Rice (Dieterich, 1979; Rice and Ruina, 1983; Ruina, 1983) has been widely used to understand the processes underlying the spectrum of fault-slip modes exhibited by tectonic faults, i.e., from aseismic to seismic slip, including slow slip behavior (e.g., Ikari, 2019; Leeman et al., 2016, among others). More recently, RSF laws have been applied to a growing number of geomechanical settings and conditions, including catastrophic failure of landslides (e.g., Agliardi et al., 2020; Handwerger et al., 2016; Helmstetter et al., 2004) and glaciers (Hudson et al., 2023; McCarthy et al., 2017; Zoet et al., 2013), as well as a broad range of non-geological materials (Shroff et al., 2014).

To determine the RSF parameters, the most widely used approach consists of running laboratory friction velocity-stepping tests, i.e., friction experiments in which time- and slip-dependent changes in macroscale frictional resistance are monitored through step changes in load-point velocity (Marone, 1998; Scholz, 2019). During the velocity stepping tests, initial steady-state (SS) friction conditions, i.e., constant friction resistance over time and slip preceding the friction velocity steps (VSs), are perturbed, which leads to a step change in frictional resistance, followed by a transient evolution in frictional resistance and consequent attainment of new SS friction conditions at the new slip velocity (Marone, 1998). RSF-dependent constitutive equations specify how frictional resistance varies as a nonlinear function of slip velocity and a set of state variables that describe the material’s slip history (see the next section for the functional form of these descriptions). In this framework, RSF parameter values are determined by identifying SS friction conditions before and after VSs, and nonlinear least-squares fitting of the data (Blanpied et al., 1998; Reinen and Weeks, 1993). Before VSs, SS is implicitly assumed (e.g., Giorgetti et al., 2015; Ruggieri et al., 2021), and after VSs, identifying SS can be made more difficult by frictional transients with slip.

A further challenge arises from the superposition of an overall slip dependency in friction (weakening or hardening), which is commonly observed in experimental data sets (e.g., Ito and Ikari, 2015; Kilgore et al., 1993). This type of trend needs to be removed prior to running the inversion scheme. To do this, the most common approach is to select a portion of the friction curve assumed to be at SS conditions, and then use its associated slope to eliminate any average slip-dependent linear trend in friction. This process is known as “linear detrending.” Some nonlinear detrending has also been employed (e.g., quadratic detrending; Noda and Shimamoto, 2009), although when using these methods, it is difficult to be certain that SS conditions are achieved both before and after VSs, as changes in trend of the friction-versus-displacement curves are typically associated with microstructural changes within the sheared material (e.g., Frank, 1965; Marone et al., 1990). In this study, we only focus on linear detrending. In the vast majority of software programs for determining RSF parameters, points to fit a detrending line are chosen after VSs.

Current practices often rely on the user “eyeballing” when new SS conditions are achieved to perform the detrend procedure or assuming a “standard” value of displacement following the VSs after which SS is assumed. The determination of RSF parameters in such methods is sensitive to how users determine (1) where new SS is reached and (2) the length of the slip window after this point so a linear detrend can be performed to average out experimental noise. This is especially the case in materials characterized by long frictional evolution to the new SS (e.g., Ito and Ikari, 2015). Given such considerations, differing RSF parameters can potentially result from the inversion scheme for the same experimental data set. For example, Skarbek and Savage (2019) pointed out that even small variations in the detrending slope can alter the estimated RSF parameter values. Appropriate identification of SS friction and the correct use of any linear detrend are thus key for reliable determination of RSF parameters and the application of models that depend on these data to characterize earthquake nucleation and fault slip (Dieterich, 1992; Gu et al., 1984; Rubin, 2008). This holds true especially for DRS, which can vary widely upon an incorrect SS assumption and motivates the need to establish consistent and reliable detrending methods for identifying new SS friction conditions and fitting experimental data.

To address the issues identified above, we present a new methodology for determining steady-state conditions that is written as a MATLAB-based routine called steadystate. This routine automates the comparison of the slopes estimated by linear regression before and after VSs (termed S1 and S2), with the latter calculated at progressively larger displacements from VSs until new SS conditions are identified, which is defined as when S2 ≈ S1. In particular, the method proposed in this study addresses the issue of finding new SS in VSs characterized by (1) a wide range of noise levels and number of points per unit slip and (2) different average values of slip dependencies in friction superimposed on VSs. To ensure consistency in new SS determination independently of (2), in our routine, estimation of S2 at increasing displacement is preceded by the removal of any average linear trend relative to S1 (provided that S1 ≠ 0), i.e., before VSs. This contrasts from the approach followed in the software currently employed for RSF determination, in which friction data are detrended linearly after VSs.

Previous attempts by the authors to identify SS conditions involved the systematic calculation of the weighted average second derivative of friction with respect to slip, using a tapered, moving slip window, with convergence occurring when the second derivative ≈0 within a tolerance. However, this approach proved less effective in determining the first new SS point after VSs than the method employed in this study, especially in VSs characterized by long frictional evolution to the new SS.

In the following sections, we give the background details on the RSF framework (section 2), then describe the RSF analysis of synthetic velocity steps including the methodology employed in steadystate (section 3). This is followed by applications of the SS analysis to both synthetic and experimental VSs (section 4), and the workflow illustrating the main architecture of the code (section 5). For the complete list of symbols and abbreviations used, refer to Table 1. steadystate is available on Github (https://github.com/pgiacomel/steadystate).

In laboratory friction data, the frictional resistance (µ) is determined as the ratio between the macroscopic shear resistance (τ) at the sliding interface or within the volume of powdered material, and the effective normal stress acting on the sample (σ'n), which is commonly quantified using the effective stress law σ'n = σnpf (Terzaghi, 1925), where σn and Pf are the applied normal stress and pore fluid pressure, respectively. Empirical RSF laws postulate that the change in frictional resistance (µ) to a sudden change in load point velocity can be described as a function of slip velocity (V) and a set of state variables (θi), whose number (i) is dependent on the friction data being fit.

Equation 1 represents the two-state variable form of the Dieterich-Ruina constitutive equation (Dieterich, 1979; Ruina, 1983):

formula

where a, b1, and b2 are dimensionless, empirically derived constants (with b = b1 + b2), and V0 and V are the load point velocities before and after the VS, respectively. We prefer the use of DRS to describe the characteristic slip distance rather than Dc as often denoted in RSF studies to avoid any ambiguity with the critical slip distances used for linear slip-weakening friction laws (Dal Zilio et al., 2022; Erickson et al., 2023; Lambert et al., 2021). The first term on the right side of Equation 1, µ0, is the reference SS friction coefficient determined at V0; the second term represents the direct effect, i.e., the direct friction response to the step change in load point velocity until the achievement of a peak value, which is scaled by a. In an infinitely stiff apparatus, the local maximum (or minimum) in friction is reached immediately following the VS. The third and fourth terms represent two separate contributions to the evolution effect, which describes the evolution of the friction coefficient with slip to a new SS following a VS. The evolution terms in Equation 1 are scaled by bi, and evolve with characteristic slip distances, DRSi. Because the characteristic slip distance relates to the slip necessary to effect a change from one SS friction value to another after a slip-velocity perturbation, its magnitude clearly affects the cumulative shear displacement at which new SS conditions are reached. In the laboratory, materials typically display DRS of ~10 µm, but in a few cases, such as in clay-rich materials, DRS can reach values of ~100 µm or greater (e.g., Ikari, 2019; Ikari and Hüpers, 2021; Ito and Ikari, 2015; Niemeijer et al., 2010; Pozzi et al., 2023; Sawai et al., 2016). DRSi scales the evolution of the state variables, θi, which reflect slip-history effects with units of time, and are often conceptualized in terms of the average lifetime of load-bearing contacts during frictional sliding (Dieterich, 1979, 1981; Dieterich and Conrad, 1984) or as contact strength (Goldsby et al., 2004; Nakatani, 2001; Thom et al., 2023).

The two most widely employed mathematical descriptions of the time- (and slip-) dependent evolution of the state variables upon a step change in load point velocity are the aging law (Equation 2A; Dieterich, 1979) and the slip law (Equation 2B; Ruina, 1983):

formula

formula

Note that at SS friction conditions, dθi/dt = 0, and consequently, both state evolution laws yield θssi = DRSi/V. Substituting this into Equation 1 quantifies the velocity dependence of friction under SS sliding conditions:

formula

where Δµss is the SS change in friction following a step change in slip velocity from V0 to V. Within this framework, positive values of a – b define the SS friction-velocity strengthening behavior associated with stable frictional sliding, whereas if a – b < 0, friction is SS velocity weakening, which is a requirement for the onset of unstable slip.

To evaluate whether materials exhibiting SS velocity weakening display dynamic slip, a single degree of freedom spring-slider system is introduced (e.g., Cook, 1981; Gu et al., 1984; Ruina, 1983). Although it is a simplified representation, the 1-D spring-slider system reasonably describes in the laboratory the interplay between the frictional properties of the surface/gouge layer and the elastic stiffness of the surroundings in terms of time derivative of the shear resistance τ, such that:

formula

where vlp is the load point velocity, v is the frictional sliding velocity resolved at the surface interface/gouge layer, and k is the elastic loading stiffness of the system (measured in unit stress/unit length). Assuming constant effective normal stress, σ'n throughout the test, Equation 4 can be rewritten as the first derivative in time of the friction coefficient, with graphic (given in 1/unit of length).

Friction stability analysis (Gu et al., 1984) and the estimation of the minimum patch size for the nucleation of frictional instabilities (Rubin, 2008; Rubin and Ampuero, 2005) rely on the correct determination of the friction velocity dependence (a–b) and of DRS. To extract the constitutive parameters a, bi, and DRSi from friction data sets, the coupled Equations 1, 2A or 2B, and 4 are solved, and a nonlinear least squares routine is employed to find the optimum fit of the RSF parameters. Cases in which RSF fit converges with DRS1 ≈ DRS2 or b2 = 0 indicate that a single-state variable constitutive law is sufficient to adequately model the data (e.g., Marone and Cox, 1994).

A detailed explanation of the iterative least squares method employed to solve the inverse problem is described in the papers of Bhattacharya et al. (2015), Blanpied et al. (1998), and Reinen and Weeks (1993). Such inversion schemes are incorporated into the software programs that are commonly used to determine RSF parameters, which include, among others, (1) xlook (https://github.com/PennStateRockandSedimentMechanics/xlook); (2) versi ons developed in Jupiter Notebook/python’s environment based on the xlook tool, such as rawPy (https://github.com/marcoscuderi/rawPy); (3) the codes of Bhattacharya et al. (2015) using Fortran, which are available at https://zenodo.org/record/2631455; and (4) the latest release using a MATLAB-based GUI, RSFit3000, which can be found on Github (https://github.com/rmskarbek/RSFit3000) and in Skarbek and Savage (2019). As the stiffness of the system can vary by >50% during shearing in a friction test (Leeman et al., 2016; Scuderi et al., 2017), k is often used as a fitting parameter during the numerical inversions (Noda and Shimamoto, 2009).

The determination of the fitting parameters is preceded by the removal of any linear slip-dependence in friction at SS conditions following the velocity steps:

formula

in which µ′ and μ are the linear detrended and experimental (i.e., undetrended) friction values pertaining to the selected VSs, δref is the displacement reference value where the velocity jump occurs, and δ(t) is the slip accumulated with time during the velocity step test. m refers to the slope calculated via linear regression within a slip window of an arbitrary length, and represents the magnitude of the average linear friction slip dependence that is removed from the friction curve µ to obtain µ′.

To investigate the effects of incorrect assumptions of SS friction conditions on the determination of RSF parameters and to develop an unbiased methodology for determining new SS, synthetic VSs from V1 = 0.3 to V2 = 3 µm/s were generated using a single degree of freedom spring-slider system (Equation 4), via the RSF constitutive law described by two sets of state variables (Equation 1) coupled with the “slip” evolution law (Equation 2B). In this study, synthetic VSs mimic experimental data sampled at a frequency of f = 100 Hz and with a random Gaussian noise of standard deviation SD = 0.0005 added to each data point to simulate the electrical noise (e.g., Chartrand, 2011; Skarbek and Savage, 2019). Friction is at SS conditions before VSs, and no slip dependencies in friction have been superimposed. The slope associated with the friction data preceding VSs (S1) is ≈0 when using a slip window length of 200 µm for the linear regression (−1.38 × 10−8 μm−1). S1 ≈ 0 is a sine qua non for applying the procedures described for determining steady-state, as it removes one variable from the system. Note that S1 = 0 occurs only in unrealistic zero-noise scenarios.

VSs were produced for DRS2 = 0 – 50 – 100 – 200 – 500 µm, and run for step lengths of 3 mm in the DRS2 = 500 µm step test and 2 mm in the remainder of the tests. Such step lengths are greater than those usually documented in the geomechanical literature to ensure that SS at the new slip velocity is reached. To mimic the elastic loading stiffness for biaxial (e.g., Collettini et al., 2014; Leeman et al., 2016; Scuderi et al., 2017) and triaxial deformation apparatuses (e.g., Bedford et al., 2021; Faulkner et al., 2018), k′ = 0.01 and 0.001 µm−1 were used, respectively. The sets of parameter values used to generate synthetic VSs are included in Table 2. Their relative VSs are displayed in Figure S11 and stored in the externally hosted files.2

For all VSs, the first new SS candidate point was chosen at 200 µm after the VS to account for the effects that could occur in compliant deformation apparatuses following the velocity jump, such as upstep-induced friction oscillations (e.g., Gu et al., 1984) or a long direct effect. Then, a slip window with a constant length of 100 µm was used from the selected point onwards to linearly detrend the friction data (Equation 5) using the slope S2. All friction points following the slip window used for detrending were neglected. Following the detrend, the RSF fitted parameters were determined using one-state and two-state variable fits and compared with the intrinsic parameters used to generate VSs. The same operation was repeated but with a slip window that was shifted forward by 100 µm until the end of the data was reached. During each RSF fit, the detrending slope S2 was systematically compared with S1.

Figure 1 illustrates the VS with intrinsic DRS2 = 500 µm and k′ = 0.01 µm−1, and highlights two examples of SS analysis carried out at a small displacement (linear detrend between 400 µm and 500 µm: panel I, red dashed rectangle, Figs. 1A and 1B) and at a large displacement (linear detrend between 2300 µm and 2400 µm: panel II, green dashed rectangle, Figs. 1A and 1C) following the VS. The complete collection of RSF values obtained by detrending friction data at increasing displacements is reported in Figure 2 and stored in the externally hosted supplemental files (see footnote 2). In Figures 2A and 2B, fitted RSF parameters are normalized by the intrinsic RSF values, so that when normalized data approach unity, fitted and intrinsic RSF parameters essentially coincide. Data are plotted as a function of their associated detrending slopes |S2| (shown by the color bar in Fig. 2). The friction velocity dependence a–Σb in Figure 2C corresponds to a–b and a–b1–b2 when friction data are modeled using a one-state and two-state variable fit, respectively; analogously, modeled DRS in Figure 2D refers to DRS for the one-state and DRS2 for the two-state variable case. For a direct comparison with fitted data, the intrinsic RSF parameters relative to the friction velocity dependence and DRS are also shown (light blue lines in Figs. 2C and 2D).

Overall, results from Figures 1 and 2 elucidate significant differences in RSF parameters resulting from incorrectly choosing SS in different parts of the datafile, as well as assuming a one-state variable fit to data when a two-state variable fit is more appropriate. Non-SS conditions are associated with a discernible slope before the VS following the detrending operation, so that S1 becomes different than zero (Fig. 1D). This is a feature commonly seen in experimental data sets. Conversely, when |S2| ≈ S1 ≈ 0, SS conditions are reached, and modeled RSF values approach the intrinsic RSF parameters (Figs. 1E, 2A, and 2B).

Our data notably show that the earlier new SS is assumed within the non-SS stages of the frictional transients following a VS, the more the detrending operation can produce misleading interpretations of the friction velocity dependence, as in Figure 1B. Under such circumstances, apparent velocity strengthening can be observed even in inherently strong velocity-weakening materials, for both one- and two-state variable fits (Fig. 2C), with DRS2 significantly underestimated (Fig. 2D). This effect is exacerbated at higher intrinsic DRS2 (Figs. 3A, 3B, and S2; see footnote 1), which testifies that the correct determination of SS becomes a more prominent issue the longer the frictional transients with slip are.

It is interesting to note that when friction data are fitted using a two-state variable constitutive law, the modeled direct effect, a, and the parameters related to the first-state variable (b1 and DRS1) map onto the intrinsic RSF parameters used to generate the synthetic VSs regardless of |S2|, and thereby, of the displacement at which new SS is assumed (Fig. 2A). Conversely, the parameters related to the second-state variable (b2 and DRS2) are very sensitive to variations in the new SS choice and tend to be their true values when new SS is assumed at larger displacements, where |S2| ≈ S1 ≈ 0 (Figs. 2C and 2D). This is due to the progressive incorporation of the longer evolution effect relative to the second-state variable while fitting the friction event. Equivalent conclusions are also drawn for the VS with k = 0.001 µm–1 (Figs. S1F, S3, and S4; see footnote 1), which illustrates that the same approach can be followed for a wide range of machine stiffnesses employed for assessing the frictional stability in the laboratory.

Building on the above observations, the systematic determination of |S2| and its comparison with S1 along VSs can be employed in automated routines as a method to ensure the correct determination of a new SS condition, which is associated with |S2| approaching 0 during shearing in all VSs characterized by S1 ≈ 0. However, since this condition is met asymptotically (Figs. 2C and 2D), it is useful to define a threshold for |S2| below which new SS is assumed in the routine steadystate. Called TH10, this threshold is defined from the synthetic VS depicted in Figure 1 (i.e., with DRS2 = 500 um) as |S2| = 5 × 10−7 µm−1, and results in less than ±10% error in the recovered parameters compared with the intrinsic RSF values (shaded gray areas in Figs. 2 and 3A). Figures 2C and 2D notably highlight the (a–b1–b2) and DRS2 parameter values from the modeled data used to constrain TH10 (red-bordered triangle).

Modeling friction data obeying a two-state variable constitutive law with a single set of state variables may, in some cases, return erroneous RSF parameters even when new SS is correctly determined (Fig. 2B). This holds true especially for modeled DRS values (Fig. 2D) in steps with long evolution effects (Fig. 3B), although at SS, the friction velocity dependence (a–b) is still well described (Figs. 2B, 2C, and 3A). Therefore, our data suggest that, after determining the first point at new SS conditions and detrending the data, it is good practice to begin the inversion scheme with a two-state variable model and consider a single-state variable fit only as a follow-up stage.

When applying the approach in section 3, determination of SS relies on the accuracy of the slope estimation before and after the VSs via linear regression analyses (i.e., S1 and S2, Figs. 1 and S3). Slope calculations can be influenced by the combination of the level of noise of the friction data and the number of datapoints used for the linear regression.

The closest approximation to the noise characterizing electric signals is the random Gaussian noise, which is quantified by its standard deviation (SD). Several factors contribute to the average noise in friction experimental data, including the apparatus, the applied effective normal stress, and the environment (e.g., presence of electronic devices interfering with the electric signal from transducers). The number of points per unit slip, Nnorm, relates the sampling frequency, f, to the slip velocity, V, as Nnorm = f/V.

To investigate the effect of noise and number of datapoints used to determine SS, the methodology described in section 3 was automated through the MATLAB routine steadystate. To investigate how we could minimize the effect of noise and sampling frequency, we also varied the length of the window (LS2) in units of displacement from 50 µm to 500 µm with increments of 50 µm. A key difference between the manual and automated methodology is that in the latter, SS is chosen on the basis of |S2| falling below the threshold formerly defined (see section 3) as TH10 = 5 × 10−7 μm−1 rather than identifying δTH10, which is the displacement at which the error in the returned rate and state parameters falls below ±10%. Clearly, when real experimental data are used, the intrinsic rate and state parameters are unknown, and hence δTH10 cannot be determined, so this alternate approach is necessary in the automated routine steadystate.

steadystate was applied to the same synthetic VS as in Figure 1, but with a broader range of random Gaussian noise and sampling frequencies applied. Three levels of noise were analyzed with various SDs: (1) 0.00025 (low noise), (2) 0.0005 (intermediate noise), and (3) 0.001 (high noise). The sampling frequency was varied from 1 Hz to 1000 Hz, with f increasing log-linearly from one test to another. The wide spectrum of SD and Nnorm used in this study encompasses, to the best of our knowledge, the vast majority of the values that can be found in velocity-stepping experiments within laboratories worldwide.

S1 was determined over a 200 µm slip window before the velocity jump occurred, which ensures S1 ≈ 0 over the investigated range of SD and Nnorm. As with the analysis in the previous section, the first new SS candidate point is chosen at 200 µm slip after the VS, and the slip window of length LS2 is moved to higher displacements by 100 µm increments until the new SS condition is met. The first new SS points compiled, obtained using the criterion |S2| ≤ 5 × 10−7 μm−1, are plotted against the slip-window length used (LS2; Fig. 4), and collected in the externally hosted supplemental files (see footnote 2).

Also included in Figure 4 for comparison is the displacement at which new SS was chosen via the criterion δTH10 described in section 3 (which produces modeled RSF parameters within ±10% of the intrinsic values). The two analyses, based on the criteria |S2| ≤ 5 × 10−7 μm–1 and δTH10, produce comparable results in terms of first new SS points (red-edged triangle versus δTH10, Fig. 4B). The displacements associated with SS that ensure fitted RSF parameters within ± 20%, 40%, and 65% of the intrinsic values are also reported (i.e., δTH20, δTH40, and δTH65, black dashed lines in Fig. 4). The error lines can be a proxy for the goodness of SS determination as a function of data noise, number of points per unit slip, and S2 window length LS2, in VSs processed with steadystate.

Overall, our data in Figure 4 confirm that combinations of noisy data (large SD) and a small number of data points (resulting from small Nnorm and/or small LS2) produce poor linear regressions for estimating S2, and thus erratic new SS determinations. This can be ascribed to the presence of influential points introduced by the superimposed noise, which becomes more relevant the fewer the number of points. Consequently, new SS outputs characterized by significant dependencies in LS2 (e.g., purple squares at LS2 ≤ 200 μm in Figs. 4B and 4C), in addition to those falling below the region delimited by δTH20 (e.g., purple squares at LS2 < 200 μm in Fig. 4A), are likely affected by poor S2 determination, and are therefore unreliable.

A solution to these issues (at a given noise level, sampling rate, and slip velocity) would be to increase LS2, so that the points inducing such issues become less influential in the slope estimation, leading to a more accurate new SS determination (Hastie et al., 2009; James et al., 2013). If a sufficiently large LS2 is chosen, whose length depends on SD (e.g., comparison of green triangles in Figs. 4A4C) and Nnorm (e.g., Fig. 4A), the analysis generates approximately the same choice of displacement value related to new SS for progressively larger LS2 (Figs. 4 and S5B; see footnote 1).

For efficiency reasons, and because slip windows are commonly short in experimental data, the routine identifies the first new SS point with the smallest LS2 associated with the constant portion of the plot (i.e., “optimum” output, red-bordered dot, Figs. 5 and 6I; see section 5.4, phase 3, for details).

Examples of the application of steadystate to experimental VSs are reported in Figure 5. These case studies illustrate the broad applicability of the routine steadystate to friction data from different laboratories, characterized by different machine stiffnesses, noise levels, average friction slip dependencies, and sheared materials with various RSF properties (Figs. 5A and 5B). Materials characterized by long frictional transients with slip, hence large DRS (e.g., clay-rich sepiolite gouge from Sánchez-Roa et al., 2016, 2017; Figs. 5B, 5D, and S5A), require larger displacements to reach new SS conditions when compared to the step length of 500 μm that is typically used in velocity-stepping tests (Figs. 5F and S5B). Conversely, in materials with short frictional evolution with slip (e.g., unaltered basalt gouge from Giacomel et al., 2021; Figs. 5A and 5C), a VS length of 500 μm is typically sufficient to attain new SS (Fig. 5E).

In this section, we outline the workflow used in steadystate. The routine follows the work protocol described in section 3, and finds its basis for the new SS determination in the key findings of the previous sections.

5.1 Phase 0: Slicing Velocity Steps and Selecting Points

The first preliminary operation is to slice the velocity-stepping experiment into its single VSs (phase 0, Fig. 6). One way to do so is via slicing_velsteps, i.e., a standalone script containing step-by-step instructions to follow in the MATLAB Command Window, which allows the whole spectrum of mechanical datafiles from worldwide rock deformation laboratories to be read and sliced. In doing so, friction-versus-displacement and friction-versus-time curves relative to each sliced velocity-stepping friction test are automatically displayed, and their corresponding data, including the effective normal stress, are stored in .txt format. Once the VSs are properly sliced, the main routine steadystate is ready to be run. In the first stage of steadystate, for each VS the user needs to select, in consecutive order, the three main points required for the new SS assessment, namely: (1) the first point of the analysis before the VS where SS has been established; (2) the reference point, ref, at which the VS occurs; and (3) the end of the datafile after slicing (i.e., the last point where the slope S2 can be estimated). Figure 6B illustrates an example of data-point selection preceding the SS analysis.

5.2 Phase 1: Preliminary Removal of S1 Slope

The automatization of the protocol outlined in section 3 requires S1 to be quantified before S2 is estimated, so the average slope before the VSs can be removed to begin the SS assessment with S1 ≈ 0 as a starting condition (phase 1, Fig. 6). This operation allows determination of the new SS in a consistent manner, as the same threshold value for slope S2 (i.e., TH10, section 3; Fig. 2) would be used for all VSs. Since the condition S1 = 0 cannot be reached in experimental data and the criterion for new SS is |S2| ≈ S1 (Fig. 1C) with |S2| ≤ TH10 = 5 × 10−7 µm−1 when S1 → 0 (Fig. 2), the value at which S1 can be reasonably considered ≈0 has been set to the same value previously determined for S2, i.e., |S1| ≤ TH10. Therefore, in the event |S1| > 5 × 10−7 µm−1 (Fig. 6D), the average linear friction slip dependence before the VS is removed from the entire friction event (Equation 5; Figs. 6E and 6F). Note that this procedure is applied only internally in the routine, without modifying the output friction data that appear in MATLAB workspace after the analysis.

S1 is calculated from the linear regression of friction-versus-displacement data between points 1 and 2 (Fig. 6B). To deal with the issues of noise and/or small numbers of data points for S1 determination, we determine the average noise level and Nnorm from the experimental friction data between the selected points 1 and 2 (Fig. 6B) and compare these values with those obtained from a database of synthetic friction data with no overall slip dependencies in friction. This database is contained within steadystate to identify where linear regression may return an erratic S1 slope from the experiment, based on the combinations of noise level and Nnorm that in the database return gradients |S1| > 5 × 10−7 µm−1.

In the synthetic database, the average random noise is defined by the user in terms of random Gaussian noise with SD, while in experimental data sets the average noise levels are estimated in steadystate through the root mean square error (RMSE). RMSE is defined as the SD of the residuals, and is computed as:

formula

with yi and σ'n denoting the i-th observed and predicted response, that in a friction-versus-displacement data set represent the i-th friction data point and the corresponding i-th friction value from the regression line, respectively, with n, the number of data points. This parameter is slope-independent an

formula

in data sets lacking any overall slip dependencies like those contained in the above database, given that in such cases the predicted responses (graphic) approach the estimated mean (y). Satisfying the condition SD ≈ RMSE in our synthetic database thereby allows the comparison of the average noise values contained therein with those from the experimental data. Based on such analyses, an error ends the routine when the combination of the low number of data points and high noise level would probably produce inaccurate S1 estimates, and prompts the user to choose a wider slip window (LS1). In the event that increasing LS1 may include data that are not considered steady state, additional velocity-stepping tests may be necessary with longer shear run-in stages or longer step lengths from the second velocity step onward.

To further address the possible presence of outliers introduced by the experimental noise, we determined the studentized residuals from a linear regression model. The studentized residual method is employed in statistics to eliminate outliers that influence the linear regression model to such an extent that the estimated regression function is “pulled” toward the potential outlier so that it is not flagged as an outlier using the classic standardized residual criterion. Such outliers have the potential to skew the regression line and influence the value of S1. A studentized residual is calculated by dividing the residual by an estimate of its standard deviation. More details can be found in Hastie et al. (2009), James et al. (2013), and Pardoe (2021). In steadystate, we define “influential outliers” as all of the studentized residuals with values outside of the ±3 bounds and run a linear regression devoid of these points to improve the accuracy of S1 determination (Fig. 6C).

5.3 Phase 2: Estimation of S2 Slopes

Phase 2 consists of multiple SS analyses that begin at 200 µm past the velocity jump (Fig. 6F), where the slope S2 is estimated within an evaluation slip window of a given fixed length, LS2 (Figs. 6F and 6J), which spans from the minimum LS2 (min_LS2) to the maximum LS2 (max_LS2) with step increments, ΔLS2, from a new SS assessment to the next one (Figs. 6G and 6L).

In the framework of a single SS analysis, during each iteration, the absolute value of S2 is compared with TH10 = 5 × 10−7 µm−1 (Fig. 2). If |S2| > TH10, a new linear regression is run, adding on net slip increments, Δδ, to produce a new starting point to determine S2 (Figs. 6F and 6H). The iterations stop once |S2| ≤ TH10 is satisfied, or until the end of the datafile is reached if such a criterion is not met. The displacement at which iterations stop is recorded as the first new SS point associated with a given LS2 (Fig. 6K; panel I, left side, Fig. 6). When the last datapoint is reached, a warning is given; at the second warning of the same type, an error ends the routine with the suggestion that a repeat test with longer VSs may be required to address the issue.

To conduct phase 2, four input parameters (in microns) are required in steadystate in the following order: min_LS2, max_LS2, ΔLS2, and Δδ. In the first instance, the parameters should be left blank or null (i.e., []). In this case, the routine uses its default settings that for min_LS2 and ΔLS2 amount to 50 µm, and calculates Δδ, which is numerically equivalent to the slip distance between two datapoints from point 2 + 200 µm onwards. If the net slip increment, Δδ, is overridden by the user, a warning sign appears if the entered value is lower than the displacement between datapoints, resulting in the replacement of the former value with the default argument to continue the analysis. max_LS2 is a function of the total window length (TWL), which is defined as between point 2 + 200 µm and point 3 (Figs. 6B and 6F). When TWL ≤ 1000 µm, max_LS2 equals 50% of TWL being approximated to the closest multiple of 100, whereas if TWL > 1000 µm, max_LS2 is limited to 500 µm. This max_LS2 upper boundary is aimed at limiting the proportion of frictional transients in S2 calculations, which may result in new SS detection at smaller displacements than true SS.

When the three arguments associated with LS2 are inputted by the user, the routine first approximates ΔLS2 to the closest multiple of 10 and according to the input arguments, it rounds min_LS2 and max_LS2 to the closest multiple of the approximated ΔLS2. Moreover, in the MATLAB command window, steadystate suggests an optimal value for max_LS2 following the same criteria used for the default settings.

5.4 Phase 3: Collection and Checking of Outputs: Optimum New Steady State

At the end of the routine, the collection of displacement values corresponding to new SS, namely the first new SS points, is plotted against their corresponding S2 window lengths (LS2), where the optimum SS output is highlighted. The optimum choice of the first new SS point is the displacement value coupled with the smallest LS2, pertaining to the approximately constant portion of the plot (Fig. 6, left side, panel I). These values should then be used for linear detrending of the VS (Fig. 6, right side, panel I). The routine also provides the user with the opportunity to save outputs in Excel format.

Following on the output generation from the SS analysis, although unlikely, there is a slight possibility that the routine returns no optimum outputs from the SS analysis. This may occur as a result of highly scattered first new SS values as a function of LS2, in data sets characterized by significant noise and a low-number of data points, despite adopting a wide range of LS2. In this case, the routine quits with an associated error, as proceeding with the linear detrend and the inversion for the RSF parameters is not recommended. Such an error is followed by the suggestion to either attempt repeating the SS analysis with a larger max_LS2 (of the amount ΔLS2), or to run another velocity-stepping test with a larger number of points per unit slip.

We present the case for ensuring proper identification of steady-state conditions when processing laboratory friction data to estimate RSF parameters. We show from the processing of synthetic data sets that serious discrepancies between the modeled and intrinsic RSF parameters can result from the processing of data where steady-state is incorrectly identified. To address these issues, we introduce a methodology to identify steady-state conditions based on the achievement of similar slopes before and after a velocity step.

The methodology has been developed into a MATLAB-based routine called steadystate to produce unbiased, consistent, and reliable determination of new steady-state conditions under a broad range of experimental conditions and sheared materials. We hope that this will promote progress toward consistent parameter estimation techniques and reduce reliance on arbitrary user-dependent choices of steady-state during the processing of laboratory friction data that can lead to notable discrepancies in RSF parameter estimation from the same data set.

Our systematic analysis of synthetic data sets emphasizes that slope estimations and, consequently, new steady-state determinations, can be negatively affected by high-noise levels and a low-number of points per unit slip of friction data, given the same intrinsic RSF properties of the material (Fig. 4). We show how these problems can be circumvented by systematically evaluating a range of slip-window lengths (LS2) when assessing new steady-state conditions. The routine not only identifies the first displacement at which steady-state conditions are achieved, it provides the minimum evaluation window length (LS2) it used to determine this point. These two values can then be used to linearly detrend the friction data, and thus, to accurately determine the RSF parameters.

Following the identification of an accurate first steady-state point, coupled with a sufficiently large LS2 to produce a good linear detrend, it is good practice to first perform a two-state variable fit, and replace it with a one-state variable model only if RSF parameter values degenerate into a one-state variable fit (i.e., b2 ≈ 0 or DRS1DRS2). This protocol is motivated by the observation that a one-state variable fit fails to describe the friction evolution with slip if the velocity step obeys a two-state variable constitutive law (Figs. 2B and 2D).

This study confirms that materials characterized by long frictional transients with slip after velocity steps, such as clays (Fig. 5D), often require larger displacements to attain new steady-state conditions when compared to the step length of 500 µm that is typically used in velocity-stepping tests (Fig. 5F). Under such circumstances, modeled RSF parameter values would significantly depart from the true values if the intrinsic DRS2 were commensurate with the whole displacement window that contains the velocity step. Conversely, in materials with short frictional transients with slip following the velocity jumps (Fig. 5C), running velocity steps with 500 µm step length is normally sufficient to determine new steady-state (Fig. 5E). Furthermore, in these materials, since RSF estimates are less sensitive to the accuracy of steady-state determination (Fig. 3), operating the detrend “by eye,” as has been done routinely so far, would still likely produce good RSF fits.

These observations suggest that further investigation is warranted to better decipher the nature of the two-state variable, which is currently used to produce better RSF fits in long-evolution velocity steps, but whose underlying physical significance is still largely unknown.

Finally, to date, steadystate has been purposely designed to work in conjunction with an RSF modeling technique. However, its applicability in the future could be expanded to any situation involving a general slope estimation (e.g., calculating the elastic moduli, any linear slip dependence in friction, etc.).

1Supplemental Material. Figures S1–S5. Please visit https://doi.org/10.1130/GEOS.S.25456102 to access the supplemental material, and contact editing@geosociety.org with any questions.
2Supplemental Material (externally hosted). Dataset and MATLAB scripts to generate the figures related to this article. Please visit https://doi.org/10.5285/0a6b0b0c-b84e-44fd-a331-7c1b4ff9165e in the BGS National Geoscience Data Centre under the Open Government Licence (ID 184041).
Science Editor: Andrea Hampel

P. Giacomel and D.R. Faulkner gratefully acknowledge Natural Environment Research Council grant NE/V011804/1 for funding this work. Reviewers R. Skarbek and M. Ikari are warmly thanked for their insightful comments and suggestions that further improved this work. R. Skarbek is greatly thanked for providing the MATLAB spring-slider model for generation of synthetic velocity steps. M. Dalla Pria is gratefully acknowledged for his suggestions on inferential statistics, and M. Battaglia is thanked for his advice on coding. M. Ikari and A.M. Eijsink are warmly thanked for providing experimental data from MARUM Laboratory to test the code. This work also benefited from insightful discussions with E. Mariani, I. Ashman, and L. Xi.

Gold Open Access: This paper is published under the terms of the CC-BY license.