Finite Source Modeling of Great Earthquakes

Douglas Dreger


The $M_w$ 7.9 Denali fault, Alaska, earthquake, which occurred on November 03, 2002 at 22:12:41.0 UTC initiated on the north dipping Susitna Glacier fault with a reverse sense of motion that apparently triggered primarily strike-slip faulting over approximately 300 km of the Denali-Totschunda fault system (Eberhart-Phillips et al., 2003). This complex fault geometry is similar in many respects to that of the San Andreas fault system, and therefore the study of this great earthquake provides needed insight into processes by which earthquakes may grow into large, multi-segmented ruptures involving complex fault geometries. Central to such study is a detailed analysis of the kinematic source process, and investigation of dynamic rupture models to obtain insight into the complex fault interaction. Teaming with Professor Oglesby (UCR), and Drs. Harris (USGS Menlo Park), Ratchkovoski (UA, Fairbanks), and Hansen (UA, Fairbanks) we investigate the connection between the Denali rupture kinematics and the dynamic source process. This work has been submitted for publication in GRL (Dreger et al., 2003).

Kinematic Source Modeling

Three data sets were employed to constrain the kinematic rupture process: observed surface offsets reported by the Denali Earthquake Field Geology Team (Eberhart-Phillips et al., 2003), continuous and campaign-mode GPS data (Hreinsdóittir et al., 2003), and local and regional distance seismic waveforms. The seismic data from 8 recorders, operated by several groups (e.g. US Geological Survey, the University of Alaska, Fairbanks, and the Alyeska Corporation) provide control on the spatial and temporal pattern of mainshock slip. The seismic station coverage is sparse necessitating the use of GPS data from 38 stations from the continuous Alaska Deformation Array and campaign-mode collections.

The 5-segment fault geometry used for the inversion accounts for the arcuate structure of the Denali-Totschunda fault system (Figure 13.1a), as well as the north-dipping structure of the Susitna Glacier fault, which was identified as the initial rupture plane based on a oblique-reverse first motion focal mechanism obtained by the Alaska Earthquake Information Center. The strike and dip of segment 1, representing the Susitna Glacier fault, was 262 and 48 degrees, respectively. The strikes of segments 2-5 were 83, 102, 117, and 143 all with 90 degree dip since relocated aftershocks do not clearly delineate a non-vertical dip. Each fault segment of the model has a width of 30 km, and the lengths of segments 1-5 are 30, 30, 50, 136, and 100 km, respectively. Segments 1 and 2 share a common hypocenter located at 63.520 N, 147.530 W, and a depth of 7.5 km, and they spatially overlap.

The inversion method that we used is based on the multiple time window approach of Hartzell and Heaton (1983). Each of the fault segments was discretized by 3.75 km by 3.75 km subfaults. We applied slip positivity, moment minimization, and Laplacian smoothing constraints (e.g. Kaverina et al., 2002), and allow for spatially variable rake. The slip rise time was taken to be 2 seconds. Three time windows, each overlapping by 1 second, were employed allowing the slip rise time to vary between 2 to 4 seconds.

The seismic data were corrected for instrument response, integrated to displacement, and bandpass filtered using a two-pass, fourth-order Butterworth filter with corners of 0.01 and 0.5 Hz, and then resampled to 2 sps. GPS displacements processed and distributed by Hreinsdóittir et al. (2003) were used.

Seismic Green's functions were computed using a velocity model from Ratchkovski and Hansen (2002), which is appropriate for the region using a frequency-wavenumber integration approach. Green's functions for the GPS data were computed for a half-space elastic structure using the method of Okada (1985), assuming a shear elastic stiffness of 3.52*$10^{10}$ Pa, which corresponds to the value at 8 km depth in the velocity model used to compute the seismic Green's functions.

