Probabilistic model of waiting times between large failures in sheared media

Using a probabilistic approximation of a mean-field mechanistic model of sheared systems, we analytically calculate the statistical properties of large failures under slow shear loading. For general shear $F\left(t\right)$, the distribution of waiting times between large system-spanning failures is a generalized exponential distribution, ${\rho }_T\left(t\right)=\lambda (F\left(t\right))P(F\left(t\right))exp\left[-{\int }_0^{t}d\tau \lambda (F\left(\tau \right))P(F\left(\tau \right))\right]$, where $\lambda (F\left(t\right))$ is the rate of small event occurrences at stress $F\left(t\right)$ and $P(F\left(t\right))$ is the probability that a small event triggers a large failure. We study the behavior of this distribution as a function of fault properties, such as heterogeneity or shear rate. Because the probabilistic model accommodates any stress loading $F\left(t\right)$, it is particularly useful for modeling experiments designed to understand how different forms of shear loading or stress perturbations impact the waiting-time statistics of large failures. As examples, we study how periodic perturbations or fluctuations on top of a linear shear stress increase impact the waiting-time distribution.

Figure 1

Plot of the probability $P\left(F\right)$ that a small earthquake triggers a large system-spanning earthquake [Eq. (4)] for low and high disorder $b=0.1,1.0$, respectively. The (dimensionless) critical small earthquake size necessary for triggering a large earthquake is $S_c=100$. The stress axis is in units of $F_c$, the mean-field value of the critical stress, for which the earthquake occurs at a well-defined stress. In our finite-fault approximation, the actual values of stress at which the fault fails tend to be slightly larger than the mean-field value $F_c$.

Figure 2

Plot of the distribution of waiting times between large system-spanning earthquakes, ${\rho }_T\left(t\right)$, for an external stress $F\left(t\right)=f+\mathrm{{\rm Y}}t$, for low disorder ($b=0.1$) and high disorder ($b=1.0$). We set $f/F_c=0.73$ and $\mathrm{{\rm Y}}/F_c=0.015{\lambda }_0$, and the critical small earthquake size is $S_c=100$. We measure time in units of the inverse small earthquake rate ${\lambda }_0^{-1}$ (effectively, ${\lambda }_0=1$). The density drops off quickly for large waiting times $t$ and is similarly small for very short waiting times. Though the mean waiting times are similar for the two disorders shown, the variances are quite different.

Figure 3

Mean, variance, and skewness [Eqs. (7, 8, 9)] of the large quake waiting-time distribution as functions of the arrest stress $f$ and inverse shear rate $1/\mathrm{{\rm Y}}$, for low disorder ($b=0.1$; left column) and high disorder ($b=1.0$; right column). The arrest stress was swept over a range $f/F_c=0.5$ to $0.95$, and the slow shear rate was swept over a range $\mathrm{{\rm Y}}/{\lambda }_0{F}_c=0.001$ to $0.03$. Time is again measured in units of ${\lambda }_0^{-1}$. The mean scales roughly as $-f/\mathrm{{\rm Y}}$ and is weakly dependent on disorder. The variance is relatively insensitive to the arrest stress at low disorder, but increases with decreasing shear rate and increasing disorder. At high disorder, the arrest stress decreases the variance for $f/F_c$ close to 1. At low disorder, the skewness is insensitive to arrest stress, but changes sign as the shear rate varies. High disorder increases the range of skewness, becoming more positively skewed for $f/F_c\lesssim 1$, and is relatively unskewed for large shear rates.

Figure 4

Low disorder fits: Fitting gamma, Weibull, and Brownian passage-time (BPT) distributions to a histogram generated by our distribution, Eq. (1), for a low disorder condition of $b=0.1$. Other parameters for the probabilistic model were $f/F_c=0.73, \mathrm{{\rm Y}}/{\lambda }_0{F}_c=0.015$, and $S_c=100$. All distributions capture the qualitative trends well. Fitting was performed with matlab's curve-fitting toolbox, which performed nonlinear regression on the simulation data. Fit parameters for these examples are given in Table 2. (Top row) Data points represent observation frequencies of binned waiting times simulated from our probabilistic model; solid curves represent fits of the gamma (left), Weibull (middle), or BPT (right) distributions to the simulated data. (Bottom row) Data points represent residual difference between simulation data and fits in the top row; solid lines represent the line of zero residuals.

