Statistical physics of fracture, friction, and earthquakes

The present status of research and understanding regarding the dynamics and the statistical properties of earthquakes is reviewed, mainly from a statistical physical viewpoint. Emphasis is put both on the physics of friction and fracture, which provides a microscopic basis for our understanding of an earthquake instability, and on the statistical physical modelling of earthquakes, which provides macroscopic aspects of such phenomena. Recent numerical results from several representative models are reviewed, with attention to both their critical and their characteristic properties. Some of the relevant notions and related issues are highlighted, including the origin of power laws often observed in statistical properties of earthquakes, apparently contrasting features of characteristic earthquakes or asperities, the nature of precursory phenomena and nucleation processes, and the origin of slow earthquakes, etc.

Figure 1

A portion of a plate (of thickness $w$) under tensile stress $\sigma$ (model I loading) containing a linear crack of length $2l$. For a further growth of the crack length by $2dl$, the elastic energy released from the annular region must be sufficient to provide the surface energy $4\Gamma wdl$ (extra elastic energy must be released for finite velocity of crack propagation).

Figure 2

Schematic variation of failure probability $F(\sigma )$ with stress $\sigma$ for a disordered solid with volume $L_1^d$ or $L_2^d$ ($L_2>L_1$).

Figure 3

The average velocity of the crack front is plotted against external force ($f=G-\Gamma$, where $G$ is the mechanical energy release rate and $\Gamma$ is the fracture energy). For $T=0$ a depinning transition is seen. For finite temperature the subcritical creep is shown. From 215.

Figure 4

Variation of average crack-front velocity against the mechanical energy release rate is shown for Botucatu sandstone. The subcritical creep region and supercritical power- law variations are shown in the top left and bottom right insets, respectively. For the subcritical regime, the data are fitted with a function $v\sim e^{-C/(G-⟨\Gamma ⟩{)}^{\mu }}$ for $\mu =0.60$ and $⟨\Gamma ⟩=65\mathrm{J}{\mathrm{m}}^{-2}$. For the power-law variation [$v\sim (G-G_c{)}^{\theta }$] in the supercritical region, $G_c=140\mathrm{J}{\mathrm{m}}^{-2}$ and $\theta =0.80$. From 214.

Figure 5

The variations of the maximum static friction with the amplitude ($K_I$) of the interchain potential for different values of the lower-chain stiffness ($K_s$). The limit $K_s\to \infty$ corresponds to the Frenkel-Kontorova model. But it is clearly seen that even for finite $K_s$ (i.e., when the lower chain can relax) there is a finite value of interchain potential amplitude up to which the static friction is practically zero and after which it increases , signifying the Aubry transition. From 162.

Figure 6

Schematic representation of the two-chain version of the Tomlinson model with (a) no disorder, (b) Cantor set disorder, and (c) the effective substrate potential. From 76.

Figure 7

The overlap distribution for $s=9$, $r=8$ is shown. The dotted curve shows the average distribution with random offset and the continuous curve is that without random offset. The distribution is qualitatively different from the Gaussian distribution expected for random disorder. From 76.

Figure 8

The Burridge-Knopoff (BK) model in one dimension.

Figure 9

The magnitude distribution of earthquake events of the 1D BK model with nearest-neighbor interaction for various values of the friction parameter $\alpha$: (a) for larger $\alpha =1$, 2, 3, 5, and 10, and (b) for smaller $\alpha =0.25$, 0.5, 0.75, and 1. The parameters $l$ and $\delta$ are fixed to be $l=3$ and $\delta =0.01$. The system size is $N=800$. From 177.

Figure 10

The local recurrence-time distribution of the 1D BK model with nearest-neighbor interaction for various values of the frictional parameter $\alpha$. Large events with $M>M_c=3$ or 4 are considered. The parameters $l$ and $\delta$ are $l=3$ and $\delta =0.01$. The recurrence time $T$ is normalized by its mean $\overline{T}$. The total number of blocks is $N=800$. The inset represents semilogarithmic plots including the tail part of the distribution. From 176.

