The Slip History of the 1994 Northridge, California, Earthquake Determined from Strong Ground Motion, Teleseismic, GPS, and Leveling Data


David J. Wald
U.S.G.S., 525 S. Wilson Ave, Pasadena CA, 91106
Thomas H. Heaton
Seismological Lab
Caltech, Pasadena, CA 91125
Kenneth W. Hudnut
U.S.G.S., 525 S. Wilson Ave, Pasadena CA, 91106

Bull. Seism. Soc. Am., 86, S49-S70

*(To View (or download) a file with the Slip Model Click Here)*


ABSTRACT

We present a rupture model of the Northridge earthquake, determined from the joint inversion of near-source strong ground motion recordings, P and SH teleseismic bodywaves, Global Positioning System (GPS) displacement vectors, and permanent uplift measured along leveling lines. The fault is defined to strike 122 degrees and dip 40 degrees to the south-southwest. The average rake vector is determined to be 101 degrees and average slip is 1.3 meters; the peak slip reaches about 3 meters. Our estimate of the seismic moment is 1.3 +- 0.2 x 1026 dyne-cm (potency of 0.4 cubic km). The rupture area is small relative to the overall aftershock dimensions and is approximately 15 km along strike, nearly 20 km in the dip direction, and there is no indication of slip shallower than about 5-6 km. The up-dip, strong-motion velocity waveforms are dominated by large S-wave pulses attributed to source directivity and are comprised of at least 2-3 distinct arrivals (a few seconds apart). Stations at southern azimuths indicate two main S-wave arrivals separated longer in time (about 4-5 sec). These observations are best modeled with a complex distribution of subevents: The initial S wave arrival comes from an asperity that begins at the hypocenter and extends up-dip and to the north where a second, larger subevent is centered (about 12 km away). The secondary S arrivals at southern azimuths are best fit with additional energy radiation from another high slip region at a depth of 19 km, 8 km west of the hypocenter. The resolving power of the individual data sets is examined by predicting the geodetic (GPS and leveling) displacements with the dislocation model determined from the waveform data, and vice-versa, and also by analyzing how well the teleseismic solution predicts the recorded strong motions. The general features of the geodetic displacements are not well predicted from the model determined independently from the strong-motion data; likewise, the slip model determined from geodetic data does not adequately reproduce the strong-motion characteristics. Whereas a particularly smooth slip pattern is sufficient to satisfy the geodetic data, the strong-motion and teleseismic data require a more heterogeneous slip distribution in order to reproduce the velocity amplitudes and frequency content. Although the teleseismic model can adequately reproduce the overall amplitude and frequency content of the strong-motion velocity recordings, it does a poor job of predicting the geodetic data. Consequently, a robust representation of the slip history and heterogeneity requires a combined analysis of these data sets.

INTRODUCTION

The January 17, 1994 Northridge (M_w=6.7) earthquake produced the largest ground motions ever recorded in an urban environment and caused the greatest damage in the United States since the great 1906 San Francisco earthquake (U.S.G.S and S.C.E.C., 1994). Peak acceleration and velocity values were among the largest ever recorded in any earthquake, and the large number of strong-motion recordings is unprecedented. Extensive portable instrument deployments following the mainshock for recording aftershocks will provide calibration data for constraining the regional velocity structure and, ultimately, for better understanding the mainshock strong motions. The extent of the damage and the abundance of recorded ground motions necessitate a systematic analysis of the source-rupture process in order to better understand the nature of the ground motions and resulting damage patterns.

Wald and Heaton (1994a) examined an earlier subset of the strong-motion data described herein and determined a rupture model based on an inversion of those data. In this analysis, we add strong-motion data, GPS and leveling-line displacements, and teleseismic waveforms to the earlier, near-source strong-motion data set, and use all, both separately and in unison, to invert for the spatial and temporal variations of slip on the fault. The models resulting from individual data sets were checked against the other data in a forward modeling sense to analyze their resolution and predictive capacity. The results presented here supersede the earlier work since we have included additional data. Unlike the earlier studies, we use different Green's functions for strong-motion stations at soil and rock sites, rather than using a single velocity structure, and use a layered earth model for the geodetic displacement calculations, as opposed to the homogeneous half-space approximation used in earlier work. Generally, the results of our earlier modeling are consistent with the improved and updated results presented here.

Our fault parameterization involves a variable-slip, finite-fault model that treats the diverse data sets in a self-consistent manner, allowing them to be inverted jointly or independently. By representing slip on the fault with numerous subfaults and slip on each subfault by the summation of many point sources over the subfault area, we can generate near-source static, strong-motion, or teleseismic synthetic Green's functions with identical fault rupture models.

There are several important advantages in combining the multiple data sets. First, neither the GPS nor the waveform stations uniformly cover the near-source region. As a combined data set, the spatial sampling is enhanced. Second, the range of frequencies covered is from DC to 1.0 Hz, allowing comparisons between slip models which sample only coseismic slip (waveform data) with those that include slip from early aftershocks as well as aseismic slip (geodetics). Finally, since the geodetically determined slip pattern is completely independent of the rupture timing, requiring the final slip in the waveform inversions to fit the static data provides an independent constraint on any a priori timing assumptions made in the waveform modeling. This has a great advantage over band-limited waveform studies alone, where there is commonly a trade-off between the rupture timing and the slip location.

