Serdar Özalaybey, Martha K. Savage*, Anne F. Sheehan**, John N. Louie, and James N. Brune
Seismological Laboratory, University of Nevada at Reno, Reno, NV
*also at Institute of Geophysics, Research School of Earth Sciences, Victoria University of Wellington, New Zealand}
**Cooperative Institute for Research in Environmental Sciences and Department of Geological Sciences, University of Colorado, Boulder, CO
A new method based on the joint inversion of receiver functions and surface-wave phase velocities results in well-determined shear-velocity structures that are consistent with the compressional-wave structure, gravity, heat flow, and elevation data in the northern Basin and Range. This new inversion method takes advantage of average velocity information present in the surface-wave method and differential velocity information contained in the receiver function method, thus minimizing the non-uniqueness problem that results from the velocity-depth trade-off.
An unusually thick (38 km) and relatively faster crust and upper mantle are found in central and eastern Nevada compared to the thin (28 to 34 km) and slower crust and upper mantle of the western Basin and Range. We interpret the regions of thicker and faster crust and upper mantle as zones that have undergone less Cenozoic extension relative to the surrounding regions to the west and north. The thick crust and consequently greater depth to the dense mantle material is consistent with the gravity low pattern in central and eastern Nevada. Simple gravity modeling shows both local and regional isostatic compensation occur within 40 km of the surface, indicating a near classical Airy type of compensation in the province.
We analyze in detail the shear-wave (S-wave) velocity model derived from the receiver functions at station BMN and compressional-wave (P-wave) velocity models derived from the 1986 PASSCAL experiment in northwestern Nevada. The most interesting feature of these models is the presence of negative-velocity gradients in the S-wave model with no corresponding velocity decrease in the P-wave models between 10 and $24 km depth. This combined velocity model may be explained by high pore fluid pressures at these depths. This model favors layered fluid porosity model proposed in the literature to explain extensive middle to lower crust continental seismic reflections and high electrical conductivity. An upper mantle, gradational low-velocity zone is present between 32 and 38 km in the S-wave model. This upper mantle shear-wave low-velocity zone is consistent with partial melt, which may be the source material for magmatic underplating in this region.
Receiver function modeling is a valuable means of obtaining local shear velocities in the crust and upper mantle beneath three component, broadband, digital stations (e.g., Langston, 1979; Owens and Zandt, 1985; Ammon and Zandt, 1993; Kind et al., 1995; Cassidy, 1995). It is well-known, however, that the receiver function inversions for shear-wave velocity structure are non-unique since there is very little absolute velocity information contained in the receiver functions. This lack of information causes the non-uniqueness problem known as the velocity-depth trade-off (Ammon et al., 1990). The velocity-depth trade-off occurs because the receiver function technique is a differential method containing no absolute travel time information; a slow, thin layer will produce the same average differential arrival time as a fast, thick layer. This velocity-depth trade-off has been shown by Ammon et al. (1990) in the form of the initial velocity model dependency of the inverted final velocity models. Using a series of initial velocity models, a group of final velocity models with differing average velocities can be found which fit the observed receiver function equally well, thus resulting in differing structural images. Based on this analysis, Ammon et al. (1990) urged that all future receiver function studies perform non-uniqueness analysis to exclude velocity models that do not fit the other a priori geophysical constraints obtained from refraction, reflection, earthquake travel-time, and surface-wave data where available.
We address the non-uniqueness problem by combining receiver function inversion with surface-wave phase velocity data. Surface-wave (fundamental-mode) phase velocities are sensitive to the average shear-velocity structure of the material within the depth ranges to which they penetrate. The phase velocities provide information on the long-wavelength vertical averages of the shear structure between any given station pair. This information is essentially absent in the receiver functions, which are sensitive to velocity contrasts between materials. Thus, at depths sampled by both methods, a joint inversion of surface-wave phase velocities and receiver functions will result in a better-determined shear-wave structure than inversion from either method alone.
Compressional-structure studies also use a similar combination of two data sets to obtain well-determined velocity models of the crust. These data sets are from wide-angle refraction/reflection and near-vertical incidence reflection profiles (e.g., Holbrook et al., 1991; Catchings and Mooney, 1991). While seismic refraction data constrain absolute velocities, reflection data constrain a high resolution structural image.
Little is known about the detailed shear-wave structure of the northern Basin and Range Province (NBRP) other than existing gross shear-wave structure of the crust and upper mantle from the regional surface-wave phase velocity data inversions (e.g., Priestley and Brune, 1978; Priestley, Orcutt, and Brune, 1980; Priestly and Brune; 1982; Taylor and Patton, 1986) and from the attenuation of teleseismic and surface waves (e.g., Taylor et al., 1986; Patton and Taylor, 1984). On the contrary, there is a large volume of work that has produced detailed compressional-wave structure of the crust from modern reflection and refraction profiles (e.g., Klemperer et al., 1986; Allmendinger et al., 1987; Potter et al., 1987; Holbrook, 1990; Catchings and Mooney, 1991).
One-dimensional shear-wave structures obtained from the inversion of receiver functions have the potential to offer vertical and horizontal resolution comparable to the compressional-wave structures obtained from refraction studies (Zandt and Owens, 1986). At a given frequency, shear waves travel at shorter wavelengths than compressional waves and thus maintain higher resolution; a shear waveform with a dominant frequency of 1 Hz is expected to be influenced by layers as thin as 1 km (Owens, 1987). Shear waveforms also display stronger anomalies than the compressional waveforms for the presence of seismic anisotropy, fluids and fractures in crustal rocks, and partial melt in the upper mantle. With high quality receiver functions, all these potentially possible effects can be studied in detail.
In the following sections, we show that the joint inversion of receiver functions and surface-wave phase velocity data is a feasible, robust, and valuable method for obtaining well-constrained shear-wave velocity models. We show a sample synthetic test of this method and then present the results of using this method on real data to obtain one-dimensional crustal shear-wave structures beneath nine stations located in the northern Basin and Range Province (NBRP). We summarize our shear-wave structures with three parameters that represent local averages, namely crustal thickness, mean crustal and upper mantle shear-wave velocities for each station. These parameters are then interpreted together with compressional-wave structure, gravity, heat flow, and elevation data to detect regional variations in the structure.
We model shear-wave velocity structure beneath five stations (DNY, BMN, WHR, KVN, and MNA) from the permanent broadband digital network of the University of Nevada, Reno (UNR) and four portable stations (MLC, RTS, WCP, and GAR) from a recently conducted field experiment (Sheehan et al., 1995) (Figure 1).
The five UNR stations are equipped with broadband three-component Geotech sensors (models 7505Z, 8700H, SL-210Z, and SL-220H, 15 s natural period) (Peppin and Nicks, 1992). The portable stations were deployed with CMG3-ESP broadband sensors for a period of ten months. These stations are capable of recording seismic signals with a flat velocity response in the frequency range 0.3 to 20 Hz. We have analyzed the calibration pulses to make sure that the vertical and horizontal component sensors at each station have nearly the same amplitude and phase response within this frequency range. In addition to these stations, UNR broadband station WCN, and short-period stations RYN and MIL (1 s natural period, vertical component) are used to increase path coverage for the surface-wave phase velocity measurements.
Receiver functions are extracted from the three-component broadband recordings of teleseismic P waveforms at epicentral distances ranging from 30 to 85^o. For this range, P waves are steeply incident and energetic, and are recorded dominantly on the vertical component while P to S converted shear waves are recorded dominantly on the horizontal component seismogram. The extraction procedure is described by Langston (1979) and consists of a simple deconvolution of vertical component seismograms from the radial and transverse components. The receiver functions obtained this way are assumed to be free of source, mantle-path, and instrument-response effects. The deconvolution is performed in the frequency domain (Owens et al., 1984) using a Gaussian filter and spectral trough filler. We have used a Gaussian filter width of 2.5 and trough filler value of 0.01, which produced stable deconvolutions. True amplitudes of receiver functions are preserved in the deconvolution processes. The use of true amplitudes preserves information about the shallow velocity structure (Ammon, 1991) and helps avoid inaccuracies due to the presence of dipping layers (Cassidy, 1992). The receiver functions have been extracted from many events clustered in distance and back-azimuth ranges. These receiver functions are stacked to decrease noise.
Figure 2 shows the stacked receiver functions for the UNR permanent stations from three back-azimuth ranges grouped as southeast (SE), northwest (NW), and southwest (SW). Northeast receiver functions are not shown because of lack of events in this quadrant. In general, the SW receiver functions are noisier and associated with the highest transverse component arrivals. The transverse component arrivals may be due to scattering and/or anisotropy and dipping layers. The transverse component energy is, therefore, a measure of the deviation from isotropic and planar structure (Ammon and Zandt, 1993). The SW receiver functions are extracted from a distance range 65 to 85^o and therefore are associated with steep angle of incidences at the base of the crust and low signal strength, generating weaker P to S conversions than closer events. Because of the low quality of the SW receiver functions, we do not use them in the inversions. The SE and NW receiver functions have strong P to S (Ps) converted arrivals (Moho Ps phase following the direct P-wave arrival) (Figure 2).
We model the radial component receiver functions from the NW and SE back-azimuths assuming isotropic and plane layered crustal structure. The transverse component receiver functions can be modeled for the effects of anisotropy and/or dipping structure when a high number of events are used in stacking and/or data are available from an array of stations. McNamara and Owens (1993) show an example of such a study for the effect of anisotropy analyzing more than 100 Moho Ps conversions from a large aperture station array in west-central Nevada (near stations BMN and WHR) (Figure 1). They find a fast azimuth of anisotropy oriented NW-SE with average delay times of 0.2 s (due to shear-wave splitting in the Moho Ps phases), consistent with their observations of minimal transverse component energy excitation associated with energetic radial Moho Ps arrivals from the NW and SE back-azimuth events. They suggest that this anisotropic fabric was caused by preferred orientation of seismically anisotropic minerals in the middle to lower crust. The seismic reflection transects throughout the province have also shown a relatively flat Moho and subhorizontally layered middle and lower crust (e.g., Klemperer et al., 1986), but a highly heterogeneous, complex upper crustal structure interpreted from the surface structural geology correlations (Catchings, 1992).
On the basis of above results, the effect of anisotropy is expected to be minimal on our NW and SE receiver functions because these azimuths are along the inferred fast direction of anisotropy, which should produce no significant transverse component energy. According to the reflection data, the effects of dipping and laterally heterogeneous structure from the middle and lower crust must also be small in generating significant transverse component energy. Thus, we believe that the transverse component receiver function arrivals we observe are probably the result of complex near surface scattering in the upper crust, which may not be modeled in any simple way.
Figure 3 shows the stacked receiver functions from the SE back-azimuth at three distance ranges for stations BMN, WHR, and MNA. We group these receiver functions from distances 28 to 45^o as Range 1, 50 to 60^o, as Range 2, and greater than 65^o as Range 3. The receiver functions corresponding to Range 1 and 2 show coherent primary and reverberated P to S converted arrivals while the receiver functions corresponding to the most distant events from Range 3 do not generate the same arrivals, possibly due to their steep angle of incidence and low signal strength (Figure 3).
The variations of amplitude and phase among different ranges, corresponding to different sampling of ray parameters, are important in constraining the velocity structure. The moveout of the multiples that travel as both P and S waves in the crust can provide average velocity information on both compressional- and shear-wave velocities. We take advantage of this by jointly inverting receiver functions from different distance ranges.
Figure 4 shows the stacked receiver functions only from the SE and NW back-azimuth quadrants for the portable stations WCP, GAR, RTS, and MLC. The quality of the receiver functions from these stations are lower than that of the UNR permanent stations because of the short period of deployment resulting in fewer stackable events. However, we have at least one well-determined back-azimuth quadrant suitable for the inversion at each station.
We have measured Rayleigh-wave phase velocities for three paths using a regional event from Oregon and a teleseismic event from the Tonga Islands (Figure 5). While the Oregon event provided data for the measurement of short period (10 s) Rayleigh-wave fundamental-mode phase velocities near station MNA (Path 3), the Tonga event provided data for longer period (20-25 s) Rayleigh-wave fundamental-mode phase velocity measurements along Path 1 and Path 2 (Figure 5). Both events have a high signal-to-noise ratio. The reason for observing these periods only is due to the energy produced by an event of a given size at a given distance and the instrument response. Interstation Rayleigh-wave phase velocities were measured using the well-known peak and trough method (Brune et al., 1960) after correcting for instrument phase response differences between stations. The peak and trough method is based on the time domain phase correlation of surface-wave peaks and troughs of nearly the same periods recorded at a station pair lying on the same great circle path.
Measured phase velocities are shown in Figure 5. The accuracy of these phase velocities is estimated to be within +/- 0.1 km/s, which is nearly equal to the scatter in the measured phase velocities between different station pairs. These phase velocities also deviate by only +/- 0.1 km/s from the Great Basin model dispersion curve of Priestley and Brune (1978).
Receiver functions are modeled by a 24 plane-layered structure extending from surface to a depth of 40 km. The layer thicknesses in our model parameterization are 1 km for the top 6 km depth and 2 km below this depth. The inversion procedure of Ammon et al. (1990) consists of a linearized-iterative, least squares waveform fitting. An initial shear-wave velocity model solution is taken to start the iterative inversions of matrix systems representing the minimization of the L2 norm of the misfit residual vector between the observed and modeled receiver functions. The smoothness constraints are applied to minimize the energy in the second difference of the model parameters. We use a small amount of smoothness during inversions (smoothness parameter of 0.1-0.3). To compute modeled receiver functions for a given velocity structure, a technique based on propagator matrix method is used (Kennet, 1983). The compressional-wave velocities (assuming a Poisson's ratio of 0.25) and densities are adjusted by the relationships (Vp=sqrt{3}.Vs and rho=0.32Vp+0.77) given by Owens (1987). The inversions quickly converge within five iterations.
In addition to the receiver function inversion alone, we use two-station surface-wave phase velocity measurements to constrain the final velocity models. The surface-wave phase velocity measurements provide a priori constraints that any velocity model fitting the receiver function must satisfy. To implement this constraint, we have modified the inversion method of Ammon et al. (1990).
We solve a system of matrix equations in the form
where D is a matrix containing partial derivatives of receiver function with respect to unknown velocity model m, mo is the initial velocity model, F is a matrix of constraints made up from partial derivatives of surface-wave phase velocities with respect to m, Delta is the smoothness constraint matrix forming the second difference of the model m, sigma is the adjustable parameter controlling the trade-off between fitting the waveform and smoothness of the model (Ammon et al., 1990), and r_r and r_c are vectors containing receiver function waveform and phase velocity residuals, respectively. The partial derivatives in the constraint matrix F contain localized averaging functions of the model parameters associated with different periods of phase velocities. We show this behavior of F in the next section on a synthetic example. This additional averaging information helps to resolve the non-uniqueness problem caused by the velocity-depth trade-off.
In this section, we present a sample synthetic test to validate the constrained inversion method formulated above. We take the shear-wave velocity model of Mangino et al. (1993) and generate a synthetic radial component receiver function representing our observed true-model receiver function data. We contaminate this receiver function by adding a real transverse component receiver function representing noise. The maximum amplitude of the noise equals 10% amplitude of the radial receiver function. Rayleigh-wave phase velocities at 10 and 20 s periods are used as the surface-wave phase velocity data. We choose the Great Basin shear-wave structure of Priestley and Brune (1978) as our initial model. Three velocity models are inverted using the receiver function and three pairs of phase velocity data. These pairs are obtained by adding and subtracting 0.2 km/s from the true-model phase velocity pair, respectively.
Figure 6 shows these velocity models together with the theoretical dispersion curves for each model, a plot of Rayleigh-wave partial derivative-depth curves for the phase velocities, and receiver function waveforms computed for these models. The Rayleigh-wave partial derivatives of phase velocities with respect to model parameters and theoretical dispersion curves for a given velocity model are computed using the method of Rayleigh's variation principle and stress-motion propagator matrices, respectively (Aki and Richards, 1980). All surface-wave computations in this paper are made with a computer program package written by Zeng (1995), who implemented these methods into discrete computer codes (personal communication).
These models differ the most in their average velocities between 10 and 35 km depth range; the dotted model represents a high-velocity crust, dashed model a low-velocity crust, and solid model a crust with velocities between the other two models (Figure 6a). The latter model is found by the constrained inversion of the receiver function and the true-model phase velocity pair (3.06 km/s at 10 s and 3.51 km/s at 20 s) (Figure 6c) and is nearly identical to the true-model velocity structure. The high-velocity and low-velocity crust models are found by using the highest and lowest phase velocity pairs, respectively. All three models fit the true-model receiver function with a waveform correlation coefficient of 0.993 or higher (Figure 6d).
The Rayleigh-wave partial derivative-depth curves shown in Figure 6b make up the constraint matrix F in the matrix equation of (1). Clearly, the 10 s partial derivative curve is a localized averaging function for a depth range between 8 and 20 km and 20 s partial derivative curve is a localized averaging function for depths of 15 km and deeper. It is this localized averaging behavior of the partial derivatives which leads to a unique velocity model as shown in this test. Thus, a pair of phase velocity data sampling different regions of the model space (corresponding to different periods of phase velocity measurements) is sufficient to provide the information needed to remove the non-uniqueness in the solutions. This is why we refer to the F matrix as the matrix of constraints because we do not need a full joint inversion of phase-velocity data with the receiver functions, but only a pair of well-determined phase velocity points. This avoids measurement of a full dispersion curve. In this way, we let the long-wavelength features of the velocity model, i.e. its average velocity, be resolved with the phase velocity data and the short-wavelength details be resolved by the receiver function data.
Another way of looking at these partial derivatives is that they construct long-wavelength smoothness constraints. No weighting, i.e. unit weighting, between the receiver function and phase velocity data was needed in solving the matrix equation of (1) for a good fit of both data sets. Further, although this study includes only fundamental-mode Rayleigh-wave phase velocity derivatives, in principal the F matrix can include different types of phase velocity derivative functions, which may be computed for travel-times of earthquake S phases and for fundamental- and/or higher-mode surface waves of Rayleigh and/or Love type. Thus, any constraint that has relevant shear-wave velocity information can be used in the $F$ matrix. We conclude this section by stating that the constrained-inversion method we describe here is a feasible, robust, and valuable method for obtaining well-constrained velocity models.
Two initial models have been chosen to invert receiver functions. For the northwestern Basin and Range, we have taken the compressional-wave reflection velocity model of Benz et al. (1990) (N-S profile from SP8) and converted to shear-wave velocity structure (assuming a Poisson's ratio of 0.25) for our model parameterization. This model is used as the initial model for stations DNY, BMN, WHR, RTS, WCP, and GAR. For stations MNA, KVN, and MLC, we use the final shear-wave velocity model of Mangino et al. (1993) obtained from the receiver function inversion at station MNV, located within 350 m of our station MNA.
We generate 45 different initial models from each initial model by perturbing them with a cubic perturbation of 1.0 km/s and 20% random component (Ammon et al., 1990). We repeat inversions for each perturbed initial model to find joint models that fit both the receiver function and phase velocity data. The inversion of many initial models is needed because inversions starting with models that do not predict phase velocity data well result in divergence. This eliminates possible solutions that would fit the receiver function but not the phase velocity data.
The following sections present the results of receiver function inversions. This section will be presented in the following format: stations are grouped by their geographical locations as central (MNA, KVN, and MLC), eastern (RTS, WCP, and GAR), and northern (DNY, BMN, and WHR). The surface-wave phase velocity measurements shown in Figure 5 are used to constrain each inversion as follows; the long-period measurements from the Tonga event from Path 2 for the central and eastern and from Path 1 for the northern stations, the short-period measurements from the Oregon event from Path 3 for the central stations (MNA and KVN only). We present an inversion summary plot (Figures 7-9) showing the final shear-wave velocity structure, the phase velocity data, and receiver function waveform fit for stations MNA, MLC, and BMN (for other stations see \"{O}zalaybey, 1996). We identify the Moho with an arrow in these plots. Figures 10 and 11 show the final shear-wave velocity models for all the stations. Our Moho picks are based on the nature of crust-mantle transition; i.e. the first increase in velocity either sharply or with a gradient at the expected velocities for the Moho (> 4.3 km/s). Where we have gradational boundaries, we pick the Moho at depths where shear-wave velocities first exceed approximately 4.3 km/s.
MNA receiver functions from the NW show more scattering than the SE back-azimuth counterpart (Figure 2). We focus on the SE back-azimuth receiver functions because they are simpler and display the same arrivals from all three ranges on the radial receiver functions (Figure 3) The transverse receiver functions are less than 1/3 amplitude of the radial receiver function arrivals. The arrivals for the first $6$ s, however, seem to be coherent from all three ranges on both the radial and transverse receiver functions. These receiver functions are nearly identical to the receiver functions of Mangino et al. (1993), who inverted shear-wave velocity structure at station MNV (only 350 m away from our station MNA). We inverted all three receiver functions simultaneously. Figure 7 shows the inversion results. The final velocity model (Figure 7a) fits both the receiver function and observed phase velocity data considerably well (Figures 7b and c). The main features of the velocity model are a sharp velocity increase in the upper 5 km (1.8 to 3.5 km/s) followed by a 2 km thick low-velocity layer between 6 and 8 km depth, a nearly constant velocity crust (3.5 km/s), and a crust-mantle transition between depths of 32 and 34 km. The shallow low-velocity layer is a required feature of this model. The removal of this layer increases the misfit to Moho Ps phase by nearly 40%. Mangino et al. (1993) interpreted a similar feature as possible artifacts of scattered energy from a three-dimensional structure.
Station KVN has one of the lowest quality receiver functions from all back-azimuth quadrants. It exhibits large scattered arrivals for the SE quadrant (Figure 2). The NW quadrant receiver function seems to be the simplest for inversion. The final velocity model for this quadrant is shown in Figure 10. We avoid fitting the details of this receiver function by using a high value of smoothness parameter (sigma=0.3) and by smoothing several alternating high and low velocity layers that provided no significant fit. The high value of smoothness constraint for this inversion caused a poor fit to the receiver function but a good fit to dispersion data. This model indicates a crustal thickness of 32 km (Figure 10).
For station MLC, we invert the NW receiver function because it displays much smaller transverse component and better defined radial component arrivals with tighter standard deviation bounds than the SE receiver function (Figure 4). The broad Moho Ps phase clearly seen arriving near 5 s implies a thick crust at this station. Figure 8 shows the inversion results. The most important feature of the velocity model is the gradational crust-mantle transition. The upper mantle shear velocities are reached at about 38 km confirming a very thick crust beneath this station. An upper and lower crustal low-velocity zone (LVZ) are also present in this model.
At station BMN, we inverted only the SE receiver functions because other back-azimuth quadrants are not characterized yet. The receiver functions from this station are one of the highest quality set and show a clear Moho Ps phase (at 3.8 s) and its multiple PpPms near 12 s from Range 1 and Range 2 (Figure 3). The amplitude of the transverse component receiver functions is half the size of the radial components for these arrivals. We inverted the receiver functions from these two ranges simultaneously to obtain a joint shear-wave velocity model. Figure 9 shows the resulting final velocity model. This model is smooth and shows a sharp crust-mantle boundary at $28$ km. An LVZ is present between depths of 32 and 38 km. We will analyze this model further in the following sections.
WHR receiver functions from the SE back-azimuth look similar at each range, but the receiver function from Range 2 has a better standard deviation and smaller transverse component arrivals (Figure 3). Note that there is a strong, coherent transverse arrival near 7 s indicating deviations from laterally homogeneous, isotropic structure. The NW receiver function is also well-determined and displays similar arrivals (Figures 2 and 3). We inverted the receiver functions from the SE Range 2 and NW jointly. The final velocity model shows a gradational crust-mantle transition starting from 28 km and reaching upper mantle shear-wave velocities at 34 to 36 km. The large negative trough seen near 1.7 s in the NW receiver function requires a sharp velocity increase beneath the surface layer (Figure 10).
The receiver functions from station DNY are the lowest quality in our data set. Both the SE and NW receiver functions have large standard deviations and transverse component arrivals. Large secondary arrivals after the Moho Ps arrival near 4.1 s dominate both receiver functions and are probably associated with large amount of scattering (Figure 2). We have only inverted the NW receiver function because it is somewhat better than the SE one. The inversion of this receiver function was problematic because the arrivals around the Moho Ps phase and other late secondary arrivals required a model with alternating high and low velocity layers. We smoothed out these layers and repeated the inversions with a high value of smoothness parameter (sigma=0.3) until a satisfactory fit with a smooth model is obtained. The waveform fit was poor at the expense of the model smoothness constraints we employed. A sharp increase in velocity from 28 to 30 km marks the crust-mantle boundary. An immediate LVZ is present beneath this boundary, which may be the result of scattered energy (Figure 10).
At station RTS, the Moho Ps phase arrives near 3.9 s in the NW receiver function (Figure 4) but about 0.5 s later in the SE receiver function. This may indicate thicker crust to the SE. But, the NW receiver function has more events in the stack and smaller transverse component arrivals. Therefore, we inverted the NW receiver function using a high value of smoothness parameter (sigma=0.3). This caused a relatively low fit to the negative polarity arrivals around the Moho Ps phase, which may be indicative of LVZ's beneath and above the Moho and/or complex crustal scattering. More data are needed to model the details of the receiver function at this station. However, the Moho Ps phase and phase velocity data were fit well. A sharp crust-mantle boundary is evident at a depth of 32 km at this station (Figure 11).
GAR receiver functions from both back-azimuth quadrants display common arrivals for the first 6 s (Figure 4) and have been inverted jointly. The Moho Ps phase arrives late (near 5 s) indicating a thick crust. The transverse component arrivals are extremely small especially for the NW back-azimuth. Both the phase velocity data (from Path 2 RTS-WCP pair) and receiver functions were fit well. A gradational crust-mantle transition with a crustal thickness around 38 km is present in the velocity model (Figure 11).
The WCP receiver function from the SE back-azimuth is dominated by large transverse component arrivals. The NW receiver function shows a strong and well-determined Moho Ps phase (Figure 4). We inverted this receiver function (grouped into two distance ranges not shown in Figure 4). The velocity model from this inversion indicates a sharp velocity gradient at the top of crust-mantle transition with a 32 km thick crust (Figure 11).
The joint inversion method eliminates many non-unique solutions, which would otherwise be obtained from the inversion of receiver functions alone (Ammon et al., 1990). This is because the velocity-depth trade-off problem inherent in receiver function inversions is largely minimized by the additional surface-wave phase velocity constraints imposed during the inversions. This constrained inversion method results in only 1-3 accepted final velocity model solutions (out of 45 inversions with different initial models) with essentially the same velocity features and average velocities for each specific inversion at a given station. In other words, there are no other velocity models significantly different from our final velocity models that could fit both the receiver functions and and surface-wave phase velocities in the whole one-dimensional solution space. This should also be clear from the synthetic test we have presented previously. The sensitivity analysis is aimed to show the possible range of velocities for a specific velocity feature in a given final velocity model due to the uncertainties in our data.
Synthetic receiver functions and dispersion curves calculated from the final velocity models presented above provide good fits to our data sets. In order to quantify how well various features are constrained in these models, we perform a series of sensitivity tests on the data fits. We believe that the sensitivity analysis is a physical and instructive way of displaying the possible range of velocities for a given distinct velocity feature.
We consider a maximum +/- 0.1 km/s phase velocity measurement error and one standard deviation bounds of the receiver functions as allowable errors in our data fits. In general, the surface-wave phase velocities are sensitive to the absolute average velocities and the receiver functions are sensitive to the change in velocity delta beta. As mentioned earlier, the receiver functions resolve the shape of the velocity structure and the surface waves resolve average velocities. Thus, we have a good control on average velocities and velocity contrasts. Our sensitivity analysis is performed in two steps; in the first step, we add positive and negative velocity perturbations delta beta to the given velocity feature in question and keep the rest of the model the same. This results in two end-member models. In the second step, we compute predicted phase velocity and receiver functions for these models and compare these predicted data with the real data. We repeat this process until the velocity perturbations are so large that we cannot accept the misfit between the predicted and observed data. This provides us with confidence bounds for each velocity feature. The details of the sensitivity analysis can be found in the work by Özalaybey, 1996.
Figures 10 and 11 present summary plots of the velocity models for all the stations with confidence bounds obtained from the sensitivity analysis and our Moho picks marked with thick black arrows. Based on these confidence levels, the low-velocity layers less than 6 km thick can have maximum velocity perturbations of delta beta = +/- 0.2 km/s. Stations KVN, MNA, MLC, GAR show upper crustal low-velocity layers. The middle crustal and upper mantle velocities can have perturbations up to delta beta = +/- 0.1 km/s. We compute average crustal velocity (beta_c) and upper mantle (beta_m) using our Moho picks. From the confidence bounds, we obtain uncertainties in beta_c and beta_m as +/- 0.05 km/s and +/- 0.1 km/s, respectively. The receiver functions and Rayleigh-wave phase velocities are very insensitive to the changes in the compressional-wave velocities. To test this, from the final shear-wave velocity we use a Poisson's ratio of 0.25 to compute the compressional-wave velocities. Then, we perturb compressional-velocity by +/- 0.3 km/s while leaving the shear-velocity constant, corresponding to a variation in Poisson's ratio of +/- 0.035. The compressional velocities within these limits all provided fits to the data within the accepted error bounds.
We perform a more detailed sensitivity analysis for the BMN velocity model because it has two distinct features that are of interest in our discussion. These features are the negative-velocity gradients between 10 and 24 km and 32 and 38 km (Figures 9 and 10). The removal of these negative-gradients (seen in regions of 2, 3, and 4 in Figure 10) by applying positive velocity perturbations one at a time and keeping the rest of the model the same as before, resulted in synthetic receiver functions that had considerable misfits to both observed bounding receiver functions and the surface-wave phase velocity data. In particular, the removal of each negative-velocity gradient caused a timing and amplitude mismatch to direct Ps arrivals from the 18 km discontinuity arriving near 2.5 s and from the Moho near 3.8 s, and the Moho multiple PpPms near $12$ s (Figure 9). These mismatchs were 0.05 to 0.1 s early in timing and 20 to 40% in amplitude. The timing of the Moho multiple was early by nearly twice that of the direct Ps arrivals because this multiple traveled twice the distance through the model. The removal of the deeper low-velocity layer resulted in overprediction of the Moho Ps amplitude by 40% with almost no effect on the Moho multiple. The overprediction of amplitude occurs because the negative-velocity gradient layer in the mantle generates a negative polarity Ps arrival, which reduces the positive polarity arrival from the Moho. This also resulted in phase velocity values that exceeded our accepted error bounds. The estimated velocity uncertainty is +/- 0.1 km/s for this feature.
We also test if the negative-velocity gradient layers can be traded-off by a model that contains smaller velocity jumps followed by smooth positive-velocity gradients. We compute synthetic waveforms for such a model that has half the velocity jumps at 8 and 18 km discontinuities followed by positive velocity gradients while leaving the rest of the model unchanged. The largest misfit for this model was in the timing of direct Ps arrivals from the Moho and the 18 km discontinuity, which was early by 0.15 s. This model also predicted other direct arrivals that were not present in the data.
Our sensitivity analysis indicates that average velocities and crustal thicknesses are the best resolved parameters. Therefore, we choose to interpret these parameters together with compressional-wave structure, gravity, and heat flow data, to detect regional variations in the crustal structure. Figure 1 shows these parameters plotted on the Bouguer gravity anomaly map of the NBRP.
Regional crustal thickness variations derived from the receiver function and surface-wave phase velocity data correlate well with the heat flow, gravity, and elevation data. The thinnest crust is found beneath stations BMN and DNY (28-30 km) in the northwestern Basin and Range. This region is characterized by the highest heat flow (> 2.5 hfu, Battle Mountain high anomaly) (Lachenbruch and Sass, 1978), higher Bouguer gravity (~ -170 mgal), lower elevations (~ 1.5 km), and younger volcanism ( < 17 Ma) than other regions. The thickest crust is found beneath stations MLC and GAR (38 km) in the central and eastern part of the province. These stations lie within the central Nevada highland with elevations reaching ~ 2.5 km, a broad region of the lowest gravity values (~ -220 mgal), and relatively older volcanism (17-34 Ma) (Blackwell, 1978). Station MLC is within the Eureka heat flow low ( < 1.5 hfu) anomaly of Lachenbruch and Sass (1978). Other stations have a crustal thickness 32 to 34 km and are located within transitional regions of elevations, heat flow, and gravity anomalies (Figure 1).
In general, average crustal and mantle velocities increase from west to east; stations located in the western part of the province (DNY, BMN, WHR, KVN, and MNA) with average $\beta_{c}$= $3.49$ km/s and $\beta_{m}$= $4.42$ km/s in contrast to the stations located in the central and eastern part of the province (MLC, GAR, RTS, and WCP) with average $\beta_{c}$= $3.55$ km/s and $\beta_{m}$= $4.58$ km/s. Our average shear velocities are in good agreement with the Great Basin model of Priestley and Brune (1978) which has an average crustal velocity of 3.54 km/s and mantle lid shear velocity of $4.5$ km/s. They also agree with other shear velocity estimates which range from $4.44$ and $4.62$ km/s (for a Poisson's ratio of 0.25) predicted from previously reported province-wide upper mantle Pn velocities ($7.7$-$8.0$ km/s) (e.g., Prodehl, 1979; Pakiser, 1989). In a more recent study, Hearn et al. (1991) found low Pn velocities ($\sim$ $7.8$ km/s) beneath the Basin and Range Province and suggested that these low velocities may result from hotter uppermost mantle in the province. Higher upper mantle velocities in central Nevada relative to western Nevada are also supported by the teleseismic P-wave residual, travel-time inversion results of Humphreys and Dueker (1994), which showed a similar contrast in mantle velocity between these regions. Average crustal shear velocities for the western stations correlate well with the average shear velocities of $3.4$-$3.5$ km/s reported for the crust by Hawman et al. (1990) in northwestern Nevada. The northwestern Basin and Range, compressional-wave, wide angle reflection/refraction, $1986$ Program for Array Seismic Studies of the Continental Lithosphere (PASSCAL) experiment (profiles shown as solid, white, straight lines in Figure 1) provides crustal thickness estimates we can compare directly with our crustal thickness estimates in northwestern and central Nevada. Several investigators have analyzed the data from this experiment using one- and two-dimensional methods and inferred similar crustal structure, mostly differing in the details of lower crust and upper mantle (e.g., Benz et al., 1990; Holbrook, 1990; Catchings and Mooney, 1991). The PASSCAL experiment determined that the crustal thickness varies laterally from $28$ to $36$ km in this region. The crustal thickness between SP1 (i.e. shot point 1) and SP4 (Figure 1) varies by $\pm 2$ km from the average thickness of $32$ km (Catchings and Mooney, 1991; Holbrook, 1990). The crustal thickness of $34$ km found beneath WHR and $30$ km beneath DNY is in good agreement with the PASSCAL experiment result. Between SP4 and SP8 (Figure 1), the PASSCAL experiment indicates a Moho with crust-mantle transition occurring at depths of $28$ and $33$ km. (e.g., Benz et al., 1990; Holbrook, 1990; Catchings and Mooney, 1991). Hearn et al. (1991) estimated a depth of $29$ km to the Moho from the station static of their tomographic study at BMN. Our velocity model at BMN, which is $\sim$ $40$ km east of SP9, displays features that warrant a detailed comparison with the PASSCAL experiment velocity models. For this comparison, we converted the two-dimensional compressional velocity models of Holbrook (1990) and Catchings and Mooney (1991) into one-dimensional P-wave models. These P-wave models are obtained by sampling their models near SP9 according to our model parameterization (the sampling of their models according to our model parameterization shifted a few interfaces at most $1$ km). Figure 12 shows these models together with our S-wave model. The structural complexity found in the upper crust with compressional-wave velocities ranging from $2.5$ to $6.0$ km/s (Catchings, 1992) between basins and ranges can easily account for the small differences present at shallow depths ($< 8$ km) between the models. Also, our surface-wave phase velocity data are too long period ($20$ s) to constrain shear-velocities at these shallow depths. Between $8$ and $28$ km, the velocity discontinuities in the P-wave model of Catchings and Mooney (1991) shows nearly one to one correspondence to the velocity discontinuities in the S-wave model. Holbrook's (1990) model, which tried to produce the simplest possible velocity variations to fit the data, shows a lesser degree of correlation with our model in this range. We note that in the S-wave model each sudden increase in velocity with depth is followed by a smooth, negative-velocity gradient (e.g., discontinuities at $8$, $18$, and $28$ km). These gradients are not present in the P-wave models (Figure 12). Below $28$ km, at the crust-mantle transition, both P-wave models have an upper mantle velocity of $7.9$-$8.0$ km/s with a sharp jump at $30$ or $32$ km whereas the S-wave model displays a gradational LVZ between $32$ and $38$ km. The S-wave velocity within the LVZ is $4.0$ km/s, while just above and below it is $4.3$ km/s. Interestingly, the velocity model of Holbrook (1990) displays a P-wave LVZ between $38$ and $42$ km with a P-wave velocity of $7.5$ $\pm$ $0.2$ km/s (Figure 12). This LVZ extends from SP4 to about $30$ km north of SP9 (Figure 1). Holbrook's (1990) LVZ (Figure 12), although located deeper, was also modeled as a gradational boundary due to lack of upper mantle reflections on vertical-incidence reflection data in this region (e.g., Holbrook, 1990; Klemperer, 1986). Holbrook (1990), however, considered the details of his LVZ as tentative because of the uncertainties in modeling LVZ's with the surface refraction method and the evidence for the LVZ from only SP8 records. Conversely, the LVZ in the S-wave model is a required feature with only a $\pm 0.1$ km/s velocity uncertainty based on our sensitivity analysis. Thus our S-wave model combined with Holbrook's (1990) interpretation provide evidence for the presence of an LVZ at these depths in this region. Smooth negative-velocity gradients present in the S-wave velocity model deserve further attention because these gradients occur in the zones of high reflectivity in the crust determined from the PASSCAL and COCORP reflection data (e.g., Holbrook et al., 1991; Catchings, 1992; Klemperer et al., 1986). The first negative-velocity gradient layer between $8$ and $18$ km is nearly coincident with the silicic and intermediate compositional layer and the second negative-velocity gradient layer between $18$ and $24$ km is coincident with the more mafic crustal composition layer both described in Holbrook et al. (1991). These layers are highly reflective (Figure 12). The mafic lower crustal layer between $24$ and $28$ km is largely non-reflective whereas the $4$ km thick $7.4$ km/s layer just below it is also reflective in the velocity model of Catchings and Mooney (1991) (Catchings, 1992). This reflective layer is interpreted as material magmatically underplated onto the crust. The non-reflective layer is not coincident with a negative-velocity gradient layer whereas the lowermost reflective layer is partly coincident with the beginning of the LVZ in the S-wave model (Figure 12). Holbrook et al. (1991) proposed a model in which the middle crustal reflectivity is caused by ductile shear strain fabric and the lower crustal reflectivity is caused by the mafic and ultramafic intrusions. They excluded the presence of a fluid phase as the possible cause of reflectivity because of high temperature gradients in the Basin and Range. In response to this model, Hyndman et al. (1991) argued that layered fluid porosity was a viable hypothesis for the origin of reflections. Below, we provide evidence which indicates that the origin of the reflectivity may indeed be the presence of pore fluids in the middle crust. Spencer and Nur (1976) reported on the behavior of the P- and S-wave velocities at high temperature ($\leq$ $500$$^o$C), and under high confining pressure ($1$ kbar), and independently controlled pore water pressure using samples from Westerly granite. Their experiments showed that increasing temperature of a sample that contains pore water at low pressure ($100$ bar) causes a continuous decrease in the bulk modulus (k) and thus a decrease in the P-wave velocity (Vp) while the shear modulus $\mu$ and the S-wave velocity (Vs) are unaffected. This results in a decrease in Vp and nearly no change in Vs and thus causes a decrease in Poisson's ratio. This result was used by Holbrook et al. (1988) to support the observation of a P-wave LVZ that was not present in the S-wave data in southwest Germany. Under the same pressure and temperature conditions but high pore fluid pressure ($1$ kbar), Spencer and Nur (1976) showed, however, an opposite effect which resulted in a decrease in $\mu$ and Vs, while Vp remained relatively unchanged. Although the depth range where the negative-velocity gradients occur in our model, would be under higher pressures ($3$-$6$ kbar) and temperatures ($300$-$600$$^o$C) than tested for by Spencer and Nur (1976), we believe that their results may be taken as representative behavior of Vp and Vs for higher confining pressures. \subsection{Pore Fluids} Pore fluids within the continental crust may be present due to progressive metamorphism resulting in reactions of dehydration and mantle devolatilization in regions of extension at great depths (e.g., Christensen, 1989; Hyndman and Klemperer, 1989; Hyndman and Shearer, 1989). The depth and lateral extent of pore fluids in the continental crust are essentially unknown (Christensen, 1989), however, the conductivity models for the continental crust in various regions, particularly in geologically younger areas, indicate the presence of conductive lower crustal layers which are coincident with deep seismic reflections. The conductive lower crustal layers are best explained by including a few percent free saline water in interconnected pores (e.g., Hyndman and Shearer, 1989; Klemperer, 1987; Jones, 1987). Based on these findings, we propose a mechanism in which continuously increasing pore fluid pressure with depth nearly balances the lithostatic pressure along with the combined effects of pressure and temperature to create the conditions reported by Spencer and Nur (1976) resulting in a negative-velocity gradient in Vs and no change in Vp in a given crustal layer. If this mechanism is true then it may explain the negative-velocity gradients in the S-wave model with no corresponding change in the P-wave models (Figure 12). Therefore, this mechanism may also explain the origin of the upper highly reflective layer coincident with the negative S-wave velocity gradients in our model (assuming a layered, fluid-filled, porosity model, e.g., Hyndman and Shearer, 1989). The non-reflective lower crustal layer ($24$-$28$ km) may also be explained by this mechanism because this layer is coincident with a normal S-wave velocity layer where pore fluids and porosity may be significantly low. The LVZ that is present in the uppermost mantle in the S-wave model is partially coincident with the magmatically underplated layer of Catchings and Mooney (1991) and the LVZ of Holbrook (1990). Thus, our LVZ may represent evidence for the presence of a partial melt that could be considered the source material for magmatic underplating. Humphreys and Dueker (1994) also interpreted relatively low P-wave velocity upper mantle as hot and partially molten ($1$-$3$\% partial melt) in this region. Jarchow (1991) identified a magma body referred to as the 'Buena Vista Magma Body' at the base of the crust beneath Buena Vista Valley near SP9 (Figure 1) from P to S converted reflections. According to his model, the Buena Vista Magma Body is small in areal extent ($< 1.8$ km) and in thickness ($< 200$ m). He determined an age of no older than 500,000 years and a melt fraction that is greater than a few tens of percent. Given the small dimensions of the magma body, it is not possible to establish any relationship with the LVZ in our model. However, these results provide evidence for the presence of magma bodies that have larger dimensions in the northwestern Basin and Range. More detailed modeling of our LVZ is needed for a better understanding of its composition and thermal state. \subsection{Thick and Fast Crust} In central Nevada, the crustal thickness increases from $32$ km to $36$-$37$ km between SP6 and SP7 (Catchings, 1992) (Figure 1). In the velocity model of Catchings and Mooney (1991), this increase is due to the thickening of upper crustal and consequent deepening of middle and lower crustal layers. Holbrook's (1990) velocity model achieves this thickening by including a lower crustal, high-velocity layer ($7.0$-$7.2$ km/s), which thickens to the east. This region of thicker crust also produces highly attenuated, weak Pn arrivals (Holbrook, 1990; Catchings and Mooney, 1991). Both Holbrook (1990) and Catchings and Mooney (1991) have also shown that the large $95$ mgal decrease in Bouguer gravity anomaly from west to east can be explained by the eastward dip of the Moho in this region (Figure 1). MLC ($\sim$ $52$ km south of SP7) is located in this region of crustal thickening and our S-wave velocity model, with a $38$ km thick crust confirms this thickening remarkably well. The upper crustal LVZ (5-12 km) in the MLC velocity model (Figure 8) also shows a good depth correlation with the LVZ present between SP6 and SP7 (3-8 km) in the velocity model of Holbrook (1990). The westernmost extent of Precambrian North American crust is inferred to lie in central Nevada on the basis of isotopic studies by Farmer and De Paolo (1983). The dashed-white line in Figure 1 is the approximate location of the ratio of $^8$$^7$Sr/$^8$$^6$Sr= $0.708$ and $\epsilon_{ND}$ = $-7$ from the isotopic study of Farmer and De Paolo (1983). The region east of this line is underlain by Precambrian sialic normal continental material whereas the region to the west is underlain by transitional or oceanic material. Farmer and De Paolo (1983) pointed out that this boundary could mark an increase in overall crustal thickness and a compositional change in central Nevada. Catchings and Mooney (1991) and Catchings (1992) interpreted the thick crust they observed in central Nevada as being consistent with that of normal continental areas. In addition to this interpretation, Holbrook (1990) pointed out that it may be difficult to distinguish between a normal continental crust and magmatically intruded and underplated crust during Cenozoic extension based on seismic evidence alone. The distinct features of the MLC velocity model are (1) higher crustal and upper mantle velocities, (2) a thick crust (38 km), (3) highly gradational crust-mantle transition. The highly gradational nature of the crust-mantle transition may also account for the weak Pn arrivals observed on the PASSCAL lines in this region. Thus, we believe that our velocity structure together with the PASSCAL results, the lowest gravity and heat flow anomalies, and relatively older volcanism in central Nevada present a body of geophysical evidence that the crust in this region has been less affected by Cenozoic extension than the region to the north and west. The velocity structure at GAR is very similar to that found at MLC. It has nearly the same average crustal and mantle velocities, and crustal thickness ($38$ km). To the north of MLC and GAR, at RTS, the crust is only $32$ km thick. This station is nearly coincident with the reversed, seismic refraction line between the Nevada Test Site (NTS) and Boise, Idaho analyzed by Hill and Pakiser (1967). This line runs in a north-south direction and is shown with a white solid line (Figure 1). From the reversed refraction shots, between Eureka, Nevada ($\sim$ $40$ km southeast of RTS) and Mountain City, Nevada (at the northern end of our map), Hill and Pakiser (1967) found a $32$ to $33$ km thick crust with a Pn velocity of $7.9$ km/s in this region. This is in excellent agreement with our estimate of $32$ km thickness and upper mantle shear velocity of $4.58$ km/s, which predicts a $7.93$ km/s Pn velocity if we assume a Poisson's ratio of 0.25, at RTS. This result also implies that the velocity structure at RTS may be more reliable than we thought despite the relatively low quality of the receiver function. The analysis of Hill and Pakiser (1967) also included nuclear and chemical explosions from NTS. From these unreversed explosions, they found a local delay of Pn arrivals ($0.2$ to $0.5$ s) in the region between just $40$ and $90$ km south of RTS. We also noticed earlier a $0.5$ s delay in the arrival time of the Moho Ps phase from the SE relative to the NW back-azimuth, indicating a $4$ km thicker crust to the south. Hill and Pakiser (1967) offered two alternative explanations for the cause of Pn delays, either an LVZ in the upper $20$ km of the crust or a deep depression in the Moho reaching $36$ km. Although they preferred the former explanation, we prefer the latter because this is consistent with possibly thicker crust to the SE of RTS and with the inferred crustal thickening at nearby MLC (Figure 1). \subsection{Comparison with Gravity} MLC, GAR, and RTS lie in a broad region of gravity low and high elevations with a bilateral symmetry in regional gravity and topography. This bilateral symmetry with a north trending axis (running at $\sim$ long $115.5$$^o$W) displays a butterfly gravity pattern as described by Eaton et al. (1978). The approximate outline of this pattern is shown with a black solid line in Figure 1. MLC and GAR are located within opposite sides of the Butterfly gravity pattern and RTS is near the axis of symmetry. We believe that the nearly identical average crustal and upper mantle velocity structure, and thus the same subsurface mass distribution, beneath MLC and GAR is directly reflected in this gravity pattern as the paired thickening of the crust. The symmetry of the crustal structure at MLC and GAR leads to our earlier interpretation that the crust beneath GAR is also less affected by Cenozoic extension at depth. The crustal structure at RTS is less well constrained but it would be consistent with the low gravity anomaly of this region if the possible LVZ's indicated by the NW receiver function and the possibly thicker crust to the SE were taken into consideration. Outside the boundaries of the Butterfly gravity anomaly, WCP to the north and KVN to the east, indicate that the crust thins to $32$ km with decreasing average crustal velocities, Bouguer gravity anomalies ($-177$ mgal) and elevations. Various refraction studies have found a $31$ km thick crust in the vicinity of WCP (Smith et al., 1989). The PASSCAL experiment indicates a $32$ to $33$ km thick crust near KVN at SP11 (Figure 1) (e.g., Holbrook, 1990; Catchings and Mooney, 1991). The refraction and receiver function results again agree well in both regions. Further south, at MNA, the crust is only $2$ km thicker than at KVN. The velocity model at this station is well-constrained by the high-quality receiver functions and short- and long-period surface-wave phase velocities. The velocity model of Mangino et al. (1993) is similar to this model, except that our model has a sharper crust-mantle transition and no velocity gradients. Mangino et al. (1993) interpreted the relatively thicker crust in this region as representative of a transition zone between the typical Basin and Range Province crust to the east and Sierran crust to the west. This interpretation is supported by our results and consistent with the gravity anomalies in this region. This region is lower in elevation ($\sim$ $1.5$ km) and associated with a moderate gravity anomaly ($-190$ mgal) bounded by the two gravity lows ($ < -220$ mgal ), high elevations, and the thick crust of the Sierra-Nevada to the west and central Nevada to the east (Figure 1). The crustal structures derived from the shear-wave velocity models have shown a good spatial correlation with the variations in gravity and elevation throughout the province, indicating both regional and local isostatic compensation. To further quantify this apparent compensation, we perform a simple, rough check to see the extent to which these models are in isostatic balance. We follow the simple approach described by Hill (1978) and compute a gravity anomaly for each velocity model using densities predicted from our average crustal and upper mantle shear velocities using the same velocity-density relationships as in the inversions. In this over-simplified case, the gravity anomaly can be represented as \begin{equation} \Delta g = 2 \pi G [M - \bar{M}], \hspace{0.2in} M=\rho_{c} H_{c} + \rho_{m} (D_{c} - H_{c})\ \end{equation} where $\Delta g$ is the predicted gravity value due to local mass anomalies in each model, G is the gravitational constant ($6.67$x10$^-$$^8$ cgs), $\bar{M}$ is the average mass per unit area for all models, and M is the average mass per unit area in a specific model and is computed from the knowledge of average crustal density $\rho_{c}$, mantle density $\rho_{m}$, the crustal thickness $H_{c}$, and the depth of compensation $D_{c}$. We take $D_{c}$=$40$ km, which is just deeper than the base of the thickest crust in all velocity models. If all the models are in perfect isostatic balance then M=$\bar{M}$ and $\Delta g$ will be zero for each model. Hill (1978) pointed out that an error of $\pm 0.1$ gr/cm$^3$ in the estimated density over a $10$ km depth interval in a specific model will result in a gravity anomaly of $\pm 42$ mgal. Figure 13 shows a plot of the predicted gravity anomalies for each model. The predicted values deviate from zero by $\pm 52$ mgal, which is only $10$ mgal higher than the expected error given above, indicating that the models are in near isostatic balance. We also plotted the mean removed observed Bouguer gravity values at each station for comparison with the predicted values. We use the mean removed Bouguer gravity values because our purpose here is to compare relative variations in gravity. The maximum difference between the observed and predicted anomalies is less than $30$ mgal for all stations except for RTS, where a $72$ mgal anomaly mismatch occurs (Figure 13). As mentioned earlier, the NW RTS receiver function displays evidence for LVZ's in the crust and upper mantle which remain relatively unmodeled. This may result in an overestimation of densities leading to the mismatch we observe. Further, the SE receiver function indicates a $36$ km thick crust at RTS, which can also account for the mismatch. Thus, our models can account for the variations in the observed Bouguer gravity in the NBRP, implying that the isostatic compensation occurs within $40$ km depth of the surface (i.e. confined to the crust), and is close to a classical Airy type of compensation. \section{Conclusions} Our results lead to the following conclusions: 1) The constrained inversion of receiver functions by surface wave data largely minimizes the non-uniqness problem caused by the velocity-depth trade-off. 2) This method takes advantage of average velocity information present in the surface-wave method and differential velocity information in the receiver function method. Long and short period surface wave phase velocity data are essential constraints to obtain smooth and reliable shear velocity models. 3) Using this method results in well-determined estimates of the crustal and upper mantle shear-velocity structures which are in excellent agreement with the wide angle refraction and reflection and earthquake tomographic compressional-wave structure, gravity, heat flow, and elevation data in the NBRP. The crustal thickness, average crustal and upper mantle velocities are the most robust parameters determined from the joint inversions. 4) These parameters show excellent correlations with the 1986 PASSCAL, wide angle reflection and refraction experiment results in northwestern and central Nevada, and other refraction and earthquake tomographic experiment results that have been used to infer the crustal and upper mantle structure of the province. 5) An unusually thick ($38$ km) and relatively faster crust and upper mantle are found in central and eastern Nevada compared to the thin ($28$ to $34$ km) and slower crust and upper mantle of the western Basin and Range. We interpret the regions of thicker and faster crust and upper mantle as zones that have undergone less Cenozoic extension relative to the surrounding regions to the west and north. 6) A simple gravity model shows both local and regional isostatic compensation occur within $40$ km of the surface, indicating a near classical Airy type of compensation in the province. 7) The analysis of the S-wave velocity model derived from the receiver functions at station BMN and P-wave velocity models derived from the 1986 PASSCAL experiment in northwestern Nevada shows the presence of negative-velocity gradients in the S-wave model with no corresponding velocity decrease in the P-wave models between $10$ and $24$ km depth range. This combined velocity model may be explained by high pore fluid pressures at these depths. This model favors layered fluid porosity models proposed to explain extensive middle to lower crust continental seismic reflections and high electrical conductivity. 8) An upper mantle, shear-wave LVZ is found between $32$ and $38$ km depths in northwestern Nevada. This upper mantle shear-wave LVZ is consistent with partial melt, which may be considered as the source material for magmatic underplating in this region. \acknowledgments This work was funded by National Science Foundation grant EAR-9319007 and EAR-9319008. The first author acknowledges financial support by the Turkish Ministry of Education. We thank K. D. Smith for providing substantial help with using ERMapper software to generate Bouguer gravity map and also thank Y. Zeng for making his surface-wave computation codes available to us. Parts of the inversion code were provided by C. J. Ammon. Discussions with G. P. Biasi, K. F. Priestley, G. Shields, and Y. Zeng helped to improve our method and interpretations. C. Jones helped in field and provided an early review. A. \c{C}etinta\d{s} provided voluntary help during the portable field experiment. \begin{references} \reference Aki, K., and P. G. Richards, Quantitative Seismology: Theory and Methods, W. H. Freeman and Company, New York, pp. 913, 1980. \reference Allmendinger, R. W., T. A. Hauge, E. C. Hauser, C. J. Potter, S. L. Klemperer, K. D. Nelson, P. Kneupfer, and J. Oliver, Overview of the COCORP $40^o$N transect, western U.S.A.: The fabric of an orogenic belt, Geol. Soc. Am. Bull., 98, 308-319, 1987. \reference Ammon, C. J., G. E. Randall, and G. Zandt, On the Nonuniqueness of Receiver Functions, J. Geophys. Res., 95, 15,303-15,318, 1990. \reference Ammon, C. J., The Isolation of Receiver Effects From Teleseismic P Waveforms, Bull. Seis. Soc. Am., 81, 2504-2510, 1991. \reference Ammon, C. J., and G. Zandt, Receiver Structure Beneath the Southern Mojave Block, California, Bull. Seis. Soc. Am., 83, 737-755, 1993. \reference Benz, M., R. B. Smith, and W. D. Mooney, Crustal Structure of the Northwestern Basin and Range Province From the 1986 Program for Array Seismic Studies of the Continental Lithosphere Seismic Experiment, J. Geophys. Res., 95, 21,823-21,842, 1990. \reference Blackwell, D. D., Heat flow and energy loss in the western United States, in: R. B. Smith and G. P. Eaton (Editors), Cenozoic Tectonics and Regional Geophysics of the Western Cordillera, Mem. Geol. Soc. Am., 152, 175-208, 1978. \reference Brune, J. N., J. Nafe, and J. Oliver, A simplified method for the analysis and synthesis of dispersed wave trains, J. Geophys. Res., 65, 287-304, 1960. \reference Cassidy, J. F., Numerical experiments in broadband receiver function analysis, Bull. Seis. Soc. Am., 82, 1453-1474, 1992. \reference Cassidy, J. F., A comparison of the receiver structure beneath stations of the Canadian National Seismogram Network, Can. J. Earth Sci., 32, 938-951, 1995. \reference Catchings R. D. and W. D. Mooney, Basin and Range Crustal and Upper Mantle Structure of Northwest to Central Nevada, J. Geophys. Res., 96, 6247-6267, 1991. \reference Catchings R. D., A relation among geology, tectonics, and velocity structure, western to central Nevada Basin and Range, Geol. Soc. Am. Bull., 104, 1178-1192, 1992. \reference Christensen, N. I., Pore pressure, seismic velocities, and crustal structure, in: L. C. Pakiser and W. D. Mooney (Editors), Geophysical Framework of the Continental United States, Mem. Geol. Soc. Am., 172, 783-798, 1989. \reference Eaton, G. P., R. R. Wahl, H. J. Prostka, D. R. Mabey, and M. D. Kleinkof, Regional gravity and tectonic patterns: Their relation to late Cenozoic epirogeny and lateral spreading in the western Cordillera, in: R. B. Smith and G. P. Eaton (Editors), Cenozoic Tectonics and Regional Geophysics of the Western Cordillera, Mem. Geol. Soc. Am., 152, 51-91, 1978. \reference Farmer, G. L., and D. J. DePaolo, Origin of Mesozoic and Tertiary Granite in the Western United States and Implications for Pre-Mesozoic Crustal Structure, 1. Nd and Sr Isotopic Studies in the Geocline of the Northern Great Basin, J. Geophys. Res., 88, 3379-3401, 1983. \reference Hawman R. B., R. H. Colburn, D. A. Walker, and S. B. Smithson, Processing and Inversion of Refraction and Wide-Angle Reflection Data from the 1986 Nevada PASSCAL Experiment, J. Geophys. Res., 95, 4657-4691, 1990. \reference Hearn, T. M., N. Beghoul, and M. Barazangi, Tomography of the Western United States From Regional Arrival Times, J. Geophys. Res., 96, 16,369-16,381, 1991. \reference Hill, D. P., and L. C. Pakiser, Seismic-Refraction Study of Crustal Structure between the Nevada Test Site and Boise, Idaho, Geol. Soc. Am. Bull., 78, 685-704, 1967. \reference Hill, D. P., Seismic evidence for the structure and Cenozoic tectonics of the Pacific Coast States, in: R. B. Smith and G. P. Eaton (Editors), Cenozoic Tectonics and Regional Geophysics of the Western Cordillera, Mem. Geol. Soc. Am., 152, 145-174, 1978. \reference Holbrook, W. S., D. Gajewski, A. Krammer, and C. Prodehl, An Interpretation of Wide-Angle Compressional and Shear Wave Data in Southwest Germany: Poisson's Ratio and Petrological Implications, J. Geophys. Res., 93, 12,081-12,106, 1988. \reference Holbrook, W. S., The Crustal Structure of the Northwestern Basin and Range Province, Nevada, from Wide-Angle Seismic Data, J. Geophys. Res., 95, 21,843-21,869, 1990. \reference Holbrook, W. S., R. D. Catchings, and C. M. Jarchow, Origin of deep crustal reflections: Implications of coincident seismic refraction and reflection data in Nevada, Geology, 19, 175-179, 1991. \reference Humphreys E. D., and K. G. Dueker, Physical state of the western U. S. upper mantle, J. Geophys. Res., 99, 9635-9650, 1994. \reference Hyndman R. D., and S. L. Klemperer, Lower-crustal porosity from electrical measurements and inferences about composition from seismic velocities, Geophys. Res. Lett., 16, 255-258, 1989. \reference Hyndman R. D., and P. M. Shearer, Water in the lower continental crust: modeling magnetotelluric and seismic reflection results, Geophys. J. Int., 98, 343-365, 1989. \reference Hyndman R. D., T. J. Lewis, and G. Markus, Comment on origin of deep crustal reflections: implications of coincident seismic refraction and reflection data in Nevada, Geology, 19, 1243, 1991. \reference Jarchow, C. M., and G. A. Thompson, The Nature of Mohorovicic Discontinuity, Ann. Rev. Earth Planet. Sci., 17, 475-506, 1989. \reference Jarchow, C. M., Investigations of magmatic underplating beneath the northwestern Basin and Range province, Nevada; Seismic data acquisition and tectonic problems of the Columbia Plateau, Washington; and the nature of the Mohorovicic discontinuity worldwide, Ph.D. thesis, Stanford University, Palo Alto, California, pp. 258, 1991. \reference Jones, A. G., MT and reflection: an essential combination, Geophys. J. R. Astr. Soc., 89, 7-18, 1987. \reference Kind, R., G. L. Kosarev, and N. V. Petersen, Receiver functions at the stations of the German Regional Seismic Network (GRSN), Geophys. J. Int., 121, 191-202, 1995. \reference Kennet, B. L. N., Seismic wave propagation in stratified media, Cambridge University Press, New York, pp. 342, 1983. \reference Klemperer, S. L., T. A. Hauge, E. C. Hauser, J. E. Oliver, and C. J. Potter, The Moho in the northern Basin and Range province, Nevada, along the COCORP $40^o$N seismic-reflection transect, Geol. Soc. Am. Bull., 97, 603-618, 1986. \reference Klemperer, S. L., A relation between continental heat flow and the seismic reflectivity of the lower crust, J. Geophys., 61, 1-11, 1987. \reference Lachenbruch, A. H., and J. H. Sass, Models of an extending lithosphere and heat flow in the Basin and Range province, in: R. B. Smith and G. P. Eaton (Editors), Cenozoic Tectonics and Regional Geophysics of the Western Cordillera, Mem. Geol. Soc. Am., 152, 209-250, 1978. \reference Langston, C. A., Structure under Mount Rainer, Washington, inferred from teleseismic body waves, J. Geophys. Res., 84, 4749-4762, 1979. \reference Mangino, S. G., G. Zandt, and C. J. Ammon, The Receiver Structure Beneath Mina, Nevada, Bull. Seis. Soc. Am., 83, 542-560, 1993. \reference McNamara, D. E., and T. J. Owens, Azimuthal Shear Wave Velocity Anisotropy in the Basin and Range Province Using Moho Ps Converted Phases, 98, 12,003-12,017, 1995. \reference Owens, T. J., G. Zandt, and S. R. Taylor, Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee, J. Geophys. Res., 89, 7783-7795, 1984. \reference Owens, T. J., and G. Zandt, The Response of the Continental Crust-Mantle Boundary Observed on Broadband Teleseismic Receiver Functions, Geophys. Res. Lett., 12, 705-708, 1985. \reference Owens, T. J., Crustal Structure of the Adirondacks Determined From Broadband Teleseismic Waveform Modeling, \reference \"{O}zalaybey, S., Seismic velocity structure in the western U. S. from: shear-wave splitting and receiver functions of teleseismic earthquakes, Ph.D. thesis, University of Nevada, Reno, NV, pp. 171, 1996. \reference Pakiser, L. C., Geophysics of the Intermontane system, in: L. C. Pakiser and W. D. Mooney (Editors), Geophysical Framework of the Continental United States, Mem. Geol. Soc. Am., 172, 57-95, 1989. \reference Patton, H. J., and S. R. Taylor, Q structure of the Basin and Range from surface waves, J. Geophys. Res., 89, 6929-6940, 1984. \reference Peppin, W. A., and W. F. Nicks, Real-Time Analog and Digital Data Acquisition through CUSP, Seismol. Res. Lett., 63, 181-189, 1992. \reference Potter, C. J., C. Liu, J. Huang, L. Zheng, T. A. Hauge, E. C. Hauser, R. W. Allmendinger, J. E. Oliver, S. Kaufman, and L. Brown, Crustal structure of north-central Nevada: Results from COCORP deep seismic profiling, Geol. Soc. Am. Bull., 98, 330-337, 1987. \reference Priestley, K. F., and J. N. Brune, Surface waves and the structure of the Great Basin of Nevada and western Utah, J. Geophys. Res., 83, 2265-2272, 1978. \reference Priestley, K. F., J. A. Orcutt, and J. N. Brune, Higher mode surface waves and the structure of the Great Basin of Nevada and western Utah, J. Geophys. Res., 85, 7166-7174, 1980. \reference Priestley, K. F., and J. N. Brune, Shear Wave Structure of the Southern Volcanic Plateau of Oregon and Idaho and the Northern Great Basin of Nevada From Surface Wave Dispersion, J. Geophys. Res., 87, 2671-2675, 1982. \reference Prodehl, C., Crustal structure of the western United States, U. S. Geol. Surv. Prof. Pap., 1034, pp. 74, 1979. \reference Sheehan A. F., C. H. Jones, K. G. Dueker, M. K. Savage, and S. \"{O}zalaybey, Colorado Plateau - Great Basin PASSCAL Seismic Deployment 1994-1995, in: IRIS - 2000, Section II, Scientific Contributions from the Proposal: A Science Facility for Studying Dynamic of the Solid Earth submitted to the National Science Foundation by IRIS, 1995. \reference Smith, R. B., W. C. Nagy, K. A. Julander, J. Viveiros, C. A. Barker, and D. G. Gants, Geophysical and of tectonic framework of the eastern Basin and Range-Colorado Plateau-Rocky Mountain transition, in : L. C. Pakiser and W. D. Mooney (Editors), Geophysical Framework of the Continental United States, Mem. Geol. Soc. Am., 172, 205-233, 1989. \reference Spencer, J. W., and A. M. Nur, The Effects of Pressure, Temperature, and Pore Water on Velocities in Westerly Granite, J. Geophys. Res., 81, 899-904, 1976. \reference Taylor S. R., and H. J. Patton, Shear-velocity structure from regionalized surface-wave dispersion in the Basin and Range, Geophys. Res. Lett., 13, 30-33, 1986. \reference Taylor S. R., B. P. Bonner, and G. Zandt, Attenuation and Scattering of Broadband P and S Waves Across North America, J. Geophys. Res., 91, 7309-7325, 1986. \reference Zandt, G., and T. J. Owens, Comparison of crustal velocity profiles determined by seismic refraction and teleseismic methods, Tectonophysics, 128, 155-161, 1986. \section{Figure Captions} \noindent{\bf Figure 1.} Bouguer gravity anomaly map of the northern Basin and Range Province plotted with crustal thickness (km), mean crustal shear velocity (km/s), and upper mantle shear velocity (km/s), respectively, shown under each station symbol. See text for interpretations and explanations. \noindent{\bf Figure 2.} Stacked radial (upper traces) and transverse (lower traces) receiver functions for three back-azimuth ranges from stations DNY, WHR, KVN, and MNA. The shaded region indicates one standard deviation bounds obtained from the stacking process. The mean receiver functions are shown with solid lines. Ray parameter p, and the number of events in computing each stack are indicated. All are plotted in the same scale. \noindent{\bf Figure 3.} As Figure 2, but for southeast back-azimuth at three distance ranges for stations BMN, WHR, and MNA. Primary and multiple arrivals are less energetic for the most distant range. \noindent{\bf Figure 4.} As Figure 2, but for portable stations. \noindent{\bf Figure 5.} Rayleigh-wave phase velocity measurements along three paths: Paths 1 and 2 are from a teleseismic event in the Tonga Islands and Path 3 from a regional event in Oregon. \noindent{\bf Figure 6.} (a) Velocity models from the inversion constrained by the three pairs of phase velocity (shown as observed) joined with curves of the corresponding line type shown in (c). (b) Rayleigh wave partial derivative curves computed for the initial model. (c) Dispersion curves corresponding to each velocity model. (d) Computed synthetic waveforms for each model shown above. \noindent{\bf Figure 7.} Inversion of MNA receiver functions from Ranges 1, 2, and 3 and Rayleigh wave phase velocities. (a) Final shear wave velocity model. (b) Rayleigh wave dispersion fit; observed phase velocities are from Path 3 ($10$ s period) and Path 2 of MNA-MLC station pair. Theoretical dispersion curve for the velocity model shown on the left fits this data well. (c) Receiver function waveform fit; the shaded region represents one standard deviation bounds for the stacked receiver functions. The solid lines indicate modeled receiver functions computed for the final velocity model. The arrow indicates the Moho. \noindent{\bf Figure 8.} As Figure 7, but for station MLC. Phase velocity data are the same as shown at MNA. \noindent{\bf Figure 9.} As Figure 7, but for station BMN. Receiver functions are from Ranges 1 and 2. Observed phase velocities are from Path 1 and stations WCN, WHR, and BMN. \noindent{\bf Figure 10.} Shear-wave velocity models for the UNR permanent stations. The shaded regions indicate confidence bounds obtained from the sensitivity analysis. A single labeled region identified as regions of 1, 2, 3, or 4 between thin arrows can be positively or negatively perturbed within these bounds without affecting the data fits. The bold arrow indicates the Moho. \noindent{\bf Figure 11.} As Figure 10, but for portable stations. \noindent{\bf Figure 12.} Comparison of P-wave velocity models derived from the 1986 PASSCAL experiment (near SP9) with S-wave velocity model derived from the receiver functions at station BMN. P-wave models are from Catchings and Mooney (1991) (solid) and Holbrook (1990) (dashed). Note excellent agreement in the location of velocity discontinuities down to the crust-mantle transition (CMT) between the S-wave and P-wave models and negative velocity gradients between $8$ and $24$ km in the S-wave model. An S- and P-wave low-velocity zone (LVZ) are present in the upper mantle. The arrow indicates the Moho. \noindent{\bf Figure 13.} A gravity anomaly is computed for each density structure shown above assuming a $40$ km depth of compensation. The density values (in cgs units) are computed from the average shear velocities at each station. The predicted and observed values show a good correlation except at RTS.