Intermittent Granular Dynamics at a Seismogenic Plate Boundary

TL;DR

This paper models the complex behavior of earthquakes at plate boundaries as a granular medium near jamming, using GPS data to identify scale-invariant fluctuations and proposing a new understanding of fault dynamics.

Contribution

It introduces a granular medium model for seismogenic plate boundaries, linking GPS velocity fluctuations to fault system rearrangements and a characteristic length scale.

Findings

GPS data shows heavy-tailed scaling behavior.

Plate boundary behaves like a granular medium near jamming.

Characteristic length scale of fault structures is approximately 91 km.

Abstract

Earthquakes at seismogenic plate boundaries are a response to the differential motions of tectonic blocks embedded within a geometrically complex network of branching and coalescing faults. Elastic strain is accumulated at a slow strain rate of the order of $10^{-15}$ s$^{-1}$, and released intermittently at intervals $>100$ years, in the form of rapid (seconds to minutes) coseismic ruptures. The development of macroscopic models of quasi-static planar tectonic dynamics at these plate boundaries has remained challenging due to uncertainty with regard to the spatial and kinematic complexity of fault system behaviors. In particular, the characteristic length scale of kinematically distinct tectonic structures is poorly constrained. Here we analyze fluctuations in GPS recordings of interseismic velocities from the southern California plate boundary, identifying heavy-tailed scaling behavior. This suggests that the plate boundary can be understood as a densely packed granular medium near the jamming transition, with a characteristic length scale of $91 \pm 20$ km. In this picture fault and block systems may rapidly rearrange the distribution of forces within them, driving a mixture of transient and intermittent fault slip behaviors over tectonic time scales.

Tables

Table 1. Table 1: Observed and predicted  Ono  et al. ( 2003 ) ratio of earthquake cycle timing τ dn / τ up subscript 𝜏 dn subscript 𝜏 up \tau_{\mathrm{dn}}/\tau_{\mathrm{up}} , where τ dn subscript 𝜏 dn \tau_{\mathrm{dn}} is the duration of an earthquake, equivalent to the time of a rearrangement event in a granular system, and τ up subscript 𝜏 up \tau_{\mathrm{up}} is the time between earthquakes where the elastic strain is accumulated, equivalent to the time in between rearrangement events. The predicted ratio is within the range of observed values. Also shown are observed values for the shear strain rate ϵ ˙ ˙ italic-ϵ \dot{\epsilon} , the crack tip evolution τ d subscript 𝜏 d \tau_{\mathrm{d}} , and the normalized shear strain rate γ ˙ ˙ 𝛾 \dot{\gamma} , substituted in Eq.  5 to yield the predicted ratio of earthquake cycle timing. Further details in the main text.

Parameter	Description	Value
${\tau }_{\mathrm{dn}}$	duration of earthquake	10-600 s (Scholz, 2002)
${\tau }_{\mathrm{up}}$	earthquake recurrence interval	${10}^8-{10}^{11.5}$ s (Scholz, 2002)
$\dot{ϵ}$	shear strain rate	${10}^{-15}$ $s^{-1}$ (Shen et al., 2007)
${\tau }_d$	crack tip evolution	${10}^{-9}$ s (Ben-Zion, 2008)
$\dot{\gamma }=\dot{ϵ}{\tau }_d$	normalized shear strain rate	$8\times {10}^{-24}$
${\tau }_{\mathrm{dn}}/{\tau }_{\mathrm{up}}$	observed ratio	$6\times {10}^{-6}-3\times {10}^{-11}$
$t_{\mathrm{dn}}^{pred}/t_{\mathrm{up}}^{pred}$	predicted ratio (eq. 5)	$5\times {10}^{-9}$

Equations

$$\delta {\mathbf{v}}_{\text{i}}={\mathbf{v}}_{\text{i}}-⟨\mathbf{v}(y_{\text{i}})⟩,$$

$$C(r)=\frac{⟨\delta \mathbf{v}(r)\cdot \delta \mathbf{v}(0)⟩}{⟨\delta {\mathbf{v}}^2⟩}.$$

$${\alpha }_n(x)=\frac{⟨x^{2n}⟩}{C_n⟨x^2{⟩}^n}-1,$$

