Compositional Attributes of the Deep Continental Crust Inferred From Geochemical and Geophysical Data

This study provides a global assessment of the abundance of the major oxides in the deep continental crust. The combination of geochemistry and seismology better constrains the composition of the middle and lower continental crust better than either discipline can achieve alone. The inaccessible nature of the deep crust (typically >15 km) forces reliance on analog samples and modeling results to interpret its bulk composition, evolution, and physical properties. A common practice relates major oxide compositions of small‐ to medium‐scale samples (e.g., medium to high metamorphic grade terrains and xenoliths) to large scale measurements of seismic velocities (Vp, Vs, Vp/Vs) to determine the composition of the deep crust. We provide a framework for building crustal models with multidisciplinary constraints on composition. We present a global deep crustal model that documents compositional changes with depth and accounts for uncertainties in Moho depth, temperature, and physical and chemical properties. Our 3D compositional model of the deep crust uses the USGS Global Seismic Structure Catalog (Mooney, 2015) and a compilation of geochemical analyses on amphibolite and granulite facies lithologies (Sammon & McDonough, 2021, https://doi.org/10.1029/2021JB022791). We find a SiO2 gradient from 61.2 ± 7.3 to 53.3 ± 4.8 wt.% from the middle to the base of the crust, with the equivalent lithological gradient ranging from quartz monzonite to gabbronorite. In addition, we calculate trace element abundances as a function of depth from their correlations with major oxides. From here, other lithospheric properties, such as Moho heat flux ( 21.6−5.6+16.0 ${21.6}_{-5.6}^{+16.0}$ mW/m2), are derived.

The calculation adopts a subdivision of the global continental crust into 13 tectonic provinces (Figures 1 and 2) to speed calculations and extrapolate results to areas with lower data coverage. Individual tile resolution of this global model is set to 1°latitude × 1°longitude × 3 km depth as a default, but can be changed in the model to suit user needs. We chose this default resolution for our global model based on the resolution of our crustal categories (each 1° × 1° of crust was assigned a tectonic province based on models such as CRUST1.0, Litho1.0, and modifications discussed further in Section 2.1), and the resolution of our crustal thickness and temperature data, the ramifications of which are discussed further in the Results section. For considering higher resolution, regional scale data, the same methods can be used. Instead of simplifying the crust into tectonic provinces, calculations are run for individual seismic velocity profiles, so that if there are, for example, 34 seismic velocity profiles as inputs, there will be 34 locations for which compositional profiles are generated.
We calculated the overlapping probability between measured seismic velocities and the Perple_X-derived velocities for amphibolites and granulites equilibrated at middle and lower crustal pressures and temperatures. A reference average crustal density of 2,900 ± 200 kg/m 3 is used to calculate lithostatic pressure , cf., Christensen and Mooney (1995)). Integrating the area under both curves, the area shown as magenta in Figure 3, for a sample of composition X yields the total probability of sample X producing the observed seismic signal. Repeating this technique for a multitude of sample compositions at various depths and temperatures yields a final Monte Carlo model for a deep crustal composition. Probability distributions are generated for Vp, Vs, and Vp/Vs and then multiplied together to further constrain the final probability.