Other recent studies have shown the benefit of combining geodetic and waveform data in source inversions. ( Wald and Heaton (1994b) found that the addition of the geodetic data to the strong-motion and teleseismic data in the analysis of the 1992 Landers earthquake added important constraints on the rupture evolution. Only with the addition of the geodetic data (Hudnut et al., 1994) could the temporal and spatial evolution be imaged with confidence. In a study of the historic data from the 1923 Kanto, Japan earthquake, ( Wald and Somerville (1995) constrained the slip on the subducting fault plane with the available geodetic data (Matsu'ura et al., 1980) and placed constraints on the rupture timing with teleseismic body-waveform data.

Unfortunately, the Northridge inversion is more limited than the Landers study for several reasons. First, with the Northridge earthquake, there is ambiguity in assuming the location and geometry of the fault rupture surface(s) at depth. Geometrical fault complexity at depth is difficult to interpret from the aftershock locations alone. In contrast, the multiple-fault segments of the Landers rupture were exposed at the ground surface. Although complex, they were clearly defined and known to be nearly vertical in the down-dip direction. Second, the greater average depth extent of the slip in the Northridge earthquake reduces the effective resolution of the geodetic observations, particularly for deeper slip. Conversely, the geodetic monuments for the Landers earthquake were often immediately adjacent to shallow, high-slip areas of the fault and were thus very sensitive to the location and amount of slip. Finally, the rupture dimensions and source duration are relatively short for the Northridge earthquake, so we are attempting to resolve shorter wavelength features with lower slip amplitudes than for both the Landers and Kanto earthquakes.

We first discuss the coherency and variations of the near-source recorded ground motions and display the ground motions in a map view, allowing the waveform and amplitude variations to be examined. Next, we invert the bandpassed (1-10 sec) velocity ground motions alone to determine the distribution of slip on the fault rupture plane. The teleseismic and geodetic data are then inverted separately in the same fashion, allowing a direct comparison between the waveform and static solutions. A combined inversion of all three data sets is then performed to find a dislocation model most compatible with all of the observations and the solution is discussed. Finally, the individual models are tested against the other independent data sets, and implications for the resolution of our analysis are discussed.

FAULT RUPTURE MODEL

Fault Parameterization

In order to model slip during the Northridge earthquake, we need to assume a fault geometry, so we chose a single fault plane that is most consistent with a broad range of observations. We use a strike of 122 degrees, compromising between the different solutions found from modeling teleseismic surface-waves (Harvard CMT) and body-waves (Thio and Kanamori, 1996) which indicate strikes near 130 degrees, and the first-motion mechanism (U.S.G.S. and S.C.E.C., 1994) which requires a strike between 100-110 degrees. Further, vertical cross-sections of the aftershock distribution ( Mori et al. , 1995) present the simplest planar structure when projected perpendicular to roughly a 120 degrees strike Figure 1. The fault plane dips 40 degrees and passes directly through the relatively simple, planar aftershock distribution (Fig. 1). We did not attempt to model the data with fault locations that violate the aftershock observations.

The depth to the top of our assumed fault is 5.0 km and the bottom fault depth is 20.4 km, giving a down-dip width of 24 km. The along-strike fault length is 18 km. We discretized the fault plane into a total of 196 subfaults, each 1.29 km wide and 1.71 km long down-dip, in order to represent variable slip along the fault. The fault parameterization and modeling procedure we employ is fully described by Hartzell and Heaton (1983) and is summarized only briefly below.

Synthetic Green's Functions

Each subfault's motion is obtained by summing the responses of 25 point sources uniformly distributed over each subfault. Each point source is lagged appropriately in time to include the travel-time difference due to the varying source-to-station positions and to simulate the propagation of the rupture front across each subfault. Thus, all subfaults separately include the correct effects of directivity. The complete point-source responses for the strong-motion synthetics and the geodetic static displacements are computed for a layered velocity model (Table 1) with the discrete-wavenumber, finite-element (DWFE) scheme (Olson et al., 1984) for frequencies up to 3.0 Hz. In practice, we calculate a master set of synthetics for 1 km increments in depth from 4.0 to 22.0 km and for ranges between 0 and 60 km, to allow for the closest and farthest possible subfault-station combinations. Then for each point-source station pair, the required response is derived by a linear interpolation of the closest Green's functions available in the master set. The linear interpolation of adjacent Green's functions is performed by aligning the waveforms according to their shear-wave travel times.

The source-region velocity model used to compute the strong-motion Green's functions at rock sites (Table 2) and all the static displacements (GPS and leveling), given in Table 1a, is modified from Langston (1978, Model C). We have added a thin (0.5 km), slower layer to Langston's model to better approximate elastic properties just beneath the strong-motion rock-site stations. Minor variations on this model have been used extensively (e.g., Dreger and Helmberger, 1990) for modeling many local and regional waveforms in Southern California. For soil-site strong-motion stations, we replace the top 0.3 km of the rock-site velocity model with slower P and S wave velocities as shown in Table 1b.

Source Time Function And Rupture Velocity

The dislocation time history for each subfault is represented by the integral of an isosceles triangle with a duration of 0.6 sec. Each subfault is also allowed to slip in any of three identical 0.6-sec time windows following the passage of the rupture front, with the initiations of each window separated by 0.4 sec. Since the windows overlap in time, they can provide a smooth overall slip history lasting up to 1.4 sec, if necessary. Examples of the resulting subfault displacement time functions are shown in a later section. With multiple time windows, we can approximate both spatial variations in slip duration and rupture velocity perturbations from the assumed uniform velocity.

The rupture velocity is assumed to be a constant 3.0 km/sec, or about 85% of the shear-wave velocity in the source region (Table 1). We iterated through a range of values from 2.7 to 3.3 km/sec but found that rupture velocities in the range of 2.8 to 3.0 km/sec provides the best fit to the strong-motion data. The faster rupture velocity (3.0 km/sec) however, provided a better match when the geodetic and waveform data are inverted jointly.

Rupture Initiation

Evidence from the strong-motion data indicates that the initial rupture was rather subdued, reminiscent of delayed initial growth of the Loma Prieta ( Wald et al., 1991) and Landers earthquakes (Abercrombie and Mori, 1994). However, there is no evidence for the unusual long-period radiation that was observed in the beginning of the Loma Prieta earthquake. That is, the beginning of the Northridge mainshock was indistinguishable from the beginning of aftershocks in the hypocentral region (Abercrombie, 1994) The strong-motion trigger times, when available, indicate that the triggering P wave arrived at least 0.5 seconds later than expected (Wald and Heaton, 1994a) given the predicted travel time from the hypocentral parameters determined from the high-gain short-period records from the Southern California Seismographic Network (SCSN). Further analysis by Ellsworth and Beroza (1994) suggests that a small nucleation phase of the rupture was followed by a secondary, larger rupture episode beginning near the hypocenter approximately 0.5 seconds later, consistent with the delayed strong-motion trigger times. We used the origin time (12:30:55.2 GMT), epicenter (34.211 North Latitude, 118.546 West Longitude), and hypocentral depth (17.5 km) determined by relocating SCSN network phase data (J. Mori, written communication, 1994). Based on the above observations, we initiate the rupture model 0.5 sec after the hypocentral time. We thus chose to ignore the foreshock or initial rupture and began modeling at the time of the first significant rupture episode. We assumed that the main (secondary) rupture began at or near the network hypocentral location and then allowed the rupture to propagate radially outward from that point.

Inversion Method

A constrained, damped, linear, least-squares inversion procedure is used to obtain the subfault dislocation weights that give the best fit to the velocity waveforms and/or geodetic displacements. The inversion is constrained by requiring that the slip is everywhere positive, and it is damped by minimizing the difference in dislocation values between adjacent subfaults. These constraints have been previously discussed in detail by Hartzell and Heaton (1983). Solving for the amplitude of slip on each subfault, given the strong-motion observations and subfault synthetic seismograms, is posed as an overdetermined system of linear equations.

The rake vector is allowed to vary within the range of 55 degrees (left-lateral thrust) to 145 degrees (right-lateral thrust), requiring two relative slip components for each subfault to be determined. Additionally, since the three-time window parameterization allows variable rake in each window, each subfault requires a total of six slip variables. Hence, for 196 subfaults, the total of unknown variables is 1176. In the case of the geodetic inversion alone, the number of data is smaller than the number of unknown slip parameters, so the inversion is formally underdetermined; however, the addition of the smoothing constraints make the inverse problem overdetermined.

STRONG-MOTION INVERSION

Strong-Motion Data and Preliminary Analysis

We use strong-motion accelerograms from the California Division of Mines and Geology (CDMG; Shakal et al., 1994), the U.S. Geological Survey (USGS; Porcella et al., 1994), the Los Angeles Department of Power and Water (LADWP), Southern California Edison (SCE), the University of California at Santa Barbara (UCSB), and the University of Southern California (USC). Table 2 lists the station abbreviations and locations as well as other site specifications, and the distribution of stations is displayed in Figure 2. Where two or more stations were located in close proximity to each other, we chose a representative location for our analysis.

The variability of the ground motions in the Northridge earthquake is examined in map view in Figure 2. Figure 2a shows the unfiltered acceleration recordings and Figure 2b. shows the unfiltered velocity waveforms. At each station, the component shown has been rotated to maximize the recorded peak ground velocity; the rotation angles and peak velocities are tabulated in Table 2. Since the integration from acceleration to velocity enhances lower frequencies, comparison of Figure 2a and 2b allows the spatial waveform and amplitude variations to be visualized as the frequency bandpass is shifted from higher frequencies to longer periods. Effectively, this allows us to separate some of the effects of wave propagation and site response (which are most profound in the accelerations) from the more obvious contributions of rupture directivity and radiation pattern, which dominate the longer-period velocity waveforms.

The variability of the observed ground motions can be attributed to a number of factors in addition to source distance. The ground velocities north of the epicenter are dominated by simple, large-amplitude pulses indicative of northward, up-dip source directivity. In the region 10-25 km north-northwest to north-northeast of the epicenter, where we would expect the combined effects of radiation pattern and directivity to be maximized for this fault geometry, large ground velocities are observed. In fact, the recorded peak horizontal ground velocity at the free-field site near the county hospital in Sylmar (stations SYL, 15 km north-northeast of the epicenter) was about 130 cm/sec; the peak velocity was over 180 cm/sec at the Los Department of Water and Power Rinaldi Receiving station (RRS) several kilometers to the south. The Rinaldi ground velocity is the largest recorded to date from any earthquake. Another impressive ground-motion recording is at station U56 (Fig. 2) which lies above the northwest up-dip corner of the modeled fault. Though the acceleration recording at U56 (Fig. 2) is not particularly large, the peak velocity is nearly 115 cm/sec, and the width of the velocity pulse is nearly 2.0 sec (Fig. 2).

These large, up-dip recorded ground velocities are significant, since peak ground velocity is a better measure of damage potential for large structures than is peak ground acceleration (EERI, 1994). It is important to note that much of the up-dip region, where directivity effects dominate, is not as densely urbanized as regions to the south (Fig. 3). For buildings of four or more stories, Figure 3 shows almost all were located outside the region that experienced the strongest ground velocities. Among the large structures located up-dip, however, were several notable freeway bridges, including two that collapsed at the Interstate 5, State Route 14 interchange and at the Interstate 5 Gavin Canyon undercrossing (EERI, 1994). Within the southern San Fernando valley, the ground velocities were much more moderate than to the north, even though this region is directly above the rupture surface. Few large man-made structures experienced the full force of the Northridge earthquake; of those that did, many were severely damaged.

The effect of directivity is less obvious in the peak acceleration data (Fig. 2). Although the ground velocities are clearly large north of the epicenter, the concentration of the largest accelerations is nearly above the fault plane. Further, several of the larger peak accelerations (e.g., station SMC) were located south of the epicenter where the large amplitudes were likely dominated by propagation and site effects rather than source radiation alone. As in other earthquakes, soft soils, basin structures, and topographic features may have produced higher ground motions locally for the Northridge earthquake (e.g., Graves, 1995; Spudich et al., 1996).

Modeling

When the trigger time was available, synthetic and observed waveforms were aligned in absolute time and only minor corrections were made for static station delays or timing errors. Inversion of only the strong-motion data from stations with absolute timing indicates a large slip at the hypocenter. Hence, for all other stations, the synthetic S waves from the subfault containing the hypocenter were aligned with the initial S wave in the data. However, a primary factor limiting the ultimate resolution of our modeling is the lack of absolute timing for many of the strong-motion recordings.

All ground velocity observations are scaled to a unit amplitude in the inversion in order to insure equal importance of smaller amplitude stations and to down-weight possible site effects. Examination of the ground-motion recordings shows that, at adjacent stations, more variability was found in the vertical components, suggesting more contamination from site and path effects in the vertical data. For this reason, the vertical components were down-weighted by a factor of two with respect to the horizontal components. With the exception of CAS and VSQZ (Fig. 2), all 30 stations used in the inversion have horizontal distances of less than 26 km from the center of the fault. We avoided more distant stations and those within the Los Angeles Basin since many of the aftershock recordings at these locations indicate waveform modifications caused by wave propagation through complex structures.

The accelerograms were bandpass filtered between 0.1 and 1.0 Hz with a zero-phase, third-order Butterworth filter and were then integrated to obtain ground velocity. This bandpass was chosen to avoid long-period integration noise and to avoid inadequacies of the theoretical Green's functions at higher frequencies. The use of velocity rather than acceleration waveforms further emphasizes longer-period characteristics of the strong motions. Our inability to adequately estimate strong-motion Green's functions at frequencies higher than 1.0 Hz is limited in part by our lack of knowledge of the crustal velocity structure, particularly at more distant stations and, again, by the lack of absolute time at many strong-motion sites. We modeled between 15 to 20 sec of the strong motions on each component, depending on the duration at individual stations. We did not rotate the stations to fault normal and parallel, since for stations in the near-source region, particularly above the fault, rotation is ambiguous. However, in order to facilitate waveform comparisons, all horizontal components are rotated to north and east, if not so recorded.

Results

Inversion of the strong-motion data alone results in the slip distribution displayed in Figure 4. The slip maximum is 273 cm and the total seismic moment is 1.10 x 1026 dyne-cm (Table 3). A comparison of the observed and synthetic strong motions is given in Figure 5. For each station waveform, all six traces (observed and synthetic north, east, and vertical components) are scaled to the peak value of the largest component (given in cm/sec). Most of the waveforms are well matched, both in amplitude and phase, by the synthetic ground motions.

The overall slip pattern is quite similar to the strong-motion model in our earlier article (Wald and Heaton, 1994a, Fig. 8). The overall slip is slightly smaller in this study, in part, because we have reduced the velocities in the near-surface layers at soil sites, increasing the impedance contrast and the resulting synthetic ground-motion amplitudes. The model is also in notable agreement with most of the features found in the slip pattern determined by Dreger (1994) using an empirical Green's function deconvolution of regional waveform data and with the strong-motion model of Zeng and Anderson (1996). All three strong-motion models require a substantial amount of slip west of, and at a comparable depth to, the hypocenter and another asperity up-dip and due north of the epicenter. As will be discussed later, however, the strong-motion model does not adequately fit the geodetic observations.

TELESEISMIC INVERSION

Teleseismic Data and Modeling

The 13 teleseismic station locations for the broadband data used in this study are listed in Table 4 and their azimuthal distribution with respect to the epicenter is shown in Figure 6. These stations provide a well distributed azimuthal coverage of the source. The instrument responses were deconvolved from the original recordings to obtain ground velocities (Fig. 7, left) and displacements (Fig. 7, right) and the data were low-pass filtered at 2 Hz. Since faulting was predominantly thrusting on a moderate dipping plane, all direct P waves are compressional. The initial portions of the P waveforms indicate at least 3 subevents.

For the teleseismic data, we modeled the first 25 sec of both the P and SH wavetrains. The synthetic arrival from the hypocenter was aligned with the observed velocity waveforms, since the initial arrival is more impulsive in velocity than in displacement. Further, since the source dimensions of the Northridge rupture are relatively small to be well resolved teleseismically, we inverted the higher-frequency velocity waveforms rather than the displacement waveforms in order to try to resolve finer scale rupture details.

Results

Inversion of the teleseismic data alone results in the slip distribution displayed in Figure 4b. The slip maximum is 312 cm and the total seismic moment is 1.31 x 1026 dyne-cm (Table 3). Although the total slip, maximum slip and the degree of slip heterogeneity are similar to those of the strong-motion model (Fig. 4b), overall, the teleseismic slip model is quite different. However, both data sets do share the requirement for significant slip at the hypocenter, but with the largest slips occurring up-dip and to the west. Furthermore, they both indicate another patch of deep slip about 8 km west of the hypocenter. A comparison of the observed and synthetic teleseismic waveforms is shown in Figure 7. Although we inverted the velocity data, we show comparisons of both the velocity and displacement data, since displacement is most often used in teleseismic waveform inversions (e.g., Thio and Kanamori, 1996). Note that large differences in the fits to the velocity data correspond to less dramatic differences in the match to the displacement data. Most of the details, including the initial subevent arrival in the P waves, are well modeled.

GEODETIC INVERSION

GPS Data

The GPS data we use consist of horizontal and vertical displacements at 46 monuments from the analysis of Hudnut et al. (1996). In their study, measurements of static ground displacements from a time period of two years prior to the earthquake and one month after the earthquake were determined from with GPS data. Coseismic displacement vectors were determined by correcting the displacements at each monument for secular variations by using the secular velocity estimates given by Feigl et al. (1993). A detailed description of the data reduction and processing plus a tabulated listing of the GPS station parameters is given by Hudnut et al. (1996). The station locations used in this study and the observed vertical uplift and horizontal displacement vectors for the closer stations are shown in Figure 8.

Leveling-Line Data

The leveling data used in this study are the result of a preliminary analysis of data collected before and after the Northridge earthquake. The 210 leveling-line station locations are shown in Figure 9. The routes include sections of Interstates 405 and 5, and State Routes 101, 118, and 126. The data are shown along profiles in Figure 10. The data were collected and adjusted by National Geodetic Survey (NGS) and California Division of Transportation (Caltrans). We took the pre- and post-earthquake adjustments and simply differenced the values at each benchmark. Any obvious outliers were then eliminated.

Pre-earthquake data were collected jointly by the NGS and Caltrans. All of the NGS data and part of the Caltrans data were adjusted using standard least-squares methods by the NGS in obtaining their North American Vertical Datum 1988 (NAVD88). The Caltrans pre-earthquake data that were not included in that adjustment were separately adjusted holding the elevations fixed at those few stations with NAVD88 heights, typically at endpoints of level-line segments. In other words, the Caltrans adjustments were done assuming line endpoint elevations (in NAVD88) to be correct. Post-earthquake data used here were collected by the NGS on contract to the USGS, in a project supported by the Federal Emergency Management Agency (FEMA). These data were collected during March-September, 1994.

A preliminary adjustment of these post-earthquake data, using standard least-squares methods, was performed by Caltrans and provided to the USGS, along with the pre- earthquake adjustments of the Caltrans data (J. Satalich, Caltrans, pers. comm., 1994). Because the data from the various leveling projects comprising the pre-earthquake data set are all adjusted to a common datum, yet the actual surveys span different times (1989-1993), some unestimated error resulting from differential monument instability may be included in these data. Also, benchmarks may have been disturbed by shaking during the earthquake, by contributions to the static displacement field from aftershocks, or by interseismic and post-seismic deformation. A more thorough analysis of these data will be carried out during 1995 by the USGS in collaboration with NGS and Caltrans, in part to refine these preliminary results and also to better analyze potential effects of non-tectonic and/or non-coseismic motions.

Modeling

The geodetic displacements are calculated for the same one-dimensional layered velocity structure used in computing the strong-motion waveforms at rock-site stations (Table 1). Since the DWFE Green's functions computed for the strong-motion modeling are complete waveforms, including near-field and static motions, we can take advantage of these calculations and simply preserve the permanent static displacements. The subfault point-source summation for the static calculation uses the same algorithm as the strong-motion Green's function, except that the synthetic time series are replaced with a single displacement value at each given station location.

This approach to computing the static displacements is advantageous in that we have a common velocity structure in the source region for the teleseismic, strong-motion, and geodetic models, so the seismic moment estimations from the different models are directly comparable. The earth structure model has more realistic attributes than a half-space approximation, including depth-dependent velocities and Poisson's ratios. In weighting the geodetic data, we weight each station proportional to its relative amplitude, unlike the strong-motion modeling where each waveform is normalized in the inversion. In this sense, monuments with larger amplitude displacements play a more important role in the geodetic inversion. We also down-weighted the leveling data by a factor of three relative to the GPS data to compensate for the large ratio of the number of leveling stations to GPS monuments.

Since we are interested in recovering the coseismic slip on the fault plane from the geodetic data, factors other that coseismic slip contributing to the movement of the geodetic monuments were examined. The dominant sources of contamination of the geodetic data are most likely 1) permanent displacement due to aftershocks, 2) local, non-tectonic movement of monuments due to intense shaking, and 3) fluid or gas withdrawal resulting in subsidence. A few leveling stations on each of the leveling lines indicate substantially smaller uplifts relative to their neighbors (Fig. 10), indicating local subsidence, either long-term or shaking induced. These sites were not weighted in the inversions. We note that all of the anomalous sites are located on alluvium (Fig. 9 and Fig. 10).

