The composition of the crust is one of the most uncertain and controversial components of continental estimates due to (1) limited direct access and (2) inconsistent indirect assessments. Here we show that by combining high-resolution shear velocity (Vs) models with newly measured with newly measured ratio of compressional wave velocity (Vp) and Vs, or Vp/Vs ratio, for the crystalline crust, a 3-D composition (SiO2 wt%) model of the continental crust can be derived with uncertainty estimates. Comparing the model with local xenolith data shows consistency at mid and lower crustal depths. The spatial patterns in bulk SiO2 content correlate with major geological provinces, including the footprints of Cenozoic and Mesozoic mafic volcanism in the western U.S., and offer new insight into the composition and evolution of the continental U.S.

The bulk silica content of the deep continental crust and its spatial variation play a foundational role in studying the growth and evolution of the continents (Condie, 1993; Gao et al., 1998; Rudnick and Fountain, 1995). Since crustal rocks with more felsic compositions bear higher heat-producing elements and other compatible elements (Hasterok et al., 2018), better compositional constraints will benefit (1) mass-balance geochemical calculations (Rudnick and Gao, 2014), (2) inferences about the thermal structure, and (3) identification of possible source regions of geoneutrino signals (Huang et al., 2013). Silica content also influences lithospheric strength and thus seismicity distributions (Lowry and Smith 1995). However, even the average concentration of this most abundant oxide is still under debate due to the lack of direct access to the deep crust (Hacker et al., 2015). Previous global average estimates, for example, span a range for the bulk crust from mafic-intermediate (i.e., ~57 wt% SiO2) to intermediate-felsic (~64 wt% SiO2) (McLennan et al., 2005), with the lower crust being the most uncertain.

The difficulties of indirectly determining silica content and its heterogeneity in the deep crust are reduced but not fully resolved using seismic properties. As a physical attribute influenced by composition, P velocity (Vp) of the crust has been used to place constraints on the chemical composition (Rudnick and Fountain, 1995; Christensen and Mooney, 1995; Christensen, 1996). However, this method suffers from either limited spatial coverage when using local active source seismic properties or limited resolution when using a global seismic model (e.g., CRUST 1.0, Laske et al., 2013). The ratio between compressional velocity (Vp) and shear velocity (Vs), or Vp/Vs ratio, a quantity that is directly related to Poisson’s ratio, has also been widely used to infer the silica content qualitatively (e.g., Lowry and Pérez-Gussinyé, 2011; Ma and Lowry, 2017) as quartz has particularly low Vp/Vs. A more recent effort used the crustal Vs and thermodynamically calculated petrophysical properties to quantify deep crustal composition, but it is only applied regionally to the southwestern U.S. (Sammon et al., 2020). Using Vp or Vp/Vs alone to quantify composition is difficult since neither provides deterministic constraints (Hacker et al., 2015). Additionally, combining Vp and Vp/Vs to infer the composition is also challenging in the arcs near subduction zones (or hot lower continental crust, as discussed later in this paper), where the pressure and temperature lie near the α-β quartz-phase transition zone (Shillington et al., 2013; Jagoutz and Behn, 2013). As a result, quantifying the silica content of the deep crust using seismic properties, or even mapping it at a scale in which tectonic conclusions can be drawn, remains a challenge.

In this work, we exploit the observation that SiO2 content for major crustal rocks (including amphibolites and mafic and felsic granulites, which are the main constituents of the deep crust) presents a monotonically varying spectrum in Vs-Vp/Vs space (Lowry and Pérez-Gussinyé, 2011; Fig. 1). Additionally, we measure Vp/Vs for the crystalline crust across the continental U.S., combine the result with recently published Vs models to provide quantified silica content and build a 3-D compositional model including quantitative uncertainties. Notably, we show that the model reveals a generally mafic deep crustal composition for the continental U.S., especially in the tectonically stable central and eastern parts. This observation is in accordance with deep crustal xenoliths collected locally and globally. Our findings put the chemical composition of deep crust on the more mafic end of the spectrum of previously estimated global models. We also provide a comprehensive evaluation of the sources of random and systematic errors in this analysis. From this, we identify that the uncertainty in the mineral phase equilibrium may be the major source of errors in this type of analysis, and can be as high as 5%–7% for some cases.

The manuscript is organized as follows: in Section 2, we first introduce the data and methods. Particularly in this part, we introduce how the Vs-Vp/Vs-SiO2 wt% interpolation surfaces are created from lab-measured/thermodynamically calculated rock properties (Section 2.1), and how the Vs is corrected to the same physical condition of the lab measurements/thermodynamic calculations (Section 2.2); The methods introduced in this section include: (1) an updated sequential H-κ (H for thickness and κ for Vp/Vs ratio) stacking method for receiver function (Section 2.3); (2) quantifying SiO2 wt% from Vs and Vp/Vs and a synthetic test to show its feasibility to distinguish crust with different compositions (Section 2.4). In Section 3, we present the major results of this work. These include a new map of the Vp/Vs for the crystalline crust across the continental U.S. (Section 3.1) and a 3-D compositional model (Section 3.2). Horizontal maps and vertical transects of the composition are exhibited in this section. Section 4 mainly discusses the errors in the results, including nonsystematic (Section 4.1) and systematic errors (Section 4.2). In Section 5, the resulting 3-D compositional model is further compared with local xenoliths (Section 5.1) and regional volcanism (Section 5.2). In particular, we discuss bulk composition differences between the western and central/eastern U.S. and deep crustal composition of the Archean and Proterozoic U.S. (Section 5.3). Then we discuss some caveats of the current method (Section 5.4). Finally, we summarize our conclusions and provide an outlook for further research in Section 6.

In this work, we use the crustal Vs and Vp/Vs ratio to constrain the SiO2 wt% for the continental U.S. using the lab-measured composition-seismic relationships based on rock properties from lab measurements (Christensen, 1996; Hacker et al., 2015) or thermodynamic calculations (Hacker et al., 2015). The relationships allow us to use the crustal Vs and Vp/Vs ratio to constrain the SiO2 wt% for the continental U.S. We use the crustal shear velocity model from Shen and Ritzwoller (2016; S&R 2016 hereafter), which is constrained by ambient noise- and teleseismic earthquake-derived Rayleigh waves, receiver functions, and earthquake-based Rayleigh wave ellipticity (horizontal to vertical) or H/V) ratio) ratio measurements. The usage of ambient noise and H/V ratio help constrain the Vs at an unprecedented resolution, especially at crustal depths. The published model is at 1 Hz in frequency (Shen and Ritzwoller, 2016), and is also influenced by other factors in addition to composition (i.e., temperature, pressure; Rudnick and Fountain, 1995). As a result, it cannot be compared directly to the seismic properties of crustal rocks, and has to be corrected to the condition under which the relationship is obtained. In addition to the details of how the S&R 2016 model is corrected, we also present how Vp/Vs ratios are measured through a 2-layer H-κ stacking method. Finally, we provide technical details of how SiO2 wt% is quantified from the seismic properties with uncertainties. In particular, we present synthetic tests to demonstrate the feasibility of constraining composition using Vs and Vp/Vs.