Model Inputs
A global model of Vp, Vs, and Vp/Vs was generated from a compilation of over 5000 (Vp) and 1000 (Vs) 1-D seismic velocity profiles obtained from the Global Seismic Structure Catalog (GSSC) database (Mooney, 2015). Both controlled and passive source seismic velocity profiles were included to increase data coverage. We included profiles with Vp and Vs data that had been sampled at a minimum of 5 depth intervals within the crust. Figure 2 shows our tectonic province designations and the location of each seismic velocity profile used. We used global Moho depths from Litho1.0, except on the continental margins, where we reference Szwillus et al. (2019) Moho values. In comparison to Litho1.0, Szwillus et al. (2019) incorporated a larger data set on the continental margins (∼1600 profiles) and did not average depths across the continent-ocean transition. Global Moho temperatures were generated from the TC15 global temperature model of Artemieva (2006). We assumed a linear temperature gradient within the continental crust, though we address the contributions from crustal heat production in a later section of this paper. The Szwillus et al. (2019) and Artemieva (2006) models are based on the same USGS-maintained database as Mooney (2015). In most areas, the middle and lower crust are not in thermal 3 of 20 equilibrium with the conditions of amphibolite and granulite facies, respectively. Temperature profiles suggest that the deep crust has cooled off since reaching peak metamorphic grade. To account for this in our composition calculations, we adopt the assumptions of Hacker et al. (2015), which uses an equilibration temperature of 500°C for amphibolite facies in the middle crust and 750°C for the granulite facies in the lower crust. We performed a temperature correction of 4 × 10 −4 km s −1 C −1 (Christensen & Mooney, 1995;Rudnick & Fountain, 1995) to account for crustal cooling. We also perform a pressure correction of −2 × 10 −4 km s −1 MPa −1 on both facies (Christensen & Mooney, 1995;Rudnick & Fountain, 1995).
The foundation of the tectonic provinces chosen for this global model are the classifications of crust provided by the Crust family of models (Mooney et al., 1998). To further identify tectonic provinces and group together geophysically similar crust, we incorporated crustal thickness, seismic velocity (Vp, Vs), gravity anomaly, sediment thickness, crust elevation, and surface heat flux observations in a tSNE test (t-distributed stochastic neighbor embedding). Results generally favored grouping the continental crust into 8-12 regimes, mostly matching the designations already given in Crust1.0. We augmented these regimes with additional groupings, such as "Thinner Himalyan" crust, when it became clear that the seismic velocity structure of the perimeter of the Himalayas differed from the thickest Himalayan crust. Areas with sparse seismic coverage such as central South America and northern Africa, rely heavily on extrapolation of measurements from similar tectonic provinces. Average Vp and Vs profiles for most tectonic provinces were created from a distribution of tens to hundreds of individual measurements. A notable exception is the "Continental Margins" province, which was represented by >1,600 profiles. Highly localized regimes, such as Andean or Himalayan crust, tended to have <100 profiles due to the uniqueness of their crustal profiles. Figure 1 and Table 1 show the proportion of different crustal provinces by surface area coverage. These tectonic provinces consider only crust exposed at the surface, so that provinces such as "Platform" have underlying crystalline crust that may be Proterozoic or Archean in age. The Proterozoic crust covers the largest fraction (32%) of the continental crust, followed by continental margins (17%). The thicknesses of the crust and lithospheric mantle in Table 1 were calculated from the Litho1.0 Moho and lithosphere-asthenosphere boundary depths. The masses and volumes of the lithospheric mantle are from Wipperfurth et al. (2020)'s physical properties associated with Litho1.0. Sammon and McDonough (2021) serves as our geochemical constraint on the deep (middle and lower) continental crust. That study represents a compilation of major and trace element abundances for amphibolite and granulite facies rocks. We modeled amphibolite facies lithologies for the middle third of the crust and granulite facies lithologies for the bottom third, in agreement with the depth assignment of Rudnick and Gao (2014). We cannot confidently determine a priori which portions of the deep crust are more appropriately represented by amphibolite versus granulite facies data with our current model. In theory, one facies-classified lithology would have greater overall overlap with the seismic velocity profile(s), thus determining which is the more accurate rock type to use. In practice, however, amphibolite and granulite facies lithologies of the same SiO 2 abundance tend to have similar seismic velocities (see Section 3.1). As such, we have assumed a metamorphic Figure 1. The weighted area proportion of crustal types, or "tectonic provinces," used for our model as (a) a fraction of total crust and (b) a fraction of continental crust. Proterozoic crust is most abundant (32% of the continental crust), followed by continental margins (17%) and Archean crust (10%). Modern and paleo-orogens, including arcs, make up a combined 15% of the continental crust in our model. Though trace elements do not participate in thermodynamic calculations, we were able to estimate trace element abundances based on a joint probability analysis with the mineral-forming major oxides. Samples were placed into bins based on the abundance of the oxide and trace element of interest (e.g., SiO 2 and U). Bin width was selected using Sturge's rule (N bins = log 2 (N) + 1). For each major oxide composition bin, there was then a correlated trace element abundance distribution. The mapped distribution of our 13 crust types and (b) the seismic velocity profile data distribution from the USGS Global Seismic Structure Catalog. Data coverage is greatest in the northern hemisphere while places with less coverage, like Africa and Antarctica, rely more heavily on extrapolation of crust type.

Model Uncertainties
Errors in seismic and/or geochemical inputs will skew results. It is imperative to understand the uncertainties in the input datasets, if we want an accurate picture of the uncertainty of our crust compositional models.
The program also will not assess modeling errors stemming from foundational assumptions about what types of lithologies should be used as geochemical inputs and the tectonic provinces assigned to the global crust. These two assumptions are expected to control the systematic error of the model, which is why we made the program flexible and modular. Our approach facilitates testing different crustal models and highlights the projected differences in crust composition.
The primary sources of model error stem from uncertainty in the assumed crustal temperature gradient and Moho depth. Again, these are parameters that can be set by the user. For our preferred model, the uncertainty on Moho depth is on the order of 10% or less in most areas of the global model (Szwillus et al., 2019). The temperature uncertainty is much greater. Global Moho temperatures are taken from Artemieva (2006), which reports no uncertainties. Therefore, uncertainty is taken as the standard deviation of all temperatures found within a given crustal tectonic province (discussed below), and the model runs a Monte Carlo simulation of at least 10,000 iterations to produce a distribution of Moho depths and temperatures. Future results could be improved with Moho temperature models that quantify uncertainty more directly.
We have also attempted to mitigate the bias introduced by oversampling of particular geochemical compositions. An over-sampled composition, such as 100 input compositions with nearly identical major oxide content  artificially inflates the probability of that composition in our final combined model. However, we do consider the reporting of compositions to be at least somewhat reflective of the proportion of rock types present in the deep crust, that is, if the distribution of reported compositions is bimodal, the rocks in the deep crust are likely bimodal in composition. Therefore, we only considered a sample redundant if its oxide content differed from another's by <3 wt.% (9 major oxides, using the distance between vectors formula d = √ 2 1 + 2 2 + ⋯ + 2 , where x n is the difference in wt.% of an oxide between two samples), and its Perple_X generated values for Vp, Vs, and Vp/Vs were within uncertainty of each other.
The internal error contributed by calculational uncertainty is minimal. The overlap between of seismic velocity measurements and Perple_X-derived seismic velocities is calculated via trapezoidal numerical integration at intervals determined by the uncertainty in the seismological data. When this interval is too large to use for the integration, the program reduces the interval by half. The precision errors of Perple_X are generally negligible compared to the uncertainty on our other inputs (Connolly, 2005). The systematic uncertainty of Perple_X, however, remains difficult to assess. As the software is built upon empirical measurements and thermodynamic equations, the empirical database for crustal minerals is incomplete and relies on ideal thermodynamic relationships to compute expected physical properties. The empirical databases themselves also vary with systematic bias.