Figure 11

The event frequency preceding the large event of $M>M_c=3$ vs the distance from the epicenter of the upcoming mainshock of the 1D BK model with nearest-neighbor interaction. The parameters $\alpha$, $l$, and $\delta$ are $\alpha =1$, $l=3$, and $\delta =0.01$. The data are shown for several time periods before the mainshock. The inset represents similar plots with longer time intervals. The system size is $N=800$. From 177.

Figure 12

The magnitude distribution $R(M)$ of the 2D BK model with nearest-neighbor interaction for various values of the friction parameter $\alpha$. The other parameters are $l=3$ and $\delta =0.01$. (a) $R(M)$ for smaller values of the friction parameter $0\le \alpha \le 3$, (b) $R(M)$ for larger values of the friction parameter $3\le \alpha \le \infty$. The system size is $60\times 60$. From 178.

Figure 13

Phase diagram of the 2D BK model with nearest-neighbor interaction in the plane of the friction parameter $\alpha$ versus the elastic parameter $l$. The parameter $\delta$ is $\delta =0.01$. From 178.

Figure 14

The local magnitude distribution preceding a mainshock of $M>M_c=5$ of the 2D BK model with nearest-neighbor interaction. The parameters are $\alpha =1$, $l=3$, and $\delta =0.01$. The data are shown for several time periods before the mainshock. The system size is $60\times 60$. From 178.

Figure 15

The magnitude distribution $R(M)$ of the 2D BK model with long-range interaction for various values of the friction parameter $\alpha$. The other parameters are $l=3$ and $\delta =0.01$. (a) $R(M)$ for smaller values of the frictional parameter $0\le \alpha \le 10$. (b) $R(M)$ for larger values of the frictional parameter $10\le \alpha \le \infty$. The system size is $60\times 60$. From 179.

Figure 16

The phase diagram of the 2D BK models with long-range interaction in the plane of the friction parameter $\alpha$ vs elastic parameter $l$, compared with that of the 2D BK model with short-range interaction. The parameter $\delta$ is set to $\delta =0.01$. From 179.

Figure 17

The magnitude dependence of the mean displacement $\Delta \overline{u}$ at each seismic event of the 2D BK model with long-range interaction. In the main panel, the friction parameter $\alpha$ is varied with a fixed system size $60\times 60$, while in the inset the system size $N$ is varied for the case of $\alpha =30$. The parameters $l$ and $\delta$ are fixed to $l=3$ and $\delta =0.01$. From 179.

Figure 18

The magnitude distribution $R(M)$ of earthquake events of the 1D viscous BK model ($\eta =0.02$) with $\delta =0.01$. The dimensionless grid spacing $d$ is varied in the range $1\ge d\ge 1/32$. (a) and (b) represent the cases of $\alpha =1$ and 3, respectively. The system size is $L=dN=200$. From 180.

Figure 19

The event frequency in the time period $t\nu =0–0.01$ immediately before the mainshock of $M>M_c=2$ of the 1D viscous BK model ($\eta =0.02$) with $\alpha =1$ plotted vs $r$, the distance from the epicenter of the upcoming mainshock. The dimensionless grid spacing $d$ is varied in the range $1/4\ge d\ge 1/32$. The parameter $\delta$ is fixed to $\delta =0.01$. The system size is $L=dN=200$. The inset represents the peak position of the event frequency, corresponding to the range of the doughnutlike quiescence, as a function of the dimensionless grid spacing $d$. The doughnutlike quiescence vanishes in the continuum limit $d\to 0$. From 180.

Figure 20

The magnitude distribution of the 1D BK model with the RSF law, for the cases of (a) a strong frictional instability $b>b_c$ and (b) a weak frictional instability $b<b_c$, with $b_c=2l^2+1$. The parameter values are $a=0$, $c=1000$, $\nu ={10}^{-8}$, $v^{*}=1$, and $l=3$ in (a) and $a=1$, $b=5$, $c=1000$, $v^{*}=1$, and $l=5$ in (b). The borderline $b$ value is $b_c=19$ in (a) and $b_c=51$ in (b). The system size is $N=800$ in (a) and $N=1200$ in (b). (a) From 193. (b) From 181.