2.1. Creating the Vs-Vp/Vs-SiO2 wt% Interpolation Surface

Shown in Figure 1, we compile the rock properties (i.e., Vs, Vp/Vs, and SiO2 wt%) from 937 lab measurements (Christensen, 1996; Hacker et al., 2015) and 5754 thermodynamic calculations (Hacker et al., 2015) and correct the seismic velocities to 600 MPa and room temperature. Both lab measurements and thermodynamic calculations show a similar trend: the rocks that have high Vp/Vs and Vs are mafic, and felsic rocks have low Vp/Vs and Vs in general.

Based on this trend of rock properties, we build an interpolation surface of rock Vs, Vp/Vs, and SiO2 wt%. The rock (Vs, Vp/Vs, SiO2) triples are binned and interpolated/extrapolated via the Generic Mapping Tool algorithms (Wessel et al., 2019; see Supplemental Material1). Two smoothed interpolation surfaces describing functions of SiO2 (Vs, Vp/Vs) on regular grids are derived from lab measurements (Fig. 2A) and thermodynamics calculations (Fig. 2B) separately. Both the interpolation surfaces capture the “felsic-intermediate-mafic” trend from low Vs and Vp/Vs to high, which can be potentially used to distinguish the compositions of the Earth’s crust. Misfits between the surface-predicted SiO2 wt% and real rock SiO2 wt% are shown as the misfit distributions in Figures 2C and 2D. We observe that the misfits on average are ~−0.2 wt% with a standard deviation of ~5 wt%. These observations indicate that if accurate in situ measurements of crustal Vp/Vs and Vs can be taken, the SiO2 wt% mapped by the interpolation surfaces will likely identify the chemical composition of silica content with precision small enough to distinguish mafic and felsic rocks.

Beyond showing similar trends, we also notice that surfaces from the calculated rock seismic properties exhibit a different distribution in Vs-Vp/Vs space compared to the lab measurements (corrected to the same pressure-temperature [P-T] condition). Particularly, the calculations result in relatively higher Vs (>3.8 km/s) and lower Vp/Vs (<1.8) compared to the lab measurements. Several possible reasons may contribute to this difference: (1) The accuracy of lab-measured seismic properties might suffer from fractures and/or volumetrically minor alteration phases. (2) The seismic properties from thermodynamic calculations are sensitive to the assumption of a given water content, which might be different from the rocks used in the laboratory. (3) The lab measurements include different igneous and metamorphic rocks formed at a variety of crustal pressure-temperature conditions, and the thermodynamics calculations are for lower crust rocks equilibrated at 650 °C (amphibolites) or 750 °C (granulites) and 1.0 GPa (Hacker et al., 2015). Thus, the surface from lab data is essentially an averaged interpolation surface for the P-T conditions across the entire crust, and the surface from thermodynamics data is more representative for a given physical condition (e.g., the lower crust for the dots shown in Fig. 1). Considering the range of observed Vs-Vp/Vs space that will be shown in sections 2.2 and 3.1, we only show the compositional results based on the interpolation surface from lab measurements in the following sections. The caveat of not using the thermodynamics surface to quantify lower crust SiO2 wt% is discussed in Section 5.4.

2.2. Correcting the Shear Velocity Model

We use the geotherm profiles in the U.S. Geological Survey crustal thermal model (Boyd, 2020) to make the temperature corrections on shear velocities. These geotherm profiles are constrained by geothermal heat flux, surface temperature, and Moho temperature derived from the uppermost mantle Vp tomography with consideration of crustal heat generation. Figure 3A shows the temperature profiles we used for two example locations. A uniform temperature derivative (−2 × 10−4 km s–1 °C–1, Rudnick and Fountain, 1995) is used to correct the Vs to room temperature. For pressure correction, we use a uniform depth-pressure relationship (i.e., 27 MPa/km) and pressure derivative (1 × 10−4 km s–1 MPa–1, Rudnick and Fountain, 1995). Frequency-dependent attenuation also affects seismic speed. In order to compare the Vs model to the lab-measured seismic properties, which are usually measured in MHz (Kern et al., 1996), we correct all Vs to infinite frequency following Equation 1 (Minster and Anderson, 1981):
formula

In which V(∞) represents the Vs at an infinite frequency that can be compared with lab-measured Vs directly, and V(ω) is the velocity at frequency ω (for the model we used, ω = 1 Hz). During the correction, we assume that the crustal attenuation (1/Q) is 1/600 across the study region, and the frequency exponent (α) is 0.15 (Kennett et al., 1995). Further in Section 4.2.3, we discuss possible systematic errors introduced due to choices of the thermal model, depth-pressure relationship, and Q model. Figure 3B shows example Vs profiles before and after the temperature, pressure, and attenuation corrections, and Figure 3C shows the corrected shear velocity of the continental U.S. at 20 km. Using different temperature/pressure/attenuation models may change the corrected seismic speed, but tests show that the change is relatively small.

2.3. Measuring the Crustal Vp/Vs

Crustal Vp/Vs has been measured mainly through receiver functions (Zhu and Kanamori, 2000; Ma and Lowry, 2017). In this work, we take advantage of sequential H-κ stacking and apply it to the USArray Transportable Array (U.S.-TA), which provides a uniform sampling of the crystalline crustal Vp/Vs. Notably, we used P-wave waveforms from events (M > 5.5, distance between 30–103 degrees) collected at 1708 U.S.-TA stations during their deployment periods (~2 years on average). Figure 3C shows the stations we used, which cover the study region regularly with a lateral space of ~70 km.

Receiver functions (RFs), computed from three-component seismograms, show the relative response of Earth structure near the receiver. The H-κ stacking method is a widely used RF-stacking method to reduce the random uncertainties for a better estimation of Moho depth and crustal Vp/Vs (Zhu and Kanamori, 2000). It involves stacking of the amplitudes of P wave converted phases (including Moho-converted P-to-s (Ps) phase and multiples, see Fig. 4A top) at discontinuities (i.e., sediment bottom and Moho) within the crust. Sediments usually have higher Vp/Vs than the crystalline crust (Brocher 2005). We adopt a 2-layer sequential H-κ stacking method to calculate the crystalline crust Vp/Vs (Yeck et al., 2013). This method allows Moho-converted phases to separate from sediment conversions and reduces bias in the resulting Vp/Vs of the crystalline crust.

Before 2-layer H-κ stacking, two sets of receiver function waveforms with different Gaussian parameters were calculated separately: a Gaussian parameter of 5 (pulse width ~1.0 s) for the sedimentary layer, and a Gaussian parameter of 2.5 (pulse width ~0.75 s) for the crystalline crust layer. A larger Gaussian parameter will generate receiver function waveforms with a higher frequency, which are more sensitive to the shallow structure. After separating sediment-generated and Moho-generated phases, we stack their corresponding energies to obtain the thickness and Vp/Vs for each layer (sediment and crystalline crust). Shown in Figure 4B, stacked energy maps for station M26A show a sedimentary layer thickness of ~3 km and Vp/Vs of ~2.5. We also tested the stacking for individual waveforms, and the result indicates that earthquakes from different back azimuths and distances provide similar results (Fig. 4B). The crystalline crust beneath M26A has a lower boundary (i.e., Moho) at 46.8 km depth and a Vp/Vs ratio of 1.83 ± 0.06 (Fig. 4C), a relatively high value compared with typical continental crust (~1.75). We adopt a constant Vp/Vs within each layer. The temperature and pressure corrections to Vp/Vs are relatively small compared to the measurement uncertainties (Hacker et al., 2015) unless the phase change of quartz from alpha to beta polymorphs occurs. Later in Section 4.2.4, we will show that this phase change should not occur systematically across our study region except in a few areas.