Quality, Expense, and Time: Global Versus Local Models
In numerical modeling, there is often a tradeoff between computation time and model resolution. For a global view of the continental crust, breadth and total model coverage may be more valuable than high data resolution, especially if results can be averaged over large areas. A large-scale, globe-encompassing model, however, comes with the choice of either short computation time and low resolution or longer computation time and higher resolution. Alternatively, those interested in a more in-depth analysis of a localized region may be able to accommodate higher resolution models. We suggest considering the following when determining whether to use a global or local scale model: data resolution (especially in seismic velocity profiles), data coverage, and model application.
Those with data resolution on the scale of >0.5° × 0.5° should consider using the global version of the script. Data sources with higher resolution, such as that provided by the Earthscope USArray, the AUSArray, or the J-ARRAY, should use the regional scale model. For the remainder of this study, we will analyze global model results. Sammon et al. (2020) presents an example of a local-scale composition analysis using a nascent version of this method.

Empirical Composition-Velocity Trends
For granulite and amphibolite facies lithologies seismic velocities most strongly correlate with SiO 2 content because of its high abundance, as compared to all other oxides. Perple_X-calculated Vp and Vs values at given pressure-temperature conditions show a quadratic relationship between SiO 2 and velocity (Figures 4 and 5). The coefficients of the quadratic best fit curve vary for different pressures and temperatures, and are ultimately correlated to the empirical mineral physics datasets used in the Perple_X Gibbs free energy minimization. Amphibolite and granulite facies lithologies span similar Vp and Vs values, though the shapes of their distributions are  Granulite results are shown for 600°C and 0.85 GPa. The color of the data points indicates percent data point density, with the brighter colors indicating more data points. There is more scatter between SiO 2 and Vs than and Vp, but together they can be combined for a tighter constraint on composition than either compressional or shear velocity alone. 7 of 20 marginally different. This is because the velocities of many amphibolite and granulite facies minerals (e.g., garnet, pyroxene, plagioclase) have similar Vp (∼7 km/s) and Vs (∼3.6 km/s) values. Despite considerable scatter in the Vs data, when paired with Vp, a clear trend emerges: increasing SiO 2 abundance leads to decreasing velocities.
Higher Vp values correlate to lower silica content (Figures 6a and 6b). Higher Vp/Vs ratios also have lower silica content, though for a given SiO 2 , there is roughly a 10% spread in Vp/Vs. A slight curve in the amphibolite facies data becomes more pronounced in the granulites, developing an arcuate shape in the Vp/Vs vs. Vp plot. The same trends appear when analyzing Vp/Vs vs. Vs (Figures 6c and 6d), though the data are more acutely curved. For both amphibolite and granulite lithologies, increasing Vs can lead to either an increase or a decrease in Vp/Vs ratio. The maximum Vp/Vs for amphibolite facies lithologies at typical middle crustal P-T conditions, is expected at a Vs of about 3.5-3.8 km/s, a Vp of 6.5-7 km/s, and SiO 2 of 55 wt.%. For granulite, this maximum is expected at compositions closer to 60-63 wt.% SiO 2 . Interestingly, the maximum Vp/Vs in granulite lithologies corresponds to the lowest Vs rather than the highest Vp, suggesting that Vs variations exert a stronger control on Vp/Vs ratios.

