Data acquisition and treatment. Geochemical data of bulk volcanic rocks from 16 recent arcs (Supplementary Data S1) collected from the Georoc database (http://georoc.mpch-mainz.gwdg.de/georoc/) were filtered and treated according to the method described by ref. 9 (Methods). Median values of 143Nd/144Nd for bins of 0.5 wt.% MgO display distinct evolutionary trends with typical differentiation indices, such as MgO and Co, for arcs with different crustal thicknesses (Fig. 1; see Figs. S1 and S2 for the plots of all individual arcs). The choice of Co was determined by the relatively well constrainable bulk coefficients of this compatible element for modelling purposes and by the fact that, because Co is not very chalcophile, sulfide fractionation occurring in thick arcs9 is unlikely to significantly affect this element, as it would affect other compatible trace elements (e.g., Ni). Nonetheless, the correlations between 143Nd/144Nd and Co reported here are also found using other compatible and incompatible elements instead of Co as it should be expected due to the fact that strong correlations occur between 143Nd/144Nd and MgO (Fig. 2a) and that MgO is a universal hallmark of magma differentiation that shows strong (positive or negative) correlations with nearly all major and trace elements.
In thin arcs, 143Nd/144Nd displays a flat trend for decreasing MgO and Co, whereas in increasingly thick arcs 143Nd/144Nd values steadily decrease with MgO and Co along slopes that increase with increasing thickness of the arc (Figs. 1 and 2; Figs. S1-S2). Indeed, slopes of the linear trends of 143Nd/144Nd-MgO and 143Nd/144Nd-Co of arc magmas ranging between MgO < 1 and > 10 wt.% correlate significantly with crustal thickness (Fig. 2; see Tables S1 and S2 for the full parameters of the linear correlations). Whereas the significance (p < 0.0001) of these correlations (R = 0.921 for the 143Nd/144Nd-MgO slope and R = 0.914 for the 143Nd/144Nd-Co slope) might seem to be strongly enhanced by the value of the Andean Central Volcanic Zone (CVZ), the correlations are statistically significant also disregarding CVZ (R = 0.796 with p = 0.0001 for the 143Nd/144Nd-MgO slope and R = 0.778 with p = 0.0002 for the 143Nd/144Nd-Co slope). The intercept values of the trends at MgO = 0 wt.% or Co = 0 ppm, i.e., a proxy for the crustal assimilant, also display significant correlations with crustal thickness (Fig. 2c). The correlations for the intercept values using the 143Nd/144Nd-MgO slopes have R = 0.921 with p < 0.0001 including CVZ and R = 0.787 with p = 0.0002 excluding CVZ, whereas correlations using the intercept values of the 143Nd/144Nd-Co slopes have R = 0.922 with p < 0.0001 including CVZ and R = 0.791 with p = 0.0002 excluding CVZ. The two sets of intercept values strongly correlate with each other along a slope of ~ 1 (143Nd/144Nd = 1.0403 + 0.0207, R = 0.995) which indicates that the intercept values for each arc are nearly the same using both 143Nd/144Nd-MgO and 143Nd/144Nd-Co regressions. The Nd model ages of the intercept values calculated for the Sm/Nd value of the bulk continental crust (0.19518), which obviously also display a correlation with crustal thickness, range in age between ~ 200 and ~ 1400 Ma (Fig. 2c).
The values of 143Nd/144NdMgO12 (and corresponding Nd model ages) calculated from the regression trends above for 12 wt.% MgO, a plausible proxy of primary magmas in subduction zones, display less significant correlations with crustal thickness (R = 0.735 and p = 0.0005 with CVZ), especially if CVZ is excluded (R = 0.505, p = 0.0387) (Fig. 2d).
Correlations were also checked for Pb and Sr isotopes. Whereas Pb isotopes did not show any significant correlations, Sr isotopes display a similar behavior to Nd isotopes (Fig. S3; Supplementary Data S2). However, the Sr isotope variation with arc thickness is less significant than that of Nd isotopes resulting in a lower resolving power among arcs, except for the very thick CVZ arc (Fig. S3). Additionally, the Sr partition coefficient during magmatic differentiation may be very different in thick versus thin arcs depending on plagioclase stability with pressure11, adding complexity to the modeling discussed below. In contrast Nd and Co have more straightforward behaviors during magmatic differentiation. Therefore, the discussion is limited to Nd isotopes.
Modelling of Nd isotope systematics of arc magmas. A large number of petrographic, geochemical and isotopic studies have shown that arc magmas are the result of various processes including fractional crystallization, mixing (e.g., through recharge of mafic magmas), and assimilation of host rocks7,8,19. The combined use of isotopic compositions of magmas and incompatible or compatible trace elements, whose behavior can be modeled through their bulk partition coefficients between crystallizing minerals and residual melt, provides first order information on the prevalent magma evolutionary processes occurring in specific volcanic edifices20. The application of such type of modeling to entire arc segments, although rocks may have different ages, may be sourced from more or less different mantle domains and may interact with more or less different crustal rocks, is informative of first order differences in large-scale processes related to crustal architecture and geodynamic setting, including crustal thickness9–14.
Several recent papers have highlighted the role that the crustal thickness of the overriding plate plays in the geochemistry of both parental arc magmas13,14,21 and their derivatives9–12,22. Because this study focuses on geochemical and isotopic data of arc magmas for MgO values between < 1 and > 10 wt.% it necessarily bears on the control that crustal thickness exerts on the intracrustal evolution of magmas. Nonetheless, intracrustal evolutionary paths may also be partly or significantly controlled by the initial geochemistry of parental magmas, which is in turn controlled by mantle source processes13,14. Deconvolving these two processes is a difficult and debated task22.
The control that crustal thickness may play on the differentiation processes of parental arc basalts is intuitive. Parental basalts intruded into thin crust cool rapidly and evolve toward derivative terms dominantly through fractional crystallization processes23. Because of the low temperature of the host rocks, due to intrusion of parental basalts at overall shallow levels in a thin crust, assimilation should be, intuitively, minor. The bulk of the geochemical differences in thin arc magma sequences should be therefore controlled by fractional crystallization and by mixing between more or less fractionated terms. In this situation, radiogenic isotope compositions of the whole evolutionary sequence are expected to remain those of the parental magma with limited variations reflecting minor assimilation of a thin oceanic crust that is anyway likely to be isotopically similar to the parental magma. In a plot of compatible elements versus 143Nd/144Nd, such evolution corresponds to a subhorizontal trend in which 143Nd/144Nd does not vary significantly whereas the compatible element decreases continuously due to its sequestration into fractionating minerals. This behavior is shown by the evolution of thin arc magmas in Fig. 1.
In contrast, parental magmas ascending from the mantle wedge through a thicker crust encounter hotter crustal rocks in the deeper parts of the arc and should interact more significantly with these rocks23, at least because they travel through a thicker warm lower crustal section. Under these conditions, interaction processes, such as assimilation and mixing with partial melts of the crust, are plausible24, in addition to the processes of recharge from continuously incoming basalt from the mantle, the latter feature also occurring in thin arcs. This has been described as MASH processes25 or hot zones19 occurring in the lower crust of thick arcs and might result in systematic radiogenic isotope changes in derivative magmas should the crust into which the parental basalt intrudes be isotopically different from the basalts themselves. Indeed, the trends of intermediate to thick arcs (Cascades, Northern Volcanic Zone NVZ, Mexico) to very thick ones (CVZ) show a continuous, albeit different, decrease in 143Nd/144Nd values with decreasing MgO and Co, which is indicative of such a process (Fig. 1).
Table 1
Ranges of values for the parameters used in the equations describing the AFC model. The range of parent Nd and Co concentrations are from the range of oceanic and continental basalts of ref. 26. The ranges of Nd and Co concentrations of the assimilant are the same as that of the parent for thin arcs and are from ref. 18 for thick arcs and CVZ.
Parameter | Range for thin arcs | Range for thick arcs | Range for CVZ |
Remaining melt | 0-100 | 0-100 | 0-100 |
R | 0–1 | 0.05-1 | 0.05-1 |
Nd parent (ppm) | 10–15 | 14–15 | 14–15 |
Nd assimilant (ppm) | 10–15 | 11–25 | 11–25 |
143Nd/144Ndparent | 0.5130–0.5132 | 0.51293–0.51294 | 0.51293–0.51294 |
143Nd/144Ndassimilant | 0.5080–0.5139 | 0.5114–0.5132 | 0.5115–0.5126 |
Co parent (ppm) | 40–42 | 40–42 | 40–42 |
Co assimilant (ppm) | 40–42 | 20–30 | 20–30 |
DNd | 0.1–0.2 | 0.1–0.3 | 0.1–0.3 |
DCo | 1–6 | 1–6 | 1–6 |
To quantify these processes, assimilation fractional crystallization (AFC) modeling was carried out in the 143Nd/144Nd versus Co space using DePaolo’s27 equations with input parameters of the elemental (Nd, Co) concentrations in the parent and assimilant, the Nd isotope composition of the parent and assimilant, the bulk partition coefficients for Nd and Co, the fraction of remaining melt and the ratio (r) of the mass of the assimilant to the mass of crystallized melt. In particular, the latter value is informative of the thermal status of the crustal level at which the AFC process is occurring, with r = 0 corresponding to pure fractional crystallization, low r values (< 0.3) being typical of the upper crust, high r values (> 0.5) being typical of high assimilation rates in the lower crust and r nearing 1 virtually corresponding to pure mixing27. A Monte Carlo approach was used due to the loosely constrained nature of most of the parameters used in the model (see Methods). Geologically significant ranges used for each of the parameters are reported in Table 1.
Dependency of crust assimilation by arc magmas on overriding plate thickness. The trends with slopes changing for different arc thicknesses in the 143Nd/144Nd-MgO and 143Nd/144Nd-Co spaces suggest that, whereas arc magmas differentiating in thin crust arcs are consistent with a pure fractional crystallization process, mantle-derived arc basalts differentiating in increasingly thick crust require assimilation of crustal rocks, with increasingly lower 143Nd/144Nd values as the crustal thickness increases.
To evaluate the output of the Monte Carlo-based model of AFC, the model results have been filtered to fit the different arc trends in the 143Nd/144Nd-Co space. Figures 3a, 3b, and 3c show the forced fitting of the model to the data available for thin-intermediate arcs (Kermadec and Mariana as examples: Fig. 3a), thick arcs (Cascades and NVZ as examples: Fig. 3b) and the very thick CVZ arc (Fig. 3c). Figures 3d, 3e, and 3f show the density scatterplots returning the most probable values of r and 143Nd/144Ndassimilant in successfully reproducing the 143Nd/144Nd-Co trends of Figs. 3a, 3b, and 3c. The results of r for the thin arcs must be considered maximum values because the subhorizontal trends of all these arcs are consistent with a pure fractional crystallization process without any assimilation, a situation that cannot be modeled in these plots because not having an assimilant would result in no 143Nd/144Ndassimilant values returned by the model and therefore in the impossibility of drawing Fig. 3d.
The results suggest that the most probable solutions for the r values of the AFC process (highest density areas in Figs. 3d, e, f) steadily increase with arc thickness, passing from < 0.25 (most likely close to 0 for the above discussion) in thin-intermediate arcs to 0.25–0.50 for thick arcs and to > 0.75 for the very thick CVZ arc. The 143Nd/144Ndassimilant in contrast steadily decreases from ~ 0.51290 through ~ 0.51265 to ~ 0.51235 for thin-intermediate, thick and very thick arcs, respectively. Solutions for differential assimilation (different r values) of a similar low radiogenic assimilant (e.g., ≤ 0.5123) in both thick and very thick arcs are less probable, according to the model (see density contours in Figs. 3d, 3e, 3f).
Although the absolute r values of the model output have to be considered with caution, due to the uncertainties in the parameters used (Table 1), the model results strongly support the intuitive hypothesis that arc magma differentiation in increasingly thicker arcs occurs in average at deeper crustal levels11,23 and is therefore characterized by an increasing assimilation rate (r value) of the continental crust until a nearly pure mixing trend between crustal melts and mantle-derived components in the CVZ (r approaching 1).
Systematic variations in incompatible elements and their ratios in primitive magmas from arcs of different crustal thicknesses suggest that crustal thickness controls the degree of element enrichment also by modulating the partial melting of the mantle13,14. It has been suggested that under thick continental arcs (and therefore under a thick lithosphere) the degree of partial melting of the mantle wedge is lower than under thin arcs because of a systematic change in the thermal structure of the mantle wedge with crustal thickness of the overriding plate17. In contrast, radiogenic isotopes, including those of Nd, do not seem to be significantly affected by this process13,14, although local correlations of 143Nd/144Nd of primitive arc rocks with crustal thickness have been reported and attributed to systematic along-arc changes in ambient mantle29. The results of the present study show a poor correlation of 143Nd/144NdMgO12 of putative parent arc magmas with crustal thickness, which becomes more significant if the CVZ data are included (Fig. 2d). A possible mechanism explaining the much less radiogenic composition of putative parental magmas of the CVZ could be fertilization of the mantle wedge by subduction erosion of the lower portion of the overriding plate crust, which has been suggested to occur in the CVZ since ~ 19 Myr, with the highest rates during the last 7 Myr30,31. Subduction erosion likely occurs also in other arcs32–34. The poor correlation of Fig. 2d leaves open the possibility of subduction erosion or any other process progressively enriching the mantle wedge source of arc magmas through time (i.e., decreasing its 143Nd/144Nd value) as the arc thickness increases.
Nonetheless, even considering such a potential decrease in the 143Nd/144Nd values of putative primitive magmas of arcs with crustal thickness (Fig. 2d), the comparison between Figs. 2a and 2b shows that the Nd isotope systematics of derivative magmas of arcs increasingly thicker are controlled to a larger extent by the isotopic composition of the assimilant and therefore by intracrustal processes rather than by primary mantle-related differences. This lends strong support to the important role of arc thickness in modulating the elemental and isotopic compositions of derivative arc magmas through intracrustal processes that significantly overprint systematic geochemical trends caused by variable arc thicknesses in primitive magmas, especially in thicker arcs.
Implications for continental crust growth processes. Another interesting and not straightforward outcome of the above data analysis is that the thickness of the arc crust correlates with the Nd isotopic composition of the assimilant (Figs. 2a, 3d, 3e, 3f), suggesting that an increasingly thicker crust is characterized by overall increasingly lower 143Nd/144Nd values. This suggests thickening of the continental crust in arcs by continuous reworking of previously formed (“isotopically older”) crust by subsequent episodes of “isotopically younger” arc magmatism23,28. Although the model ages returned by the assimilants through the regression of the slopes in the 143Nd/144Nd-MgO and 143Nd/144Nd-Co spaces do not have an absolute geological meaning, they roughly represent the time-integrated “average” model ages of the reworked continental crust in the investigated arcs, suggesting reworking of progressively older crust as the arc thickness increases. In its turn, crustal thickness controls arc magma chemistry both at the source and in the crust through a feedback loop process. In fact, increasing crustal thickness, achieved through arc magmatism reworking of previous thinner crust segments, results in primitive magmas with increasingly higher contents of incompatible elements13,14 and perhaps an increasingly crustal Nd isotope composition (i.e., lower 143Nd/144Nd) of the ambient mantle wedge (Fig. 2d), e.g., through subduction erosion or other processes. Parental magmas, during intracrustal evolution, become increasingly contaminated by their interaction with a thicker, older and more felsic crust (Fig. 2c) driving it into an even thicker, more felsic and isotopically evolved crust. The thicker and older the crust the parental arc magmas are interacting with, the stronger the divergence of parental and derivative magmas from parental and especially derivative magmas of previous stages of arc crust build-up (Fig. 4), highlighting the significant role of magmatic reworking in crustal growth processes.
The results of this study, therefore, in addition to vindicating a role of crustal thickness on the chemistry and isotopic composition of arc magmas both at their source and during crustal differentiation, are consistent with views of crustal growth as a continuous process characterized by increasing thickness, SiO2 content and isotopically evolved compositions through time6,28, of which thin-intermediate and variably thick recent arcs could be a present-day snapshot.