It was possible to fit the GPS and seismic data sets quite well with respective variance reductions of 99.4% and 44.3%. The preferred slip model is shown in Figure 13.1b. The slip on Susitna Glacier fault (S1) is confined to depths less than the intersection depth (7.5 km) with the vertical Denali fault. The slip on this portion of the fault sums to a seismic moment of 3.12*$10^{19}$ Nm equivalent to $M_w$ 7.0. Following the initial reverse event very little slip accumulated on the Denali fault for a distance of about 50 km to the east (S2 & S3) consistent with surface slip reports (Eberhart-Phillips et al., 2003). The main Denali segment (S4) of the rupture model has 3 shallow asperities, which are principally strike-slip. The bulk of the slip is shallow (less than 10 km) consistent with the depth of relocated aftershocks (e.g. Ratchkovski et al., 2003). The kinematic model presented in Figure 13.1b includes a 15 km jump in the rupture front at the Denali/Totschunda junction that is an important feature of the dynamic rupture models described below. The jump is not constrained by the seismic data however it qualitatively improves the model in the sense that the amplitude of deep slip is reduced. Slip on the Totschunda fault (S5) in the kinematic model has relatively low amplitude, and is unfortunately poorly constrained by both the GPS and seismic data sets. The total scalar seismic moment was found to be 8.45*$10^{20}$ Nm, and the average and peak slips were 2.14 and 9.94 m.

A relatively high rupture velocity of 3.3 km/s was required to place the recovered strike-slip asperities far enough east to correlate with surface slip observations, which is consistent with the findings of Frankel et al. (2002). An average rupture velocity of 3.3 km/s also maximizes the level of fit to the seismic waveform data. Direct measurement of the S-wave arrival time at the closest recording station, site PS10 only 3 km from the fault, requires a rupture velocity of at least 3.3 km/s. It is evident from those records that a typical rupture velocity of 0.8$beta$ is much too slow. A 3.3 km/s rupture velocity represents 98% of the shear wave velocity at 8 km depth with respect to the velocity model we used. Infrasound recordings also suggest a relatively high rupture velocity (approximately 3.3 km/s) over the length of the Denali fault from just west of the trans-Alaska oil pipeline to the Denali/Totschunda junction in the east (Olson et al., 2003).

Interestingly, the western region of the fault that requires a possibly super-shear rupture velocity also has very low levels of slip suggesting that this portion of the fault may be relatively weak. The teleseismic model of Ozacar et al. (2003) also shows that this region had little slip, and they report that a low Bouguer gravity anomaly in this region indicates that the low density rocks in this region are weak; a claim also supported by the lack of aftershocks (e.g. Ratchkovski et al., 2003). The dynamic modeling also required this section of the Denali fault to be weak in order to allow slip to transfer from the Susitna Glacier to the Denali fault.

Static stress drop was computed from the kinematic slip model using the relation, $\Delta \sigma=2.44*M_0/A^{1.5}$, which assumes a circular fault model for each subfault. Moment ($M_0$) and area (A) were obtained for each of the 5 fault segments by integrating only those subfaults with non-zero slip. The stress drop varies considerably over the fault with values of 2.9, 1.0, 3.2, 5.2 and 2.4 MPa for segments 1-5, respectively.

Figure: . A) 5-segment fault model, the first motion (FM) and Harvard CMT solutions for the mainshock, and regional faults (blue lines). B) 3D projection of the kinematic slip model. C) Map view of fault geometry (red) and finite element mesh used in dynamic models. The mesh is embedded in a much larger buffer mesh to eliminate spurious reflections from the model boundaries. The imposed epicenter is shown by the green star. D) Final slip distribution for the preferred dynamic model. Color version of this figure is found in Chapter IV
\epsfig{file=dreger03_1_1.eps, width=17cm}\end{center}\end{figure*}

Dynamic Source Modeling

To dynamically model the Denali event, we use the 3-D finite element method (Oglesby, 1999; Whirley and Engelmann, 1993) and a slip-weakening friction law. The material and computational parameters used in our models are; Vp 6.25 km/s, Vs 3.55 km/s, density 2790 kg/m3, shear stress 8.46 Mpa, normal stress 27.63 Mpa, static friction 0.6, dynamic friction 0.12, and critical slip-weakening distance 0.624 m. The fault geometry and finite element mesh are shown in Figure 13.1c, and the discretization on the fault is 2 km by 2 km. The fault geometry is a slight simplification of that used in the kinematic models, with the Susitna Glacier fault having the same strike as the Denali fault, and meeting the Denali fault at its deepest extent. Unfortunately, the large size of the fault system requires a coarse discretization of 2-km along the fault, which in turn means that we cannot resolve the slip weakening distance of 0.624 m in our models. However, as will be argued below, this lack of resolution does not significantly affect the major features of our models. For simplicity, we use a homogeneous, infinite-Q half-space in the dynamic models, corresponding to an intermediate depth in the velocity mode used in the kinematic modeling.