2.4. Quantitatively Constraining SiO2 wt% from Vs and Vp/Vs

Using published Vs and thermal models as well as the receiver function H-κ stacking, we have corrected Vs and measured Vp/Vs across the crystalline crust. This part describes how SiO2 wt% is drawn from those two profiles and how uncertainties are estimated. As shown in Figure 1, Vs and Vp/Vs of crustal rocks generally follow a smooth trend. Therefore, a given pair of Vs and Vp/Vs values will determine a SiO2 wt% value based on this trend. Applying this to each depth of the Vs and Vp/Vs profiles, a crustal SiO2 wt% profile is then constructed. In this process, we assume that Vp/Vs does not change with depth and in Section 4.2.5 we assess the possible error contribution to the result from this assumption. After a SiO2 wt% profile is generated, the bulk SiO2 wt% for the crystalline crust can be calculated from the depth integration of this profile.

Since both the Vs and Vp/Vs have random errors (Shen and Ritzwoller 2016), we use a bootstrapping method to quantify the random errors in the resulting SiO2 wt% versus depth profiles. In this approach, an ensemble of Vs profiles and Vp/Vs ratios is randomly selected out of Gaussian distributions whose widths are determined by their respective measurements (2% for Vs, Shen and Ritzwoller 2016; Fig. 4C for example for Vp/Vs). This produces an ensemble of SiO2 wt% versus depth profiles whose mean at each depth interval defines the result and a standard deviation that represents the uncertainties of the profile. The distribution of depth-integrated values from the ensemble of SiO2 wt% profiles (which can be shown as histograms as in Figs. 5D and 5H) defines the probability of the bulk silica content for the crystalline crust. Mean and standard deviation are drawn to represent our final assessment of the average composition with uncertainty.

To show that this simple method is able to distinguish crustal compositional variations, we perform synthetic tests on rocks with known compositions and seismic properties. We show that the results match the input within uncertainties. Since seismic properties of real crustal rocks can either be measured in the lab or be calculated using thermodynamic and petrophysical calculations (e.g., Perple_X, Connolly and Petrini, 2002; Excel worksheets from Abers and Hacker, 2016; A&H hereafter), we present the two cases using each method, respectively. We adopted two 40-km-thick crustal profiles in each test, one felsic and one mafic, with a mid-crustal discontinuity separating it into upper and lower crust. The composition we adopted for each test case is shown in Table 1, and the details of the rock sample and Perple_X parameters we used are summarized in the Supplemental Material (see footnote 1). The seismic speed is calculated at a 15 °C/km geotherm.

For the first case, depth-dependent Vs and Vp/Vs are calculated from the lab-measured P-T derivatives and reference velocities (Kern et al., 1996, 1999). The calculated Vs and Vp/Vs profiles are shown in Figures 5A and 5B. Both the felsic and mafic profiles have distinct Vs and Vp/Vs predictions: the felsic profile, on average, has a lower Vs and Vp/Vs across all depths compared with the mafic profile. Using these Vs predictions and the average Vp/Vs value of the crust and assuming proper uncertainties (2% for Vs and 0.07 for Vp/Vs), we construct two SiO2 wt% profiles following the method presented earlier. Results for this test are shown in Figures 5C and 5D and summarized in Table 1. As seen in Figure 5C, the upper-lower crustal composition differences in both profiles are reproduced. A shift in the average SiO2 wt% for the mafic profile is observed: the resulting composition is 3 wt% more intermediate, mostly due to the fact that the Vs-Vp/Vs-SiO2 wt% relationship we adopt has a misfit to the particular rock composition we choose in this test. The bulk composition distributions are shown in Figure 5D, and we find that the resulting estimates of the average crustal composition are still separated: when the input felsic profile has ~15 wt% more SiO2 than the mafic profile, we capture 75% of the original difference (~11 wt%) in average composition. The misfit to the input models (~1–4 wt%) is substantially smaller than the estimated uncertainties (~5 wt%).

For the second case, we take a different route to calculate the seismic speed in the lower crust. Rocks with the same chemical compositions may have different modal assemblages at different pressure and temperature conditions and may thus predict different seismic speeds. The depth-dependent phase equilibria may introduce bias to our results when using a uniform Vs-Vp/Vs-SiO2 relationship to all depths and temperatures. For example, the stable mineral assemblage of deep crustal felsic rocks might include garnet under cold geothermal conditions, which might predict a higher Vs and possibly lead to estimating a mafic composition. To test this scenario, we calculate the stable mineral assemblage at 30 km depth and 600 °C for the reported composition of the lower crustal samples used in test case 1 (Kern et al., 1996, 1999). We use the thermodynamic modeling software Perple_X (version 6.9.1, Connolly and Petrini, 2002; thermodynamic database hpha02ver.dat, Holland and Powell, 2003; Connolly and Kerrick, 2002). For the composition of the sample HT1, this calculation shows ~12.5% stable garnet at 30 km depth (calculations made between 20 and 40 km depths and 400–1000 °C produce ~4.5%–14.0% garnet). Based on the resulting phase partition, seismic properties are then calculated by extracting the modes and mineral chemistry from Perple_X output for use in A&H Excel worksheets (Abers and Hacker, 2016; details of the calculation can be found in the Supplemental Material). For this garnet-bearing felsic profile, a higher Vs (3.8 km/s) is found in the lower crust compared with the non-garnet-bearing original sample (3.65 km/s), as fast as the amphibole-bearing mafic profile (Fig. 5E). However, the calculated Vp/Vs remains lower (~1.65) than the mafic profile (~1.76). As a result, the changes in calculated seismic properties do not significantly alter the resulting quantification of lower crust SiO2 wt% (Fig. 5G). Notably, changing the forward calculation to Perple_X in the lower crust also alters the bulk crustal Vp/Vs which is the input for the compositional quantification. It thus leads to a slightly different composition profile for the upper crust as well. Overall, the resulting bulk SiO2 wt% distributions are still separated, and the misfits are still less than one standard deviation (Fig. 5H). We also performed synthetic tests on other samples reported by Kern et al. (1996, 1999) and found similar results. In conclusion, the method using Vs and Vp/Vs appears a viable way to quantify the SiO2 wt% of the continental crust. In the next chapter, we report the application of this method to over 1400 stations of the USArray and present the first continental-scale model of crustal composition based on in situ seismic measurements.

This chapter presents the primary results for crystalline crustal Vp/Vs across the contiguous U.S. and the 3-D compositional model built based on that map. The discussion on the uncertainties of the resulting model is presented in the following section.

3.1. Crustal Vp/Vs of the Continental U.S.

