Crustal Thickening of the Northern Central Andean Plateau Inferred From Trace Elements in Zircon

The timing of crustal thickening in the northern Central Andean Plateau (CAP), at 13–20°S, and its relationship to surface uplift is debated. Zircon qualitatively records crustal thickness as its trace element chemistry is controlled by the growth of cogenetic minerals and relative uptake of light and heavy Rare Earth Elements. Jurassic to Neogene zircons from volcanic rocks, sandstones, and river sediments reveal shifts in trace element ratios suggesting major crustal thickening at 80–55 Ma and 35–0 Ma, coincident with high‐flux magmatism. An intervening magmatic lull due to shallow subduction obscures the magmatic record from 55 to 35 Ma during which thickening continued via crustal shortening. Protracted thickening since the Late Cretaceous correlates with early elevation gain of the CAP western margin, but contrasts with Miocene establishment of near modern elevation in the northern CAP and the onset of hyperaridity along the Pacific coast, highlighting their complex spatial and temporal relationship.

SUNDELL ET AL. 10.1029/2021GL096443 2 of 9 south which instead shows a coupling between shortening, crustal thickening, and uplift for the last ∼40 My .
The temporal relationship between crustal thickening and surface uplift is central to the debate on what process(es) best explain the generation of high topography and development of orogenic plateaus in Cordilleran orogens. The greater Andean orogen is the archetype Cordilleran orogenic system on Earth in the modern, and determining the crustal growth history is valuable to testing wholistic models of plateau development and explaining the relationships between crustal shortening, magmatism, basin development, surface uplift, and climate. Although earlier studies suggest a shortening deficit at most latitudes of the Central Andes (Hindle et al., 2005;Kley & Monaldi, 1998), more recent studies show that shortening can account for most crustal thickening (e.g., Hernandez et al., 2021). Still, the timing of crustal thickening remains an outstanding question in the northern CAP because the region is mantled by Cenozoic volcanic rocks, thus precluding detailed structural analysis as has been done farther to the south (e.g., Anderson et al., 2017).
Trace element geochemistry provides information on the history of continental crustal growth (and decay) by exploiting the relative uptake of light and heavy Rare Earth Elements (LREE and HREE) in cogenetic mineral phases. Specifically, trace elemental ratios can be used to track relative changes in crustal thickness because they reflect the presence or absence of specific mineral phases formed from residual magmatic melt (Chiaradia, 2015;Heaman et al., 1990;Mantle & Collins, 2008;Profeta et al., 2015;Turner & Langmuir, 2015). Similarly, zircon (ZrSiO 4 ) can be used to infer the chemistry of SiO 2 rich melts because it is refractory and resistant to post-crystallization diffusion (Cherniak et al., 1997), and it forms over a wide range of conditions (Watson & Harrison, 1983). Zircon preferentially incorporates HREEs via xenotime coupled substitution wherein REE 3+ + P 5+ = Zr 4+ + Si 4+ (Finch et al., 2001), and is thus sensitive to other cogenetic minerals that incorporate HREEs such as amphibole and garnet. Zircon with low LREE/HREE can be interpreted to have formed in thin crust where there is an abundance of HREE and limited LREE in the residual melt, as the latter are preferentially incorporated into plagioclase and the former are generally incompatible. Conversely, zircon with elevated LREE/ HREE can be interpreted to have formed in thick crust where HREE are preferentially incorporated into garnet and amphibole and there is an abundance of incompatible LREE. Zircon Eu/Eu* ( ∕ √ × ) is also a function of crustal thickness due to the preferential incorporation of Eu 2+ into plagioclase, and due to elevated Fe 3+ / ΣFe by garnet sequestration resulting in increased oxidation of Eu 2+ to more compatible Eu 3+ (Tang et al., 2021). Collectively, the combination of trace element proxies with U-Th-Pb geochronology (e.g., Gehrels et al., 2008;Schoene, 2014) makes zircon is a useful petrochronometer (Kylander-Clark, 2017) that can be applied to produce time-resolved records of geochemical change in igneous melt chemistry and used to infer qualitative changes in crustal thickness (e.g., Balica et al., 2020;McKenzie et al., 2018;Tang et al., 2021).
We analyzed Jurassic to Neogene zircons extracted from a suite of volcanic rocks, sandstones, and modern river sediments in the northern CAP ( Figure 1). Identifying crustal growth trends works best for igneous rocks and minerals that formed in intermediate melts in order to avoid deep mafic melts that undergo major diversification in the melting-assimilation-storage-homogenization (MASH) zone and highly fractionated melts or melting of metasedimentary rocks in the upper crust (Chapman et al., 2016;Profeta et al., 2015). We interpret geochemical results using a multivariate predictive model to infer zircon melt chemistry following the methods and training data outlined in Belousova et al. (2002). Temporal trends in zircon trace element data are then placed in the context of published regional igneous whole rock data, surface uplift estimates, basin stratigraphy, and paleoclimate records from the Late Cretaceous to Cenozoic. Finally, results are synthesized into a regional tectonic model for the northern CAP applicable to Cordilleran orogens in Earth's deep past.