We investigated the nature of the monument movements due to aftershock activity by summing the computed displacement contributions at each monument for every aftershock with a moment magnitude of 4.0 or larger that occurred within the time span of the surveys. Smaller magnitude aftershocks do not contribute measurable signals. Aftershocks were represented by displacements at depth with the locations, seismic moments, and mechanisms determined by Thio and Kanamori (1996) from local surface-wave source inversions. One of the two largest aftershocks occurred one minute after the mainshock, and while its location and local magnitude (M_l=5.9) are known (Hauksson et al., 1996), the mechanism and depth are not. Since this event is potentially a large source of immediate post-seismic surface displacement (depending on the source depth) we included this event with the assumption that it occurred on the mainshock source plane with the mainshock mechanism. However, the lack of more definitive information on one of the two largest aftershock adds uncertainty to this calculation.

The total seismic moment of the aftershocks considered is approximately 2 x 1025 dyne-cm, or about 15% of the estimated mainshock moment. However, the mechanisms and locations of the aftershocks are diverse, so their cumulative contribution to the post-seismic deformation field may be less than if their slip contributions were all the same orientation as the mainshock and on the mainshock fault plane. The total displacements at the GPS stations due to the aftershocks are in excess of 5% of the total motion only at sites RESE, CHRN, and NEWH (Fig. 8), and their orientations with respect to the total displacements at not consistent at these sites. However, these postseismic displacements due to the aftershocks may add significantly to the total slip estimated by modeling the geodetic data.