Deep Crustal Density
We calculated deep crustal density by tracking the Vp and Vs values from Perple_X that overlapped with our seismological database back to their input samples. Then, we report the Perple_X-derived density of those input samples. We found that the depth-dependent variations of deep crustal densities among the different tectonic provinces appear more similar when normalized to crustal thickness (Figure 7). Uncertainties in density are ±3% for each tectonic province, a result which is likely controlled by the ∼3% uncertainty on velocities. Deep crustal density ranges from 2,700 to 2,780 kg/m 3 at 13 km depth to 3,290-3,340 kg/m 3 at the Moho.
We note that, in order to calculate deep crustal pressure, and thus mineralogy and composition, we already assumed a bulk crustal density of 2,900 kg/m 3 . This initial assumption, though, does not significantly affect our composition results because there is, at most, a calculated pressure difference of <15% caused by using the 2,900 kg/m 3 a-priori density versus our model-generated density. This <15% pressure difference does not greatly change the stable mineral assemblages or velocities calculated by Perple_X for the deep crust.  The average densities and total masses of each tectonic province listed in Table 1 are calculated by this study. The average thicknesses, total surface area, and total volume of each tectonic province are calculated from Litho1.0's Moho and LAB depths, whose uncertainties are ∼12% (c.f. Pasyanos et al. (2014); Huang et al. (2013); Dziewonski and Anderson (1981)). Continental Margin Moho depths are an exception and are from Szwillus et al. (2019) (Section 2.1). The densities of the oceanic crust and lithospheric mantle, not calculated by this study, were assumed to be 2,900 ± 116 kg/m 3 and 3,400 ± 136 kg/m 3 , respectively (Carlson & Herrick, 1990;Huang et al., 2013). The total volume of continental crust is 9.2 × 10 9 km 3 . The continental crust, including submerged continental margins, comprises 44% of Earth's surface area. Multiplying the density of each tectonic province by its volume yields the mass of crust predicted in each tectonic province. We assumed one-third of the crustal column was upper crust with a density of 2,650 kg/m 3 (Christensen & Mooney, 1995). The total mass of the continental crust, 2.32 × 10 22 kg, is 4% higher than Cogley (1984); Albarède (1998); Rudnick (1995); Wipperfurth et al. (2020).

Composition
Our main analysis focuses on SiO 2 abundance and its uncertainties because of its strong correlation to seismic velocities, particularly Vp. The SiO 2 content at typical middle and lower crustal depth intervals ( Figure 8) is given in Table 2. All 9 major oxide inputs (SiO 2 , TiO 2 , Al 2 O 3 , CaO, MgO, FeO T , MnO, K 2 O, Na 2 O) can be found in Table 3 and corresponding maps in Supporting Information S1. We use the notation "M x% ," where x% is the percent distance to the Moho (M) from the surface. This notation normalizes comparisons among tectonic provinces of various thicknesses. We found that when we normalize depth in this manner, the composition versus depth profiles for each province become more similar. The top of the deep crust starts at an intermediate composition, globally ranging from 58 to 68 wt.% SiO 2 , and the composition gradually decreases to 50-55 wt.% SiO 2 as it approaches the Moho. Global scale SiO 2 composition of the continental crust mostly decreases (or remains steadily mafic) with increasing depth for all tectonic provinces (Figure 9). Uncertainty in global SiO 2 also decreases with increasing depth due to fewer samples fitting the seismic signal at depth. However, the Andean and Himalayan tectonic provinces have larger uncertainties at depth because of the variation in geochemical data fitting the seismic signal and the sparsity of deep seismological profiles.
The CaO content of the deep crust is also of interest due to its absolute abundance and significance as a contributor to sedimentary deposits, though only siliciclastic rocks and not carbonates were considered viable deep crustal components (Hartmann et al., 2012;Wilkinson et al., 2009). In our model, Ca is mostly contained in plagioclase, pyroxene, and garnet. The CaO abundance tends to increase with depth because of the increasingly mafic nature of the deep crust, and therefore regions of low SiO 2 correlate with regions of high CaO. Globally, the median CaO at crustal depths of M 85% is 9.1 ± 3.1 wt.% ( Figure 10).
We can also derive the global distribution of a trace element if that trace element has a quantifiable relationship to one of the thermodynamic components (major oxides) used in our model. We used a geochemical database of samples with both major and trace element concentrations (Sammon & McDonough, 2021) to generate trace element maps as a function of major oxide abundances. We used a bivariate probability analysis to generate trace element distributions from a major oxide abundance, such as SiO 2 , at a specific depth or location. Although we suggest using regional analyses for high resolution interpretations of trace element abundances, we present here global predictions and uncertainties for Sr ( Figure 11) and U ( Figure 12) content based on their relationships with CaO and SiO 2 , respectively, as examples. Global average Sr increases with increasing CaO until plagioclase, whether due to compositional changes or increased pressures, is no longer the dominant Ca-bearing mineral. Uncertainties on global U concentration span an order of magnitude because the abundance of U in a given metamorphic sample set ranges from a few hundreds of ppb to a few ppm. U and SiO 2 abundances, however, are positively correlated, with median U increasing as median SiO 2 increases. Figure 9 shows steady or decreasing SiO 2 with increasing depth in the deep crust. Figure 8 also makes it apparent, though, that the absolute SiO 2 at a given depth is not equal across different crustal types. For example, "Extended" crust appears mafic at 30 km depth while the "Thick Himalayan" crust is felsic at that depth, and "Proterozoic" crust falls in between. However, a more laterally consistent trend appears when comparing percent of the crustal column traversed rather than absolute depth (Figure 8). Most regions show a 5-10 wt.% decrease in median SiO 2 through the deep crust regardless of crustal thickness, so that SiO 2 decreases more rapidly in areas of thin crust than in areas of thick crust. We predict the global median SiO 2 at 50% above the Moho (or, alternatively, 50% crustal column thickness) to be 61.2 ± 7.3 wt.% SiO 2 with CIPW normative mineralogy of <10 wt.% alkali feldspar and <15 wt.% quartz. The average composition of middle continental crust is therefore expected to resemble that of a quartz monzonite; the lower crust average composition, with 53.3 ± 4.8 wt.% SiO 2 and 9.1 ± 3.1 wt.% CaO, is expected to resemble a gabbronorite. It must be noted, however, that these compositions can be the result of seismic velocity profiles averaging felsic and mafic endmembers of the deep crust, and need not be produced by a monzonite or norite.