Figure 5

High disorder fits: Same as Fig. 4, but for a high disorder condition of $b=1.0$.

Figure 6

Mean, variance, and skewness of the large quake waiting-time distribution in the presence of periodic perturbations to linear shear, normalized by their respective values in the absence of periodic perturbations (denoted by a subscript 0). We plot these quantities as functions of angular frequency $\omega$ and amplitude $F_0$, at low disorder ($b=0.1$; left column) and high disorder ($b=1.0$; right column). The frequencies were swept over a range $\omega {\tau }_0/2\pi =0.001$ to 10, where ${\tau }_0$ is the mean waiting time in the absence of periodic perturbations. We focus on small amplitudes $F_0/F_c=0$ to $0.05$. We fixed the phase of the stress at $\varphi =\pi /4$ (see comment in text). The arrest stress was set to $f/F_c=0.73$ and the slow shear rate was set to $\mathrm{{\rm Y}}/{\lambda }_0{F}_c=0.015$. At low disorder, the mean waiting time, variance, and skewness are approximately equal to their unperturbed values (${\tau }_0, {\sigma }_0^2$, and ${\gamma }_0$) near small integer multiples of the “natural” failure frequency $2\pi /{\tau }_0$, but vary significantly at intermediate frequencies. The effect diminishes as frequency increases. At high disorders, similar effects are observed, but they diminish much more rapidly as frequency increases. Note that the unperturbed values ${\tau }_0, {\sigma }_0^2$, and ${\gamma }_0$ are disorder-dependent (as shown in Fig. 3). White contours or bands represent regions where $\mathrm{mean}/{\tau }_0\approx 1, \mathrm{variance}/{\sigma }_0^2\approx 1$, and black contours or bands represent regions where $\mathrm{skewness}/{\gamma }_0\approx 0$.

Figure 7

Normalized histograms of waiting times between large earthquakes generated by simulations of the probabilistic model subject to white noise stress fluctuations. Histograms have all been normalized so that the total area is 1, giving probability density functions for the waiting times (i.e., $y$-axes are normalized counts). The standard deviation of the fluctuations is $\kappa$, given in units of $F_c=1$. At high disorder, the theoretical prediction is shown, along with the $\kappa =0$ prediction to emphasize the shift in the distribution caused by the stress fluctuations. For values of $\kappa$ smaller than $0.01$ the difference between the $\kappa \ne 0$ and $\kappa =0$ curves are very small, so we omit these results. At low disorder, the high slope of the large earthquake rate invalidates Eq. (14), so a theoretical curve is shown only for the no noise case ($\kappa =0$, top row). We use a linear stress increase $F\left(t\right)=f+\mathrm{{\rm Y}}t$ with parameter values $f/F_c=0.73, \mathrm{{\rm Y}}/{\lambda }_0{F}_c=0.015, S_c=100$, and ${\lambda }_0\delta \tau =0.1$. Low disorder corresponds to a value of $b=0.1$ and high disorder corresponds to $b=1.0$. In plots where only simulation data is shown, the value of $\kappa$ is given in the legend. In plots with theoretical curves, the value of $\kappa$ used for the simulation and theoretical prediction corresponds to legend value for the solid blue curve, if shown, or $\kappa =0$ otherwise.

Figure 8

Numerical solutions for the steady-state phase distribution ${\varrho }_{\mathrm{\Phi }}^{*}\left(\phi \right)$ (solid curves) for different frequencies and amplitudes. We generated histograms by drawing random earthquake waiting times from Eq. (2) and computing the $n\mathrm{th}$ phase via ${\varphi }_n=\omega t_n+{\varphi }_{n-1}$. We ran the simulation for several thousand events to allow the simulation to reach a steady state before collecting phases to form the histograms. (a) Low frequency histogram, $\omega \lesssim \mathrm{{\rm Y}}/F_0$. (b) High frequency histogram, $\omega \gg \mathrm{{\rm Y}}/F_0$. In both cases, as the amplitude $F_0$ increases, the distributions are increasingly sinusoidal and less uniform.