The key ingredient in our dynamic models is the fault stress. Ratchkovski and Hansen (2002) show that the maximum principal stress is oriented approximately 85 degrees from the strike of the main Denali fault segment. However, we find it difficult to produce rupture with correct timing on all fault segments using such a simple stress field. A stress field that rotates with the fault strike is one means of solving this problem, but for simplicity we choose a principal compressive stress oriented 65° from the strike of the main Denali segment. The stress field in the western part of the fault, at the intersection between the Susitna Glacier and Denali faults, is likely to be very complicated, due to the complex fault geometry in this region. For computational reasons a simpler fault geometry in this area is used. For both these reasons, we do not attempt to use a single regional stress field to produce fault stresses in this region. Rather, for simplicity, we impose shear and normal stresses that give a stress drop of approximately 30 bars on the Susitna Glacier fault and 15 bars on the western 50 km of the Denali Fault. This difference in stress drop is consistent with the kinematic observations. To better match the slip distribution of the earthquake we reduce the initial shear and normal stress by 50% between 82 and 130 km along the main Denali Segment, reduce the shear stress by 50% on the Totschunda fault starting 20 km from the branching point and continuing to the end of the segment, and reduce the shear stress to zero in the top 2 km of the Totschunda and eastern Denali faults.

The slip for our preferred dynamic model is shown in Figure 13.1d. The seismic moment in this model is 7.28*$10^{20}$ Nm, corresponding to $M_w$ 7.9. While we do not match the fine details of the kinematic model, we match many of the main features, including patches of high slip at both ends of the main Denali branch. The modeled surface slip matches the mapped surface slip quite well in both pattern and amplitude. Most importantly, rupture propagates to the Totschunda fault, and abandons the Denali fault at the branch in the east, although the rupture propagates (with relatively low slip) for about 40 km. In our model, this propagation pattern is a direct result of the more favorable orientation of the Totschunda fault with respect to the discussed regional stress field, coupled with the dynamic effect that slip on the Totschunda fault sends the eastern Denali segment into a stress shadow (e.g. Kame et al. (2003); Bhat et al. (2002)) Another interesting feature of the dynamic model is that rupture propagation is discontinuous: Due to the more favorable orientation of the Totschunda fault near the western branch, stress waves from the Denali fault trigger slip on the Totschunda segment approximately 14 km ahead of the primary rupture front. Experiments with finer fault discretization in the branching region show similar jumping patterns, indicating that the discontinuous rupture is not an artifact of our low spatial resolution. Minor changes to the stress field produce slightly different jumping lengths, but the main effect remains. Interestingly, stress fields that do not produce jumping rupture also tend to allow significant slip on the eastern Denali segment, in conflict with observations.


Regional seismic waveforms, continuous and campaign-mode GPS data, and surface slip measurements were used to obtain a kinematic model of the rupture process of the November 3, 2002 $M_w$ 7.9 Denali, Alaska, earthquake. The event initiated as a $M_w$ 7.0 reverse slip event on the north-dipping Susitna Glacier fault with subsequent right-lateral slip distributed over approximately 300 km of the Denali fault system. Near-shear rupture velocity was inferred from direct measurement of S-wave arrival time and the kinematic modeling. The average and maximum slips were found to be 2.14 m and 9.94 m. Static stress drop varies from 1 to 5.2 MPa over the 5-segment fault model. The dynamic models point toward the possibility of discontinuous rupture propagation in this event. Such discontinuous rupture propagation has been implied in kinematic models of earlier events (e.g., the 1984 Morgan Hill event (Beroza and Spudich, 1988)), but in this case the dynamic source of the effect is identified in terms of the fault geometry rather than heterogeneity in strength on a planar fault. In addition to qualitatively improving the kinematic slip model for this event, discontinuous rupture propagation may help to explain some of the apparent rapid rupture propagation, and aid in the abandonment of the Denali segment, in agreement with observations.