Results

The dislocation model resulting from the geodetic inversion is displayed in Figure 4c. The total seismic moment is 1.42 x 1026 dyne-cm and the maximum slip is 256 cm (Table 3). Overall, the slip is smoother than that of the strong-motion model (Fig. 4a; both the slip amplitude and the rake angles are more uniform. The slip is confined to a single, central asperity that is consistent with the central lobes of slip in the strong-motion model. However, the hypocentral and deep asperities found in the waveform models are not as prominent. A comparison of the observed and predicted geodetic-based model displacements is given in Figure 8 and Figure 10. The overall pattern of horizontal displacement and vertical uplift is well matched.

Qualitatively, our geodetic slip model (Figure 4c) is similar to the variable slip models presented by Hudnut et al. (1996) and Shen et al. (1996); the location of the slip maximum and depth extent of the slip are comparable. However, our model fault plane is required to pass directly through the plane of aftershocks (Fig. 1), whereas Hudnut et al. (1996) and Shen et al. (1996) use a plane 3.5 km shallower than that defined by the aftershock distribution. Although those studies prefer the shallower plane because it can better fit the GPS data, we find that our modeling approach is both consistent with the aftershock locations and fits the GPS and the leveling data reasonably well.

COMBINED INVERSION

Modeling

We now present an inversion that combines the strong-motion, teleseismic, and geodetic data sets. Since the fault model parameterization and inversion remain unchanged for the variety of data sets, the combined inversion is a natural extension of the prior inversions. Although the number of unknowns remains fixed, the total number of data is increased. The main difficulty encountered was determining the relative weighting factors for each data set, so that each would be represented in the joint inversion. This was accomplished by perturbing the relative weights on a trial-and-error basis to insure that the fit to each data set was not strongly degraded.

Results

The dislocation model determined from inverting the combined geodetic and waveform data is shown in Figure 4d. The seismic moment of the combined solution is 1.40 x 1026 dyne-cm and the peak slip is 319 cm (Table 3). The solution is clearly a compromise between the slip patterns determined for the strong-motion and the GPS data independently (Fig. 4a and 4c, respectively); the teleseismic data plays a lesser role (Fig. 4)b. While the central asperity in the geodetic slip model remains, so does the deep slip 12 km west of the hypocenter present in the strong-motion and teleseismic model. Likewise, the slip at the hypocenter, required to fit the waveform data, is preserved (i.e., the deep slip patches are neither resolved nor rejected by the geodetic data).

For the combined model, the comparisons of the observed and predicted leveling data, GPS data, strong-motion data, and teleseismic data are shown in Figures 11, 12, 13, and 14, respectively.

The slip contours for the combined model are also shown projected on a map view in Figure 1. As has been noted with other earthquakes, there is a slight tendency for the aftershocks to concentrate at the edges of high slip concentrations, indicating a redistribution of stresses in these areas (e.g., Mendoza and Hartzell, 1988).

It is often difficult to estimate stress drop for earthquakes since one must normally make assumptions concerning the relationship of the known rupture duration with the unknown rupture area. However, our finite fault modeling allows us the advantage of determining both the amount of slip and the area over which it occurred. Even so, the stress drop calculation is only approximate; since it is difficult to determine where the slip goes to zero, defining the boundaries of the rupture area is ambiguous. For our slip model, the stress drop expression of Eshelby (1957) for a circular fault is appropriate, Delta sigma = ( 7 pi mu overline{u} )/(16a), where mu is the rigidity, overline{u} is the average dislocation, and a is the radius. Using mu = 3.6 x 1011 dyne-cm^2, overline{u} = 140 cm, and a = 9.4 km, we obtain a stress drop of about 74 bars.

DISCUSSION

Analysis of Combined model

The combined dislocation model is used for further discussion of the source process since it provides the best overall fit to the available data. First, we view the variations of slip duration and image the rupture propagation. We then examine the source complexity and describe its contribution to the recorded ground velocity waveforms.

The displacement time histories for selected subfaults are shown in Figure 15. The displacement plots correspond to subfaults A, B, and C labeled in the cross-section of the combined rupture model shown. The rise times at these locations are 1.4, 1.0, and 0.6 sec, respectively. For these regions of relatively high slip, the particle velocity---assuming symmetric motion on both sides of the fault---is nearly 1 m/sec for region A and B, but it is closer to 2 m/sec near the hypocenter. For most regions of lower slip values, the rise times are 0.6 sec, with corresponding particle velocities of less than 1 m/sec. These velocities agree well the average values tabulated by Heaton (1990) for a number of earthquakes for both the average particle velocity over the entire fault, and the value for asperities alone.

Figure 16 shows the slipping portion of the fault and the amount of slip during 1-sec "slices" in time; hence, the images depict the slip velocity. Local slip durations are less than about 1.5 sec, whereas the total rupture duration is about 7 sec. Figure 16 indicates that only a portion of the rupture surface is slipping at one time (Heaton, 1990). We performed tests fixing longer duration time functions, but we could not match the velocity waveforms as well.

The general pattern of the strong-motion duration and waveform complexity can be partially explained by the relative position of individual stations with respect to the regions of concentrated slip shown in Figure 4d. Those stations up-dip and to the north, show simple, large amplitude, short duration shear-wave pulses since the rupture propagated towards them (Figure 2b). At other azimuths the shear-wave arrivals are both separated farther in time, and smaller in amplitude. At those locations the source complexity is more apparent, and separate subevents can be identified.

Figure 17 shows how four regions of the fault rupture contribute to the makeup of the synthetic waveforms. We considered the contributions to the ground motion from the shaded regions above each column of Figure 17. with the fourth region comprised of the hypocentral asperity. In the column below each shaded fault depiction we show the observed velocity components with the corresponding ground motion contribution from that fault area. The last column represents the complete simulation from the entire fault.

Column one in Figure 17 shows that the slip concentration that dominates the geodetic slip model (Fig. 4c) also contributes the largest ground motions in the Northridge earthquake. The large amplitudes at up-dip station PAR (as well as JFP, SYL, NHL) are predominantly from the integrated effect of the central asperity. The fourth column indicates that the slip in the hypocentral region contributes the impulsive, initial shear-wave arrivals at most stations. Although this slip is absent from the geodetic model Fig. 4c), it is clearly required to reproduce the initial portion of the strong-motion observations, and the geodetic data permit it.

