Crustal Deformation and Fault Mechanics

Numerical Simulation of Slow Slip and Dynamic Rupture in the Cascadia Subduction Zone

One of the most exciting discoveries in recent decades has been the recognition that many subduction zones undergo transient slip events at depths below the locked mega-thrust zone. These slow slip events, or silent earthquakes, have been detected by GPS networks in Cascadia, southwest Japan, Mexico, New Zealand, and Alaska. We have also discovered slow slip events beneath Kilauea volcano in Hawaii. Slow slip events in Cascadia propagate along strike at rupture speeds of roughly 10 km/day, while others such as the Tokai slow event slip for several years in largely the same locality.

Understanding the physics of slow slip events and how they differ from normal, fast earthquakes is one of the most pressing current challenges in seismology. Understanding what causes some events to slip slowly compared to others that slip rapidly, radiating damaging seismic waves would provide a fundamental insight into how earthquakes work. The societal importance is that every slow slip event adds stress to the adjacent locked megathrust zone bringing it closer to failure. For this reason we are very interested in the triggered earthquakes that occur in Hawaii, Japan, and New Zealand.

We have been exploring the possibility that dilatant strengthening may provide a reasonable explanation for slow slip events. Dilatancy is the tendency for the pore space in granular materials to expand during shear. The increased pore-space decreases the fluid pressure within the fault zone, increasing fault strength according to the effective stress principal, unless pore-fluid can flow into the fault zone fast enough to keep up with the rate of dilatancy. Our hypothesis is that frictional weakening allows slip to nucleate, but that as the slip-rate increases fluid flow becomes increasingly unable to keep up with the rate of dilatancy. Depending on the constitutive parameters and the ambient effective normal stress, dilatancy can quench the instability resulting in a slow slip event.

Our approach has been to use numerical modeling and theoretical analysis to understand the conditions when dilatancy stabilizes fault slip. This involves coupling the equations of elasticity, rate-and state dependent friction laws, to fluid transport equations, generating a large system of coupled partial differential equations. We take depth-variable frictional properties, based in part on lab experiments, and assume low effective stress, presumed to be driven by dehydration reactions, in the depth range of slow slip events. Simulations reveal generic behavior: dynamic events (DE) repeat every few hundred years, and between each DE is a quiescent period and then a long sequence of SSE. If the width of the low effective stress region exceeds a critical dimension, the SSE penetrate up-dip with time. During this period, the SSE moment rates generally (but not monotonically) increase with time. Eventually slip speeds become high enough to induce thermal pressurization, which nucleates a DE. The predicted behavior, in terms of SSE slip, stress drop, and repeat time bear many similarities to SSE in Cascadia.

Our simulation methodology is as follows. Partial differential equations in pressure and temperature are solved on profiles normal to the fault. The diffusion equations are discretized in space using finite differences on a nonuniform mesh having greater density near the fault. The full system of equations is a semiexplicit index-1 differential algebraic equation (DAE) in slip, slip rate, state, fault zone porosity, pressure, and temperature. We integrate state, porosity, and slip explicitly; solve the stress-balance equation on the fault; and integrate pressure and temperature implicitly. For speed, we adaptively switch between this DAE time-integration technique and a pure-ODE approach; the latter is faster but less accurate. We use relative error control to calculate the time step adaptively over approximately twelve orders of magnitude. The nonlinear equations in pressure, temperature, and slip speed decouple by fault cell. Therefore, at each time stage, many relatively small nonlinear equations are solved rather than one large one. This decoupling is in itself faster and also provides opportunities for additional speedups. The primary work is solving linear systems associated with solving the nonlinear equations. For efficiency, we group fault cells by physical properties, perform only one sparse LU factorization per group, and efficiently update the LU factorization at each solve, yielding a method that is linear in the number of diffusion profile nodes. The stress-slip relationship on the fault is calculated by a boundary element method; computationally, it amounts to a matrix-vector product. We compress the matrix to an H-matrix (see below for more). These algorithms together have two important properties: the overall work per time step is almost linear in the number of fault elements, and the time step is governed only by stability considerations related to rate-state friction and accuracy of the solution.

Figure 1 shows slip rate from a simulation. The along-fault direction is the x axis; time increases nonlinearly along the y axis. Time progression is indicated on the right side. Friction properties a and b and bacground effective normal stress are shown at the bottom of the figure. Three excerpts from the simulation are shown, each excerpt separated from the others by a black horizontal line. In the first excerpt, an earthquake nucleates at the transition between the locked and ETS zones. A quiescent period follows. Then a long sequence of slow-slip events occur. These gradually increase in size. Two periods are shown in the next two excerpts. Not shown is the the next earthquake, which starts the cycle again.

Figure 1. Simulated slip rate on a 1D model fault.

To test model predictions against GPS data, we develop a pseudo-3D method that accounts for the markedly non-2D geometry of the plate interface. The approach employs 3D elastic Green's functions but assumes that slip rate is a function of depth only, as computed in the physics based model. Figure 2, top left, shows the average slip distribution for one case, with slip tapered along strike to match a typical SSE (contours at 20, 30, 40, 50 km depth). Figure 2, bottom left, shows the predicted displacements compared to the mean SSE displacement over the past decade. The good agreement strongly suggests that the depth distribution of physical properties is reasonable.

Figure 2, top and bottom right, shows the model inter-ETS slip-rate distribution; the subduction interface is fully locked to 30 km depth, with slip-rate increasing to the plate-rate at about 40km. The predicted surface displacement rates are compared to GPS inter-ETS velocities in the figure below. The fit is quite good in the eastern part of the network suggesting that the down-dip extent of fault locking is accurate. The GPS station velocities are somewhat over-predicted at coastal stations, suggesting some shallow creep may be required to fit the data. The model SSE do not accommodate the full plate motion at depths below 25 km, suggesting that significant slip deficit may exist at depth. Indeed the model earthquakes rupture through the ETS zone, with obvious implications for seismic hazards in the Seattle region.

Figure 2. Surface displacement and slip rate.

Our focus now is on understanding several correlated questions related to Cascadia: whether there is creep updip of the ETS transition; whether dynamic ruptures indeed propagate into and through the ETS zone as our models to date predict; and the extent of the slip deficit in the ETS zone.

Meanwhile, we are developing algorithms and software for fully 3D simulations. Following is a very preliminary and under-resolved example of a simulation of the ETS zone.

Figure 3. Simulated slip rate on a 2D model fault.

The down-dip direction is to the right.  No attempt has been made to obtain realistic physical phenomena at this early stage.

Our present efforts are focused on algorithm development. H-matrix compression is increasingly important at the problem sizes necessary to model a two-dimenstional fault. Recently, we submitted a paper that uses careful error analysis to derive a way to increase H-matrix compression relative to the standard method; see this paper also for a description of H-matrices. For our problems, we obtain approximately a factor of four improvement in compression.

Earlier we remarked that the work per time step is almost linear in the number of fault elements coming from the discretization. We are presently focused on developing and implementing a variable time step integrator so that we can actually do better than a linear amount of work per time step. A key component of this algorithm is the ability to perform the matrix-vector product y(r) = B(r,c)*x(c), where B is an H-matrix, very quickly even as the index sets r and c change rapidly. We have a solution to this problem and are now focused on the time integrator itself.