We apply the 2-layer H-κ receiver function stacking routine to 1708 USArray TA stations and obtain 1406 meaningful Vp/Vs measurements for the crystalline crust after quality control (Supplemental Material). The main results, depth to Moho and Vp/Vs ratio of the crystalline crust, are presented together with traditional 1-layer H-κ stacking results in Figure 6. After correcting for sediment-generated phases, our new map of Moho depth has fewer small-scale variations in regions with unconsolidated sediment (e.g., the central U.S.). More importantly, the new Vp/Vs map of the crystalline crust contains few extreme Vp/Vs values (e.g., Vp/Vs > 1.9, see Supplemental Material). The overall pattern in Vp/Vs is more closely correlated with surface geology: In the western U.S., higher Vp/Vs values are found in the regions with Cenozoic mafic magmatism such as the Snake River Plain and High Lava Plains, while lower values occupy extensional and granitic areas (i.e., the Basin and Range and Idaho Batholith). In the central/eastern U.S., higher Vp/Vs is found (>1.8), except for the Atlantic coastal plain. Factors such as crustal melt and extensive cracks seem unlikely to be pervasive at a continental scale, especially in tectonically stable areas. We thus hypothesize that these high Vp/Vs ratio measurements mainly reflect compositional variations.

3.2. A 3-D Crustal Compositional Model

Using the new Vp/Vs map across the U.S., we apply the compositional modeling to all 1406 stations and construct 1-D compositional models following the method described in Section 2.3. Based on individual 1-D compositional models with uncertainties, a 3-D SiO2 wt% model with uncertainties is then constructed. Here we only briefly describe the resulting model and leave the discussion of uncertainties to Section 4.1. The benchmark of the resulting model and its tectonic implications are presented in Section 5.

3.2.1. Individual SiO2 wt% Profiles

Figures 7A7C presents input Vs and Vp/Vs as well as resulting SiO2 wt% profiles for stations F12A and M26A. The two stations were chosen since they have distinct Vs profiles and Vp/Vs measurements (Figs. 7A and 7B). As expected, the two locations exhibit very different compositional profiles, and thus different distributions of bulk SiO2 wt% (Figs. 7C and 7D). Station F12A is in the Atlanta lobe of the Cretaceous Idaho Batholith and is largely felsic (average SiO2 wt% is ~67.0 ± 3.2), midway through the range of observed chemistry in the surface exposure of the batholith (55–75 wt% SiO2; Hyndman, 1984; Gaschnig et al., 2011). Site M26A is on the Great Plains and has a more mafic crust (~50.3 ± 2.8 wt% SiO2) (Fig. 7D). This site is covered by Phanerozoic sediment but is thought to be underlain by Proterozoic basement (e.g., Worthington et al., 2016), and our estimated composition is consistent with crustal xenoliths from the Stateline kimberlite district ~160 km to the southwest (average 48 wt% SiO2; Farmer et al., 2005). For each of the 1406 stations, SiO2 wt% profiles and bulk SiO2 wt% distributions are constructed using the same approach.

3.2.2. Lateral Variations

With 1406 bulk average SiO2 wt% distributions, we calculate their mean and standard deviation and drew the smoothed map views in Figure 8, representing the final crystalline crustal SiO2 wt% map and uncertainty for the continental U.S. In terms of horizontal variations in the model, this plot presents a first-order pattern with a division between a generally mafic central/eastern (bulk average ~54.6 wt% SiO2) and a more felsic and compositionally diverse western bulk average (~58.3 wt% SiO2) and the easternmost U.S. This observation of dichotomy in the crustal composition is qualitatively consistent with observations that Archean and Proterozoic crust is more mafic than Phanerozoic crust from investigations based on Vp/Vs ratio or Vp compilations (Rudnick and Fountain, 1995; Zandt and Ammon, 1995). At regional scales, we observe the felsic bulk crystalline crust in the: Idaho Batholith, Basin and Range, southern Basin and Range, Rio Grande Rift, the central part of the Rocky Mountains, and parts of the Colorado Plateau. These felsic regions are mostly distributed in the western U.S. The High Lava Plains, Snake River Plains, and much of the Colorado Plateau, in turn, is more mafic. In the eastern U.S., the northern Appalachians and Coastal Plains are generally more felsic than the southern Appalachians and the cratonic core. These province-related variations are also seen in vertical transects discussed in Section 3.2.4 as well. Additionally, variations in the resulting model also reveal both spatial and vertical patterns that are in line with the petrological records such as the distribution of volcanism and locally collected deep crustal xenolith samples, which will be discussed in Section 5.

3.2.3. Horizontal Sections at Upper, Middle, and Lower Crust

Figure 9 shows the vertical patterns in the resulting compositional models at upper (5 km below the bottom of the sedimentary layer as determined from 2-layer stacking), middle (mid-depth of the crystalline crust), and lower (5 km above the Moho) crustal depths. In general, the SiO2 wt% increases with depth, following the increasing trend of Vs. We note similarities between the central/eastern U.S. upper crust (~60.0 wt% SiO2, Fig. 9A) and western U.S. middle crust (~58.8 wt% SiO2, Fig. 9B), as well as the central/eastern U.S. middle crust (~54.1 wt% SiO2, Fig. 9B) and western U.S. lower crust (~54.1 wt% SiO2, Fig. 9C). The similarities suggest the possibility that the upper felsic crust of the old central/eastern cratons and lower mafic crust in the tectonically active western U.S. have been diminished by long-term erosion and delamination, respectively. In addition, we note that SiO2 wt% in the lower crust of the eastern/central U.S. is the least uncertain (~3 wt%, Fig. 9F). This is because the variations in inferred SiO2 wt% due to the uncertainties in seismic properties are relatively small when both the Vs and Vp/Vs are high (Vs ~4.1–4.2 km/sec; Vp/Vs ~1.8). An example can be found in Figure 7C, showing a relatively smaller uncertainty in lower crustal composition beneath station M26A.

3.2.4. Vertical Transects of SiO2 wt% along Specific Profiles

Figure 10 presents three vertical transects across major tectonic provinces in the study area: AA’ at 39° N latitude crosses the whole continent; BB’ and CC’ located mostly in the western U.S. cut along or across the Snake River Plain-Yellowstone hotspot track, respectively. These transects are shown with surface topography (with elevation exaggerated vertically) to highlight the tectonic provinces they intersect with. Profile AA’ crosses the entire continental U.S. from the Sierra Nevada to the eastern coastal plains. A general trend is that more felsic crust (e.g., > 60 wt% SiO2) is thicker on the western side, including the Basin and Range and the Rocky Mountains. In comparison, a thicker, more mafic crust (<60 wt% SiO2) is seen beneath the Archean and Proterozoic central/eastern U.S. Notably, the Colorado Plateau, sitting between the continental rift of the Basin and Range and the Colorado Rocky Mountains, has a relatively thick mafic deep crust. At the eastern end of AA’, the felsic upper crust is thicker beneath the coastal plains than in the cratonic core to its west, suggesting a different crustal origin or tectonic history. Transects BB’ and CC’ run along and across the Snake River Plains and Yellowstone hotspot track, respectively. BB’ shows a compositional slope changing from the Basin and Range (mostly felsic) to the Snake River Plains and Yellowstone (felsic upper crust with a relatively thick mafic deep crustal root) to the Archean and Proterozoic aged crust of the Great Plains and is mafic for most of the crystalline column. In CC’, the crustal composition exhibits variations from province to province: Idaho Batholith (mostly felsic to 25 km), Snake River Plain (mostly mafic below 15 km); Basin and Range (thin and felsic); Colorado Plateau (thick and more mafic crust); southern Rio Grande Rift (more felsic). The variations reflect the complex modification to the chemistry of the crust due to different tectonics in the western U.S. and deserve further detailed investigations for each province.

