The robustness of Sr/Y and La/Yb as proxies for crust thickness in modern arcs

Convergent-margin magmatism contributes significantly to the growth of Earth’s continental crust (Rudnick, 1995; Kelemen et al., 2003, Davidson and Arculus, 2005). An essential process involves the partial melting of upper mantle and the subsequent modification of the resulting melts to produce material with compositions that are broadly similar to bulk continental crust (the andesite model [Taylor and White, 1965] and its modifications [Gill, 1981]). Details of these processes (e.g., mass balancing of fractionated products and their cumulates) are continually being refined, but the fundamental tenets of these processes are supported by the broad similarity of major and trace element compositions between modern-arc igneous rocks and continental material, suggesting that similar processes generated both. Given that continental crust can form in arcs, a parameter to characterize crustal formation is the rate of crust addition. This requires that the volume of arc crust be estimated. Volume is a function of thickness, and the present-day crust thickness is most reliably obtained by crustal reflection and refraction techniques. But these techniques are expensive, and only a fraction of Earth’s convergent plate margins have been studied in this way. These approaches are further complicated because the subarc Moho is commonly not marked by a sharp P-wave velocity (Vp) increase to 8 km/sec. This is because this region is composed of hot and partially melted mantle and delaminated cumulates, which act to diffuse the seismically defined crust-mantle boundary and make the Moho more difficult to discern (Arndt and Goldstein, 1989; Shillington et al., 2004, 2013; Kodaira et al., 2007). Other geophysical approaches for estimating crustal thickness such as gravity modeling and receiver-function analysis complement the active-source techniques but are subject to uncertainties. Studies by Dickinson (1975) and Coulon and Thorpe (1981) have shown that variations in some arc-lava major element contents correlate with crustal thickness, providing a non-geophysical method for crustal thickness determination. Refinements of these empirical studies have expanded into correlation of Moho depth with certain trace elements (Dhuime et al., 2015; Chapman et al., 2015; Profeta et al., 2015; Hu et al., 2017). This type of crustal-thickness estimation comes with the possibility of inferring paleo–crust thicknesses. In this study, we examine trace element (TE) correlations with crustal thickness and compare these TE proxies’ crustal thickness estimates against geophysics-derived Moho depths. We test the TE depth proxies’ limitation and accuracy and extend their use by subjecting them to statistical assessments designed to quantify crust thickness variations on the scale of arc segments. As we will show, these correlation methods are not foolproof, but in the best cases they match and complement estimates based on geophysical techniques. We also use our results to examine the implications of the Alaska–Aleutian Islands (northwesternmost USA) arc system’s crust thickness variation within the context of crustal evolution models.

The utility of arc TE crustal-thickness proxies depends on the partitioning of certain TEs (Sr, Y, La, and Yb) as melts interact with minerals through the lithospheric and crustal column through which they traverse, from the region of melt generation in the mantle wedge to the near surface (we exclude TEs from slab melt, for reasons discussed below). Here we briefly outline the origin of these correlations.

New oceanic and continental arc crusts are derivatives of mafic magma generated from partial melting of the mantle wedge (Kimura, 2017). These mafic primary melts are generated by decompression melting accompanying flow in the mantle wedge and by flux melting when the subarc geotherm exceeds the wet solidus, which is in turn controlled by fluid released by the subducting oceanic slab. Primitive mafic magmas commonly underplate at the base of the crust (Annen and Sparks, 2002) where the magma’s geochemistry is further modified by processes summarized as coupled assimilation–fractional crystallization (AFC; DePaolo, 1981; Spera and Bohrson, 2001) and magma assimilation, storage, and homogenization (MASH; Hildreth and Moorbath, 1988), which drive the mafic melts toward intermediate compositions. Additionally, there may be processes such as sinking of mafic-ultramafic cumulates from the lower crust (Kay et al., 1994; Lee and Anderson, 2015; Behn and Kelemen, 2006; Jagoutz and Behn, 2013; see Discussion section below) as well as accretion and underplating of buoyant subducted material (Kelemen and Behn, 2016). Magmas evolve further as they ascend and inject into the lower and mid-crust, where they change the local geotherm, fractionate, and cause partial melting and mixing (Annen et al., 2006a, 2006b). They may then reside in one or more magma chambers and may further evolve before they are ultimately extruded from volcanoes or emplaced as plutons. Clearly, these processes make interpreting the geochemical components of a given igneous rock or lava sequence, from its inception as a mantle melt through its interactions with the crust en route to the surface, challenging and uncertain. However, the variability of these processes can be constrained and smoothed by statistics such that the complexity in their behaviors is averaged out, allowing for broad patterns of chemical changes with pressure to emerge.

Here we concentrate specifically on the TE ratios Sr/Y and La/Yb of arc lavas and their relationship to different pressures in the upper mantle and crust beneath arcs where these melts form and evolve. These relationships arise because TE concentrations are by definition rarefied in magmas, and Henry’s law describes their activities, so linear equations govern mineral-melt partitioning during partial melting and crystal fractionation. Because the stability fields of equilibrium mineral assemblages are known and their TE partition coefficients can be approximated, analysis of key TE ratios can reveal the conditions under which magma generation and differentiation occur. A bonus is that ratios of highly incompatible TEs (e.g., K/Rb, La/Nb) are conserved during magma fractionation, and therefore the initial ratios in the source region and primary melt (if not complicated by contamination and mixing) are preserved in the sampled lavas (White, 2013). Sr/Y and La/Yb ratios of arc igneous rocks are related to pressure due to the different affinity of these TEs to garnet ± amphibole and plagioclase (Moyen, 2009; Davidson et al., 2007). Garnet and plagioclase are stable at high (>1.4 GPa) and low (<0.4 GPa) pressure, respectively (amphibole is stable between 0.5 and 1.5 GPa for a typical arc geotherm), so La/Yb and Sr/Y in arc lavas reflect the pressure at which melting and fractionation occurred (middle diagrams of Figs. 1A, 1B). Partial melting or fractionation at high pressure (>∼1.4 GPa), where garnet peridotite is stable, will result in a high-Sr/Y melt due to the combined actions of (1) stable garnet ± amphibole acting as a sink for the heavy rare earth elements (HREEs) and Y, and (2) the absence of plagioclase causing Sr to behave as an incompatible element (Moyen, 2009). At lower pressure (<∼1 GPa), the absence of garnet ± amphibole causes Y to behave incompatibly at the same time that stable plagioclase absorbs Sr, leading to a melt with low Sr/Y. Similar processes control the La/Yb ratio, with garnet ± amphibole absorbing Yb relative to La and plagioclase having no such effect. If we assume that partial melting of the mantle occurs just below the base of the arc crust and that MASH and AFC processes happen in the deep hot zones in the lower and mid-crust (Annen et al., 2006a, 2006b), we can use the systematics of these two sets of TE ratios to infer fractionating assemblages and thus, indirectly, crustal thickness and depth of the Moho.

