Submitted to the Bulletin of Seismological Society of America, May, 1996; revised Sept. 1997.
The motivation behind resorting to the inversion scheme to be dealt with in this paper comes from the absence of even a one-dimensional local velocity model of the area, and also of enough geologic information to make one. Eberhart-Phillips (1993) talks about two possible approaches for setting up an inversion. The first approach is to start with a velocity model based on previous geologic interpretations and try to fit the data with the model by perturbing the latter as warranted by the data. The second approach is to start with a simple 1-D model obtained with the same data and let the inversion solution add complexity wherever needed by the data. Both the approaches have the common implicit assumption of the starting model being close to the true model as required by any linear or linearized inversion process. Eberhart-Phillips (1990, 1993) prefers using a simple minimum 1-D model derived by the program VELEST (Ellsworth 1977; Kissling 1988; Ellsworth et al., 1990) from the same data set as would be used for 3-D inversion. Kissling et al. (1994) prescribe the minimum 1-D model as the ideal initial reference model for local earthquake tomography in order to minimize the artifacts caused by an inappropriate initial model. But the "recipe" to calculate the said minimum 1-D models includes the establishment of an a prioi 1-D model from information regarding the stratification of the area as the first step (Kissling et al., 1994). But in a hitherto little studied place like Eureka Valley none of the above approaches seems adequate simply because of the absence of reliable a priori geological and geophysical information. Hence, the necessity of a method that does not depend on the initial model. Two series of methods satisfy this condition. They are iterative strategies requiring matrix manipulations (Spakman, 1993) and global optimization schemes such as simulated annealing (Rothman, 1985, 1986; Sen and Stoffa, 1991) which do not require any matrix manipulation. The sparse station distribution in the Eureka Valley area precludes the effective use of any iterative strategy. Optimization by simulated annealing is appealing in this special case because of its independence of the initial model and ability to produce solutions with very sparse data. So our method of choice for investigating structure and aftershock locations builds on a first step of initial 3-D model estimation using nonlinear optimization by simulated annealing (Pullammanappallil and Louie, 1993), followed by a fine-tuning step using 3-D linearized inversion (Thurber, 1983; Eberhart-Phillips, 1990). We also run a 3-D linearized inversion using a minimum 1-D model calculated by VELEST and compare the two results.
Fig. 1: Maps showing the location of Eureka Valley and the major geological structures
of the region, permanent and portable seismographic stations (triangles), and
preliminary event locations (diamonds). In the bottom, smooth curves are Bouguer
anomaly contours and rough curves are topographic contours. G+, G-, T+ and T-
are gravity high and low, and topographic high and low respectively.
Figure 1 shows the locations of permanent and temporary stations, along with preliminary UNR catalog locations of the events. At first glance it becomes evident that the spatial distribution of the permanent stations around the epicentral area offers only limited azimuthal coverage. To locate the aftershocks and find the local velocity structure of the epicentral area we used P and S arrivals at 20 permanent stations timed by the UNR Seismological Laboratory and at 8 portable stations that we timed ourselves.
The ratio between the number of S and P observations is several times greater for the portable stations than for the relatively more distant permanent stations. Figure 2 shows a time-distance plot for P and S arrivals from 430 events at portable and permanent stations at maximum 100 km epicentral distance, based on the preliminary estimates of origin times and hypocentral coordinates obtained at the UNR Seismological Laboratory from a one-dimensional regional velocity model and the FASTHYPO program (written by Robert Hermann of Saint Louis University). The figure excludes obvious outliers by keeping only those picks that lie within a window centered around the average straight line fit to the cloud of data points. The arrival/travel times have weights of 0 to 4 assigned depending on quality of pick, 0 being the best and 4 being the least reliable.
Fig. 2: Travel time versus preliminary epicentral distance plot of compressional first arrivals (P) and S arrivals used in the study. The upper plot has about
9000 P and S arrival times from 430 events to 28 permanent and portable stations. The lower plot has 1492 P and 768 S arrivals from the 430 events to 8 portable stations.
For comparison, we also run an exclusively linearized inversion, where the stage 1 in the above approach is replaced by a step of finding the minimum 1-D model using VELEST (Kissling et al., 1994). A brief description of the three techniques involved in the two parallel approaches follows.
We use an iterative Monte-Carlo based optimization process to solve the inverse problem. Each iteration comprises the following steps: 1) travel time computation through the current model and determination of the least-square error; 2) simultaneous random perturbation of origin time, hypocenter location, and velocity structure, and computation of new least square error; 3) acceptance of the new model if the corresponding new error is less than that of the previous iteration, plus a provision for probability-based conditional acceptance of new models if the new error is larger; 4) repetition of the previous steps until the annealing converges, where the difference in the least square error between successive models and probability of accepting new models become very small. An important component of this study is the random forcing of hypocenter location perturbations at each iteration. As we will show below, this forcing appears to allow the optimization of a reliable 3D-velocity model while relegating unacceptable errors in the final, perturbed hypocenter coordinates.
Model parameterization in this method assumes a continuous velocity field. The earth structure is represented in three dimensions by velocities defined at discrete points, and velocity at any intervening point is determined by linear interpolation among the surrounding eight grid points. Values at the velocity nodes are systematically perturbed during inversion. Except the outermost nodes, which are always fixed, every node can either be kept fixed or included in the inversion. Derivative weight sum (DWS) is a useful measure of ray density in the neighborhood of a node (Toomey and Foulger, 1989) and is used to determine which of the nodes should be included in the inversion. DWS for a node is similar to the ray hit count, but weighted by the ray-node separation and raypath length in the vicinity of the node.
The ART algorithm selects the path with the least travel time from a suite of circular arcs connecting the source and receiver. The iterative PB method fine tunes the ray path obtained by ART to approximate better the true ray path dictated by local velocity gradients (Um and Thurber, 1987). We adopt Eberhart-Phillips' (1990) method of adjusting an initial ray path rather than Um and Thurber's (1987) method of calculating a new midpoint for each segment because the former allows much easier updating of the raypath each time a source-receiver path has to be calculated.
A circular arc is a poor approximation at large distances and overestimates the path length, because a refracted wave is likely to be the first arrival at such distances. For an area like Eureka Valley, where alluvium is apparently underlain by rigid rocks, refracted waves can start arriving at even 10 km distance. Because of this, and the fact that timing error usually increases with distance, we chose a distance-weighting scheme of weight 1 for source-receiver distances of up to 30 km, and 0 for 80 km and beyond, with a linear tapering in between. The iterative pseudo-bending scheme (Eberhart-Phillips, 1990) reduces the difference between exact and ART raypaths for large source-receiver distances.
Velocity model parameterizations for the linear inversion and nonlinear optimization schemes somewhat differed from each other. In the nonlinear optimization, we used a 35 (x) by 35 (y) by 20 (z) grid of 1 km blocks. The linearized inversion used a 13 (x) by 15 (y) by 9 (z) grid of nodes with a node spacing of 2 km at the center and near surface portion of the model and increasing outwards and in depth. A minimum DWS of 7 is used as a condition for any node to be included in the inversion. We plot all our results in the central 30 by 30 by 20 km volume by linearly interpolating wherever necessary.
We relocated the hypocenters of the 430 events and modeled the velocity structure simultaneously using four different processes (see Figure 3):
In view of the widely observed phenomenon that earthquakes in the Basin and Range Province occur almost exclusively in the upper 15 km of the crust (Eaton, 1982; Vetter and Ryall, 1983; dePolo et al., 1992), we used the constraint of maximum hypocentral depth of 20 km for all four of the above processes.
Fig. 3: Chart describing the 4 different processes of linearized inversion and nonlinear optimization for hypocenter relocation and velocity model estimation.
Fig. 4: Mapview of relocated hypocenters (diamonds). Size of the diamonds proportional to the magnitude (M=1.0 to M=6.1) (a) obtained simultaneously with the minimum 1-D model by VELEST (process 1). (b) obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with minimum 1-D model as the starting model (process 2).
Fig. 5: Mapview of relocated hypocenters (diamonds). (a) obtained by simultaneous relocation and 3-D velocity estimation by simulated annealing (process 3) (b) obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with starting model estimated by simulated annealing (process 4).
The most consistent spatial distribution is the north-northwest trend of epicenters, well-defined in both the cases of simultaneous relocation and three-dimensional velocity estimation by linearized inversion (Figures 4(b) and 5(b)). Cross-sections along AB in Figures 6 and 7 reflect the same consistency, though the cross-section for the nonlinear optimization case (Figure 7(a)) shows the locations much more scattered in comparison with the other three cases. The main event has been located at depths of 11, 11.5 and 10.5 km by processes 1, 2 and 4 respectively (Figures 6(a), 6(b) and 7(b)). Process 3, whose input data consisted of arrivals to portable stations only, does not have a mainshock relocation. Deducing a unique fault plane through the hypocentral scatter in any of the cross sections is not straightforward. Relocated hypocenter locations from process 3 are much more scattered than those from the other 3 processes. Average rms residuals of location in processes 3 and 4 are 0.10 and 0.12 respectively. The mainshock has an rms residual of 0.05 in both the processes.
Fig. 6: Section of relocated hypocenters (diamonds) along AB (see Figure 4) with the inferred fault shown as a thick gray line through the hypocenters. (a) obtained simultaneously with the minimum 1-D model by VELEST (process 1). (b) obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with minimum 1-D model as the starting model (process 2).
Fig. 7: Section of relocated hypocenters (diamonds) along AB (see Figure 5) with the inferred fault shown as a thick gray line through the hypocenters. (a) obtained by simultaneous relocation and 3-D velocity estimation by simulated annealing (process 3) (b) obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with starting model estimated by simulated annealing (process 4).
We infer two different planes of roughly same strike (165 degrees) from the scatter of points in the cross sections of Figures 6(a), 6(b) and 7(b). One plane, which does not contain the relocated mainshock, dips steeply (about 75 degrees) eastwards; and the other, including the relocated mainshock, dips slightly less steeply (about 60 degrees) westwards. Loper et al. (1993) report a normal mechanism for the mainshock obtained from regional surface wave inversion. By incorporating broadband data from several institutions to the data recorded by Berkeley Digital Seismic Network, they relocated the large aftershocks relative to the mainshock using a master event technique and suggest a relatively complex rupture history for the event based on examination of mainshock and aftershock records, and identify at least two subevents. Massonet and Feigl (1995) model both east and west dips with normal mechanism in their satellite synthetic aperture radar (SAR) interferometric work, though they rule out the eastward dip as a local minimum in their modeling scheme. Their best-fitting focal mechanism has a north-northwestward strike and westward dip. Another SAR study by Peltzer and Rosen (1995) confirms a west-dipping fault with a slightly NNE strike. They reconcile the observations of north-northwest alignment of epicenters and north-northeast strike of mapped local faults to infer that the earthquake occurred on a west-dipping north-northeast-striking fault plane with the rupture propagating diagonally upward and southward. The Caltech TERRAscope moment tensor (MT) and Harvard centroid-moment tensor (CMT) solutions give gentler westward dips. Table 1 summarizes the focal mechanism parameters determined by the different researchers along with those inferred in this study (using standard sign convention of strike and dip):
Method Strike Dip Rake Reference (deg.) (deg.) (deg.) ------ ------ ------ ------ --------- TERRAscope MT -173 48 west - (Caltech website) Harvard CMT 210 30 west -93 (Dziewonski et al., 1994) SAR No. 1 173 54 west - (Massonet and Feigl, 1995) SAR No. 2 187 50 west - (Peltzer and Rosen, 1995) This study 165 60 west -
We see that the strike and dip obtained in our study match fairly well with the SAR results of Massonet and Feigl (1995). It is also not very different from the results of the other SAR study and TERRAscope except that their faults strike a little east of the north whereas ours does it north-northwestwards. Our results differ considerably from the Harvard CMT solution, the latter giving a 45 degree more eastward strike and dip only half as steep. Finally, the rather complicated distribution of hypocenters on the east-west cross-section (Figures 6(a), 6(b) and 7(b)), along with the somewhat curved distribution of epicenters (Figures 4(a), 4(b) and 5(b)), can be the expression of a concave fault involvement and the wide range of values and directions of the strike and dip of the fault found by different methods may also be due to complications brought into play by the concave shape of the fault. A complex rupture history and the existence of two subevents as suggested by Loper et al. (1993) mentioned above can also be the reason behind the observed spatial configuration of the relocated hypocenters.
Layer Velocity (km/s) Starting Depth (km) ----- --------------- ------------------- 1 4.79 0.00 2 5.13 2.00 3 5.58 4.00 4 5.84 6.00 5 5.96 8.00 6 6.09 12.00 7 6.21 16.00 8 6.50 20.00 9 7.85 30.00
This minimum 1-D model was the starting model for process 2 (Figure 3). The initial velocity model for process 4 (Figure 3) was derived by smoothing the output model of the nonlinear optimization (process 3), then constraining the low-velocity anomalies and reversals that appear in the regions not well controlled by the data, in the manner of Pullammanappallil (1994).
Fig. 8: Map view in gray scale of the velocity model obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with minimum 1-D model as the starting model (process 2). (a) at surface (depth = 0km) along with epicenters. (b) at 6km depth.
Fig. 9: Map view in gray scale of the velocity model obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with starting model estimated by simulated annealing (process 4). (a) at surface (depth=0km) along with epicenters. (b)at 6 km depth.
Figures 8(a) and 8(b) show map views at 0 and 6 km depth respectively of the final three-dimensional velocity model obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with minimum 1-D model as the starting model (process 2 preceded by process 1 - Figure 3). Figures 9(a) and 9(b) show corresponding map views for the model obtained by 3-D linearized inversion with starting model estimated by simulated annealing (process 4 preceded by process 3 - Figure 3). In general, the velocity model from processes 1+2 has a higher average velocity than that from processes 3+4 (5.97 km/s versus 5.55 km/s). Velocity at surface in the first case varies from 3.6 km/s to 6.8 km/s whereas that in the second case varies from 2.7 km/s to 6.0 km/s only. There is quite good match in the spatial distribution of low and high velocity patches in the northeastern half of the the two maps, though the absolute velocities differ. In the southwestern half, whereas the map for processes 1+2 shows a more or less uniform velocity of about 4.8 km/s, that for processes 3+4 shows a 10 by 10 km high velocity patch of 5 km/s in a background of 4.4 km/s. In both the cases, the north-northwest trending epicenter distribution is underlain by a parallel relatively low velocity area flanked by higher velocities. The maps at 6 km depth have fewer features in common. A relatively low velocity north-northwest trend is present in the northeastern part of the maps. The trend is much better-defined for processes 3+4 than the other case. Also, the high velocity patch featured in the southwestern part of the surface map persists at 6 km. Figures 10(a) and 10(b) present standard errors with respect to velocity models at 6 km depth obtained by processes 2 and 4 respectively. The northeastern half of the models show higher errors than the southwestern half. Treatment of these errors without referring to the distribution of rays can be misleading, as we will see. Figures 11(a) and 11(b) show the ray density in terms of derivative weight sum (DWS - discussed earlier) and diagonal elements of the resolution matrix at 6 km depth for process 4. The northeastern half shows much higher ray density compared with the other half. The highest resolution of 0.636 is obtained in the northern part of the model. Resolution in the southwestern part of the model is very low. If we examine the velocity maps along with the ray density and resolution maps, we can say that even though the northeastern part of the model has higher errors its velocities are more reliable than those in the southwestern half because both ray density and resolution in the latter is low. It means that the 10 by 10 km high velocity block found in the southwestern part of the model by process 4 is not as reliable as the NNE low velocity trend found by both processes 2 and 4.
Fig. 10: Map view in gray scale of the standard error at 6 km depth of the velocity models obtained by (a) VELEST followed by linearized inversion (process 2). (b) simulated annealing optimization followed by linearized inversion (process 4).
Fig. 11: Map view in gray scale of ray density and resolution at 6 km depth of the velocity model obtained by process 4. (a) ray density in terms of derivative weight sum (DWS) (b) resolution. Reliability of any part of the model is proportional to the magnitudes of resolution and DWS in that part.
Fig. 12: Cross-sectional view in gray scale of the velocity model obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with minimum 1-D model as the starting model (process 2). (a) along DD' in Figure 8(a). (b) along CC' in Figure 8(a). Hypocenters shown on the cross-sections are those located inside the volume limited by vertical planes along the dashed lines on either side of CC' and DD'.
Fig. 13: Cross-sectional view in gray scale of the velocity model obtained by simultaneous relocation and 3-D velocity estimation by linearized inversion with starting model estimated by simulated annealing (process 4). (a) along DD' in Figure 9(a). (b) along CC' in Figure 9(a). Hypocenters shown on the cross-sections are those located inside the volume limited by vertical planes along the dashed lines on either side of CC' and DD'.
Figures 12 and 13 show vertical cross-sections of models obtained by processes 2 and 4 along CC' and DD' defined in Figures 8(a) and 9(a). The hypocenters shown in the cross-sections are those located inside the volume limited by vertical planes along the dashed lines on either side of CC' and DD'. The low velocity trend observed on the map views are found as a basin in cross-section DD' of both the processes. The basin is about 8 km wide and 3-4 km deep. As observed in the case of map views, the absolute velocity of the basin in the case of process 4 is 0.8 km/s lower than that in the case of process 2. The mainshock is situated at depth under the western flank of the basin.
Because most of the hypocenters are shallower than 10 km, ray coverage at greater depths is poor and as such, the velocity cross-sections are not quite reliable at depths more than 10 km. The near-surface DWS is quite high in the eastern half of both sections CC' and DD' (not shown). Between 5 and 10 km depths the area with high DWS extends farther westwards. Therefore, the low velocity basin featuring in both map views and cross-sections for processes 2 and 4 is a structure supported by data, whereas the existence of the high velocity block found in process 4 is not equally supported by data. But its existence cannot be ruled out altogether. A better coverage of local stations in the southwest quadrant of the model would help in constraining the existence of the high velocity block (see Figure 1).
We can see a good correspondence of the gravity and topographic highs and lows with the high and low velocity areas in the velocity. The linear north-northwest low velocity trend coincides with the strike of the alluvial valley and also with a trend between two gravity lows. The 10 by 10 km high velocity block matches with a gravity high of comparable dimension.
In addition to qualitatively matching between trends in topography, gravity, and computed velocity model, we made a quantitative comparison of the fit between gravity anomalies and the two velocity models. Any gravity modeling from seismic velocity data first requires the velocities to be converted to densities. This conversion is almost always a problematic job because of the widely varying conditions under which different researchers have found their empirical relationships. Some of the formulas are derived entirely from studies on igneous and metamorphic rocks (Birch, 1961; Bateman and Eaton, 1967), while others are derived from exclusively sedimentary rock samples (Nafe and Drake, 1963; Gardner et al., 1974). Not a single one of these relationships can be used with perfect relevance to our study area. We used the following empirical formula of Gardner et al. (1974) to convert three-dimensional velocity models into three-dimensional density models:
Fig. 14: Bouguer gravity anomaly map (smooth contours) of the area superimposed on the gravity anomaly (gray scale) computed from densities estimated from velocity model obtained by (a) process 2; (b) process 4. Plus and minus signs indicate relative gravity highs and lows. The numbers on the contours indicate the Bouguer gravity anomalies according to Oliver and Robbins (1978).
Figures 14(a) and 14(b) show the computed gravity values (in gray scale) superimposed on the Bouguer gravity anomaly contour map of Oliver and Robbins (1978) for processes 2 and 4. Both the computed gravity maps show maximum contrasts twice as much as that in the Bouguer anomaly map. Whereas the maximum contrast on the Bouguer anomaly map is 25 mGal, those for processes 2 and 4 are 52 and 56 mGal respectively. Computed gravity values for process 4 are mostly negative, while those for process 2 are predominantly positive. Evidently this difference is due to the lower velocities of the process 4 model, which is on average 0.8 km/s slower than the process 2 model. The shape of the valley is modeled well by both the processes. The discrepancy between the magnitudes of gravity contrast on the Bouguer map and our modeled maps may suggest that even though outlines of the major low and high velocity regions have been obtained by the inversion, the details of lateral heterogeneity and of vertical variation of velocity have not been adequately discerned. Of course, another reason for the discrepancy lies in the inevitable difference between the physical conditions under which the empirical relationship was obtained and those prevailing in the Eureka Valley (exclusively sedimentary versus both sedimentary and, igneous and metamorphic rocks).
So far, we have seen that both processes 2 (linearized inversion with initial model obtained by minimum 1-D velocity modeling) and process 4 (linearized inversion with initial model obtained by nonlinear optimization by simulated annealing) map a few common features. The low velocity basin is mapped well by both the processes. But at depth the process 2 performs much more poorly than process 4 (Figures 8(b) and 9(b)). At 6km depth, the general appearance of the process 4 model is much more coherent and realistic. Also, the initial minimum 1-D model used for process 2 has a too high P wave velocity (4.79 km/s) at surface. Even though we attempted quite a few trial and error steps for the minimum 1-D modeling the results probably correspond to a local minimum. This is the inherent danger in any linearized inversion used for solving a nonlinear problem, especially in the absence of prior information. The relatively higher velocity obtained in the process 2 is reflection of the high velocity initial 1-D model.
We see that the simulated annealing optimization provides a velocity model that can be used as an initial model in a simultaneous linearized inversion of velocity and hypocenters, to obtain realistic velocity models and hypocentral locations. Hypocentral locations as such from the optimization alone are very much scattered. This may be an effect of not allowing the right degree of hypocentral perturbation during the optimization of successive trial models. Pullammanappallil and Louie (1993), in a similar simultaneous optimization of velocities and reflector locations, found that inappropriate constraint of reflector "hypocenters" led to poor results. These observations led us to do simultaneous linearized inversion for hypocenters and velocity using the velocity results of nonlinear optimization as the initial velocity model, producing the successful results above.
As a result of combining a linearized inversion with a nonlinear optimization in tandem, we were able to obtain a realistic velocity model and fault plane definition from hypocenters. The velocity model adequately describes broad structural features at the local scale. The nonlinear simulated annealing optimization scheme, by virtue of its capability of finding the global error minimum irrespective of any starting model, produced a velocity model for an area like Eureka Valley where little is known about the local velocity structure. The linearized inversion following the nonlinear optimization could thus start from a fairly valid initial model and preliminary hypocentral locations, and fine tune them.