The deepest slip concentration (Fig. 17, second column) plays an important role in providing the later arrivals at southern and western azimuths, here shown with stations U53, TPG, and WOO. This is also apparent at stations SCC and SCR (Fig. 2b), where the large secondary arrivals come from this deep asperity. It is not clear whether the deep slip represents the down-dip extension of the mainshock rupture plane or a separate fault rupture.

Curiously, as shown by ( Wald and Heaton (1994b), Figures 5 and 15), the large (anomalous) peak acceleration at the Santa Monica site is on the east component even though the source radiation from our modeling predicts a substantially larger amplitude on the north component. This is also true of the Tarzana accelerogram, where our source model predicts a larger north component, and indeed the nearby ENR recording and our simulations have much larger north components, yet the east component dominates the Tarzana accelerations (Spudich et al., 1996). This suggests that either local propagational complexities dominate the Tarzana and Santa Monica recordings, or that there is complex source radiation not accounted for in our study; perhaps this isolated deep patch of rupture has a significantly different focal mechanism from the main rupture fault plane.

Further analysis of Figure 17 also indicates that at near-source stations, the nearest asperity often controls the character of the ground motion. For example, the contribution from the shallow, northeast region (Fig. 17, third column) dominates the arrivals at station VNY, just above and east of this fault area. This can be attributed to both the additional distance from the further patches of concentrated slip and the favorable source radiation pattern for stations above the rupture surface.