SiO 2 and Overall Deep Crustal Composition
The processes of magmatic differentiation and density sorting could produce the compositional structure of the continental crust. The process of crustal genesis leaves mafic, restitic material at the base of the crust. Magmatic Figure 8. Global SiO 2 composition at a depth of 30 km shows regional distinctions, whereas measuring composition at a crustal depth relative to the Moho (M 85% notation = 85% of the total crustal depth) produces a view of a deep crust that is relatively more homogeneous, with SiO 2 decreasing gradually with depth. Areas of high projected SiO 2 include the Himalayas, Andes, East African rift, and some continental margins. While the apparently more felsic nature Himalayas and Andes may be a genuine compositional feature, the high SiO 2 in some rifts and continental margins are likely from model input inaccuracies. differentiation can explain the composition correlation with % depth in the crust on both the thin and the thicker scales. More buoyant, felsic material ascends to the top of the crust, producing a gradient of SiO 2 that scales with crustal thickness. Alternatively, the deep crust could be more mafic because it is simply closer to the mantle and therefore has a greater number of mafic intrusions. Our results do not indicate any need for sharp compositional boundaries in the deep crust. The M X% notation reinforces the importance of scaled, relative depth in the crust rather than absolute depth for making compositional comparisons.
Two regions that appear conspicuously more felsic than the global deep crustal median are the Andes and the Thin Himalayan crust (Figure 9). A low temperature gradient could be the cause of the compositional difference between these two tectonic provinces an the global average, but we also must consider two other possibilities, particularly around the northern and northeastern Tibetan Plateau and Himalayan ramp. The first is that thick, convergent margins, especially in the Himalayas, might have layers of upper crustal material thrust deeper within the crust, incorporating more felsic material into the deep crust. In contrast, underthrust upper crustal material is less likely to appear in the Andes, which is a continent-ocean subduction zone. Alternatively, pockets of melt and partially melted material in the Andean and Himalayan middle and lower crust could reduce the shear wave velocity (Caldwell et al., 2009;Hacker et al., 2014;Nelson et al., 1996;Regis et al., 2016;Schilling & Partzsch, 2001;Schmitz et al., 1997;Searle et al., 2009). Because our current model does not factor in melt, slower Vs speeds would be attributed to a more felsic composition instead.
Other anomalous regions in Figure 8, particularly the continental margins of Antarctica, the East African rift zone, and the Sea of Japan, are likely caused by inaccurate temperature and Moho inputs. The East African Rift could appear felsic because the model's temperature gradient for that actively rifting region is too low; a cooler felsic composition can produce the same velocities as a warmer mafic composition. On the other hand, the highly localized, extremely felsic borders around Antarctica and between Japan and China likely indicate a misclassification of crust type and/or Moho depth. Oceanic-type crust has been documented in both regions (Cho et al., 2004;Gohl, 2008;Hirata et al., 1992;McCarthy et al., 2020). Better Moho and temperature resolution of the ocean-continent transition should increase the accuracy of compositional models in these regions.
Mafic granulite lithologies reach gravitational instability in the lower 10%-20% of the average crustal column (Jagoutz et al., 2011), surpassing the upper mantle's density of 3,300-3,400 kg/m 3 . Therefore, according to Figure 7, most of the granulite facies lower crust for continental margins, Andean crust, Tethyan crust, and Phanerozoic crust should be gravitationally unstable. On the other hand, most other tectonic provinces would just reach mantle-like densities around the Moho depths. Thinner Himalayan type crust has a middle crustal density ∼9% lower than other provinces, correlating with negative seismic velocity anomalies. Arcs have the next lowest densities on average, suggesting that a potential denser lower crust beneath some arcs has already foundered (Jagoutz et al., 2011). Figure 7b shows that accreted arcs of Andean type crust in particular display a stark decrease in density interpreted as being due to delamination of the lowermost crust (Ducea, 2011;Gao et al., 2021;Kay & Kay, 1993).
Forming continental crust via island arc processes, however, would then require the deep crust to become denser over time, since most of our tectonic provinces have a lower crust that is denser than arc lower crust. This can be achieved by cooling the crust, thickening it further, intra-crustal differentiation, or by mafic igneous injections into the lower crust. Moreover, if our   There is a tradeoff between temperature and composition. Vp and Vs both carry a temperature dependency through their bulk and shear moduli, so accurate temperature estimates are imperative for modeling the crust; decreased seismic velocities can be the result of either higher temperatures or greater SiO 2 content. As a reminder, this study implements a linear temperature gradient through the crust from the TC15 global temperature model (Artemieva, 2006). Table 4 reports compositional models for the middle and lower continental crust, a practice that is required to make meaningful comparisons to previous crustal models. While we recognize the assumption of a three-layer crust as an oversimplification of the diversity of crustal compositions, it is useful for some calculations to have average composition numbers for the crust.
A few examples include mantle tomography studies which require crustal correction, crustal corrections for geoneutrino studies; models of Earth's thermal history; and planetary scale compositional model for comparison with other rocky bodies. Our middle crustal composition falls between two possible compositions given by Hacker et al. (2015): the fastest Vp endmember composition for the middle crust (62.7 wt.% SiO 2 ), and the middle crustal composition expected when the crust takes on a two compositional layer (upper and lower) structure, instead of three, (57.3 wt.% SiO 2 ). These SiO 2 estimates also overlap with the 62 wt.% SiO 2 reported by Christensen and Mooney (1995) and the 63.5 wt.% SiO 2 middle crust reported by Rudnick and Gao (2014). Our proposed lower crust composition is agrees with Rudnick and Gao (2014) and other mafic models (e.g., Hacker et al. (2015)'s fast Vp lower crust; Jagoutz and Schmidt (2012)). Models which predict a more intermediate-felsic lower crust, such as the North China craton lower crustal model of Liu et al. (2001) or the higher SiO 2 , lower Vp options listed by Hacker et al. (2015), are not consistent with our global average, though isolated regions of more felsic lower crust may exist.