High-Sr/Y, high-La/Yb lavas of the far western Aleutians (west of ∼175°E) (Kay, 1978; Yogodzinski et al., 2001; Kelemen et al., 2003; Yogodzinski et al., 2015) reveal processes likely related to highly oblique subduction and melting of the subducted slab (and subsequent interaction with the metasomatized upper mantle) producing andesites with variable silica (54% to 65% SiO₂) with high Mg# (>0.45) (Moyen, 2009). The Sr/Y and La/Yb ratios of these lavas likely do not reflect the processes we are interested in, and they do not form a statistically relevant population in the TE depth correlation studies used here; for these reasons they are excluded from this study (see the Results section for the studied region of the Alaska-Aleutian arc).

We note that there are additional processes that control subduction-related Sr/Y and La/Yb ratios, although these are probably second-order effects. We look at possible factors that could cause deviations from the TE ratios–depth correlation in the Other Explanations for Differences in Crustal Thickness subsection of the Discussion section.

Global correlations of arc Moho depth with variation in arc lava major elements (MEs) have been established by a number of studies (Coulon and Thorpe, 1981; Leeman, 1983; Arculus, 2003; Mantle and Collins, 2008; Plank and Langmuir, 1998; Turner and Langmuir, 2015a, 2015b). Coulon and Thorpe (1981) concluded that crust thickness largely controls arc lava composition. They identified a crustal-thickness threshold that separates dominantly tholeiitic volcanism on thin (<20 km) crust and dominantly calc-alkaline volcanism on thick (>20 km) crust. Leeman (1983) reported a relation between the silicic content of arc magmas and arc crustal thickness. He regressed global arc lavas’ percent of andesite-dacite-rhyolite with arc crustal thicknesses and found a logarithmic fit with high r (correlation coefficient ∼0.8) and concluded that arc magma evolution must scale with the amount of magma-crust interaction. These conclusions are corroborated by recent studies (e.g., Farner and Lee, 2017; discussed below).