Predictive Capacity of Independent Models.

Our general philosophy is that no source model is sufficient if it clearly violates independent observations. Hence, it is informative to examine the ability of each model, derived from the separate data sets, to predict the remaining independent observations. Further, when near-source ground motions and geodetic displacements are lacking, source slip patterns and ground-motion estimations are often made using independent teleseismic observations alone, so it is useful to test this approach in a case like the Northridge earthquake, where the data are sufficient to do so.

Figure 18 shows the prediction of the leveling data from the slip model derived from strong motions alone. Some of the general features of the geodetic displacements are not well predicted by the strong-motion model. The uplift on the western portion of Route 118 line is a consequence of the western-most lobe of slip in the strong-motion slip model (Fig. 4a) that is not apparent in the geodetic model (Fig. 4c). Likewise, the under-prediction of the leveling data near the peak uplift on the Route 118 and Interstate 5 lines is due to the reduced slip in the central portion of the fault in the strong-motion model compared with the geodetic solution. We do not show the comparison of the GPS data and predictions since it is difficult to see the systematic misfits with so few stations, but the GPS data misfit by the strong motion model is comparable to the corresponding leveling misfit. Clearly, the strong-motion model alone cannot adequately reproduce the observed static displacements.

Further evidence for limitations of the strong-motion data comes from a comparison of the observed and synthetic strong-motion records due the strong-motion model alone (Fig. 5) with the fits of the combined model (Fig. 19). Although these models have have some substantial differences in the slip patterns, only a small degradation to the waveform matches is apparent. This implies a limited resolution of the slip details at this scale for the strong-motion data alone.

