We present a MATLAB graphical user interface (GUI) software package for analyzing rate and state friction experiments. Called RSFit3000, the software allows users to easily determine frictional parameters by fitting velocity-step and slide-hold-slide events using the aging- and slip-law forms for state variable evolution. RSFit3000 includes features for removing strain hardening or weakening trends from the data, and provides options for using two state variables, applying a weighting function, and treating stiffness as a fitting parameter. Completed fits are conveniently saved in MATLAB structure arrays that contain fitted parameter values with their error intervals, and all information required to reproduce a given fit. The GUI makes the program simple to use, as all fitting tasks are completed via interaction with the interface. Here we describe how to use the software, and illustrate its flexibility and utility by analyzing two sets of synthetic friction data, as well as some previously published experimental data. Although descriptions of rate and state friction fitting routines have been published in the past, RSFit3000 marks the first time a software package for analyzing friction experiments has been described in the literature.
Rock friction experiments have provided the basis for the rate and state friction equations, which represent the current best description of sliding friction available to the Earth science community. Although rate and state friction theory is mostly used to describe various aspects of fault slip (e.g., Marone, 1998; Lapusta et al., 2000; Rubin, 2008), it has also had success in describing glacier sliding (e.g., Zoet et al., 2013; McCarthy et al., 2017), landslides (Handwerger et al., 2016), chemical bond–induced friction at nanoscale contacts (Li et al., 2011; Tian et al., 2017, 2018), and frictional slip for a wide range of materials (Heslot et al., 1994; Dieterich and Kilgore, 1996; Ullah et al., 2006). The equations describe the evolution of the friction coefficient as a nonlinear function of sliding velocity, and an internal state variable that can be thought of as representing the average lifetime of frictional contacts (Scholz, 1998).
The rate and state equations are empirical (although some physical interpretations have been attached to the constitutive frictional parameter; e.g., Rice et al., 2001), and due to the nonlinear nature of the equations, it is necessary to carry out a numerical inversion to determine values of the constitutive parameters for any given experiment. The details of numerical routines for carrying out such inversions have been published by Reinen and Weeks (1993), Blanpied et al. (1998), Noda et al. (2009), and Bhattacharya et al. (2015). At present, there are some inversion codes used by experimental rock friction community with varying levels of documentation. Among these, the best documented is xlook (https://github.com/PennStateRockandSedimentMechanics/xlook), many features of which are shared by RSFit3000. The inversion scheme described by Bhattacharya et al. (2015) can be found at https://zenodo.org/record/2631455.
Here we describe a MATLAB graphical user interface (GUI)–based program, called RSFit3000, that provides an easy-to-use method for determining rate and state frictional parameters from experimental data. The code is included in Supplemental File 11, and is also available on GitHub (https://github.com/rmskarbek/RSFit3000). RSFit3000 is capable of fitting velocity-step events as well as slide-hold-slide (SHS) events using both the aging- and slip-law formulations for state variable evolution. It can also perform fits using two state variables, and includes options for treating the stiffness as a fitting parameter and for producing weighted fits. A detailed user guide is included in Supplemental File 1. Here we describe the main features of the program and the workflow involved in obtaining fits to data. Finally, we suggest some best practices for fitting and reporting rate and state friction variables.
Other equations have also been proposed (see Marone, 1998; Nagata et al., 2012), but are not in common use for experimental data, and we note that recent experimental and theoretical work has shown that the slip law is capable of describing a larger range of experimental results than the aging law (Bhattacharya et al., 2015). RSFit3000 makes use of both Equations 2 and 3, with one or two state variables θi.
RSFit3000 is suitable for the analysis of any friction experiment where the normal stress remains constant. It does not take into account any geometrical or frictional effects that can occur as a result of changing normal stress, such as Linker-Dieterich effects (e.g., Linker and Dieterich, 1992; Prakash, 1998; Richardson and Marone, 1999; Hong and Marone, 2005; Kilgore et al., 2017), or experiments conducted on sliding surfaces that are not orthogonal to the loading axes (e.g., He et al., 1998).
Experimental measurements of rock friction are commonly made using two types of imposed friction events: velocity steps, and SHS events. In a velocity step, the load point is set in motion at a constant rate vinit for a sufficient displacement that a steady friction coefficient is achieved. The sliding rate of the load point is then changed as quickly as possible to a new value vf. The change in sliding rate induces a “direct effect” and an “evolution effect” in the measured value of the friction coefficient (Fig. 1B). We describe these effects for an infinite-stiffness system undergoing a velocity increase (vf > vinit); the effects are reversed if vf < vinit. Immediately after the change in sliding rate, the friction coefficient increases to a local maximum, changing by an amount aln(vf /vinit) (the change is smaller than this amount for finite stiffness); this is the direct effect. After reaching a peak value, the friction coefficient decays to a new steady-state value as sliding continues, changing by an amount −bln(vf /vinit) from the peak value; this is the evolution effect. If , the steady-state friction value at vf will be less than that at vinit, and vice versa if .
A SHS event also initiates at a steady-state value of friction at a load-point velocity vload. The load point is stopped altogether for some amount of time (the hold), and then started again at a new value vreload, usually, but not necessarily, equal to the previous sliding rate. During the hold period, the friction coefficient decays logarithmically with time. When sliding is reinitiated, the friction coefficient increases rapidly to a peak value, and then decays to a new steady state (Fig. 1C). SHS events can also be considered as two velocity steps: a decrease from vload to zero, and an increase from zero to vreload (e.g., Blanpied et al., 1998; C. Marone, 2019, personal commun.).
Due to the nonlinear nature of Equations 1–3, the frictional parameters a, bi, and dci must be determined by conducting forward-modeling fits to experimental friction events. Additionally, because a simple spring-slider model, Equation 4, is used to represent the combined elastic response of the testing machine (which depends on its geometry and materials) and the sample assembly (which likely evolves throughout a given experiment), the stiffness k should also be determined through the fitting procedure as well. Recent experimental studies on gouge materials have shown that the effects of compaction, comminution, and shear localization can cause the stiffness to vary by >50% as a function of shear strain (Leeman et al., 2015, 2016; Scuderi et al., 2016). Alternatively, one may attempt to directly measure the stiffness, for example, by performing a long SHS event periodically throughout an experiment and taking k as the slope of a linear fit to friction versus load-point displacement during the initial part of reloading (e.g., Bhattacharya et al., 2015). In either case, it is important to note that treating the stiffness as a fitting parameter would generate fits with different frictional parameter values than would otherwise occur.
The value of μ0 can also be used as a fitting parameter, but is generally well constrained through direct measurement. Reinen and Weeks (1993), Blanpied et al. (1998), and Noda et al. (2009) all provided descriptions of techniques for fitting Equations 1–4 to experimental data. In developing RSFit3000, we made use of a combination of techniques from each of these papers, although our methods are perhaps most similar to those of Noda et al. (2009).
RSFit3000 is built around a nonlinear least-squares fitting routine that uses the Levenberg-Marquardt method (e.g., Noda et al., 2009). The routine makes use of the built-in MATLAB function lsqnonlin, which is part of the MATLAB Optimization Toolbox. The program will not function properly on MATLAB installations that do not have this toolbox. We give a brief description of the algorithm. Detailed descriptions of the Levenberg-Marquardt method can be found in Levenberg (1944), Marquardt (1963), and Moré (1978). In cases where more complex state variable evolution laws are used (e.g., Nagata et al., 2012), or second-order terms are included in Equation 1 (e.g., Baumberger et al., 1999; Bhattacharya et al., 2017), it may be more appropriate to use methods other than Levenberg-Marquardt (Bhattacharya et al., 2015).
When the fitting routine is started, it computes μpred(ti; ξ0) by numerically simulating the experimental conditions using Equations 1, 2 or 3, and 4. Simulations are carried out using the built-in MATLAB function ode45, which is an explicit Runge-Kutta–based solver with an adaptive time step (Shampine and Reichelt, 1997; Ashino et al., 2000). The simulation is rerun every time the parameter values are updated.
The fitting routine updates the parameter values vector until one of several stopping criteria is met. Figure 2 schematically illustrates how the stopping criteria are evaluated. The routine will stop successfully if any of three criteria are met: (1) a local minimum in the objective function is found, (2) the relative sum of squares of the objective function (vertical distance between points A and B in Fig. 2) changes by less than a defined tolerance, or (3) the norm of the step size (horizontal distance between points A and B) changes by less than a defined tolerance. The routine stops unsuccessfully if the number of function evaluations or solver iterations exceeds defined limits, or if the current parameter values produce unstable, stick-slip events. The tolerance values and function evaluation and iteration limits can all be changed by the user (see the user guide [footnote 1] for details), but we have found that using the default values works very well. The stopping criteria are part of lsqnonlin, and more information on tolerances and stopping criteria in general can be found in the MATLAB documentation. The routine generally only stops unsuccessfully when the current set of parameter values leads to unstable slip. When this occurs, the simulated sliding velocity diverges and causes the solver routine to crash. Also, in some cases, the routine may result in a poor fit to the data depending on the initial parameter values. In the next section, and in the user guide, we discuss techniques for avoiding these outcomes.
In this section, we summarize the steps needed to produce fits and describe some of the mathematical procedures that RSFit3000 uses. The following data are required to perform fits: friction coefficient, load-point displacement, normal stress, and time. Each of these data sets must be loaded into the MATLAB workspace as vectors with identical sizes, and entered into the Experimental Data panel (Fig. 3A). Once these data are input to the program, initially identical plots of friction coefficient against load-point displacement appear in the Static and Windowing axes (Figs. 3B and 3C). To continue, the zoom and pan tools are used to isolate an event of interest in the Windowing axes, such as a single velocity step or SHS event. A red box will appear on the Static axes showing the location of the windowed data within the entire data set.
Next, detrending is performed in the Windowing axes if necessary. In real experimental friction data, there is commonly a global strain weakening or hardening trend. In other words, for sliding at a steady load-point velocity, the friction coefficient decreases (weakening) or increases (hardening) at an approximately constant rate. The common practice is to remove the global trend before fitting the data with the rate and state equations (Lockner et al., 1986; Tullis and Weeks, 1986; Marone et al., 1990; Chester, 1994; Beeler et al., 1996; Blanpied et al., 1998). Detrending is performed in RSFit3000 by fitting a line to selected portions of the experimental data that display a weakening or hardening trend, and using the slope of the fitted line to remove the trend (Fig. 3D).
To detrend data, the user selects two points (δ1, μ1) and (δ2, μ2) on the windowed data, where δ is the load-point displacement, and the program fits a line to the data between these points. Then the user selects a reference point (δr, μr) on the data, and the program removes the trend defined by the fitted line relative to the reference point, according to the equation μ′(ti) = μ(ti) + m[δr − δ(ti)], where m is the slope of the line between points (δ1, μ1) and (δ2, μ2), and μ(ti) and δ(ti) denote the friction coefficient and load-point displacement data that are visible in the Windowing axes.
After detrending, the detrended and windowed data appear in the Fitting axes (Fig. 3E). The proper event type (velocity step or SHS) must be selected in the Event Type panel (Fig. 3F). If the event is a velocity step, the detrended data are plotted against load-point displacement; if it is a SHS event, the data are plotted against time. In either case, the x-axis data in the Fitting axes are referenced to zero at the beginning of the windowed event. There are four options available for performing a fit: (1) use of a weighting function, use of (2) μ0 or (3) stiffness k as fitting parameters, and (4) use of two state variables (Fig. 3G). Any of these options may be combined with one another. The user may also select whether to use both or one of the state evolution laws.
We implement a weight function similar to those described by Reinen and Weeks (1993) and Blanpied et al. (1998). When use of the weight function is desired, the user selects the appropriate option, then sets a value for the exponent p in the Weight Parameters panel (Fig. 3H), and clicks the Set Weight Location button. The weight location is set by selecting a data point in the Fitting axes. Once selected, the corresponding values of load-point displacement and friction coefficient appear in the Weight Parameters panel, and the program generates a plot of the resulting weight values.
The effects of using a weight function following Equation 9 have been discussed by Reinen and Weeks (1993), who used a weight function, and also by Noda et al. (2009), who did not use a weight function. Reinen and Weeks (1993) and later Blanpied et al. (1998)’s motivation for using a weight function is that the fitted values of the rate and state parameters depend mostly on the evolution of the friction values after the peak value is reached following a velocity up-step (i.e., on the evolution effect), whereas the rise in friction values to the peak itself (i.e., the direct effect) is highly influenced by the value of the stiffness k. As such, neither Reinen and Weeks (1993) nor Blanpied et al. (1998) treated the stiffness as a fitting parameter. Instead, they gave k a constant value, and argued that the friction values recorded after the peak value should be weighted more heavily than other portions of a velocity up-step event. Thus they used weight functions similar to Equation 9 with the weight location set at the peak friction value. Alternatively, Noda et al. (2009) did treat the stiffness as a fitting parameter, and pointed out that using a weighting function in such a case is likely to increase uncertainty in the fitted value of k, because the portion of the up-step event that depends most on the stiffness would be weighted less than other portions of the data. We leave it to the user to decide upon the use of a weight function, but note that its use can occasionally produce better fits, even when fitting the stiffness, as illustrated in the next section.
The procedures for fitting velocity steps and SHS events are very similar. For a velocity step, the user selects the location of the step in the data on the Fitting axes. The values of load-point displacement and friction coefficient from the selected point appear in the Velocity Step Parameters panel (Fig. 3I), and the program generates a plot of the load-point displacement against time from the windowed data, showing the location of the selected step. The program also calculates and displays the load-point velocity before and after the selected step location, using the load-point displacement time series before and after the selected step. All of these values can also be manually entered.
For a SHS event, the user also selects the location of the beginning of the hold on the Fitting axes. After selecting the hold location, the values of time and friction coefficient appear in the Slide-Hold-Slide Parameters panel (Fig. 3J), but the program does not automatically compute the remaining values in the panel. The user must manually enter values for the length of time that the hold lasted (Hold Duration) and the load-point velocities before (Load Velocity) and after (Reload Velocity) the hold. As an additional characterization of the event, when either a velocity-step or SHS location is selected, the normal stress at the time of the location appears below the Velocity Step Parameters panel (Fig. 3K).
To apply the fitting procedure (Fig. 3L), the user selects whether to perform fits using the aging and/or slip laws for state evolution. The user then enters trial values for μ0, a, b1, dc1, and k, and for b2 and dc2 if using two state variables. Values for μ0 and k must be entered regardless of whether the options to constrain these two values are selected. After entering the trial values, clicking the Test button runs a spring-slider simulation using the trial values, and the results appear in the Fitting axes. Depending on the event type, the program uses the event parameters to simulate the correct behavior. For a velocity-step event, the program determines the time that the step occurred when the user selects the step location. This information, along with the initial and final load-point velocities, is passed to the spring-slider program and used to simulate the proper behavior of the load point. A similar procedure is used to simulate a SHS event. Information about when the hold initiates, how long it lasts, and the load-point velocities before and after the hold are passed to the spring-slider program. If the test lines do not closely overlie the data, any of the trial values can be changed to run another test. This process should be repeated until the test lines overlie the data as much as possible. Obtaining trial values in this manner is critical to obtaining good fitted parameter values.
After generating satisfactory initial parameter values, the user clicks the Fit button, and the program carries out the nonlinear least-squares fitting routine. Information regarding the number of solution attempts and the parameter values for each attempt is displayed in the MATLAB command window as the fitting routine is running. When the routine has completed, the results appear in the Fitting axes, and the fitted parameter values are displayed in the Fitting Parameters panel. Each fitted parameter value is displayed with its 2s error interval. These values, along with all of the information used to produce the fit, can be saved as a structure in the MATLAB workspace by entering a name in the Structure Name field and clicking the Create Structure button (Fig. 3M).
The created structure contains the load-point displacement, measured friction coefficient, detrended friction coefficient, and time values from the windowed event. It also contains the parameters that were used to detrend the friction data, as well as the initial parameter values, the fitted parameter values with the error intervals, the R2 value and covariance matrix, and the predicted friction-coefficient values generated by the fitting procedure. If the event was a velocity step, the structure contains the initial and final load-point velocities and the time and the normal stress at which the step occurred. If the event was a SHS, the structure contains the load and reload velocities, the time and normal stress at which the hold initiated, and the duration of the hold. If the weight function was used, the relevant parameters are in the structure as well. The user can window another velocity step or SHS event in the Windowing axes using the zoom and pan tools, and repeat the procedures described above to obtain more fits. When moving to a new event, all of the parameters from the previous event remain in the GUI, and can be used to detrend and fit additional events (see the user guide [footnote 1] for more details).
In this section, we present some examples that illustrate the features of RSFit3000 and techniques that can be used to fit friction events. First, we present some results from a synthetic data set, which allows us to produce events that require different techniques to obtain good fits. In contrast to real experimental data, synthetic data are generated using a defined set of parameter values, as well as a specific choice of state evolution law. However, to mimic real data, we add a significant amount of noise to our synthetic data, and the goal of our fits is not to exactly recover the original parameter values. The addition of noise allows for a more realistic demonstration of the fitting routine. Indeed, fits to synthetic data that did not have added noise would always produce the original parameter values with extremely small error estimates, and be of little instructive value. Finally, we present some results from fits to experimental data previously published by Ikari et al. (2014), which illustrate the effects that detrending can have on fit results. All of the information needed to reproduce the fits shown in this paper can be found in Supplemental File 22.
We generated two synthetic experimental data sets by running a spring-slider simulation using the slip law. The synthetic data and the original parameter values are shown in Figure 4. The data points were generated at 100 Hz. The data shown in Figure 4A use a single state variable, and those shown in Figure 4B use two state variables. Both sets were generated using an identical load-point history, consisting of four events: (1) a velocity up-step from 1 μm/s to 10 μm/s; (2) a velocity down-step from 10 μm/s to 1 μm/s; (3) a SHS of 10 s duration with identical load and reload velocities of 1 μm; and (4) a SHS of 1000 s duration with a load velocity of 1 μm/s and a reload velocity of 10 μm/s. Random Gaussian noise was added to both data sets (e.g., Noda et al., 2009) using the MATLAB function randn in the equation μ = ν × randn(N, 1) + μs, where μs is the original synthetic value of the friction coefficient, N is the number of samples, and ν = 0.005 is the variance of the added noise signal.
The fits to the events from both data sets are shown in Figures 5 and 6, along with the fitted parameter values and their error intervals. The fits demonstrate the flexibility of the features included in RSFit3000. Frequently, it’s possible to fit an event with both state evolution laws using a single set of initial parameter values ξ0. Then the aging- and slip-law fits can be generated simultaneously by selecting both laws in the Fitting Parameters window. We used this procedure to fit the velocity up-step (Figs. 5A and 6A) and the first SHS event (Figs. 5C and 6C) for both the single- and two-state-variable data. Sometimes a set of initial parameter values generates a good fit using one state evolution law, but not the other; this was the case for the velocity down-step in the one-state-variable data (Fig. 5B). To fit that event, a different set of values for ξ0 can be used for each state evolution law, and the fits generated individually by selecting the desired law.
It may also occur that the routine cannot produce a fit for one of the evolution laws; this occurred for the second SHS in both data sets (Figs. 5D and 6D), as well as the velocity down-step in the two-state-variable data (Fig. 6B). If both state evolution laws are selected, the fitting routine attempts to produce an aging-law fit first. If the aging-law attempt fails, then the routine stops and no fitting results are displayed on the Fitting axes or in the Fitting Parameters window for either law. If the slip-law attempt fails, results for a successful aging-law attempt are still displayed. The law that failed to produce a fit can be deselected, and the routine run successfully for the other law. In the case of the synthetic data, we could only successfully generate fits for these specific events while using the slip law.
The first SHS event in the two-state-variable data illustrates the sometimes ambiguous nature of the rate and state friction equations. For this event, we can obtain fits using a single state variable (Fig. 6E) that are just as good as fits obtained using two state variables (Fig. 6C). Hence, without prior knowledge of how the data were generated, one would be unable to properly determine the state evolution behavior from examining this particular SHS event. Note that when using two state variables, if dc1 ≈ dc2, then the equations are equivalent to using a single state variable with b = b1 + b2.
The second SHS event in the two-state-variable data demonstrates a case where using the weight function can produce a better fit (Fig. 6D). Both the weighted and unweighted fits for this event were produced using the slip law. The unweighted fit visibly overshoots the peak friction that is attained at reloading after the hold, and also does not quite capture the slight oscillation that occurs before friction attains a new steady-state value. The clear visual offset of the unweighted fit motivates the use of the weight function. For this event, we applied the weight function by choosing an exponent of p = 0.2 and setting the weight location at the peak friction value. The resulting weight values are shown in Figure 6F. When using the weight function, we find that it is best to start with a small exponent (e.g., p = 0 corresponds to an unweighted fit) and slowly increase the value until the desired effect is achieved. The weighted fit does a much better job at reproducing the entire frictional response during the SHS event. We note that it is entirely possible that using a different set of parameter values could have generated a good fit without use of the weighting function. The ability to apply a weight function to the fitting procedure increases the utility of RSFit3000.
In this section, we use experimental data previously published by Ikari et al. (2014) to demonstrate how the detrending procedure can influence the fitted parameter values. We use friction data that were measured on powdered gouge material collected from the Alpine fault in New Zealand during the Deep Fault Drilling Project (Experiment p3783; Ikari et al., 2014). The measurements were conducted using a double-direct shear apparatus at Rock and Sediment Mechanics Laboratory, Penn State University (Pennsylvania, USA).
A velocity up-step from 10 μm/s to 30 μm/s, conducted at 30 MPa normal stress, is shown in Figure 7A, along with four fits to the data. The up-step is identical to that shown in Figure 2C of Ikari et al. (2014). Each of the fits in Figure 7A was generated in RSFit3000 using a different detrending slope but identical detrending reference points. The fit labeled “ideal” was generated using the detrending procedure described in the Workflow Section: we picked two points on the data to fit a detrending line, finding a slope of m = 7.29 × 10−7 1/μm. The fits labeled A, B, and C were generated by manually entering different values for the detrending slope. Each fit was generated using the same location for the velocity step, and identical initial parameter values. So, the only thing that differentiates the fits is the value of the detrending slope. Figure 7A shows fits that were generated using the slip law, but we also generated fits using the aging law. The aging-law fits are visibly indistinguishable from the slip-law fits.
Because the detrending procedure is a subjective process, there is no “correct” fit. The four fits all have similar R2 values, and could be said to approximately fit the data equally well. The fitted parameter values are also similar for each of the fits. Figures 7B (slip law) and 7C (aging law) show the fitted parameter values for all of the fits, normalized by the value obtained in the “ideal” fit. Note that the scatter is smallest for the quantity for both evolution laws. We emphasize that the detrending procedure does affect the fitted parameter values, and therefore any parameters associated with detrending should be reported when publishing fits to experimental data.
In conclusion, we recommend a set of best practices to follow when analyzing rock friction experiments using RSFit3000. RSFit3000 is designed to enable easy reporting of all relevant information that is involved in producing a model fit to an experimental friction event. This information includes the fitted parameter values and their error intervals, the covariance matrix between the parameters, the R2 value of the fit, parameters describing the event, parameters used to detrend the data, and parameters used to weight the data. All of this information can be saved in a MATLAB structure array when the user produces a fit. Any published fits should include these structure arrays; we also include MATLAB scripts in Supplemental File 1 (footnote 1) for writing information contained in RSFit300 structure arrays into Microsoft Excel and text files. Finally, we suggest that the stiffness should be treated as a fitting parameter.
We thank Chris Marone and Pathikrit Bhattacharya for careful reviews that helped to improve this manuscript and the RSFit3000 code. We also thank Dan Faulkner, Christine McCarthy, Kristina Okamoto, Hannah Rabinowitz, Alex Rösner, and Caroline Seyler for code testing. RSFit3000 is maintained on GitHub, and the most up-to-date version can be found at https://github.com/rmskarbek/RSFit3000.