Methods
Trace element geochemistry and U-Th-Pb geochronology of zircon was conducted using laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) at the Arizona LaserChron Center (ALC) following methods outlined in Balica et al. (2020). We analyzed the trace element concentrations of zircons with published U-Pb ages from late Oligocene to Pleistocene volcanic rocks across southern Peru (Sundell, Saylor, Lapen, et al., 2019), detrital zircons from late Eocene to late Miocene sandstones in the Peruvian Altiplano (Sundell et al., 2018), and detrital zircons from modern river sediments in southern Peru and northern Chile (Pepper et al., 2016;Figure 1b). A small subset of previously unanalyzed zircons from the Altiplano samples were analyzed for both geochronology and geochemistry. In total 760 zircons were analyzed. Outliers were detected and removed if the individual element or ∑LREE/∑HREE results were outside median ±2σ standard deviation for the complete data set, and/or if chondrite-normalized Lu was less than 1,000 ppm, the latter as a proxy for analysis of non-zircon (Belousova et al., 2002). Data filtering reduced the data set size to n = 710, including n = 326 zircons from N = 58 volcanic rock samples, and n = 384 detrital zircons from N = 7 sandstone samples and N = 5 modern river sediment samples. We used a model-based approach to infer igneous melt chemistry from zircon trace element concentrations. Although early work suggested zircon trace element chemistry is not useful for determining whole rock chemistry (Hoskin & Ireland, 2000), Belousova et al. (2002) demonstrated that application of multivariate analysis can be used as a robust predictor of igneous rock type. We apply the classification and regression tree (CART) model developed in Belousova et al. (2002). CART models fall under the broader umbrella of multivariate predictive models in supervised machine learning. Here, multivariate data with known responses are used to "train" a predictive model which can then be applied to new data to make predictions (Kuhn & Johnson, 2013). Belousova et al. (2002) developed a CART based on igneous rock types with known zircon trace element concentrations based on n = 330 zircon analyses from N = 18 different rock types to train a predictive CART (Figure 2a). Passing zircon trace element data through this model produces a predicted rock type at confidence levels above 75% (Belousova et al., 2002). We first evaluate results in comparison to Ti-in-zircon crystallization temperature calculated following methods of Watson and Harrison (2005), then temporally. CART analysis, as applied here, effectively acts as a filter for interpreting temporal trends in zircon melt chemistry by excluding zircons inferred to be from highly fractionated, high SiO 2 melts and zircons formed in deep melts formed during MASH zone diversification. Results were also compared to data compiled in the Central Andes Geochemical Database (Mamani et al., 2010) between 13°S and 20°S to compare zircon geochemistry to regional whole rock geochemical records (Figure 1a).
Because intermediate magmatic compositions and zircon saturation span such a wide range of depths within the crustal column, trace element concentrations are not necessarily representative of the true crustal thickness. Rather, we consider the range of trace element concentrations to be representative of various (and inherently unknown) depths of crystallization, with the highest LREE/HREE and Eu/Eu* ratios approaching the base of the crust. With this in mind, we seek temporal trends by binning the results into 4 Myr intervals and reporting the maximum values for each bin. We chose a 4-Myr interval because this is the smallest possible range given the limited data density of samples older than 40 Ma.
Plotting ∑LREE/∑HREE against Ti-in-zircon crystallization temperature (Figure 2d) shows that analyzed zircons from high SiO 2 melt (70%-75% and >75%) generally formed at lower temperatures and lower ∑LREE/∑HREE, whereas zircons from 65% to 70% SiO 2 melt formed over a wide range of generally higher temperatures and higher ∑LREE/∑HREE. Low ∑LREE/∑HREE is consistent with melt chemistry in thin crust where there is an abundance of HREE and a deficit of LREE, the latter due to its preferential uptake in plagioclase feldspar (Heaman et al., 1990). High ∑LREE/∑HREE results point to residual melts where there is a depletion of HREEs due to competition with increased abundance of cogenetic, high lithostatic pressure (>1 GPa) minerals such as garnet and amphibole (e.g., Mantle & Collins, 2008).
High SiO 2 melts undergo fractionation at shallower crustal levels, and thus do not provide reliable proxies of crustal thickness. Temporally, this can be seen in ∑LREE/∑HREE and Eu/Eu* where there is a clear separation from zircons derived from 65% to 70% SiO 2 melt after 35 Ma; Eu/Eu* separation likely reflects the combined effect of decreased uptake of Eu 2+ in plagioclase and oxidation of Eu 2+ to Eu 3+ (Figures 3a and 3b). Similar trends can be seen in Ti-in-zircon crystallization temperature and Th/U, both of which show an increase in range after 35 Ma corresponding to increased range in crustal thicknesses and a diversification of magma. For the latter, lower Th/U likely corresponds to Th removal by accessory minerals that saturate before zircon such as monazite (Bea, 1996), rather than the presence of metamorphic fluids as only few zircon analyses show log(Th/U) values < −1 (Figure 2b), further supporting the idea that the crustal column became extremely thick after 35 Ma. In the absence of an assumed geothermal gradient, Ti-in-zircon temperatures can be interpreted as qualitative proxies of melt depth; results show a steady increase after 80 Ma, with a pronounced decrease in melt depths between 55 and 35 Ma (Figure 3c). Additionally, La/Yb, Sm/Yb, Ce/Y, and Dy/Yb zircon ratios are useful proxies for crustal thickening, with higher ratios indicating thicker crust, where garnet and amphibole are stable. Zircon data inferred to be from 65% to 70% SiO 2 granitoid melt compositions show an increase in all of these ratios between 80 and 55 Ma and 35 Ma to present (Figures 3e-3h). Available igneous whole rock data earlier than Neogene are limited; however, similar trends of increased whole rock trace element ratios after 35 Ma are interpreted as the result of increased crustal thickness, which appears to have continued to present day (Figures 3e-3h).