Similarly, the geodetic model does a poor job at predicting the recorded strong-motion velocities (Fig. 20). In order to predict the strong-motion data from the geodetic model, we assumed a rupture velocity of 3 km/sec and a rise time of 1.0 sec to be consistent with findings from other earthquakes. The large simple slip concentration that makes up the geodetic model produces amplitudes at up-dip strong-motion stations (e.g, NHL, JFP, PAR) that are too large (by a factor of two or three). Furthermore, the pulse widths of the large velocity arrivals are substantially longer period than those observed and the waveforms are simpler than the observations. In a predictive sense, and in terms of earthquake engineering concerns, the spectral response of these synthetics are a poor representation of the actual recordings. Whereas a particularly smooth slip pattern is sufficient to satisfy the geodetic data, the strong-motion data require a more heterogeneous slip distribution in order to reproduce the velocity amplitudes, frequency content, and waveform complexity.

Finally, the large initial arrivals at many stations (e.g., SHR, SCR, and ENR) attributed to slip near the hypocenter in the strong-motion model (Fig. 17) are not generated by the geodetic slip model. Since the ability of the geodetic data to resolve details of the slip pattern diminishes greatly with depth, some of the deep slip variations, which explain portions of the ground-motion waveforms, are not imaged. Considerably different details in the slip distribution at depth (compare Figs. 4a and 4d) give rise to only small differences to the geodetic predictions (compare Fig. 8 and Fig. 12), also indicating that the geodetic data alone has limited resolution of the deeper slip patches that are required by the waveform data.

Finally, we compare the observed strong motions with those predicted by the teleseismically-derived dislocation model (Figs 4b) in Figure 20. Clearly, the overall slip pattern and degree of slip heterogeneity of the teleseismic model was (with some exceptions) sufficient to produce reasonable estimations of the amplitude and frequency content of observed strong motions. This suggests that estimating near-source strong motions based on teleseismic modeling alone can be a useful endeavor.

CONCLUSIONS

A summary of the preferred fault parameters determined from our study is presented in Table 5. The rupture began at a depth of 17.5 km and propagated predominantly in the up-dip direction approximately along the direction of the average rake vector of 101 degrees. Slip terminated at a depth of about 6 km. Rupture occurred over an area approximately 15 km along strike (west-northwest from the hypocenter) and nearly 20 km up-dip. Our estimate of the seismic moment is 1.3 +- 0.2 x 1026 dyne-cm (potency of 0.4 cubic km) with an average slip about 1.3 meters over the rupture area, yielding an average stress-drop of about 74 bars.

The peak slip value is about 3 meters. The rupture velocity is estimated to be 3.0 km/sec, though slightly slower rupture velocities give comparable solutions. The rise time is best approximated with durations less than or equal to about 1.4 sec.

The up-dip, near-source strong-motion velocity waveforms are dominated by large S-wave pulses (Fig. 2b) attributed to source directivity and consist of at least 2-3 distinct arrivals (a few sec apart). Stations at southern azimuths (e.g., SCC, SCR, an SHR) indicate two main S-wave arrivals further separated in time (about 4-5 sec). These observations are best modeled with a complex distribution of subevents. The initial S wave arrival comes from an asperity that begins at the hypocenter and extends up-dip. The largest subevent is centered (about 12 km away) up-dip to the north. Secondary S arrivals at southern azimuths are best fit with additional energy radiation from another high-slip region near a depth of 19 km, 8 km west-northwest of the hypocenter. The correspondence to secondary arrivals observed at more distant stations to the south (e.g., Santa Monica) is more tenuous, since several of the aftershocks recorded there also indicate later arrivals as well, but, clearly, a secondary source contribution is expected at the time of the peak acceleration (Fig. 2a) based on our model of the near-source stations.

Both the observations and simulations indicate that the strongest long-period (1 to 3 sec) ground motions occurred up-dip from the rupture surface where source directivity is greatest. We note that much of this same region is sparsely populated and has few, if any, larger steel- or concrete-frame structures. Consequently, the engineering problems discovered in these structures after the earthquake occurred at relatively modest levels of ground motion relative to the ground motions experienced north of the epicenter.

After comparing the observed and predicted strong motions and geodetic displacements generated from the separate and combined slip models, we conclude that, given these specific data and source dimensions of the Northridge earthquake, neither the strong-motion nor the geodetic data alone have the resolving power to adequately recover the detailed source slip heterogeneity, at least well enough to independently predict the other observations. Further, even though the overall slip and degree of heterogeneity derived from the teleseismic data was sufficient to predict the character of the observed strong-motion velocities, it could not adequately constrain the slip pattern in detail or reproduce the static displacements. Consequently, an adequate representation of the source requires the combined analysis of these data sets.

ACKNOWLEDGMENTS

G. Hawkins of Southern California Edison, M. Trifunac of the University of Southern California, J. Steidl of the University of California at Santa Barbara, and R Tognazinni, C. Davis, and P. Lahr of the Los Angeles Department of Power and Water generously provided information and data from their strong-motion stations. R. Stein of the U.S.G.S. graciously allowed use of the post-earthquake leveling data collected under his direction, with funding support from FEMA. J. Satalich of Caltrans provided all of the preliminary leveling data we required for this analysis. Figure 3 was provided by N. Blaise and Kevin Miller of FEMA. Discussions with D. Dreger, R. Graves, N. Smith, and H. Thio were helpful. Reviews by P. Spudich and D. Jackson improved the original manuscript.

REFERENCES