$$F(r)=e^{-(r/\xi {)}^{\beta }},$$

$${\tau }_{\mathrm{dn}}/{\tau }_{\mathrm{up}}=7.9{\dot{\gamma }}^{0.4},$$

Full text

Fault slip at seismogenic plate boundaries is rarely continuous, instead occurring intermittently during short duration earthquakes which last from seconds Abercrombie (1995) to minutes Ishii et al. (2005) during which fault slip may reach up to 50 meters Ozawa et al. (2011). Throughout the interseismic interval between large earthquakes, plate boundaries are not dormant but instead slowly accumulate the elastic strain that will be released in future earthquakes Lawson and Reid (1908), while the faults themselves are frictionally locked. In addition to this temporally bi-modal behavior, fault slip at plate boundaries is not localized along a single fault but rather spread across fault systems where faults are bounded by tectonic blocks to form anastomosing fault systems Plesch et al. (2007); Jennings (1975); Wallace et al. (2004).

Dynamic models of activity at seismogenic plate boundaries have generally focused on the physics within narrow ($<1$m) individual fault shear zones Fisher et al. (1997); Dahmen et al. (1998); Aharonov and Sparks (1999, 2004); Bar-Sinai et al. (2013). However while large faults (e.g. the San Andreas fault in California) accommodate the majority of deformation between two plates (e.g. the Pacific and North American plates), they are embedded within a geometrically complex network of branching and coalescing faults fau (2014); Plesch et al. (2007), as can be appreciated in Fig. 1, and the interactive map in California Department of Conservation (2010).

Here we investigate the macroscopic dynamics of planar tectonics within such fault and block systems occurring on scales $>100$km, while considering the suggestion Anderson (2007) that sheared granular systems may provide an analog. Both systems exhibit intermittent dynamics manifested in a gradual increase in elastic energy released by a precipitous event; in slowly sheared granular media these are rearrangement events of single granules, while in seismogenic plate boundaries these are earthquakes. Moreover in both systems the motion at the boundaries is defined and results from the differential motion of the granules or tectonic blocks adjoining a potential slip surface. Here we consider the extent to which scaling laws derived from numerical simulations of slowly sheared 2D deformable granular systems (foam) Ono et al. (2003); Durian (1995, 1997), may describe earthquake cycle activity at tectonic plate boundaries, including spatial variabilities of nominally interseismic GPS velocities across the Pacific-North America plate boundary in California, as well as characteristic earthquake cycle time scales.

We consider $1,106$ nominally interseismic GPS velocity observations, recorded over a period of $\sim 10$ years on average, in the area sheared between the Pacific and North American plates in California (Kreemer et al., 2014). The area is defined as running from the Pacific plate at the bottom, where we define $y=0$, to the North American plate at the top where $y=L=565$ km, while the $x$-axis runs along $0<x<1,200$ km, approximately parallel to the San Andreas fault. The velocities $\mathbf{v}_{\textmd{i}}$ at each GPS reading i are calculated relative to the North American plate, and are shown in Fig. 1a, represented by a vector field. Relative velocities exhibit a clear trend, with lower absolute values closer to the North American plate (here considered fixed), and increasing values closer to the moving Pacific plate, with a maximum of $\sim 50$ mm/year. We calculate the average velocity profile along the $y$-axis, $\langle\mathbf{v}(y)\rangle$, by discretizing the $y$-axis into $N=65$ bins, with $\textmd{d}y=8.7$ km. The discretized velocity profile is then calculated as the average velocity of readings within the same bin. For simplicity in what follows we will adopt the continuous notation $\langle\mathbf{v}(y)\rangle$. The average velocity profile is shown in Fig. 1b, exhibiting non-linearity, i.e. shear bands, as also observed experimentally in different sheared granular systems  Lauridsen et al. (2002). The shear strain rate, defined as the relative velocity of the two plates divided by the distance between them, i.e. $\dot{\epsilon}=(\langle v_{x}(L)\rangle-\langle v_{x}(0)\rangle)/L$, yields $\dot{\epsilon}=7.66\mathrm{\times}10^{-8}$ (1/yr). We define the fluctuation of the velocity of a reading from the average velocity profile as:

[TABLE]

where $\mathbf{v}_{\textmd{i}}$ is the velocity of reading i, and $\langle\mathbf{v}(y_{\textmd{i}})\rangle$ is the average velocity at position $y_{\textmd{i}}$ associated with reading i. To characterize the kinematics of the system as manifested by the velocity fluctuations, we calculate the distribution of velocity fluctuations of the GPS readings $P(\delta\mathbf{v})$, and the equal-time spatial correlation function

[TABLE]

We consider the analogy between the intermittent earthquake events, which release energy accumulated due to the shearing of tectonic plates along a fault, and the intermittent rearrangement events in sheared 2D granular systems, and compare our results to those of numerical simulations of slowly sheared 2D deformable granular medium (foam) Ono et al. (2003). For slow shearing rates the distribution of velocity fluctuations has been found to deviate from a Gaussian (expected for fluid-like flow at high shearing rates), and spatial correlations follow a stretched exponential (deviating from the exponential expected for fluid behavior).

The distribution of velocity fluctuations of the GPS readings in the x and y directions , $P(\delta v_{x})$ and $P(\delta v_{y})$, are plotted in Fig. 2, and their fit to a Gaussian distribution is added for comparison. The distributions deviate significantly from a Gaussian distribution, exhibiting heavy tails. This is further confirmed by calculating the non-Gaussian parameter (NGP) Rahman (1964) for different moment ratios:

[TABLE]

where $C_{n}$ is a known constant. By definition, $\alpha_{n}$ equals zero for Gaussian distributions. Here we calculate $\alpha_{2}$ and $\alpha_{3}$, where in 2D $C_{2}=3$ and $C_{3}=15$. The NGP calculated for the fluctuations in the x direction yield $\alpha_{2}=0.25$ and $\alpha_{3}=0.47$, and in the $y$ direction $\alpha_{2}=0.62$ and $\alpha_{3}=2.18$. For comparison, we calculate these values for a numerically generated Gaussian distribution with an identical population size, yielding $\alpha_{2}=-0.05$ and $\alpha_{3}=-0.18$, highlighting that the NGP values calculated for the original distributions clearly indicate non-Gaussianity.

We note that a part of the San Andreas fault deviates from the horizontal orientation of the plates within the area analyzed here, as depicted in Fig 1a. In order to rule out the possibility that the broad distribution of velocity fluctuations may be an artifact due to this slope, we carry out the analysis for the left third of the plate boundary, where the fault is relatively horizontal (details in the Supporting Information - SI). We find that the heavy tails are even more pronounced (see Fig. S2 in the SI), showing the independence of our results on the geometric form of the San Andreas fault. Moreover, substantial fluctuations are not limited to the San Andreas fault, but rather are exhibited in an area of $\sim 100$km across the fault (see Fig. S1 in the SI).

To calculate the spatial distribution of velocity fluctuations we limit ourselves to $\sim 200$km around the fault (roughly $100\mathrm{km}<y<300\mathrm{km}$), thus avoiding recordings associated with the plates themselves which would mask the relevant correlation. In Fig. 3 we plot the calculated spatial correlation of velocities fluctuations according to Eq. 2. We fit this to the general form of a stretched exponential:

[TABLE]

yielding the values $\beta=0.75$, and $\xi=92$ km, with $R^{2}=0.92$ of the least squares fit. Simulations for 2D sheared foam at low shearing rates Ono et al. (2003) exhibit values of $\beta<1$, deviating from a regular exponential behavior with $\beta=1$, and $\xi$ is found to be associated with the average granule or bubble size.