Crustal Thickening, Surface Uplift, and Paleoclimate
Zircon trace element data provide a nearly continuous record of crustal thickness change through time (Figures 3  and 4e). Although data are limited between 200 and 80 Ma, they provide a baseline of relatively thin crust. The overlap in zircons inferred to be from 65% to 70% and 70%-75% SiO 2 melts between 200 and 80 Ma followed by a separation after ∼80 Ma supports a later diversification in the crustal column (Figures 3a-3d). Results support contrasting hypotheses based on two apparent pulses of crustal thickening at 80-55 Ma and 35-0 Ma (Figures 3  and 4e). In one scenario, the crust thickened in two separate stages, with an intervening period of crustal thinning (removal). Alternatively, the crust thickened continuously, albeit with subdued magmatism during shallow subduction, thus obscuring any thickening trends that would otherwise be seen in zircon trace element data. We prefer the latter hypothesis based on (a) an inboard stepping of the magmatic arc after ∼45 Ma (Figure 4f; Mamani et al., 2010); (b) a decrease in Ti-in-zircon temperatures at 55-35 Ma, consistent with limited high-T (i.e., deeper) melt production; (c) nearly continuous advancement of the retroarc thrust front since ∼80 Ma (Figure 4g; Data are separated into groups of zircons inferred to be from 65% to 70%, 70%-75%, and >75% SiO 2 melts based on CART analysis (see Figure 2a). Parts e-h include igneous whole rock data (yellow squares) from the Central Andes Geochemical Database (see Figure 1 for spatial range; Mamani et al., 2010). Horton, 2018); (d) Permo-Triassic rocks of the Eastern Cordillera that reached granulite facies at ∼45 Ma ; and (e) increased retroarc foreland basin subsidence rates at 55-35 Ma, suggesting a coupled plate margin (Figure 4i; Horton, 2018); all of which is consistent with the proposed timing and location of shallow subduction and tectono-magmatic model for the northern CAP of Sandeman et al. (1995). This interpretation suggests that magmatism, in addition to crustal shortening, is a dominant driver of crustal thickening (Haschke & Günther, 2003) during normal (i.e., not shallow) subduction.
The absence of proxy data precludes robust estimation of paleoelevation in the late Cretaceous-early Cenozoic. However, there are clues to this early paleoelevation history that can be reconciled with new qualitative estimates of crustal thickening (Figures 4e and 4h). The Altiplano and western Andean margin appear to have different surface uplift histories (Figure 4h; Boschman, 2021). Before 80 Ma, the northern CAP was likely a neutral plate margin undergoing orogenic stasis and sediment bypass with low elevation, as evidenced by the lack of retroarc foreland basin stratigraphy (Figures 4a and 4i; Horton, 2018). Localized surface uplift of the CAP western margin starts at ∼80 Ma (Figure 4h), consistent with 80-55 Ma zircon trace element data which show an increased range of crustal depths associated with different parts of the crustal magmatic column (Figures 3 and 4b). Although marine rocks now in the orogenic hinterland require that at least parts of the northern CAP Altiplano be at sea level in the latest Cretaceous (Figure 4h; Sempere et al., 1997), retroarc foreland basin stratigraphy suggests the development of a Western Cordilleran crustal load of thickened crust and moderate elevation between 60 and 25 Ma resulting in increased sediment accumulation rates associated with the development of a Paleogene retroarc foreland basin system (Figure 4i; DeCelles and Giles, 1996;Horton et al., 2002;Sundell et al., 2018). Late Cenozoic surface uplift in the northern CAP Altiplano may have been partly decoupled from crustal thickening, in contrast to the earlier elevation gain in the western CAP margin (Figure 4h) and Puna region of the southern CAP (Canavan et al., 2014;Quade et al., 2015). For example, phylochronologic analysis of potato parasite nematodes suggests the maximum paleoelevation of the northern CAP between 25 and 18 Ma was no more than half of its current elevation (Picard et al., 2008). Gregory-Wodzicki (2000) similarly required a ∼2 km elevation maximum by early Miocene time based on leaf physiognomy, although this study disregards climatic effects. Garzione et al. (2008) presented δ 18 O isotopic data that allowed for maximum elevations of 2.2 km between 25 and 10 Ma, which was bolstered by later paleoelevation estimates documented by δ 2 H analysis of volcanic glass, δ 18 O and Δ47 clumped isotope analysis of soil carbonates, pollen assemblages, and leaf wax lipids which collectively produced a maximum paleoelevation of 2.8 km between 18 and 7 Ma (Kar et al., 2016). Regional δ 2 H stable isotopic analysis of volcanic glass yield a maximum elevation of ∼2 km between 25 and 19 Ma (Sundell, Saylor, Lapen, et al., 2019). Finally, Monte Carlo based river incision modeling of the Rio Ocoña in southern Peru provides a maximum paleoelevation estimate of ∼3 km at 20 Ma (Jeffery et al., 2013).
The transition to hyperarid conditions occurred long after the establishment of moderate to high elevations along the western CAP margin and the attainment of near modern crustal thickness (Figure 4e). Miocene palustrine and lacustrine carbonates and fossil soils place the onset of extreme aridity in the Atacama Desert at 12-10 Ma (Rech et al., 2019), after the attainment of extreme elevation in the Altiplano (Figure 4h). Similar evidence is seen slightly earlier in hinterland basin stratigraphy of the Peruvian Western Cordillera where a nearly bimodal switch from dominantly fluvial to gypsum-bearing lacustrine-evaporitic deposition takes place at ∼18 Ma (Figure 4i; Carlotto, 2013). This is in contrast to paleoclimate models which predict that the opposite took place: that the CAP became wetter, not drier, due to surface uplift (Chandan & Peltier, 2017;Ehlers & Poulsen, 2009), and is inconsistent with earlier climate modeling work (e.g., Garreaud et al., 2010) as well as recently documented fossil records of pollen, leaves, and wood (Martínez et al., 2020).
Despite the protracted history of crustal thickening since ∼80 Ma and maximum zircon trace element ratios reaching near modern values by ∼30 Ma (Figures 4e and 4g), the northern CAP Altiplano did not reach near modern elevations until the middle to late Miocene (Figure 4h), implying that these processes were in part decoupled in the late Cenozoic. Alternatively, if our preferred hypothesis is incorrect and results indeed suggest that the crust thickened in two separate stages, with a period of crustal thinning in between, results might indicate several periods of uplift and subsidence. In either scenario, the apparent inconsistences in paleoelevation we observe may be related to the complex effect of global climate change, at least in Miocene time, and warrants further research on stable isotopes used as paleoelevation proxies (e.g., Carrapa et al., 2019).

Data Availability Statement
The geochemical data used in the study are available from the authors upon request and at Zenodo via https://doi. org/10.5281/zenodo.5800271.