We thank Jeff Freymueller and Sigrun Hreinsdóittir for providing their GPS solutions for both the continuous and campaign-mode sites, and Mark Murray for his help with GPS solutions at the initial stages of this study.


Beroza, G.C., and P. Spudich, Linearized inversion for fault rupture behavior: application to the 1984 Morgan Hill, California, earthquake, Bull. Seism. Soc. Am., 93, 6275-6296, 1988.

Bhat, H.S., R. Dmowska, J.R. Rice, and N. Kame, Testing theory for fault branching: Denali to Totschunda, Alaska, November 3, 2002, EOS Trans. AGU, 83, Fall Meet. Suppl., Abstract S72F-1364, 2002.

Dreger, D. S., D. Oglesby, R. Harris, N. Ratchkovski, and R. Hansen, Kinematic and dynamic rupture models of the November 3, 2002 $M_w$ 7.9 Denali, Alaska, earthquake, Geophys. Res. Lett., submitted, 2003.

Eberhart-Phillips, D., et al., The 2002 Denali fault earthquake, Alaska: A large magnitude, slip-partitioned event, Science, 300, pp. 1113-1118, 2003 DOI 10.1126/science.1082703.

Frankel, A. D., N. N. Biswas, A. H. Martirosyan, U. Dutta, and D. E. McNamara, Rupture Process of the M7.9 Denali Fault, Alaska, Earthquake Determined from Strong-Motion Recordings, EOS Trans. AGU, 83, Fall 2002 AGU meeting S72F-1340, 2002.

Hartzell, S. H., and T. H. Heaton, Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull. Seism. Soc. Am., 73, 1553-1583, 1983.

Hreinsdóittir, S., J. T. Freymueller, H. J. Fletcher, C. F. Larsen, and R. Bürgmann, Coseismic slip distribution of the 2002 MW7.9 Denali fault earthquake, Alaska, determined from GPS measurements, Geophys. Res. Lett., 30, 1670, doi:10.1029/2003GL017447, 2003.

Kame, N., J.R. Rice, and R. Dmowska, Effects of pre-stress state and rupture velocity on dynamic fault branching, J. Geophys. Res., 108, 2265, doi: 10.1029/2002JB002189, 2003.

Kaverina, A., D. Dreger, and E. Price, The combined inversion of seismic and geodetic data for the source process of the 16 October, 1999, $M_w$ 7.1 Hector Mine, California, earthquake, Bull. Seism. Soc. Am., 92, 1266-1280, 2002.

Okada, Y., Surface deformation due to shear and tensile faulting in a halfspace, Bull. Seism. Soc. Am., 75, 1135-1154, 1985.

Oglesby, D.D., Earthquake dynamics on dip-slip faults, Ph. D. thesis, University of California, Santa Barbara, 1999.

Olson, J.V., C.R. Wilson, and R.A. Hansen, Infrasound associated with the 3 November, 2002 Denali Fault, Alaska earthquake, Geophys. Res. Lett., submitted, 2003.

Ozakar, A.A., S.L. Beck, and D.H. Christensen, Source process of the 3 November 2002 Denali fault earthquake (central Alaska) from teleseismic observations, Geophys. Res. Lett., 30, 1638, doi:10.1029/2003GL017272, 2003.

Ratchkovski, N. A., and R. A. Hansen, New constraints on tectonics of interior Alaska: earthquake locations, source mechanisms, and stress regime, Bull. Seism. Soc. Am., 92, 998-1014, 2002.

Ratchkovski, N.A., A. Hansen, J.C. Stachnik, T. Cox, O. Fox, L. Rao, E. Clark, M. Lafevers, S. Estes, J.B. MacCormack, T. Williams, and K.R. Kore, Aftershock sequence of the MW 7.9 Denali, Alaska, earthquake of 3 November, 2002 from the regional seismic network data, Seism. Res. Lett., in press, 2003.

Whirley, R.G., and B.E. Engelmann, DYNA3D: A Nonlinear, Explicit, Three-Dimensional Finite Element Code for Solid and Structural Mechanics - User Manual, University of California, Lawrence Livermore National Laboratory, 1993.

Berkeley Seismological Laboratory
215 McCone Hall, UC Berkeley, Berkeley, CA 94720-4760
Questions or comments? Send e-mail:
© 2004, The Regents of the University of California