Abercrombie, R. (1994). Stress Drops of the Northridge mainshock and aftershocks, Program for Northridge abstracts, Seismo. Soc. of Am. Annual Meeting, Pasadena.

Abercrombie, R. E. and J. Mori (1994). Local observations of the onset of a large earthquake: 28 June 1992 Landers, California, Bull. Seism. Soc. Am., 84, 725-734.

Dreger, D. S. (1994). Empirical Green's function study of the January 17, 1994 Northridge Mainshock (Mw 6.7), Geophys. Res. Lett., 21, 2633-2636.

Dreger, D. S. and D. V. Helmberger (1990). Broadband modeling of local earthquakes, Bull. Seism. Soc. Am., 80, 1162-1179.

Earthquake Engineering Research Institute (1992). EERI Special Earthquake report, August 1992, EERI Newsletter, 26, 1-12.

Ellsworth, W. L. and G. Beroza (1994). Nucleation and intial growth of the Northridge, California, earthquake, Program for Northridge abstracts, Seismo. Soc. of Am. Annual Meeting, Pasadena.

Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion and related problems, Proc. Roy. Soc. London, Series A, 241, 376-396.

Fiegl, K., D. Agnew, Y. Bock, D. Dong, A. Donnellan, B. Hager, T. Herring, D. Jackson, T. Jordan, R. King, S. Larsen, K. Larson, M. Murray, Z. Shen, and F. Webb (1993). Space geodetic measurements of the velocity field of central and southern California, 1984-1992, J. Geophys. Res. , 98, 21,677-21,712.

Graves, R. W. (1995). Preliminary analsis of long-period basin response in the Los Angeles region from the 1994 Northridge earthquake, Geophys. Res. Lett., 22, 101-104.

Hartzell, S. H. and T. H. Heaton (1983). 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.

Hauksson, E., L. M. Jones, and K. Hutton (1996). The 1994 Northridge earthquake sequence in California: seismological and tectonic aspects, J. Geophys. Res. , in press.

Heaton, T. H. (1990). Evidence for and implications of self-healing pulses of slip in earthquake rupture, Phys. Earth Planet. Inter., 64, 1-20.

Hudnut, K. W., Y. Bock, M. Cline, P. Fang, J. Freymueller, K. Gross, D. Jackson, S. Larson, M. Lisowski, Z. Shen, and J. Svarc (1994). Coseismic displacements in the Landers sequence, Bull. Seism. Soc. Am., 84, 625-645.

Hudnut, K. W., Z. Shen, M. Murray, S. McClusky, R. King, T. Herring, B. Hager, Y. Feng, P. Fang, A. Donnellan, and Y. Bock, (1996). Co-seismic displacements of the 1994 Northridge, California, earthquake, Bull. Seism. Soc. Am., 86, this issue.

Langston, C. A. (1978). The February 9, 1971 San Fernando earthquake: a study of source finiteness in teleseismic bodywaves, Bull. Seism. Soc. Am., 68, 1-29.

Matsu'ura, M., T. Iwasaki, Y. Suzuki, and R. Sato (1980). Statical and dynamical study on the faulting mechanism of the 1923 Kanto earthquake, J. Phys. Earth, 28, 119-143.

Mendoza, C. and S. H. Hartzell (1988). Aftershock patterns and main shock faulting, Bull. Seism. Soc. Am., 78, 1438-1449.

Mori, J., D. J. Wald and R. L. Wesson (1995). Overlapping fault planes of the 1971 San Fernando and 1994 Northridge, California earthquakes, Geophys. Res. Lett., 22, 1033-1036.

Olson A. H., J. Orcutt and G. Frazier (1984). The discrete wavenumber/finite element method for synthetic seismograms, Geophys. J. R. Astr. Soc. 77, 421-460.

Porcella, R. L., E. C. Etheredge, R. P. Maley, and A. V. Acosta (1994). Accelerograms recorded at USGS national strong-motion network stations during the $M_S = 6.6$ Northridge, California earthquake of January 17, 1994, U.S.G.S. Open-file Report 94-141, 100 pp.

Shakal, A., M. Huang, R. Darragh, T. Cao, R. Sherburne, P. Malhotra, C. Cramer, R. Sydnor, V. Graizer, G. Maldonado, C. Peterson, and J. Wampole (1994). CSMIP strong motion records from the Northridge, California earthquake of 17 January, 1994. Report No. OSMS 94-07, California Strong Motion Instrumentation Program, 308 pp.

Shen, Z., B. X. Ge, D. D. Jackson, D. Potter, M. Cline, and L. Sung (1996). Northridge earthquake rupture models based on the Global Positioning Systems measurements, Bull. Seism. Soc. Am., 86, this volume.

Spudich, P., M. Hellweg, and W. K. Lee (1996) Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northridge, California, earthquake: implications for mainshock motions, Bull. Seism. Soc. Am., 86, this issue.

Thio, H. K. and H. Kanamori (1996). Source complexity of the 1994 Northridge earthquake and its relation to aftershock mechanisms, Bull. Seism. Soc. Am., 86, this issue.

U.S. Geological Survey and the Southern California Earthquake Center (1994). The magnitude 6.7 Northridge, California, earthquake of January 17, 1994, Science, 266, 389-397.

Wald, D. J. and T. H. Heaton (1994a). A dislocation model of the 1994 Northridge, California, earthquake determined from strong ground motions, U.S.G.S. Open-file Report 94-278, 52 pp.

Wald, D. J. and T. H. Heaton (1994b). Spatial and temporal distribution of slip for the 1992 Landers, California, earthquake, Bull. Seism. Soc. Am., 84, 668-691.

Wald, D. J., T. H. Heaton, and D. V. Helmberger (1991). Rupture model of the 1989 Loma Prieta earthquake from the inversion of strong motion and broadband teleseismic data, Bull. Seism. Soc. Am., 81, 1540-1572.

Wald, D. J. and P. G. Somerville (1995). Variable-slip rupture model of the great 1923 Kanto, Japan earthquake: geodetic and body-waveform analysis, Bull. Seism. Soc. Am., 85, 159-177.

Zeng, Y., and J. G. Anderson (1996). A composite source modeling of the 1994 Northridge earthquake using genetic algorithm, Bull. Seism. Soc. Am., 86, this issue.


Return to Dave Wald's Homepage