Model errors include systematic and nonsystematic errors. The systematic errors come from the fact that the assumptions we made during the compositional model building may not hold completely true for real Earth. The nonsystematic errors come from the uncertainties in the seismic data (Vs, Vp/Vs). This section discusses the nonsystematic errors first in Section 4.1 and discusses factors that contribute to the systematic errors in Section 4.2. The result discussed here is based on the surface derived from the lab-measured seismic properties only.

4.1. Nonsystematic Errors

The nonsystematic errors in the SiO2 wt% come from the random errors in the Vs model and crystalline crust Vp/Vs model. The uncertainties in the Vs model have been discussed extensively by Shen and Ritzwoller (2016). To simplify the discussion, here we conclude that the uncertainty in the crustal Vs is generally 2%. Traditionally, the uncertainty for crustal Vp/Vs is calculated during the H-κ stacking, using either the second-order derivative at the maximal energy (Zhu and Kanamori, 2000) or the bootstrapping uncertainty bounds method (Efron and Tibshirani, 1991). In this study, we not only stacked waveform energy of all event summation, but also calculated Vp/Vs from each individual event using H-κ stacking. The standard deviation of the resulting Vp/Vs ratios is then taken as the uncertainty estimate for Vp/Vs (Fig. 6D). On average, the uncertainty on Vp/Vs estimated through this method is ~0.07, an estimate that is generally greater than the traditional methods. This uncertainty definition reflects a more conservative choice in assessing the errors in seismic properties.

Based on the uncertainties of Vs and Vp/Vs, the method discussed in Section 2.3 allows an assessment of the nonsystematic errors in the resulting SiO2 wt% profiles (as shown in Fig. 7C for individual stations) or bulk SiO2 wt% estimates. As shown in Figure 8B, the nonsystematic errors of the bulk average SiO2 wt% is estimated to be ~4.4 wt% across the continental U.S., much smaller than the regional variations in bulk SiO2 wt% (from ~45 wt% in Montana to ~68 wt% in Idaho).

4.2. Systematic Errors

Precisely quantifying total systematic errors is more challenging since multiple factors contribute to them, and they may constructively or destructively interfere with each other. These factors include the uncertainty of the Vs-Vp/Vs-SiO2 wt% relationship, choice of Vs models, assumptions in thermal, pressure, and attenuation corrections to the Vs models, lack of knowledge on possible quartz phase transitions and partial melt in the deep crust, and the assumption of a constant Vp/Vs across the crystalline crust. Here we discuss each source individually and present some tests to show how changing some of these assumptions impacts the resulting SiO2 model, either at individual stations or at a continental scale. Particular attention is paid to how our conclusions might be affected by the systematic errors, and we find these effects are relatively small compared with the uncertainties that we report.

4.2.1. Uncertainty of the Vs-Vp/Vs-SiO2 wt% Relationship

The uneven distribution of the rock samples in Vs-Vp/Vs space, as well as some SiO2 wt% spikes, bring uncertainties to the derived Vs-Vp/Vs-SiO2 wt% relationship. Applying a bootstrapping method in constructing the relationship will help us to quantify the uncertainty level. We randomly choose 90% of the rock samples to create a new interpolation surface and repeat the process 100 times. The standard deviation of the ensemble of the interpolation surfaces is calculated and regarded as the uncertainty of the relationship (Fig. 11A). We find that the overall uncertainty level is low (<1 wt%), but large uncertainty exists when one sample has different SiO2 wt% from other samples with similar Vs and Vp/Vs (e.g., Vs: 3.4–3.5 km/s and Vp/Vs: 1.70–1.75). Then, we use the ensemble of interpolation surfaces instead of one single surface to calculate the SiO2 wt% models for stations M26A and F12A. Figure 11B shows that for M26A, the SiO2 wt% profiles with/without considering the relationship uncertainty are similar. For F12A, the uncertainty in the interpolation surface increases the overall uncertainty by 1–2 wt% in the SiO2 wt% profile. Figure 11C shows the distributions of bulk average SiO2 wt% of the two example stations are still well separated. Compared to Figure 7, the differences in bulk SiO2 wt% are <1 wt% and the uncertainties defined by the standard deviations increase <1 wt%. We conclude that the uncertainty of the Vs-Vp/Vs-SiO2 wt% relationship is minor and the large-scale features in crustal composition revealed in Figure 8 should hold.

4.2.2. Choices of Different Vs Models

The Vs models that we use may contain bias from the presumptions made during their construction. Using different published crustal Vs models in the U.S. can help us further test the reliability of the crustal composition. The same approach has been applied to the crustal model made by Schmandt and Lin (Schmandt et al., 2015; S&L 2015 hereafter) to quantify the SiO2 wt%. The difference between the two Vs models and the resulting SiO2 wt% for station M26A is shown in Figure 12. There is a ~1 wt% DC shift, mostly due to the lower Vs in the upper crust, compared to the SiO2 wt% from S&R 2016. This difference is much smaller than the 1 standard deviation uncertainty. In the scale of the whole continental U.S., the pattern of more felsic western (~57.4 wt%) and more mafic central/eastern portions (~54.6 wt%) still holds, and the difference is very subtle.

4.2.3. Assumptions in Thermal/Pressure/Attenuation Corrections

The assumptions in thermal/pressure/attenuation (T/P/Q hereafter) models may bring bias to the corrections to the Vs profiles since the assumed models are not error-free. To test the sensitivity of the result to the thermal correction, we perform a test with a different thermal model based on the latest geothermal heat flux measurements of the U.S. (Blackwell et al., 2011) with a pre-set crustal heat generation model (Rudnick and Fountain, 1995). A similar approach is performed based on this new thermal model and the resulting compositional patterns persist (i.e., more felsic western (~58.7 wt% SiO2) and more mafic central/eastern U.S. (~55.0 wt% SiO2)). The influences from the pressure model and Q model are smaller than the temperature effects. Assuming Vs is 4 km/s at 40 km, the difference between the two thermal models we tested may exceed 300 °C, which results in 0.06 km/s (1.5%) difference in corrected Vs, and up to ~0.5–2 wt% of SiO2 wt% (varying from different Vp/Vs) at this depth. Using a different pressure model (e.g., 30 MPa/km rather than 27 MPa/km), the pressure difference of 120 MPa generates 0.012 km/s (0.3%) difference in corrected Vs, leading to ~0.05–0.4 wt% change in SiO2. Additionally, using a different crustal Q value (e.g., 200 rather than 600), the attenuation corrections to Vs change from ~0.3% to ~0.9%, leading to a compositional difference of ~0.1–0.7 wt%.

