Project Title: 3-D Elastic Finite-Difference Simulation of a Dynamic Rupture
Investigators: Dr. Kim B. Olsen, (805) 893 7394, kbolsen@quake.crustal.ucsb.edu
Institution: Institute for Crustal Studies, UC Santa Barbara

Problem

How can we simulate realistic spontaneous rupture due to shear cracks in heterogeneous 3-D elastic media including a free-surface?

Methodology

Our approach to solve the problem stated above is a staggered-grid velocity-stress finite-difference solution to the three-dimensional elastic wave equation. We use a new implementation of the boundary conditions on the fault which allows the use of general friction models, including slip-weakening and rate-dependent laws with a length scale (Do), for simulating spontaneous rupture propagation along an arbitrarily loaded planar fault. The crack is introduced using the appropriate mixed boundary conditions along a plane that coincides with one of the mesh planes. Our numerical method includes full elastic wave interactions as well as all the usually accepted properties of dynamic faulting, including frictional instability, rupture initiation from a finite initial patch, spontaneous rupture growth and healing both by stopping phases and rate-dependent friction.

Results
Spontaneous rupture initiation, propagation, and healing

We use the method to model rupture starting from a localized asperity on circular and rectangular faults. The shape of the fault is close to circular but depending on details of the friction law tends to become more or less elongated in the direction of the dominant traction drop (the in-plane direction). The rupture shows a strong tendency to propagate at super-shear speeds in the direction of in-plane shear, promoted by high initial stresses and small slip weakening distances. Rate-weakening friction tends to reduce super-shear rupture speeds and generally produces narrow rupture pulses (Fig. 1). Comparison of scalar and vector boundary conditions for the friction shows that slip is dominant along the direction of the pre-stress, with the largest deviations near the rupture front and the edges of the fault (Fig. 1 and 2).

Dynamic Simulation of the 1992 Landers Earthquake

In order to test the method on a realistic earthquake scenario we simulated the 1992 Landers earthquake dynamically (Olsen et al., 1997, in Science 278). The initial stress field was estimated from the slip distribution proposed by Wald and Heaton. A constant yield stress level of 120 bars and Do equal to 80 cm produces a total rupture time and final slip distribution in agreement with kinematic inversion results. Before the simulation we constrained the initial stress on the fault to values just below the specified yield level (120 bars) in order to prevent rupture from several locations. We use the same regional one-dimensional (1D) model of velocities and densities as Wald and Heaton. Rupture was forced to initiate by lowering the yield stress in a small patch of radius 1 km inside a high-stress region near the hypocenter towards the southern end of the fault, as inferred from the kinematic results (Fig. 3).

The earthquake effectively released the stress on all the most heavily loaded areas of the fault (Fig. 4A-C). Note that the stress has increased within and particularly near the edges
hin and particularly near the edges of the areas that showed little or no slip (for example, about 25-27 km and 37-48 km from the northwestern edge of the fault at 3-12 km depth) to values similar to or larger than the yield stress specified for the earthquake that just occurred. This result has to our knowledge not been reported in previous studies and suggests that these areas could be the location of nucleation for a future earthquake on the fault. The final slip distribution (maximum 6.1 m) represents a smooth version of the one used to compute the initial stress (Fig. 4D).

There are important similarities and differences between the kinematic inversion results and our dynamic reconstruction of the Landers event. Both methods find that a confined rupture pulse generally propagates unilaterally along a continuous path on the fault. A more detailed comparison of the rupture pulse is complicated by the relatively lower resolution of the kinematic inversion. The rupture velocity shows a strong variation during the earthquake for the dynamic simulation (Fig. 5A). On average, the rupture here propagates with velocities close to the local S-wave speed and terminates about 22 s after initiation, in agreement with the kinematic results. However, the rupture velocity varies strongly from sub-sonic to super-sonic during the simulation. The sub-sonic rupture velocities generally occur within and near the low-stress areas on the fault while super-sonic ones dominate within highly-stressed patches of the fault where the rupture resistance is relatively low. Super-sonic velocities are particularly dominating near the surface, allowing the shallow fault to lead the rupture. This result is not expected from the kinematic inversions where the rupture generally propagates faster at depth. These shallow super-sonic rupture velocities are likely caused by the free surface which promotes the generation of S-P converted head waves (S*). The nearly-horizontal time-distance contours in certain areas on the fault represent 'jumps' into higher-stresses areas (e.g., at 4-10 km in Fig. 5 and snapshots at 18.7-19.5 s in Fig. 3).

We have compared the synthetic seismograms generated by our model with the observed velocity seismograms at four stations located in the forward (Yermo, YRM) and backward (Desert Hot Spring, DHS) rupture directions as well as southwest (Big Bear Lake, BIG) and northeast (Amboy, AMB) of the fault. (Fig. 6). The main features of the low-frequency ground motion for both amplitude anci waveform are reproduced by the synthetic seismograms at the four stations. The best fit is obtained for the relatively stronger ground motion recorded in the forward rupture direction (YER). The fit at the other three stations is somewhat worse, in part related to the weaker ground motion in the back-rupture direction where the relative influence of the crustal structure and fault geometry are larger. Effects of deviations of the simple 1D model is furthermore increased by the large distance from the fault to AMB (65 km). Unfortunately, AMB is the only station available at northeastern azimuths from the fault.

Movies and high-resolution images of the results as well as the paper by Olsen et al. (1997) in Science can be found on the world wide web on: http://quake.crustal.ucsb.edu/~kbolsen/DYN.html

References

Madariaga, R., K.B. Olsen, and R.J. Archuleta (1996). Modeling dynamic rupture in a 3D earthquake fault model, submitted to Bull. Seism. Soc. Am.

Madariaga, R., K.B. Olsen, and R.J. Archuleta (1997). 3-D elastic finite-difference simulation of a spontaneous rupture, SRL 68, 312.

Madariaga, R., K.B. Olsen, and R.J. Archuleta (1997). Rupture dynamics of a 3D planar fault: role of friction and stress heterogeneity, EOS, Trans. Am. Geophys. Union (Supplement, fall 1997).

Olsen, K.B., R. Madariaga, and R.J. Archuleta (1997). Three-dimensional dynamic simulation of the 1992 Landers earthquake, Science 278, 834-838.