While the early studies involved MEs, more recent works noted that some arc lava TE ratios also correlate to the thickness of the underlying crust (Chapman et al., 2015; Chiaradia, 2015; Dhuime et al., 2015; Profeta et al., 2015; Farner and Lee, 2017). In particular, Chapman et al. (2015) showed that the TE ratio Sr/Y in intermediate and felsic whole-rock samples correlates linearly with crustal thicknesses to ∼70 km. Profeta et al. (2015) showed similar correlations with low-MgO calc-alkaline rocks and added a power-law correlation between La/Yb and crustal thickness (Figs. 2A, 2B). Profeta et al.’s results are empirical curves fitted to large data sets of arc lavas from global geochemical repositories (e.g., GEOROC, http://georoc.mpch-mainz.gwdg.de/georoc) with crustal thicknesses derived from the CRUST1.0 global model (Laske et al., 2013) and individual published studies (Profeta et al., 2015, and references therein). While the correlation is derived from global arcs and their median TE ratios, the Profeta et al. (2015) proposed that the relationship could be extended to the scale of individual arcs. They demonstrated this possibility by superimposing Sr/Y and La/Yb versus Moho depth for individual volcanoes from the Southern Volcanic Zone of the Andes onto their Sr/Y and La/Yb versus Moho depth regressions (red dots in Fig. 2) and noted that these fall subparallel to the regressed curve. The interpretation is that the depth proxies of individual volcanoes (from at least this arc segment) conform to the correlation for global arcs, thus extending the correlations to a finer scale. The present study expands on Profeta et al.’s results and tests these correlations to infer crustal thickness along strike for three circum-Pacific arc segments: the Central Volcanic Zone of the Andes, the Central America arc at Nicaragua and Costa Rica, and part of the Alaska-Aleutian arc. We present the results as depth-to-Moho estimates or, equivalently to within a few kilometers, crustal thicknesses for these arcs, compare them with geophysical constraints, and examine implications of our results.

While Profeta et al. (2015) interpolated their correlations from the arc scale to the scale of individual arc volcanoes, Farner and Lee (2017) investigated similar correlations at a more granular level. They argued that if arc crusts are in isostatic equilibrium, then an arc’s elevation reflects the total crust column thickness beneath it via an empirical relation derived by Lee et al. (2015), and elevations of lava samples are proxies for the crustal thickness. Then it is straightforward to construct correlation of the inferred crustal thicknesses at the sample locations with the samples’ TE ratios. Farner and Lee’s correlation is between mean TE values of samples binned into 10-km³-volume bins (see below) and their derived crustal thickness via the isostasy relation. In Figure 3 we superimposed the La/Yb versus crustal thickness curve of Farner and Lee onto the one from Profeta et al. The two regressed curves are subparallel and coincide for a large portion of the range of crustal thickness. The Profeta et al. curve is data heavy for thin crust because it is skewed to more oceanic arcs (where marine seismic experiments can be conducted) whereas the Farner and Lee curve is derived from a more evenly distributed data set (closer to being homoscedastic) because the authors could infer arc crust thicknesses merely by knowing the samples’ mean elevations above sea level. There is an offset of ∼3 km in Moho depth (or a difference of ∼6 in La/Yb value) between the two curves for thick crust, but in general, the curves are within each other’s 95% confidence value. Notably, the two studies differ in the length-scale treatment of the data: Farner and Lee’s analysis involves taking the mean TE ratios from a collection of samples of a unit volume of crust [length × width = 10 km × 10 km and binned to their elevation ∆(h) = 0.1 km] and using isostasy to obtain thicknesses, whereas Profeta et al.’s analysis relates median TE values of whole rocks and median Moho depths (from geophysics) of entire volcanic arcs. The quantitative agreement of the two studies, each using a different method and performed at different data scale, bolsters confidence that correlations of TE ratios and crustal depth are robust. We propose that, given that the two approaches for estimating crustal thickness give comparable results, the correlations are plausible. In the present study, we further test the TE–crustal thickness relationship by applying the correlation to individual samples of volcanic and plutonic rock and use the rock sample population of entire arcs to derive crustal thickness variations along arcs.

To reproduce and assess Profeta et al. (2015)’s results and to obtain a confidence level for the statistics of their approach, we reanalyzed their data (Fig. 4). We produced least-squares best-fit lines and Pearson’s coefficients (represented by r) for Sr/Y versus crust depth using Monte Carlo simulation of their data assuming Gaussian distributions for their quoted 1σ errors. We calculated 10,000 bootstrapped regressed lines (∼60 are shown in Fig. 4) and the corresponding correlation coefficients (r). These form population distributions that we use to estimate confidence intervals. Our median simulated derived slope and intercept are 0.96 and 11.8, respectively, compared to Profeta et al.’s 1.11 and 8.05. The difference between the analyses is 10%–30%, which corresponds to a depth difference that is generally <3 km (5 km maximum) for most values of Sr/Y. More importantly, our Monte Carlo derivation of r values gives their median as 0.84, with confidence interval between 0.61 and 0.92 (lower right of Fig. 4). This means that the correlation level of significance over the null hypothesis for the 22 arcs analyzed is >99% and a median goodness of fit (R²) is 0.71. We therefore accept the correlation between Sr/Y and depth, as conceptualized by Profeta et al. (2015) and Chapman et al. (2015). A similar result was obtained for the La/Yb versus depth correlation. For what follows below, we will use the correlation equations (Sr/Y and La/Yb versus crustal thickness) of Profeta et al. (2015).

Our analysis depends on the assumption that globally derived TE behavior holds at arc-segment and volcanic-sample scales. We also rely on the assumption that the observed local-scale geochemistry of rocks is part of a rational representative of the underlying whole population. We accept that lavas from a given volcano are variable and may or may not follow Profeta et al. (2015)’s correlation. We show below that when we consider all of the Sr/Y and La/Yb ratios from all available lava samples of an arc and treat these statistically, we can map meaningful thickness variations for the arc using the correlation. We do this by employing a regression and smoothing algorithm to the sample suite and bootstrapping for accuracy estimates. We show that such a procedure is equivalent to regressing the median value of each volcano for that arc; both methods give nearly the same regression curve. We repeat the procedure using lava samples from three circum-Pacific arcs and compare each with geophysics-derived Moho depths to assess the robustness of our technique as well as the validity of the TE depth proxies.

Data for volcanic samples were downloaded from EarthChem (https://www.earthchem.org), which includes data sets from GEOROC, PetDB (http://www.earthchem.org/petdb), and the U.S. Geological Survey (https://mrdata.usgs.goc/ngdb/rock). A typical download includes all major elements, all available trace elements, isotope ratios, volatiles, and location, among others. Harker-type diagrams were used to inspect the general geochemistry of the samples, and altered samples and other anomalous samples were deleted by inspection. We used robust statistics (e.g., median rather than mean) for our analysis to minimize the influence of outliers. Extreme outliers were evaluated individually for inclusion or deletion. Finally, sample TE data used for crustal thickness estimation were subjected to filters following Chapman et al. (2015) and Profeta et al. (2015) in order to utilize their crustal-thickness proxy correlation equations. One exception is that we did not apply the Thompson tau test for outlier used by Chapman et al. (2015). A brief description of the filter used on the samples for inclusion is the following: major element totals between 97 and 103 wt% and SiO₂, MgO, and TE ratio Rb/Sr ranges of 55–68 wt%, 1–6 wt%, and 0.05–0.2, respectively.

In devising methods to characterize arc geochemistry from a population of samples accurately, we strove to minimize biases and sampling effects that may skew the data set. The large number of samples gives some protection against bias from outliers. In addition, the estimators we used (median, boxplot, non-parametric regression, etc.) are less influenced by outliers than standard estimators (e.g., mean, standard deviation). The non-parametric LOWESS (locally weighted scatterplot smoothing) regression estimator described below uses polynomial regression in one of the steps that is model based, but bootstrapping (see below) the regression adds robustness to the procedure.

When locations of lava samples are plotted as a function of distance along the arc, the distribution is usually uneven because volcanoes are unevenly sampled due to accessibility, regional politics, and other factors. In this situation, geochemical characterization of a given arc will be skewed by the more densely sampled volcanoes, which are overrepresented relative to the population. We sought a method to de-bias the overrepresented volcanoes. One standard method is to “normalize” the geochemical value to the location parameter, e.g., 54 wt% SiO₂ per 10 km³ of an arc (e.g., see discussion of Farner and Lee [2017] methods in the Robustness of the TE–Moho Depth Correlations section above). The method we employed here to mitigate sampling bias is tied to the method of estimating accuracy when characterizing lava composition: the weighted bootstrap.

The LOWESS estimator (Cleveland, 1993) is a nonparametric regression devised to extract patterns in bivariate plots. It is nonparametric in that it does not rely on, for example, the sample population to be Gaussian, or any other a priori distribution. This estimator is well suited for characterizing Earth chemistry, as we do not expect a parametric control in sample chemistry distribution (but if it is present, the LOWESS regression is likely to detect it). It is not in the scope of this study to detail LOWESS operation, but a brief description is apt. A LOWESS regression fits a smoothed curve to characterize a set of bivariate data. Each particular point of this curve is calculated from data in the neighborhood of that point. The width of this neighborhood controls the smoothness of the final curve, and this width constitutes one of two input parameters for the procedure. The LOWESS regressed value at this point is the fit of a weighted polynomial regression of the data within its neighborhood. For this study, we used a degree 1 linear regression (least squares) with a tri-cubic weight function throughout. The sharpness of the weighting function is the other parameter of the LOWESS. The LOWESS regressor can be thought of as a generalized “moving average”–type regression. Its salient features are that it is nonparametric, adaptive, and outlier resistant. We will further discuss the merits of using LOWESS regression in the next section.

To estimate the accuracy of the curves, we use the weighted bootstrap technique. Bootstrapping (Efron and Tibshirani, 1993) is an accuracy-estimating procedure that synthesizes many (typically from 100 to 10,000) populations of data sets by repeatedly resampling (with replacement) from the existing data set. These synthesized populations represent the underlying (unobservable) data of which the actual data are a subset. The relevant statistics that were performed on the original data set are performed on these synthesized populations to obtain a distribution from which the variance and accuracy of the original data set are extracted. The bootstrapping distribution gives a measure (in many cases, a confidence interval) of the statistical accuracy of the samples’ measurement. Bootstrapping also allows for correcting bias in the resampling process (see Sample Location Bias section above). To do this, we assign a weight to each sample during the bootstrap resampling. A sample’s weight is based on the sample’s proximity to all other samples of the arc system: a sample with many neighbors is weighted lower than a more isolated sample (see caption of Fig. 5 for equation). This sample weight influences the probability of the sample being selected in the bootstrap resampling, with the overall effect of a more even selection of samples among all the volcanoes. The sample-biasing mitigation can be seen in data handling for the Alaska-Aleutian arc (Fig. 5), specifically mitigating the larger number of samples in the Katmai region with the much smaller number of samples in the more remote western Aleutians (Fig. 5B). In Figure 5C, a weighted bootstrap resample data set shows that the Katmai data have been down-weighted and the western Aleutian data up-weighted. (If there is a small number of anomalous samples, they will be up-weighted, but in practice, the smoothing process will tend to negate this type of leverage in producing the final, coherent curve.) This procedure results in a more balanced sampling of individual volcanoes for the whole arc, from which statistical analysis is performed. We repeat this process thousands of times to derive a distribution of the statistics leading to a less-biased data set when operated on by the LOWESS estimator.

Here we demonstrate that we can subject lava samples to LOWESS regression to construct the geochemical variation trend along an arc, and that doing so is equivalent to obtaining a trend by plotting median geochemical compositions for each volcano. The difference between using LOWESS regression of samples and taking volcano medians is nevertheless significant. This is because the LOWESS regressed curve takes values from a percentage of nearby samples to derive a value for any one point, so that the generated curve does not allow for complete independence of samples from one location from those of its neighbors. Although it may seem restrictive that the constructed curve is so constrained and not completely independent, this is in fact what we assume when we calculate, for example, the mean silica weight percent of a volcano: we assign a mean to that volcano and infer that all of the lavas generated from it are related to (or constrained to be near) that mean. In LOWESS regression, we take this idea beyond a single volcano and say that at any single locale, lavas are related to, or are influenced by, a certain percentage of lavas from other nearby locales. This means that, geochemically, the magmas originated from a process that was “simple” (i.e., partial melting of the mantle) but was acted upon by other processes that increased the variance of the magma compositions. For example, fractionation, anatexis, mixing, etc., contribute to the scattering of, say, a volcano’s sample coordinates in a Harker diagram. However, due to the statistics of large numbers, perturbations that tend to increase an oxide’s value are counteracted by others that tend to decrease it, so that LOWESS regression, which seeks the median of these scatters, can recover and better characterize the original signal. We propose that LOWESS regression can better represent the underlying geochemical signature of the relevant region by deemphasizing the processes that scatter the sample data.

To demonstrate the near equivalency of LOWESS regression and volcano medians, we apply this regression to SiO₂ content versus location along arc for a suite of Alaskan-Aleutian arc lavas. The curve generated uses LOWESS regression on weighted bootstrapped resampling (Fig. 6A) with 5% and 95% confidence curves that bound the LOWESS estimate of the median value curve. All lava samples (n = 3250) are plotted as black dots. Note that the uneven density of samples of the Alaska Peninsula volcanoes (relatively dense) versus the Aleutian Islands volcanoes (relatively sparse) are reflected in the spacing between the 5% and 95% bounding curves, with a narrower spacing corresponding to denser data. This LOWESS regressed variation curve is interpreted here as the typical silica value at a locale along the arc. It clearly shows a gradual decrease in silica from the western Aleutian Islands volcanoes to a sharp inflection at the continental shelf break (approximately the continent-ocean boundary), then a sharp increase (at ∼165°W) from the peninsula to the continental interior. Figure 6B shows boxplots for the median SiO₂ content of individual volcanoes. The LOWESS regression curve is shown for reference and to demonstrate the equivalency of the two plots. The inset shows a histogram of the residuals, the difference between the volcano median and the LOWESS curve at each longitude. The residuals histogram is close to being a Gaussian distribution with mean near zero, indicating that the prediction of the LOWESS regression curve is not systematically biased against the prediction of the individual volcano medians. A feature of the LOWESS regressed curve that is distinct from the volcano medians is that the large number of samples guards against the curve being influenced by a few outliers, whereas in the volcano-median plot, a small volcano with an outlier median value may unduly influence the overall trend. That the generated curve gives a succinct and clear graphical representation (with accuracy estimate) of the underlying data is the main reason we propose that this type of curve characterizes the geochemical value against distance along an arc more accurately than other methods.

We have shown above that bootstrapped LOWESS regression is potentially useful for elucidating lava geochemical variations along an arc such as the Alaska-Aleutian system. We now address the question: What is the likelihood that such variation comes about from pure chance? To do this, we have devised a test to evaluate the possibility that the curve-generating procedure indicates a chemical variation–versus–distance relationship when no relation exists in actuality. For this test, we use the Alaska-Aleutian arc La/Yb proxy depths and, as before, perform LOWESS regression to examine their variation against along-arc location (longitude). We pose the null hypothesis that the La/Yb variations have no dependency on longitude. If the null hypothesis is true, it means we can randomly permute the samples’ longitude values, perform the LOWESS regression on this longitude-randomized sample set, and expect a resultant curve that is equivalent to the original unpermuted data set. We repeat this process 1000 times, then examine the distribution of these curves and compare them with the curve generated with the original data. Figure 7 shows a few dozen of the location-randomized La/Yb-depth variation curves along the arc (green curves). There are statistical variations and “kinks” in the curve, where data density is low and a few samples control the curve. If we take the mean of the one thousand location-randomized curves, we obtain a near-horizontal curve (dark green) fluctuating near the mean value of the La/Yb of the population (here equal to 26 km crustal thickness), as is expected for a randomly located set of samples. Comparison of the randomized data with the actual variation of La/Yb depth versus longitude (red curve, copied from Fig. 13) shows how different the real variation is, compared to the ones made with the randomized data; the actual La/Yb depth variation with longitude is such that this curve barely intersects the field defined by the randomized curves. The LOWESS regression curve generated from the data is significantly distinct from those that come from chance alone. Thus the null hypothesis is false: the La/Yb variations along this arc do depend on longitude. The fact that there are two populations of La/Yb ratio (or Moho depth) amongst the samples is evident in the two distinct peaks in the histogram (left side of Fig. 7). How these samples (as depth proxies) are actually distributed along the arc is the subject of the Results and Discussion sections.

In general, given enough data, we can characterize lava geochemical variations along an arc using weighted bootstrap LOWESS regression as described in this section. These characterizing curves can extract trends in noisy data sets, can correct for known bias, and come with estimation of accuracy.

Here we show results of crustal thickness estimates of the three arcs: the Central Volcanic Zone of the Andes, Central America, and Alaska-Aleutian. First, we show the studied regions (Figs. 8, 9, and 10 for Central Volcanic Zone of South America, Central America, and Alaska-Aleutian, respectively) superimposed with sample locations along with the tectonic boundaries. Seismic-derived Moho depths are also shown and referenced to the associated studies. Note that our curves are for the volcanic front (where the samples are located) but seismically determined Moho depths are for a much broader region. The relevant tectonic boundaries that define crustal thicknesses are shown. In the Alaska-Aleutian case (Fig. 10), note the location of our boundaries (“west,” “transition,” and “east”): they differ from “western” and “eastern” denotations in past studies. Also, we emphasize that our Alaska-Aleutian studied area excludes the far western Aleutians (what is referred to by others as “Western Aleutian”) because we only consider the arc associated with a subduction zone as defined by the extent of the seismically detectable slab mapped by Syracuse and Abers (2006). In the Alaska-Aleutian arc, the slab is not defined west of ∼180°, where the convergence velocity tends to very low value and the tectonics are dominated by strike-slip faulting. The key TE correlations that we address do not apply there.

Crustal thicknesses for the Andes Central Volcanic Zone were calculated from the individual samples’ La/Yb and Sr/Y ratios based on the Profeta et al. (2015) correlations and regressed against arc strike (Fig. 11). Confidence-interval curves (5% and 95%) are shown bounding the median (50%) curve in the figure. Both the individual samples’ La/Yb and Sr/Y depth proxies and the derived curves plot close to each other and are in phase as they track changes along the northern (Altiplano) segment of the arc, with estimated crustal thickness of 68–72 km. South of ∼18°S near the Puna-Altiplano boundary, there is a change in slope from flat to a decrease in thicknesses (of −3 to −6 km per degree latitude for La/Yb and Sr/Y, respectively) to 21.5°S, where the Sr/Y proxy predicts a shallower Moho depth than do the La/Yb proxy and the geophysics estimates by ∼17 km [mean of (La/Yb)_n depth and geophysics depths versus Sr/Y depth]. In general, the northern portion of the geochemically derived crustal-thickness curves compares well with estimates from the McGlashan et al. (2008) teleseismic study and from CRUST1.0 (Laske et al., 2013) (we show both surveys’ data as LOWESS regressed curves made with the same method we used for the TE ratios). The density of samples along the arc, with few gaps, results in high fidelity of the TE proxy curves. The boundary between the Puna and Altiplano appears in the TE proxy curves as a difference in crustal thickness but is not as well resolved in the geophysically derived Moho depths. The TE-ratio depth proxies decouple in the Puna region, with Sr/Y depths ∼20 km shallow than La/Yb depths and geophysics depths. This reflects the complex nature of the petrogenesis of the region as documented by Kay et al. (1994) and Kay and Coira (2009), where ignimbrite, calc-alkaline, intraplate, and shonshonitic centers are superimposed in a small region (23°S–27°S). The Altiplano region is less complicated; here the two TE proxies agree (to within 10 km) with the McGlashan et al. (2008) seismic depth, but show ∼15 km greater crust thickness than CRUST1.0. The Sr/Y proxy depth departs from the La/Yb proxy depth and geophysics depths by ∼15 km south of 19.5°S but is consistent with thinning of the crust south of the Altiplano (Kley and Monaldi, 1998; Kay and Coira, 2009). In general, the TE proxies give plausible Moho depths compared to those available from geophysics. The Sr/Y proxy may reflect the effect of plagioclase in cumulates or restites in lowering Sr.

We constructed LOWESS curves of crustal thickness derived from La/Yb and Sr/Y along the Central American volcanic arc, from latitude 9°N to 15°N (Fig. 12). Here we compare our TE-based crustal thickness estimates with geophysically based estimates from CRUST1.0 (Laske et al., 2013) and MacKenzie et al. (2008). Both the Sr/Y and La/Yb Moho depth estimates decrease from the thick Chorotega block going north into the Nicaragua basin. The two geochemical proxies show different crustal thicknesses, with the La/Yb proxy predicting deeper Moho than the Sr/Y proxy by ∼20 km. The La/Yb proxy predicts thick crust (50+ km) beneath the southern Chorotega block compared to a thinner (35–40 km) estimate given by the Sr/Y proxy. The Sr/Y depth proxy predicts thinner crust beneath the center of the Chorotega block compared to the other estimates. The reason for this discrepancy may relate to excess plagioclase, or it may reflect the plume-influenced nature of magmatic processes beneath the southern Central American arc (Gazel et al., 2011). Both Sr/Y and La/Yb give slightly thinner crust beneath the Nicaragua basin as compared to the geophysically based estimates but converge with the CRUST1.0 curve north of 12°N. The two geophysical approaches show more moderate changes (compared to the TE-based proxies) in crustal thickness beneath Central America. CRUST1.0 (Laske et al., 2013) indicates a gradual decrease in Moho depth northward from 38 to 28 km within 5° of latitude, while the MacKenzie et al. (2008) study shows a similar decrease of ∼8 km (39–31 km) northward along 3° of latitude. The two geochemical proxies for Moho depth diverge in detail but consistently show modest crustal thinning from the Chorotega block to the Nicaragua terrane.

We apply our method of estimating crustal thicknesses beneath volcanic arcs using TE ratios to the Alaska-Aleutian system (Fig. 13). We note that the Alaska-Aleutian convergence system has an overriding plate that is oceanic in the west (Aleutians) and continental (Alaska Peninsula) in the east. We are especially interested to know whether the geochemical proxies show a significant crustal thickness difference for these two different crustal thickness end members. Figure 13 shows the Moho depths inferred from La/Yb and Sr/Y along the arc. Crustal thicknesses agree between the two TE proxies, and the curves are coherent and in phase. The proxies predict a crust thickness of ∼20–25 km between 174°W and 160°W, ∼30 km crust west of 174°W, and notably, a ramp-like increase in thickness from ∼23 to 35–40 km east of 155°W. This contrasts with Moho depth estimates modeled from the Alaska Seismic Experiment (a two-ship reflection and refraction seismic survey in 1994; lines A1, A2, and A3 in Fig. 10; see Holbrook et al., 1999, for details), which show Moho depths of ∼28 km to ∼35 km at 172°W and 164°W (lines A1 and A3; Holbrook et al., 1999, and Van Avendonk et al., 2004), and along-strike Moho depth of 30 ± 4 km (line A2; Fliedner and Klemperer, 1999). Moreover, the CRUST1.0 (Laske et al., 2013) LOWESS regression curve shows Moho depth increasing eastward with a near-linear slope, from ∼19 km at 180° to ∼38 km at 153°W. A recent receiver-function study by Janiszewski et al. (2013) placed the Moho at a (regressed) near-constant 39 km depth from 160°W to 177°W, almost 20 km deeper than the TE proxies’ predictions. Discrepancies between the geochemical proxies and the geophysics-derived crust thickness are further addressed next in the Discussion section.

How reliable is the TE-based method of applying global correlation to a high-spatial-resolution and geochemically variable study of an individual arc? We note that in the case of Central America and the Andes Central Volcanic Zone, when the two TE ratios are coherent, they tend to agree with geophysical estimates, but when the two TE ratios disagree, they depart from geophysical estimates. We conclude that statistically treated geochemistry-based proxies can provide useful estimates of crustal thickness and are complementary to geophysical methods when estimates for the two ratios agree. Given this caveat, we look in more detail at the Alaska-Aleutian arc result and discuss the implications of the disagreement between our coherent TE-inferred Moho depths and seismically inferred Moho depths. We examine the geochemistry of this arc as it pertains to crustal thickness, then crustal thickness as it relates to volcanic output. Finally, we combine both analyses with a possible crustal model that may reconcile the thickness difference between the results of geochemistry and geophysics.

The Aleutian-Alaskan arc erupts both tholeiitic (TH) and calc-alkaline (CA) lavas (Kay et al., 1982; Kay and Kay, 1994; George et al., 2003; Mangan et al., 2009), but does the systematic relation noted by Coulon and Thorpe (1981) between CA-TH suites and crustal thickness exist in the Alaska-Aleutian arc? If so, what are the geodynamic implications? Recently, Farner and Lee (2017) observed in their global compilation that there was a correlation between elevation (their proxy for crustal thickness; see the Robustness of the TE–Moho Depth Correlations section) and calc-alkalinity. In this section, we apply LOWESS regression to construct curves to characterize CA-TH affinity of volcanoes and rock samples and correlate these with the crustal thickness proxies derived above. Figure 14 shows our resultant LOWESS regression curves of published CA-TH indices (the tholeiitic index of Zimmer et al. [2010] and calc-alkaline–tholeiitic index of Hora et al. [2009]) generated from volcanic samples from our compiled data set. Zimmer et al. (2010)’s tholeiitic index (THI) is a measure of the Fe enrichment of a volcano expressed in the ratio of FeO at 4% MgO over FeO at 8% MgO. THI is a per-volcano measurement, whereas Hora et al. (2009)’s calc-alkaline–tholeiitic index (CATH) assigns an index to individual rocks based on Miyashiro (1974)’s separation of arc TH and CA rocks (see defining equations in Figure 14 caption).

Alaska-Aleutian igneous rocks define two distinct populations, with most plutonic samples having CA affinities (Kelemen et al., 2003; Cai et al., 2015). For this reason, we performed the regression with only volcanic samples. For the per-volcano LOWESS curve generated by the THI index of Zimmer et al. (2010), TH volcanoes are preferentially located along the transitional segment of the arc, with CA volcanoes to the west (Aleutians) and east (Alaska Peninsula). The CATH index of Hora et al. (2009) shows a similar result with a few exceptions (see Fig. 14). We conclude that Aleutian-Alaska TH arc lavas are associated with thinner crust (as inferred from our results) of the transitional segment between ∼161°W and 174°W, and that CA lavas are associated with the thicker crust of the oceanic and continental arc on either side (Fig. 13). The regression parameters we used are sensitive to long-wavelength changes, so we did not reproduce the results of Kay et al. (1982), in which CA and TH volcanoes are associated with the variable stress regime in and around rotated tectonic blocks. They proposed that TH volcanoes concentrate between rotated tectonic blocks, where magmas evolved at low pressure and ascended through extensional basins, whereas CA volcanoes are concentrated on the blocks, associated with a thicker crust and higher-pressure magmatic evolution. This CA-TH relationship associated with tectonic blocks in the Aleutians is subsumed into the longer-wavelength variation in the whole-arc treatment of the LOWESS regression of the geochemical data. In our analysis, the entire Aleutian-Alaskan arc from 176°E (Buldir Island) to 152°W (Mount Spurr) manifests a single CA-to-TH-to-CA cycle.

Likely controls on the magma flux for a mature arc are convergence velocity, subducted fluids and sediment, magma plumbing system, and thermal structure of the mantle wedge. For example, Fournelle et al. (1994), George et al. (2003), and others have pointed out a relation between Aleutian volcano volume (inferred to reflect magmatic flux) and convergence velocity. We checked their correlation with a polynomial fit (degree 1 or 2) and found that volume and velocity are related by an r-value = 0.4, with n = 34, a 95% significant correlation. Here we examined crustal thickness as a possible control of magma addition. We assumed that all arc volcanoes are about the same age (Jicha et al., 2006) and that each magmatic cell associated with an individual volcano reflects a similar proportion of intrusive to extrusive rocks. In this case, volcano volume scales with magmatic flux, and we will take the volcano volume as a proxy for flux. We plot the Alaska-Aleutian volcano volume along the arc (Fig. 14, top), then we examine the relationship between volcano volume and crustal thickness. Figure 15 shows the correlation of the arc segment volcanoes’ volume with the corresponding crust thickness beneath each. The negative correlation for the linear fit is significant (at 99% for n = 34, with similar significance for a polynomial fit). If the volcanoes are grouped by the three arc segments (western, transitional, eastern), the larger volcanoes tend to be located on the transitional segment (Fig. 15, volume histogram). The distribution of volcano volume is similar for the eastern continental and western oceanic segments even though the two have distributed crustal thickness values (Fig. 15, thickness histogram). At the same time, the transitional segment volcanoes are underlain by thin crust (21–25 km) and show a bimodal distribution of volumes. These observations reveal that the arc segments host distinct distributions of volcano volumes: the thin-crust transitional-segment volcanoes exist in two volume modes, ∼75 km³ and ∼300 km³, whereas the thicker-crust western oceanic (Aleutian) and eastern peninsula and continental (Alaskan) segments’ volcano volumes are in a positively skewed distribution with median of 40 km³. It is also noteworthy that the transitional segment contains most of the large edifices (seven out of nine volcanoes with volume ≥200 km³) but only four of 25 volcanoes with volume <200 km³.

We performed a randomization test (similar to the test described in the Null Hypothesis Check section) to see if the volcano volume–versus–segment relation arose from chance alone. The test assumes a null hypothesis, in which case randomizing the relation between volume and crustal thickness would not affect the correlation. In fact, the real correlation is >95%, significantly different than in the randomized case. Thereby the null hypothesis is rejected and the correlation is significant. There is no similar correlation of volume versus crustal thickness for the other studied arcs (Central America and Andes Central Volcanic Zone). The proposed correlation between magma addition and crustal thickness is simplistic in that it assumes all magma addition to be expressed in volcano volume, ignoring the cryptic addition of mass by mechanisms such as underplating and “relamination” of subducted material (Castro et al., 2013; Kelemen and Behn, 2016). However, the correlation exists and it is compelling. The thin-crust transitional segment of the arc hosts all of the large-volume (>200 km³) volcanoes. This result suggests that thinner arc crust has a simpler plumbing system for magma migration to the surface, or that the Alaska-Aleutian transitional arc segment is a region that is more conducive to mantle melting. We explore the latter supposition in the next section.

We use derived crustal thickness along the Alaska-Aleutian margin to suggest reasons for the discrepancies with geophysical derived Moho depths in the context of current ideas on arc crust construction. Assuming that the TE-inferred Moho depths are useful approximations of reality, we wonder what is responsible for thickening and thinning the crust between the transition segment and the continental segment of the arc. The LOWESS regressed curve of TE proxy depth predicts that the thin crust of the transition segment would thicken to more continental-like thickness as the arc extends into the Alaska Peninsula, with the thin crust coinciding with the high-volume, TH volcanoes’ location on the transition segment, and the thick crust with the low-volume, CA volcanoes of the western and continental (eastern) segments (Figs. 13, 14, 15). We relate these features to a model of oceanic arc crust evolving to continental crust, integrating observations from the Alaska-Aleutian Lava Affinities and the Alaska-Aleutian Volcanic Activity and Crustal Thickness sections above with current crustal-evolution theories.

It is known that the construction of continental crust from convergent-margin magmatism is a multistage process that transforms basaltic mantle melts to the andesitic composition of bulk continental crust. The processes of fractionation, melting, and mixing to evolve high-Mg# andesite with TE contents that match those of continental crust requires a complementary cumulate in the lower crust or upper mantle (Kay and Kay, 1993; Rudnick and Fountain, 1995; Holbrook et al., 1999; Kodaira et al., 2007). Seismic velocity profiles and fossil arc crust sections (e.g., Talkeetna arc, south-central Alaska; Greene et al., 2006) mostly do not show a mafic cumulate layer in the lower crust for the required mass balance. Jagoutz and Behn (2013) observed that the exposed Kohistan lower crustal section (northwestern Pakistan) is denser than the upper mantle and suggested repeated delamination of the mafic lower crust as a mechanism to evolve andesitic arc crust and the continental Moho. This also addresses the presence of a sharp seismic boundary at the continental Moho (where Vp, the P-wave velocity, transitions from ∼7 km/s to ∼8 km/s in a step function) that is absent for the Moho beneath magmatic arcs. Most magmatic arcs show a thick (∼10–15 km) transitional layer where Vp increases from ∼7.4 km/s to ∼7.8 km/s, which is generally assumed to represent upper mantle (Calvert and McGeary, 2012). However, these velocities are characteristic of lithologies such as pyroxenites and eclogites that are denser than the upper-mantle peridotite, and are weak and thus are likely to founder, as noted by Behn and Kelemen (2006). Similar lithologies compose the crustal section that Jagoutz and Behn (2013) modeled as negatively buoyant and the absence of which creates a sharp Vp contrast between crust and mantle that characterizes the continental Moho. Such a sharp Moho is exposed at the Talkeetna crustal section, whereas the analogous exposed Kohistan crustal section exhibits a more gradational crust-mantle boundary with a (calculated) gradual Vp ramp that may be more typical of arc crust.

Further insights into the nature of the lower crust and upper mantle beneath arcs are provided by measurements of shear-wave splitting or SKS. Shear-wave splitting reveals seismic anisotropy induced by preferred orientation of olivine in the upper mantle (Karato, 2009) or by melt-filled cracks oriented parallel to the maximum compressive stress direction (Yang et al., 1995; Nowacki et al., 2012). Both mechanisms are taken to indicate the direction of mantle flow. Measurement of seismic anisotropy beneath mid-ocean ridges shows that this is related to convection in the uppermost mantle, which is predicted to flow parallel to plate motion and perpendicular to the ridge strike. In arc settings, the fast-polarization directions are commonly oriented parallel to arc strike (see Yang et al. [1995] for seismic anisotropy beneath the Shumagin Islands within the Alaska-Aleutian transition segment), indicating paradoxically that mantle flow is perpendicular to the subduction vector. Behn et al. (2007) proposed delamination (or Rayleigh-Taylor instability foundering) of dense lower arc crust as a mechanism that induced such trench-parallel flow. They modeled gravitationally unstable diapir-like sinking masses of ∼15 km diameter (by a few kilometers thick) spaced ∼40 km apart and observed induced flows that are similarly oriented to SKS fast directions documented in arcs. We propose that the transitional segment of the Aleutian arc is shedding dense lower-crust granulite and pyroxenite (Vp of >7.4 km/s) formed as cumulates by basaltic fractionation at ∼25 km depth, as indicated by the La/Yb and Sr/Y proxies (Figs. 1 and 16). At 15 km diameter and a few kilometers thick, these sinking masses would not be resolved in tomographic images, and their Vp contributions would be averaged into the surrounding mantle.

Furthermore, delamination of these cumulate masses induces upwelling of the subarc asthenosphere, leading to enhanced decompressional melting near the base of the crust (Kay et al., 1994; Behn et al., 2007) and the observed high-volume TH volcanism of the transition segment. Behn et al. (2007) estimated that each downgoing mass could induce mantle upwelling to generate 10 km³ of melt. We propose that such processes are responsible for the enhanced volcanic activity of the transition segment. If the masses are spaced 40 km apart with a descent velocity of ∼1 cm/yr and delamination occurs every 10⁶ to 10⁷ yr (Behn et al., 2007, from arc-magma production estimates of Jicha et al., 2006), the upper mantle beneath the volcanoes of the transition segment would accumulate, in volume, ∼30% cumulate rock, mixed in with the upper-mantle peridotite, becoming a layer of mixed lithologies between lower crust and mantle, giving the high-Vp structure at 20–40 km depth seen in teleseismic studies (Fig. 1). This interpretation is consistent with the results of Shillington et al. (2013), who analyzed Vp/Vs ratio (where Vs is the S-wave velocity) in addition to Vp to determine likely lower-crustal lithologies in the Aleutian transition segment and concluded that a mixture of material is needed to explain the observed high Vp and low Vp/Vs ratio. These masses are smaller than delaminated slabs such as those proposed, say, for the Andes (Puna Plateau; Kay et al., 1994) or the Sierra Nevada (California, USA; Ducea and Saleeby, 1996; Manley et al., 2000; Farmer et al., 2002; Lee et al., 2015) where detachment of region-scaled lower crust and lithosphere occurred with the induction of prolonged region-wide uplift and magmatism. Our proposed scheme is more on the order of a dynamic exchange of material across the crust-mantle boundary, cycling at ∼106 yr; in this situation, uplift is suppressed by the balanced flux of cumulate loss and magma inflow recharge. In our interpretation, the TE-inferred Moho depths mark the transition from the lower crust (at ∼25 km) into a transitional lithological layer of melting and fractionating with a mixture of cumulate and peridotite rock, the bottom of which (at ∼40 km) is interpreted as the Moho in the seismic survey.

Here we address the discrepancy between the TE depth proxy predictions and the geophysics-based estimates (Holbrook et al., 1999; Lizarralde et al., 2002; Van Avendonk et al., 2004; Janiszewski et al., 2013; Shillington et al., 2013; CRUST1.0 model: Laske et al., 2013) without using the crustal construction model discussed in the previous section. The correlations of TE-estimated Moho depth seem robust because they are based solely on the abundance of four elements. The application of LOWESS regression to the correlations acts to average out the variable complexities of these elements in arc processes to exhibit the TE correlation to crustal thickness. This dependency on only four elements is also the correlations’ weakness, for the correlations are vulnerable to systematic variations in those elements, as seen for the decoupled TE ratios in northern Puna of the Andes Central Volcanic Zone and the Chorotega region of Central America. In those regions, the TE proxies for crustal thickness may not be applicable: the multitude of processes operating in the regions decoupled the two TE ratio correlations from showing an unambiguous and valid result. Similarly, for the crust of the Aleutian arc transition segment, there may be local processes acting to disturb the global correlation enough to invalidate its usage. The TE ratios there may reflect, for example, variable mantle wedge chemistry such as Sr differentially leached from the subducted slab. These types of effects may be large compared to the global trend responsible for the correlation to crustal thicknesses. It is also possible that the Profeta et al. (2015) correlation itself may not be applicable to discern crustal thickness within an arc at the resolution we consider here. The global variations exist but have a variance that cannot be scaled down to show thickness variation at the resolution of this study. However, one must still explain the systematically lower Sr/Y for 15° longitude along the Alaska-Aleutian arc that anti-correlates to both sediment input (Kelemen et al., 2003) and convergence velocity. There is also the possibility that geophysical techniques infer a Moho that is too deep: a crustal thickness of ∼40 km at longitude 160°W to 175°W as proposed by Janiszewski et al. (2013) implies that the Moho at longitude 170°W is <20 km from contacting the subducted slab as defined by the depth of the Wadati-Benioff zone of Syracuse and Abers (2006), a geometry that is problematic because it implies a cooler geothermal gradient in the mantle wedge than is capable of producing melt. The corresponding TE-proxy Moho depth of 20 km gives a more realistic 40 km separation between the subducted slab and the lower crust (Fig. 1). Beneath the easternmost continental volcanoes, the Wadati-Benioff surface gives what is probably the minimum separation of ∼35 km for the TE-proxy crustal thickness of 40 km. In the Central American arc, the TE proxies diverge from each other at latitude 10°N–11°N, which is at the limit of the volcano line beneath which the subducting slab ceases to have seismicity and where the influence of the Galápagos plume (Cocos Ridge) may affect arc magma compositions.

We have demonstrated a proof-of-concept method using statistical techniques applied to geochemical data to extend the usage of the Sr/Y and La/Yb TE depth proxies to resolve variations in crustal thickness within arc segments. These estimates are plausible when compared to geophysical refraction and reflection studies and thus serve as a complementary technique to model crustal thickness, especially for where geophysical coverage density is low. We have shown that when the two TE proxies are in phase, they can complement geophysics in the Andes Central Volcanic Zone and the Central America cases. Then we examined an intriguing case in the Alaska-Aleutian arc where the TE proxies–geophysics disagreement necessitated novel interpretation of arc crust construction where the Moho interface of arcs maybe more opened, with influxes of mantle-derived magmas and delamination of fractionated cumulates and restites. This interpretation implies that the Moho beneath an active magmatic arc may be more challenging to identify than generally acknowledged, with geophysical estimates being deeper than those from TE proxies.

We appreciate critical and constructive reviews by an anonymous reviewer and Jeffrey Ryan. This is The University of Texas at Dallas Geosciences contribution number 1327.