Figure 21

The magnitude distribution of the 2D BK model with the RSF law, for the case of (a) a strong frictional instability $b>b_c$, and of (b) a weak frictional instability $b<b_c$, with $b_c=4l^2+1$. The parameter values are $a=1$, $c=1000$, $\nu ={10}^{-8}$, $v^{*}=1$, and $l=2$ in (a), and $a=1$, $c=1000$, $\nu ={10}^{-8}$, $v^{*}=1$, and $l=2$ in (b). The borderline value is $b_c=17$ in both (a) and (b). The system size is $N=60\times 60$ in (a), and $N=30\times 30$ in (b). From 122.

Figure 22

The typical rupture process realized in the 1D BK model with the RSF law for (a) a strong and (b) a weak frictional instability, corresponding to (a) $b>b_c$ and (b) $b<b_c$ with $b_c=2l^2+1$. The color represents the rupture velocity. The parameter values are $a=1$, $c=1000$, $\nu ={10}^{-2}$, $v^{*}=1$, and $l=5$ for both (a) and (b), corresponding to $b_c=51$, whereas $b=60$ in (a) and $b=3$ in (b). From 181.

Figure 23

(a) Recurrence of Kamaishi earthquakes of nearly the same magnitudes and recurrence intervals. (b) Cumulative seismic moment of Kamaishi earthquakes. (c) Coseismic slip distribution of the 1995 and 2001 Kamaishi earthquakes estimated from seismic waveforms. Broken contours and blue contours denote the seismic slip of the 1995 and 2001 earthquakes, respectively. From 195.

Figure 24

Spatial distribution of cumulative slip for 30 days of afterslip of the 2003 Tokachi-oki earthquake ($M=8.0$), off Hokkaido, northern Japan, estimated from GPS data (color contours) by 169. Black contours with an 0.5 m interval show seismic slip in the 1973 Nemuro-oki (right), 1968 Tokachi-oki (left), and 2003 Tokachi-oki (center) earthquakes (294). The black star and small circles denote the epicenter and aftershocks of the 2003 earthquake.

Figure 25

Spatial distribution of interplate coupling estimated from geodetic data (colored circles) along the Sunda trench, where the Australian plate subducts beneath the island of Sumatra. Dark circles indicate that the plate interface is nearly locked and strain is accumulated during an interseismic period, and light circles indicate that continuous aseismic sliding occurs and strain is not accumulated. Contours with a 5 m interval show seismic slip in the 2004 Sumatra-Andaman ($M=9.1$) and the 2005 Nias-Simeulue ($M=8.7$) earthquakes. Lines show the approximate source are as of the 1797 and 1833 great earthquakes. From 51.

Figure 26

Snapshots of simulated slip rate $V$ on the model plate interface normalized by the relative plate velocity $V_{\mathrm{pl}}$ in a model for recurrence of great earthquakes along the Nankai trough, central Japan. Different colors show seismic slip rates, stable sliding with sliding velocity nearly equal to the plate velocity, and nearly locked, respectively. Adapted from 109.

Figure 27

The duration of slip events vs $r/r_c$ in numerical simulations using the RSF law, where a circular asperity of radius $r$ with velocity-weakening friction is embedded in a fault of velocity-strengthening friction. $r_c$ denotes the critical fault radius for unstable slip defined in the text. From 131.

Figure 28

Schematic diagram of the relation between frictional force and slip distance during the slip-weakening process, where $F_i$, $F_p$, $F_d$, and $D_c$ denote the initial force, the peak frictional force, the dynamic frictional force, and the critical slip distance, respectively. The shaded area indicates the fracture energy $G_c$.

Figure 29

The size distribution of the OFC model under open boundary conditions for various values of the transmission parameter $\alpha$. The slope of the data gives the value of $1+B$, which is shown in the figure. From 140.

Figure 30