In this work we investigate macroscopic planar tectonic dynamics across a geometrically complex network of faults and blocks accommodating the deformation between the North American and Pacific plates in California on scales $>100$km. Here we consider the similarity between fault block systems and the intermittent dynamics of granular systems sheared at low shearing rates, both of which exhibit a volatile elastic energy over time. The elastic energy increases gradually as granules slowly deform (as blocks deform during the interseismic phase of the earthquake cycle), and decreases rapidly due to intermittent granule rearrangement events (due to earthquakes). Pursuing this analogy we find that the dynamics displayed by the readings of the seismogenic plate boundary, as manifested by their velocity fluctuations, are consistent with those characterizing analogous 2D sheared deformable granular media, where the velocities $\mathbf{v}_{\textmd{i}}$ are associated with densely packed granules or bubbles constrained between two shearing plates at a low shearing rate. Namely we show that, consistent with findings for slowly sheared granular media Ono et al. (2003), the distribution of velocity fluctuations $P(\delta\mathbf{v})$ deviates from a Gaussian, exhibiting broad tails, and the correlation function decays as a stretched exponential. We note that the velocity profile deviates from the linear profile expected for a homogeneously sheared elastic body, ruling this out as an alternative mechanism.

This analogy, based on kinematic characteristics, allows to deduce the dynamics of the system, and to make predictions concerning time-scales and length-scales associated with the dynamics at a seismogenic plate boundary. One of the primary observations in fault block systems is the ratio of time needed for the coseismic release of elastic strain (earthquakes) and the time it takes for slow plate motions (10-100 mm/yr) to accumulate elastic strain. Following the analogy with sheared granular systems Ono et al. (2003), we consider the timescales of the systems, defining $\tau_{\mathrm{up}}$ as the average duration of the energy buildup in between intermittent sharp energy releases, and $\tau_{\mathrm{dn}}$ as the average duration of these energy drops. Ono et al. Ono et al. (2003) find that the ratio of these characteristic times follows

[TABLE]

where $\dot{\gamma}=\dot{\epsilon}\tau_{\mathrm{d}}$ is the normalized shear rate, and $\tau_{\mathrm{d}}$ is the characteristic time scale in the model, the duration of a rearrangement event Durian (1995); Ono et al. (2003). For the earthquake system, $\tau_{\mathrm{dn}}$ is the characteristic duration of an earthquake, and $\tau_{\mathrm{up}}$ is the earthquake recurrence interval. Observed values and their ratio $\tau_{\mathrm{dn}}/\tau_{\mathrm{up}}$ are given in Table 1. The ratio is found to be consistent with the ratio predicted by Eq. 5, thus relating these timescales to the observed sheared strain rate $\dot{\epsilon}^{\mathrm{obs}}$, and the timescale of the system $\tau_{\mathrm{d}}$, suggested to be the nano-second value for the crack tip evolution time scale Ben-Zion (2008).

Moreover, we recall that in slowly sheared granular material the value of the stretched exponential fit of the spatial correlation function in Eq. 4, $\xi$, is associated with the average granule or bubble size Ono et al. (2003), suggesting that the length scale associated with the seismogenic plate boundary is roughly $90$ km. This value is 4-5 orders of magnitude larger than granular fault gauge (Aharonov and Sparks, 2004), indicating that plate boundaries may be treated as granular systems at a macroscopic scale with the length scale $\xi$ associated with the average mesh size of the fault network. Furthermore, we show that the heavy tails of the fluctuation distributions do not depend on the form of major fault accommodating most of the deformation between the two plates - the San Andreas fault (details in the SI). This further strengthens the concept that the macroscopic tectonic dynamics are governed by the intricate fault and block system, rather than just the major fault.

Lastly, this model provides an explanation for the clustered and intermittent fault behaviors observed in the geologic record as a response to the rapid reorganization of force chains by earthquakes themselves. Geological observations of macroscopic fault system activity have revealed an array of behaviors, for example at the seismically active southern California plate boundary observations of paleo-earthquake activity range from nearly periodic (Fumal et al., 2002; Scharer et al., 2007), to clustered (Grant and Sieh, 1994), and out of phase (Dolan et al., 2007). These observations of diverse fault activity can be explained by a model where plate boundaries are considered as macroscopic granular shear zones near the jamming transition (Ben-Zion, 2008) with effective granule sizes $>10$ km. Granule sizes at this scale enable earthquakes themselves to redistribute forces within plate boundaries by creating and destroying force chain paths and producing complex time evolving fault slip rate histories.

Acknowledgments

The Authors thank Ido Regev and Yohai Bar-Sinai for helpful observations and fruitful conversations.