CaO and Sr
Bulk CaO concentration increases with depth ( Figure 10), but as a component of mafic, siliciclastic rocks, not carbonate. The lack of carbonate is due, in part, to our imposed amphibolite/granulite grade lithology restrictions on possible deep crust composition, but is reinforced by high density and Vp values observed in the deep crust. Carbonates, with densities of approximately 2,750 kg/m 3 and Vp's of 6.6-6.8 km/s at deep crustal conditions (Christensen & Mooney, 1995), cannot substantially contribute to the observed deep crustal velocities. Also, there are few carbonate-dominated granulite facies xenoliths and terrains as compared to silicate granulites (Hartmann et al., 2012;Wilkinson et al., 2009). A comparison of Figures 8b and 10 shows tight correlation globally between regions of high SiO 2 and low CaO (see Figure S26 in Supporting Information S1). Uncertainties in CaO estimates follow the same trends as SiO 2 uncertainties, though the former's uncertainty is roughly 10%-20% higher. The relationship between Vp or Vs and CaO content has more scatter than with SiO 2 , which increases the % uncertainty on CaO estimates. Sr abundances cannot be directly derived from velocity calculations, but they can be predicted from Sr's geochemical relationship with CaO or Na 2 O. In amphibolites and granulites, the CaO content tracks predictably with Sr concentration until CaO abundance reaches 5-6 wt.% ( Figure S28 in Supporting Information S1). However, once CaO exceeds 6 wt.%, Sr abundances plateau and scatter increases. This shape of the Sr versus CaO distribution Figure 9. Global SiO 2 decreases with increasing depth from the middle to the bottom of the continental crust. The middle crust M 50% ranges from 60 to 65 wt.% SiO 2 in most areas and increases at a rate of about 1 wt.% per km until reaching the base of the crust. Uncertainties can be found in Figure S19 in Supporting Information S1. correlates with Perple_X-calculated anorthite number. Sr abundance depends on both CaO and Na 2 O, and a Sr map can be predicted from either. Figure 11 shows Sr abundances derived from a bivariate probability analysis between Sr and CaO. The extremely low (dark blue) Sr abundances in thick regions of the deep crust could be a result of low CaO abundances overall and/or the production of pyroxene and garnet at the expense of Ca-rich plagioclase. A subtler drop in Sr is seen in these areas when predicting Sr from Na 2 O abundance ( Figure S33 in Supporting Information S1), possibly because regions of thick crust maintain relatively constant amounts of Na 2 O (and Na-rich plagioclase). While Figure 11 appears to show abrupt changes in median Sr content, note the 50% uncertainty. The variable Sr median abundances shown on the map are a result of the scatter in the Sr-CaO relationship. On the contrary, a Sr abundance map derived from Na 2 O ( Figure S33 in Supporting Information S1) has lower apparent resolution because of the high percent uncertainty on Na 2 O, but less overall uncertainty because of the tighter correlation between Na 2 O and Sr. These findings highlight the usefulness of being able to constrain a trace element from relationships with multiple major oxides, but also stress the need for a higher resolution understanding of crustal composition and mineralogy in unique tectonic provinces.