The time dependence of the frequency of aftershocks (a) and of foreshocks (b) of the OFC model under open boundary conditions on a log-log plot for several values of the transmission parameter $\alpha$. Mainshocks are events of size greater than $s\ge s_c=100$. The time $t$ is measured with the occurrence of a mainshock as the origin. The range parameter is $r_c=10$. From 140.

Figure 31

Log-log plots of the local recurrence-time distributions of large avalanches of their size $s\ge s_c=100$ for a fixed range parameter $r_c=10$, with variation of the transmission parameter $\alpha$. The arrow in the figure represents the expected peak position for $\alpha =0.245$ corresponding to the period $T^{*}=1-4\alpha =0.02$. The inset is a magnified view of the main peak for the case of $\alpha =0.17$. From 140.

Figure 32

Snapshots of the stress-variable distribution of the OFC model under open boundary conditions for the case of $\alpha =0.2$: (a) immediately before a large event at time $t=t_0$, (b) immediately after this event, (c) immediately before the following event which occurs at time $t=t_0+T^{*}(T^{*}=0.2)$, and (d) immediately after this second event. The two events are of sizes $s=15891$ and 15 910 on an $L=256$ lattice. The region surrounded by bold lines represents the rupture zone, while the star symbol represents an epicenter site which is located at the tip of the rupture zone. From 140.

Figure 33

The stress-variable distribution $D(f)$ of each site contained in the rupture zone of an asperity event of the OFC model under open boundary conditions, just before (a) and after (b) the asperity event. An asperity event is defined here as an event of size greater than $s\ge s_c=100$ belonging to the main peak of the local recurrence-time distribution function. The transmission parameter is $\alpha =0.17$. The inset is a magnified view of the main peak. From 140.

Figure 34

The time sequence of the stress-variable distribution $D(f)$ at the time of toppling of each site contained in the rupture zone of the asperity events. An asperity event is defined here as an event of size greater than $s\ge s_c=100$ belonging to the main peak of the local recurrence-time distribution function. The transmission parameter is $\alpha =0.17$. As the events repeat, the stress-variable distribution at the time of toppling becomes more and more concentrated on the borderline value $f_c=1$. From 140.

Figure 35

The fiber bundle consists initially of $N$ fibers attached in parallel to a fixed and rigid plate at the top and a downwardly movable platform from which a load $W$ is suspended at the bottom. In the equal-load-sharing model considered here, the platform is absolutely rigid and the load $W$ is consequently shared equally by all the intact fibers.

Figure 36

The simple model considered here assumes uniform density $\rho ({\sigma }_{th})$ of the fiber strength distribution up to a cutoff strength (normalized to unity). At any load per fiber level ${\sigma }_t$ at time $t$, the fraction ${\sigma }_t$ fails and $1-{\sigma }_t$ survives.

Figure 37

The breaking rate $R(t)$ vs the rescaled step variable $t_f/t$ for the uniform threshold distribution for a bundle of $N={10}^7$ fibers. Different symbols are for different excess stress levels $\sigma -{\sigma }_f$: 0.001 (circles), 0.003 (triangles), 0.005 (squares), and 0.007 (crosses). From 225.

Figure 38

The burst size distribution for different values of $x_0$ in the equal-load-sharing model with uniform threshold distribution. The number of fibers is $N=50000$. From 224.

Figure 39

Crossover signature in the local magnitude distribution of earthquakes in Japan. During the 100 days before the mainshock the exponent is 0.60, much smaller than the average value 0.88. From 139.

Figure 40

(a) Schematic representation of the rough Earth surface and a tectonic plate. (b) The one-dimensional projection of the surfaces, forming overlapping Cantor sets.

Figure 41

(a) Two Cantor sets along the axes $r$ and $r-r^{\prime }$. (b) The overlap $s_1(r)$ along the diagonal. (c) The corresponding density ${\rho }_1(s)$.

Figure 42

The overlap densities (probabilities) $\rho (s)$ at various generations; (a) zeroth, (b) first, (c) second, and (d) infinite generation.