Considering that a lower geotherm is usually associated with higher densities and lower attenuation (and vice versa), the biases from assumptions in T/P/Q may constructively interfere with each other. We use stations M26A and F12A to exemplify how different combinations of thermal, density, and attenuation models shift the final result. A test using a substantially colder geotherm (12 °C/km), higher density (31 MPa/km), and lower attenuation (Q = 1000) at station M26A would result in a < 0.8% difference in the Vs profile and produce a compositional model with a bulk crustal SiO2 of 50.8 ± 3.3 wt%, 0.5 wt% higher than the published model here. For station F12A, if a substantially hotter geotherm (30 °C/km, 25% higher compared with 24 °C/km in the thermal model we used), lower density (25 MPa/km), and higher attenuation (Q = 200) model is assumed, Vs would increase by up to ~2%. The resulting bulk crustal SiO2 wt% would be 66.4 ± 4.1 wt%, 0.6 wt% lower than the presented model. In conclusion, when a reasonable combination of thermal, density, and attenuation values is assumed, the bias is limited to up to 1.0 wt%.

4.2.4. Possible α-β Quartz Transition and Partial Melt

In this work, we consider chemical composition as the dominant factor that contributes to crustal Vp/Vs. In reality, there are other factors that may also affect it. These include the α-β quartz phase change (Mechie et al., 2004; Kuo-Chen et al., 2012) and locally existing partial melt, with both increasing Vp/Vs substantially, and placing challenges in constraining crustal compositions (e.g., Shillington et al., 2013; Jagoutz and Behn, 2013). Partial melt also greatly decreases Vs. Ignoring these effects on Vp/Vs might cause an underestimate in SiO2 wt% for areas with such complications. With a crustal thermal model (Boyd, 2020), we identify the stations that might be affected by partial melt if the crustal temperature is above 650 °C when the pressure reaches 500 MPa. For a possible α-β quartz phase transition, we identify the stations whose thermal structure exceeds the pressure-temperature condition using Equation 2 (Shen et al., 1993):
formula

Out of 1406 stations analyzed, 157 stations may experience either partial melt (29, 2% of the total stations) or the quartz phase transition (154, ~11% of the stations) within the crust, and their locations are shown in Figure 13, given the temperature conditions that are assumed here. We note that amongst the 154 stations with a possible quartz phase change, only 34 (~2% of total) stations might have beta quartz stable in a significant portion of the crust (e.g., thicker than 10 km). We notice that most of the identified stations are located in the western U.S., especially for some regions along the Snake River Plain and northwestern Basin and Range. For these regions, interpreting their higher Vp/Vs without considering the partial melt/quartz phase change might underestimate the resulting SiO2 wt%. For most of the eastern U.S., where the result shows a more mafic deep crust, few stations are affected. In summary, these observations indicate that the majority of the stations we analyzed are not affected by partial melt or quartz phase change, and the large-scale trend in crustal composition revealed in Figure 8 should hold. The quantitative estimate of the effects on the interpretation of crustal composition for the Snake River Plain, Yellowstone, and Basin and Range is beyond the scope of this paper but warrants further investigation in the future.

4.2.5. Depth-Varying Vp/Vs

Except for the cases discussed in Section 4.2.4, the assumption that Vp/Vs is constant across the crystalline crust may not be true in the real crust. To evaluate how this assumption influences our quantification of SiO2 wt%, here we perform a test using station M26A that removes this assumption. Instead of using a constant Vp/Vs ratio across the crystalline crust, we allow it to vary with depth, with the bulk Vp/Vs following the measured value and associated uncertainty. These synthetic randomly sampled, depth-varying Vp/Vs profiles are parameterized by seven cubic B-splines. As shown in Figure 14A, the removal of the constant Vp/Vs assumption leads to larger uncertainty at each depth, ending up with a wider corridor for the SiO2 wt% profiles (Fig. 14B). As a result, the resulting bulk SiO2 wt% shows a higher uncertainty compared with the original result (3.1 wt% compared with 2.8 wt%). Notably, the average resulting compositional profile is shifted toward a slightly more felsic composition (the bulk is ~1.4 wt% more felsic, Fig. 14C). This is because when Vp/Vs for M26A is allowed to vary with depth, its impact on composition is asymmetrical: when the Vp/Vs randomly decreases (e.g., from 1.83 to 1.73) at some depths, a more felsic composition (~7 wt% increase in SiO2) will be obtained, while when it increases (e.g., from 1.83 to 1.93), the decrease of SiO2 wt% is not as significant (only ~1.5 wt% decrease in SiO2), since the average Vp/Vs is already high. In summary, assuming that Vp/Vs stays constant with depth will cause an underestimate to the uncertainty of the resulting compositional model and introduces a bias of up to ~1.5 wt% for the station M26A. Because of the asymmetry effect, for other sites with a felsic crust, this bias would overestimate SiO2 wt%; for mafic crust, the bias underestimates SiO2 wt%. For the continental U.S., there are more areas with mafic-intermediate crust than with felsic-intermediate crust, and we thus conclude that not knowing the Vp/Vs variation in depth may introduce a bias toward underestimating the overall SiO2 wt % by ~1 wt%.

In this section, we discuss the resulting model in more detail. In particular, we interpret the model at both local and regional scales, with benchmarks to locally collected xenolith samples, surface geology, and globally aggregated deep crustal composition. We pay extra attention to the bulk composition of the deep crust of the U.S. with respect to the estimated global compositional model and its tectonic implications. Additionally, we discuss the differences between the western and central/eastern bulk crustal compositions. Finally, we discuss some caveats of the current method.

5.1. Comparison with Local Xenoliths

On a local scale, thermobarometry and bulk compositional estimates from deep crustal xenoliths exhumed by Eocene minettes in central and northern Montana (Barnhart et al., 2012; Mahan et al., 2012) were recently compiled (Schulte-Pelkum et al., 2017), and they provide another independent comparison. These data represent direct measurements of the deep crustal composition and allow us to perform an in situ benchmark with the seismologically modeled composition. As shown in Figure 15A, almost all deep crustal granulite xenoliths collected at the three sites (blue error bars) exhibit intermediate-mafic compositions at equilibrium depths between 20–54 km except for one. When compared with the composition model for a geographically close station, C19A, most of the measured SiO2 wt% are within the uncertainty estimate, as well as possible ranges at 14 other nearby stations (Fig. 15B). The striking consistency between the petrologically and seismologically obtained compositions confirms that the deep crust of this area is most likely mafic, and further demonstrates the usefulness of determining deep crustal SiO2 wt% using Vs and Vp/Vs at a regional scale. At other localities, our model is consistent with xenoliths at the Santa Lucia, California, USA (63.2 ± 5.1 versus 64.7 wt% from xenoliths, Ducea et al., 2003), Colorado-Wyoming (USA) state line (52.9 ± 4.8 versus 45–52 wt% from xenoliths, Farmer et al., 2005), and Leucite Hills sites, Montana, USA (49.8 ± 3.6 versus 42–52 wt% from xenoliths, Farmer et al., 2005).

5.2. Trend in the Bulk Crustal Composition at Regional Scales