Heat Production and Moho Heat Flux
Radioactive elements produce heat inside Earth, with K, Th, and U, the heat producing elements (HPEs), accounting for 99.5% of the total radiogenic energy output . The few remaining radioactive elements in the planet contribute <0.5% to Earth's total heat production. We can model the abundance and distribution of HPEs in the crust, and Moho heat flux, with our model and supplementary geochemical relationships.
HPE abundance data are provided in the (Sammon & McDonough, 2021) geochemical data set. Using a joint probability analysis, we derived U abundances with uncertainties from SiO 2 . The middle continental crust has a median U concentration of 1.3 ± 0.2 ppm and the lower crust has 0.22 ± 0.03 ppm. The uncertainty on these median values is reported as the standard error of the median. We do this to highlight that this is the confidence in the central value, the median. The full range of possible abundances spans five orders of magnitude ( Figure 12b) Th abundances could also be generated with a joint Th-SiO 2 probability analysis, but we use a constant Th/U mass ratio of 3.77 ± 0.1 (Sammon & McDonough, 2021;Wipperfurth et al., 2018) which provides a narrower constraint on Th concentration. K 2 O abundance is calculated directly as a thermodynamic component in our compositional model. We calculate heat production on the same 1°latitude by 1°longitude cell scale that we use for calculating major element abundances. The concentrations of K, Th, and U in each cell are multiplied by each cell's density (calculated in previous section) and volume. The resulting masses of K, Th, and U are then multiplied by the heat produced during their respective decays , resulting in a heat production value for each cell in units of W/kg. Heat production values can be converted to W/m 3 via multiplying by cell density once more. Uncertainties are calculated through standard propagation of errors on concentration and density. Figure 13 maps the total (depth-integrated) median contribution to surface heat flux from the upper, middle, and lower continental crust in mW/m 2 .
We calculated a median heat production of 0.3 nW/kg (0.8 μ W/m 3 ) for the middle third of the crust and 0.05 nW/kg (0.1 μ W/m 3 ) for the lower Figure 10. Global CaO abundance and uncertainty at 85% of the total crustal depth. Areas of low CaO correlate to areas of high SiO 2 . There does not appear to be any correlation between CaO content and uncertainty, with most regions having 3-4 wt.% uncertainty regardless of CaO abundance. Figure 11. Global Sr abundance and uncertainty derived from a joint probability analysis with CaO at 50% of the total crustal depth. third of the crust. Previous studies have also found low heat production in the deep continental crust (Fountain et al., 1987;Jaupart et al., 2016;Kukkonen et al., 1997). Our model is also consistent with local studies of HPE analyses of deep crustal xenoliths, such as Gruber et al. (2021); Pinet and Jaupart (1987); and Ashwal et al. (1987).
The uncertainties on this global model are dominated by uncertainties on U abundances, which span an order of magnitude. Even so, our uncertainty on the median or central value (standard error of the median) of HPEs or heat production is well constrained at ± 0.1%. While the range of possible heat production values span an order of magnitude, the median/average heat production value is tightly constrained. Provinces with high predicted SiO 2 content, such as the Andes, have estimated U content up to four times higher than the global lower crustal median because of the correlation between high SiO 2 samples and high U. We recommend using regional HPE data for understanding smaller scale variations and reserve this study's results for continent-or global-scale models.
We took our calculations a step further by estimating Moho heat flux. Once we determined deep crustal heat production, we calculated a subcontinental Moho heat flux by subtracting continental heat production from the global surface heat flux of Lucazeau (2019), supplemented with Shen et al. (2020)'s Antarctica surface Figure 12. Global U abundance derived from a joint probability analysis with SiO 2 at 85% of the total crustal depth. Uncertainties span orders of magnitude because of the range of possible U values, but the global median at this depth is ∼0.2 ppm U. Regions of high SiO 2 , especially the potentially inaccurate continental margin of Antarctica, correlate with high U abundances and the highest uncertainties. heat flux. We assumed a Gaschnig et al. (2016) upper crustal heat production of 1.75 μ W/m 3 to complete the initial calculation. Pairing an upper crustal composition of Gaschnig et al. (2016) with our deep crustal composition yielded Moho heat fluxes for tectonically stable regions that agree with Jaupart et al. (2007) and marginally agree with (Hacker et al., 2015)'s models of crust with slow Vp. However, Moho heat flux calculations depend substantially on the assumed HPE abundance model for the upper crust, which contributes ∼70% of the total crustal heat production in most regions. Using Gaschnig et al. (2016)'s upper crustal U and Th abundances in low heat flux, cratonic regions, results in roughly 6% (by area) of the continents having a negative heat flux across the Moho (an unreasonable condition) -or more likely, other factors, such as heat dissipation through fluid circulation in the near surface, would be needed to explain these low surface heat flux (e.g., 20-40 mW/m 2 ) regions. Alternatively, the assumed upper crustal heat production values may need to be lowered. Most low heat flux areas coincide with stable cratonic lithosphere, where low heat flux and low heat production is not a new observation (e.g., Nyblade and Pollack (1993); Kukkonen et al. (1997);Jaupart et al. (2007); Cammarano and Guerri (2017)). Estimates of upper crustal heat production in cratonic regions range from 0.6 to 1 μ W/m 3 (Gruber et al., 2021;Jaupart et al., 2014Jaupart et al., , 2016Mareschal & Jaupart, 2013;, so we approximated Archean upper crustal heat production at 0.8 μ W/m 3 , which is the maximum permissible heat production value found by Rudnick and Nyblade (1999) for the Kalahari craton and also the maximum average crustal heat production expected for crust ≥2.5 Ga (Jaupart et al., 2016). Our final model for Moho heat flux uses the data and resources listed in Table 5. The median global continental Moho heat flux, shown in Figure 14, is 21.6 +16.0 −5.6 mW/m 2 . Figure 15 summarizes the distribution of Moho heat flux values for each tectonic province.    Figure 16 compares the different models of U abundance in the bulk continental crust from Table 6 to their corresponding Moho heat flux values.
Notably, all models reported in Table 6 assume a similar compositional model for the upper crust (Gaschnig et al., 2016;Rudnick & Gao, 2014). A 1% change in upper crustal U has eight times the effect on the Moho heat flux as a 1% change in lower crustal U content. Thus, having an accurate model for the upper crust is of paramount importance to determining Moho heat flux. Only after upper crustal heat production is constrained will the uncertainties on deep crustal heat production impact Moho heat flux calculations. . Continental crustal heat production is dominated by upper crustal heat production but with significant contributions from the middle crust.

Parameter Value
Global surface heat flux Inputs from Lucazeau (2019) Antarctica surface heat flux Inputs from Shen et al. (2020) Upper crust heat production 1.75 μW/m 3 (Gaschnig et al., 2016) Upper crust heat production (cratonic) 0.8 μW/m 3 (see discussion for sources) Deep crustal densities Inputs from this study Table 5 Heat Production Calculation Parameters

Conclusions
We have constructed a global compositional model of the deep continental crust by synthesizing seismic, temperature, heat flux, and geochemical data. We predict deep crustal compositions on the global scale using major and trace element compositions from amphibolite and granulite facies lithologies, and seismic velocity profiles. Our proposed global compositional model reconciles the USGS Global Seismic Structure Catalog of crustal velocities, compositions for thousands of medium and high grade metamorphic rocks, and constraints on Moho depth (Pasyanos et al., 2014;Szwillus et al., 2019), crust temperature (Artemieva, 2006), and surface heat flux (Lucazeau, 2019;Shen et al., 2020).
Vp values, and to a lesser extent Vs and Vp/Vs values, correlate with bulk rock SiO 2 content, and SiO 2 can be used as a predictor of seismic velocities if temperature can be estimated accurately. Globally, SiO 2 concentrations tend to decrease with increasing depth, leading to a predominantly mafic to intermediate-mafic base of the crust. The decreased density and less mafic nature of the lower crust in younger and tectonically active crust, such as arcs and active mountain ranges, suggests that those tectonic provinces are hotter than our temperature model predicts, that they have undergone lower crustal delamination, or both. Global median SiO 2 concentrations for the middle and lower crust are 61.2 ± 7.31 and 53.3 ± 4.8 wt.%, respectively, though steady composition and velocity gradients in the deep crust urge us to embrace a less distinctly layered view of the crust. This mid-to-deep crustal gradient in wt.% SiO 2 is the equivalent of a lithological gradient ranging from quartz monzonite to gabbronorite. We predict the abundances of multiple major oxides, many of which are correlated to trace element abundances. This correlation allows us to derive expected heat production in the deep crust. We also predict a Moho heat flux of 21.6 +16.0 −5.6 mW/m 2 . Note. Data sources: RG13 = (Rudnick & Gao, 2014), HKB11 = (Hacker et al., 2011), HKB15 = (Hacker et al., 2015), WSM20 = , Šrámek et al. 13 = (Šrámek et al., 2013). Moho heat flux is calculated by taking the average surface heat flux for the continents 63 mW/m 2 (and the oceans 65 mW/m 2 ) (Lucazeau, 2019) and subtracting the continental crust's (and ocean) radiogenic contribution. Upper crust contribution is calculated using the upper crust composition from Rudnick and Gao (2014) relative to the total crust production; note all models assume this same upper crust composition.

Table 6
Heat Production in the Continental Crust Figure 16. There is a linear dependency between crustal U concentration and calculated Moho heat flux. When the different models of Table 6 are compared, differences in each models' predicted Moho heat flux are mostly the result of differences in upper crustal heat producing element concentrations.

Data Availability Statement
Data and modeling software can be accessed at https://doi.org/10.5281/zenodo.6967908.