On a regional scale, our result is also consistent with the geological features and division of tectonic provinces and additionally provide new insights into the tectonism of the continental U.S. Shown in Figure 15B, regions that have a mafic bulk crust, estimated by the resulting compositional model, geographically overlap with large amounts of Cenozoic and Mesozoic mafic volcanics (NAVDAT, Walker et al., 2004) and with areas where previous geophysical and petrological studies suggest that the bulk of the crustal column is mafic (e.g., McCurry and Rodgers, 2009). In contrast, extensional and granitic provinces (i.e., the Basin and Range and Idaho batholith) are generally more felsic, indicative of a crustal extensional model involving delamination of mafic lower crust rather than brittle faulting (Sammon et al., 2020). These features suggest that our model is consistent with the surface geology of the tectonically active western U.S. On the east coast, the higher SiO2 areas appear to closely correspond to terranes that originated as microcontinents with high arc-related plutonic activity prior to accretion during the Appalachian orogen (Carolina superterrane in the southeast and Avalon terrane in New England; e.g., Hatcher, 2010). The distinct high SiO2 estimate in southernmost Florida probably reflects exotic West African granite similar to that identified in well-bore cuttings by Dallmeyer et al. (1987). As much of these terranes, particularly in the south, lie beneath the coastal plain, the spatial variations in composition from our model offer additional detail that could be helpful in identifying terrane boundaries. Compared with previous studies on continental crustal composition (e.g., Rudnick and Fountain, 1995), we note that this model is constructed by sampling the continental U.S. with the regularly deployed USArray, and reveals signatures of how tectonism has introduced modifications to the crustal composition and perhaps crustal strength (Lowry and Smith, 1995) in an unprecedented way. Our resulting model contains geological implications that should be further investigated in the future.

5.3. Seismic Signatures of Deep Crustal Composition of the Archean and Proterozoic U.S.

We note a considerable difference between the western and central/eastern U.S. bulk crustal compositions, no matter which Vs model and thermal model is used. We display the reasons for the difference visually. In Figure 16, we show that the bulk average Vs and Vp/Vs for stations in western and central/eastern U.S. have different distributions in Vs-Vp/Vs space. The complex compositional patterns in the western U.S. result from relatively low Vs (<3.8 km/s) and scattered Vp/Vs. The felsic crust in the Idaho batholith and Basin and Range is similar to granites and felsic gneisses, while the mafic crust in the High Lava Plain and Snake River Plain is close to basaltic rocks. The central/eastern U.S. crust is mostly mafic because both the bulk average Vs (>3.7 km/s) and Vp/Vs (>1.7) are high. The Archean and Proterozoic crusts may consist of different types of rocks (from mafic amphibolites and granulites to gabbro).

As a result of the observation of high Vp/Vs and Vs, a thick, mafic lower crust beneath much of the stable, cratonic central/eastern U.S. is found in the resulting model (Figs. 9E and 10). This result is partly based on the relatively fast Vs in the lower crust (>3.8 km/sec, Shen and Ritzwoller, 2016). However, in areas with a thick crust, garnet may form long after orogenesis (Fischer, 2002; Blackburn et al., 2018), and garnet-bearing felsic rocks are known to have relatively high Vs (e.g., > 3.8 km/sec, Hacker et al., 2015; Williams et al., 2014). As a result, garnet-bearing felsic rocks should not be ruled out in interpreting faster seismic speeds in the Vs model (e.g., Schulte-Pelkum et al., 2017). Here we examine the likelihood that the seismic signature can be explained by a more felsic, garnet-bearing composition.

Indeed, thermodynamic calculations on other natural garnet-bearing felsic rocks at lower crustal conditions predict a faster Vs (e.g., Canadian Shield: 3.7–3.9 km/sec, Williams et al., 2014; Robinson Range, Bearspaw Mountains, Montana, USA: up to 3.8 km/sec, Barnhart et al., 2012; Sweetgrass Hills, Montana, USA: ~3.88–3.94 km/sec, Mahan et al., 2012; Schulte-Pelkum et al., 2017). However, the Vp/Vs for the same samples are generally low (Canadian Shield: 1.74; Robinson Range, Bearspaw Mountains: 1.73–1.74; Sweetgrass Hills, central Montana: ~1.70 on average), all distinguishable from mafic diabase (Vp/Vs ~1.82). Instead, garnet-bearing mafic rocks have a higher Vp/Vs (>1.8). In the synthetic test shown in Section 2, we examined a felsic lower crustal composition that may allow >10 wt% garnet stability, the thermodynamic calculation also reveals a relatively low Vp/Vs (<1.7), allowing us to reproduce a felsic composition as the input.

For the Archean and Proterozoic central U.S., a higher Vp/Vs for the average crystalline crust is found (~1.82 for much of this region, Fig. 6D). Given that Vp/Vs is likely increasing with depth, the lower crust of the central U.S. is more likely to be higher than the crustal average (i.e., > 1.82). This Vp/Vs of deep crust is not compatible with felsic, garnet-bearing rocks. Finally, the deep mafic crust is also well benchmarked by the thermobarometry result for deep crustal xenoliths in this area (as shown in Section 5.2). Given this observation, we conclude that if garnet is pervasively stable in the deep crust of the central U.S., the seismic observation (especially the Vp/Vs) suggests that the rocks are more likely to be mafic than felsic. This suggests that our result of deep mafic crust beneath the cratonic U.S. is perhaps unlikely to be strongly biased by garnet stability. However, as shown in the next section, using an alternative sample database derived from thermodynamic calculations could bias the results to more felsic compositions by a few percent, although the deep crustal composition of the central/eastern U.S. is still on the mafic side of the xenolith-terrane discrepancy (Supplemental Material).

5.4. Additional Caveats for the Method and Compositional Model

In Section 4, we have identified the major factors that may introduce bias into the resulting SiO2 wt% model. In addition to them, there are additional caveats in using the petrophysics database for future refinement in the model.

One of the most significant caveats is that the resulting compositional model depends on how the Vs-Vp/Vs-SiO2 wt% relationship is derived (i.e., lab measurements versus thermodynamic calculations), and the result can potentially be biased by the choice of database. Using the interpolation surface from thermodynamics calculations leads to different Vs-Vp/Vs-SiO2 wt% interpolation surfaces, especially for typical crustal Vs and Vp/Vs ranges (Figs. 2A and 2B). In Figure 17, we compare two interpolation surfaces corresponding to 37 km in depth derived from lab-measured data and thermodynamic calculations. The P-T condition at 37 km is close to the P-T condition, which the phase equilibrium has been calculated at thermodynamically (1 GPa, 650 °C [amphibolites] and 750 °C [granulites], Hacker et al., 2015). The contrasts in Figures 17A and 17B indicate that the usage of thermodynamics-calculated rock seismic properties produces a ~8.0 wt% difference in the absolute SiO2 wt%. One fact that needs attention is that the thermodynamic-based surface is not applicable to some stations with high Vp/Vs (>1.85). Only considering the stations that appear in both maps, the average contrast in SiO2 wt% at 37 km SiO2 wt% is ~6.7 wt% (Figs. 17C and 17D). Additionally, we note that the spatial patterns are similar, supporting that (1) the central/eastern U.S. is more mafic than the western U.S.; (2) crust east to the Appalachians is more felsic than the Archean and Proterozoic crust; (3) the Snake River Plain and High Lava Plains are more mafic than the Idaho batholith. Therefore, we conclude that the choice of Vs-Vp/Vs-SiO2 wt% relationship from a different data set of rock with seismic properties calculated based on local phase-equilibrium will induce a systematic shift in the estimated deep crust SiO2 wt%, and this shift is greater than the nonsystematic errors we discussed in Section 4. However, the major compositional variations across tectonic boundaries in the continental U.S. remain unchanged and the associated discussions in the main manuscript remain valid. However, having addressed this point, we suggest that additional investigation of the surfaces from the rock properties calculated at different P-T conditions should be performed for future refinement.

Moreover, focusing on using the interpolation surface from lab-measured rock properties, our methods include certain assumptions and simplifications. Accordingly, the second caveat is that we ignore microcracks in the upper crust when the Vs is corrected for temperature and pressure. These cracks, especially in unconsolidated sedimentary rocks, may reduce observed Vs, and ignoring their effect may introduce bias when we apply the pressure correction (Kern et al., 1996, 1999). In this work, we only construct a compositional model for the crystalline crust, but such bias may still exist in the uppermost crust.

The third caveat in this work is that the observed Vp/Vs is based on single station receiver function analysis, which is sensitive to the depth-averaged local structure near the station location within a ~20 km range. However, the Vs model we use has a lower horizontal resolution (~50–100 km) and higher vertical resolution (~10 km). This differential resolution may introduce further bias for certain locations (e.g., a spottier compositional model than the input Vs model). Additional work to reduce this differential resolution would be desirable for improving this technique. Finally, we note that the seismic model we used here is mainly from Rayleigh waves, which are sensitive to the velocity of vertically polarized shear waves. As a result, the corrected Vs used in this study may differ from the Voigt averaged speed for certain areas by up to 2% (e.g., the western U.S. where stronger crustal radial anisotropy caused by the lattice preferred orientation of anisotropic minerals such as micas and amphiboles have been identified, Moschetti et al., 2010). This will not fundamentally change the patterns we see in the composition. Future improvement thus would include using a Voigt averaged Vs model considering both azimuthal and radial anisotropy.

In this paper, we present a method that combines the Vs and Vp/Vs model to quantify the SiO2 wt% of continental crust, based on the observation that composition varies monotonically with these two seismic properties. We summarize our findings below:

  1. Lab measurements of amphibolites, mafic and felsic granulites, and other crustal rocks show that SiO2 wt% of the crustal rocks vary monotonically with Vs and Vp/Vs. Mafic rocks generally have a higher Vs and Vp/Vs, while the felsic rocks exhibit lower Vs and Vp/Vs values.

  2. Based on conclusion 1 and synthetic tests using both the lab-measured rock properties and thermodynamic calculations, we show that using a Vs profile and average crustal Vp/Vs can quantitatively distinguish between mafic and felsic crust, even when garnet is stable in the lower crust.

  3. A sequential 2-layer H-κ stacking method applied to USArray with rigorous quality control produces a smooth, high-resolution crystalline Vp/Vs map at a continental scale. Variations in Vp/Vs follow the tectonic provinces, indicating a geological origin of the seismic signature.

  4. Combining the new Vp/Vs map with published Vs models, a 3-D compositional model for the continental U.S. emerges. The model shows that crustal composition varies with tectonism and surface geology:

    • - for the tectonically active western U.S., the regions with Cenozoic active rifting, granitic batholiths are more felsic, while the regions with mafic volcanism such as the Yellowstone hotspot track and High Lava Plains show a more mafic crust;

    • - for the tectonically stable central and western U.S., a more mafic deep crust is observed;

    • - for accreted terranes along the east coast (Coastal Plains, northern Appalachians, and Piedmont regions), the crust is more felsic-intermediate, differing from the cratonic core.

  5. Benchmarks with local xenolith data collected in central Montana and their thermobarometry results show high consistency between seismologically constrained deep crustal composition and xenolith samples.

  6. A benchmark test of the crustal SiO2 wt% model with localities of Mesozoic and Cenozoic volcanism shows that regions with those volcanic rocks commonly have mafic compositions.

  7. The ignorance of complexities in mineral phase stabilities and the possibility biased choice of rock property database, differential horizontal resolution of Vs and Vp/Vs measurements, upper crustal microcracks and fluids, lack of vertical resolution in crustal Vp/Vs, and lack of consideration of the quartz phase change, partial melt, and anisotropy are the major caveats of this work. A rigorous error analysis, however, shows that some of them either only affect certain areas (e.g., quartz phase change and partial melt are not significant sources of error for central/eastern U.S.), or impose limited impacts on the resulting bulk crustal compositional model (e.g., lack of vertical Vp/Vs resolution only contributes an error of up to 1.5 wt%). The ignorance of the more complex phase stability for deep crustal conditions contributes to the potential bias up to a few percent (~6.7 wt% at 1 GPa), contributing the largest bias to the approach used in this study.

Since the development of ambient noise tomography, more accurate Vs models of the crust have been built across all of the major continents (North America, South America, E. Asia, Australia, Europe, Africa, Antarctica). These models, together with the widely applied H-κ receiver function investigation across the globe, allow a broader investigation of the crustal composition using the method described in this paper. In addition to the average crustal composition that can be further constrained, the potential compositional models will provide fruitful implications for other geochemical and geophysical subjects, and we only list a few here: (1) Since the silica content of crustal rocks is approximately correlated with heat generating elements (i.e., K, U, and Th), the new compositional model allows the derivation of heat generation of the crust. (2) As the rheology of crustal rocks is highly dependent on the chemical composition, the new model will also facilitate better quantification of rheological properties and strength of continental crust (Shinevar et al., 2015; Lowry and Pérez-Gussinyé, 2011). (3) Using seismology to constrain the composition will allow additional investigation into the geology of continents whose surface geology is not well assessed (e.g., continents covered by ice sheets, such as Antarctica and Greenland).

1Supplemental Material. Supplemental text, Figures S1–S4, and Tables S1–S3. Please visit https://doi.org/10.1130/GSAB.S.21091078 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Brad S. Singer
Associate Editor: Marcin Dabrowski

Aspects of the work were supported by National Science Foundation (NSF)-1744883, NSF-1945856, NSF-13002097, and NSF-1937343. The authors thank Brad Hacker for providing the petrological database, and Mark Behn for providing the thermodynamically calculated database. We appreciate Editor Marcin Dabrowski and the reviewers Mark Behn and Emily J. Chin for their valuable opinions. Seismic data (USArray-TA) was downloaded from Incorporated Research Institutions for Seismology (IRIS). IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope proposal of the NSF under cooperative agreement Division of Earth Sciences 1261681. The relationship surfaces and figures are generated/plotted by the Generic Mapping Tool (Wessel et al., 2019). The Cenozoic magmatism data was downloaded from the North American Volcanic and Intrusive Rock Database, part of the EarthChem initiative (Walker et al., 2004).

Gold Open Access: This paper is published under the terms of the CC-BY license.