Slow Long-Term Exhumation of the West Central Andean Plate Boundary, Chile.

We present a regional analysis of new low-temperature thermochronometer ages from the Central Andean fore arc to provide insights into the exhumation history of the western Andean margin. To derive exhumation rates over 10 million-year timescales, 38 new apatite and zircon (U-Th)/He ages were analyzed along six ~500-km long near-equal-elevation, coast parallel, transects in the Coastal Cordillera (CC) and higher-elevation Precordillera (PC) of the northern Chilean Andes between latitudes 18.5°S and 22.5°S. These transects were augmented with age-elevation profiles where possible. Results are synthesized with previously published thermochronometric data, corroborating a previously observed trenchward increase in cooling ages in Peru and northern Chile. One-dimensional thermal-kinematic modeling of all available multichronometer equal-elevation samples reveals mean exhumation rates of <0.2 km/Myr since ~50 Ma in the PC and ~100 Ma in the CC. Regression of pseudovertical age-elevation transects in the CC yields comparable rates of ~0.05 to ~0.12 km/Myr between ~40 and 80 Ma. Differences between the long-term mean 1-D rates and shorter-term age-elevation-derived rates indicate low variability in the exhumation history. Modeling results suggest similar background exhumation rates in the CC and PC; younger ages in the PC are largely a function of increased heat flow and consequently an elevated geothermal gradient near the arc. Slow exhumation rates are suggestive of semiarid conditions across the region since at least the Eocene and deformation and development of the Andean fore arc around this time.

Introduction

In mountain belts, low‐temperature thermochronometry can constrain the cooling histories of rocks that were exhumed to the surface over timescales of millions to tens of million years (e.g., Reiners & Brandon, 2006). Because of the complex surface and subsurface feedback mechanisms involved in mountain evolution (e.g., Beaumont et al., 2001; Willett, 1999; Zeitler et al., 2001), accurate and complete rock cooling records are critical for understanding climate and tectonic factors contributing to orogen deformation and topography. After the Himalayas, the Central Andean plateau is the second highest in the world and separates some of the driest and wettest locations on Earth. In the wake of these extremes and at the margin of an actively subducting plate, one could expect to see ongoing mountain‐producing deformation at the leading edge of the range. Instead, indicators of exhumation and erosion such as mineral cooling ages have documented processes an order of magnitude older in parts of the Andes than what we see in other convergent orogens (e.g., Barnes & Ehlers, 2009; Masek et al., 1994). Understanding the mechanisms behind such long‐lived slow mountain building, and decoupling the relative contributions made by tectonic and climatic processes, has long been a challenge in this field of research.

In the fore arc of the Central Andes, ~15°–25°S, the significance of thermochronometric data has been subject to various interpretations; mineral cooling ages have been used to track the timing, depth, and spatial shifts of pluton emplacement (Andriessen & Reutter, 1994; McInnes et al., 1999), rock uplift induced by subduction of oceanic ridges (Juez‐Larré et al., 2010; Wipf et al., 2008), regional erosion and exhumation (Maksaev & Zentilli, 1999; Schildgen et al., 2009), and changes in the thermal structure (Juez‐Larré et al., 2010). While rock cooling rates are intrinsically tied to the rate of erosion, as well as rock and surface uplift, the comparison of these values is highly sensitive to changes in geothermal parameters and whether an erosion‐uplift system is in an exhumational equilibrium state (England & Molnar, 1990; Ring et al., 1999; Willett & Brandon, 2002). While directly equating cooling rates with erosion and uplift may be specious, with thoughtful sampling and consideration of thermal properties in a kinematic model, mineral cooling ages can track, to a first order, the rates of rock exhumation and thus the erosion history of a region. Previous thermochronometric studies quantifying exhumation in the Central Andean fore arc have often focused on canyon incision rates as proxies for timing of surface uplift (Cooper et al., 2016; Hoke et al., 2007; Schildgen et al., 2007; Schildgen et al., 2009). Rapid and sustained river incisions have been suggested as mechanisms for inciting tectonic deformation. For example, the “tectonic aneurism” model (e.g., Zeitler et al., 2001) predicts that removal of cold crust by river incision could enhance thermal gradients and induce temperature‐dependent crustal weakening leading to focused deformation. Conversely, in arid environments such as the Atacama desert, located in the Central Andean fore arc between latitudes ~18°S and 30°S, river incision may be very local and too low to induce such feedbacks. Tectonic forces may instead play a significant role in controlling exhumation of uplifted terrains, particularly at plate corners such as the Central Andean syntaxis where the subducting plate assumes geometrically induced rigidity (Bendick & Ehlers, 2014). To decipher the origins of deformation, regional background exhumation rates on longer timescales (e.g., tens of million years) are important in constraining the full history of a region's climatic‐erosional regime.

To better constrain the impact of long‐term climatic and tectonic conditions that aided in forming the Central Andes, we calculate mean regional background exhumation rates in the Andean fore arc. To do so, we synthesize existing low‐temperature thermochronometric cooling ages with newly collected apatite and zircon (U‐Th)/He ages from the northernmost part of the Chilean coastal margin. After assessing the possibility of crystallization cooling effects and the geothermal context of the region, we interpret our ages as resulting from exhumation. Using a 1‐D transient Monte Carlo thermal‐kinematic modeling and a multichronometer approach, we synthesize a range of plausible mean exhumation rates across an ~500‐km stretch of the fore arc, rather than focusing on the cooling history of a localized area. We compare these results to new and existing age‐elevation vertical transects that provide evidence for periodic exhumation rates on shorter timescales. We address the timing and rates at different elevations and distances along the subduction margin to allow for a more holistic picture of possible spatial variations or gradients in erosion. Finally, we present our exhumation history in the context of what is known from previous studies about the climate and tectonic history of the margin.

Field Setting

Regional Geology and Tectonic Setting

The Andes have formed in a continental‐oceanic plate convergent setting. Subduction of the Nazca (and previously Farallon) Plate under the South American Plate has resulted in a mountain belt characterized by the largest trench‐to‐peak relief on the planet, frequent seismic activity (M  > 8.5), and a long‐standing history of volcanism (e.g., Armijo et al., 2015). Notably, three distinctly oriented sections of the range (Figure 1) may be the result of differential shortening recording periods of Andean deformation and uplift (Allmendinger et al., 1997; Isacks, 1988; Roperch et al., 2006), though this uniquely bent geometry may itself lead to focused deformation in the upper plate, particularly at plate corners (e.g., Bendick & Ehlers, 2014; Jordan et al., 1983). Other tectonic controls have been invoked to explain variation in Andean morphology (i.e., angle of convergence, convergence velocity, age of oceanic crust, and trench orientation), though no clear consensus exists on their relative importance (e.g., Barnes & Ehlers, 2009; Maloney et al., 2013; Oncken et al., 2006). Broad regional uplift has been tied to removal of lithospeheric material through either gradual subduction erosion or rapid lower lithosphere removal of a high‐density crustal root (e.g., delamination; Garzione et al., 2008, 2017). Subduction of aseismic ridges (e.g., Nazca and Juan‐Fernandez ridges) may also “indent” overlying topography and focus uplift (Gutscher et al., 1999; Jordan et al., 1983; Saillard et al., 2011; Wipf, 2006; Yáñez et al., 2002), though the significance of subducted bathymetric anomalies in upper plate dynamics has also been questioned (e.g., Skinner & Clayton, 2013).

Figure 1

Tectonic overview map of Andean margin from ~0 to 40° south. Topography is illustrated with the USGS GTOPO30 30‐arcsecond global DEM with labeled major segments based on strike and morphological deflections. Volcanic distribution in the Andes is closely correlated to areas of flat‐slab subduction, also visible as deviations in the 100‐km slab contour lines bending in the direction of subduction (Slab 1.0 model; Hayes et al., 2012). The field area of this study, outlined in red, sits at the northern end of Chile directly under the Arica bend within a zone of relatively consistent slab age, volcanic distribution, subduction angle, distance from trench, convergence orientation, and morphology. Location of trench/plate boundaries and oceanic crustal ages are derived from the GPLATES global plate model (Müller et al., 2008; Seton et al., 2012). Locations of aseismic ridges, labeled and outlined in gray, are derived from the GEBCO global bathymetry model (Weatherall et al., 2015). Modern relative plate motion direction and magnitude is from the MORVEL plate model (DeMets et al., 2010). Volcanic centers are from Siebert and Simkin (2002). Political boundaries are outlined in dashed black lines.

Although subduction along the Andean margin started as early as the Jurassic (~200 Ma) (James, 1971; Rutland, 1971), most of the uplift that has produced the Altiplano‐Puna Plateau (~3,800‐m mean elevation) has occurred only since the Paleocene and Eocene (~60–40 Ma; Barnes & Ehlers, 2009; McQuarrie et al., 2005; Sempere et al., 1990). The timing and mechanism of this uplift has been largely debated (e.g., Barnes & Ehlers, 2009; Ehlers & Poulsen, 2009; Garzione et al., 2008, 2017), and the timing and rate of exhumation along the margin, intrinsically tied to uplift, can thus be difficult to decouple from its tectonic controls.

The Central Andean fore arc has not generally been well integrated into Andean deformation models (e.g., Juez‐Larré et al., 2010) despite the region's potential to record the oldest history of subduction‐related uplift in the Andes. Although subduction erosion removes much of the previous history of deformation in the fore arc, in conjunction with subcrustal accretion at the leading edge of the continental plate, subduction has driven deformation in the fore arc and resulted in the formation of three distinct physiographic N‐S striking domains, which lead into the present‐day arc, the Western Cordillera (WC; e.g., Hartley et al., 2000; Kukowski & Onken, 2006). These domains are, from west to east, the Coastal Cordillera (CC), largely made of sedimentary and volcanic deposits ranging from the Paleozoic through the Paleogene with Jurassic and Cretaceous andesites, diorites, and granodiorites formed in an early arc formed at the birth of the Andes; the Central Depression (CD), a basin filled with late Eocene to early Pliocene sediments and volcanic deposits; and the Precordillera (PC), containing plutonic and volcanic remnants of a Cretaceous‐Paleogene arc rising into the WC with Mesozoic sedimentary and volcanic sequences underlying more recent Cenozoic cover (Figure 2b; Hartley & Chong, 2002; Scheuber & Reutter, 1992). The CC is crosscut by N‐S striking normal faults, E‐W striking reverse faults, and some NW‐SE striking strike‐slip faults (Figure 2a; Allmendinger et al., 2005; González et al., 2003). Separating the Central Depression from the Coastal Cordillera on the west and the Precordillera on the east are the Atacama Fault Zone and the West Andean Thrust, respectively, accommodating some of the differential uplift between the physiographic domains (Armijo et al., 2015; González et al., 2003; Muñoz & Charrier, 1996; Muñoz & Sepulveda, 1992; Victor et al., 2004). Most of northern Chilean (e.g., Isacks, 1988; Jordan et al., 2010; Wörner et al., 2002) and Peruvian (2.5° northward; e.g., Schildgen et al., 2007) Western Cordilleran uplift, however, is thought to have occurred by way of a crustal‐scale monocline of which the Precordilleran slope constitutes the limb.

Figure 2

(a) Digital elevation model (SRTM 3‐arcsecond) hillshade of field area with locations of apatite (U‐Th‐Sm)/He (AHe), apatite fission track (AFT), and zircon (U‐Th)/He (ZHe) samples. Previously published age locations in northern Chile, depicted with hollow symbols, are from Andriessen and Reutter (1994), Maksaev and Zentilli (1999), and Juez‐Larré et al. (2010). Profile lines used in Figure 3 across the northern and southern sections of the field area are depicted in red. Individual faults are mapped in black, and the Atacama Fault Zone and approximate location of the Western Andean Thrust (WAT) is delineated in white dashed lines. A solid transparent white line delineates the political boundary of Chile. (b) Simplified geologic map of field area showing location of plutonic units and sedimentary/volcanic units. Deposits younger than the Paleozoic are grouped by geologic period. Swaths used for latitudinal topographic analysis plotted in Figures 5 and 6 are outlined in dashed red lines. CC = Coastal Cordillera; CD = Central Depression; PC = Precordillera; WC = Western Cordillera. Bottom right inset map depicts extent of both (a) and (b). Fault locations and geologic units were mapped by SERNAGEOMIN (2003).

Figure 3

Thermochronometer age and topographic profiles A‐A' and B‐B′ depicted in Figure 2. Mean (green‐solid), maximum (red‐dotted), and minimum (blue‐dotted) elevation is calculated in swaths extending 100 km in both directions perpendicular to each profile line. Elevation and location of thermochronometric sample locations (projected onto the nearest profile line) from this study are plotted as white circles over the topographic profiles. Age data are plotted overhand. All data are plotted against the distance from the present‐day subduction trench, measured at the center line of each swath. Error bars represent 2σ uncertainties. Maximum and mean elevations show clear morphological variation between the northern and southern sections of the Coastal Cordillera.

Figure 5

North‐south topographic swath profile for the Coastal Cordillera outlined in Figure 2b with iso‐elevation cooling ages for sea level and ~1,000‐m bins plotted overhand. Mean (green‐solid), maximum (red‐dotted), and minimum (blue‐dotted) elevation is calculated in the longitudinal direction and plotted against latitude. Multichronometer samples used in Pecube modeling are outlined in dashed black lines. Dashed blue and red lines represent the inverse‐squared‐error‐weighted means across the sea level elevation bin for apatite and zircon (U‐Th)/He ages, respectively. Error bars represent 2σ uncertainties.

Figure 6

North‐south topographic swath profile for the Precordillera outlined in Figure 2b with iso‐elevation cooling ages for ~1,500, 2,000, 2,500, and 3,000‐m bins plotted overhand. Mean (green‐solid), maximum (red‐dotted), and minimum (blue‐dotted) elevation is calculated in the longitudinal direction and plotted against latitude. Multichronometer samples used in Pecube modeling are outlined in dashed black lines. Dashed blue and red lines represent the inverse‐squared‐error‐weighted means across the ~2,000‐m elevation bin for apatite and zircon (U‐Th)/He ages, respectively. Error bars represent 2σ uncertainties.

Regional Climate and Geomorphic Setting

Climate in the northern Chilean Central Andean fore arc is characterized by long‐lived and spatially homogenous aridity attributed to the location of the Atacama region at the eastern edge of the subtropical Pacific. At these latitudes, strong and stable atmospheric subsidence maintains the SE Pacific anticyclone and promotes upwelling of cold waters, in turn, promoting further subsidence and drying of the air at the coast (Garreaud, 2009; Garreaud et al., 2010). However, estimates of onset, origin, and intensity of Atacama aridity have been largely debated (e.g., Dunai et al., 2005; Hartley & Chong, 2002). Long‐term estimates based on evaporitic deposits suggest at the very least intermittent semiarid conditions since as early as 200 Ma (Clarke, 2005) though a decrease to arid (<200 mm/yr) and subsequently hyperarid (<20–50 mm/yr) conditions is generally agreed to have occurred since the Miocene (e.g., 19–13 Ma; Jordan et al., 2014; Rech et al., 2006). Repeated onset of hyperaridity since 14 Ma may be further responsible for ambiguous records and interpretations (Jordan et al., 2014). These transitions have been tied to (1) a coastal temperature inversion from the cold Humboldt current formed concurrent with the opening of the Drake Passage (Gregory‐Wodzicki, 2000; Lawver & Gahagan, 1998; Wörner et al., 2002), (2) an Andean rain shadow blocking off Atlantic moisture as a result of post‐Miocene uplift (e.g., Houston & Hartley, 2003; Rech et al., 2006, 2010), and (3) additional shut off from Pacific moisture blocked by the Coastal Cordillera(Rech et al., 2010). Present‐day precipitation rates are recorded as low as <1 mm/yr in the coastal cities of northern Chile (e.g., Schulz et al., 2012). The implication of arid climate in Andean mountain building mechanisms is critical as it is associated with the absence of significant and widespread erosion, and removal of surface material can compensate for tectonically induced rock uplift, weaken the thermal integrity of the crust, or initiate an isostatic response (Gregory‐Wodzicki, 2000; Jordan et al., 2010; Schildgen et al., 2009; Zeitler et al., 2001). Additionally, the lack of sediment supply to lubricate subduction in the trench has been suggested to act as a roughening mechanism, enhancing the transfer of compressional forces from the subducting plate to the overlying crust (Lamb & Davis, 2003).

Despite climatic homogeneity, geomorphic variation in the area is starkly expressed at ~19.5°S by a threshold between exorheic basins to the north and endorheic basins to the south, draining water either directly to the Pacific or through the Central Depression and out through the Rio Loa (e.g., Schlunegger et al., 2017). As a result, in the north, topography transitions smoothly from the CC at ~1,000 m above sea level (m.a.s.l.) into the Central Basin and up toward the slopes of the PC, subsequently cut through by deep canyons, which are much less frequent toward the south (Figures 3 and 4). Coudurier‐Curveur et al. (2015) attributed this change to latitudinal precipitation gradients over the WC and Altiplano, which carry monsoonal Atlantic humidity and ultimately drain toward the Pacific. The endorheic/exorheic threshold may thus coincide with a transition from aridity/semiaridity to hyperaridity (García et al., 2011) rather than one based on tectonic controls (e.g., subducted ridges and effects of enhanced subducting plate rigidity in the Arica bend).

Figure 4

Map of bedrock thermochronometry results in northern Chile. Black contour lines show elevation values in m.a.s.l. Samples from this study are indicated as white circles. Sample elevations in meters above sea level are shown in black italics. Apatite and zircon (U‐Th)/He ages with 2σ uncertainties are shown in blue and red, respectively. Previously published thermochronometer sample locations in the area are shown as hollow blue squares, blue circles, and red diamonds for apatite (U‐Th)/He, apatite fission track, and zircon (U‐Th)/He data, respectively. Samples used for vertical transects are outlined with red dashed lines. The vertical transect at Tocopilla from Juez‐Larré et al. (2010) is outlined in pink dashed lines.

Bedrock Thermochronometry

Apatite and Zircon (U‐Th)/He Dating

The (U‐Th)/He method is used to track the time since a mineral cooled below closure temperatures of ~55–90 °C for apatite and ~160–200 °C for zircon (e.g., Ehlers et al., 2005). The effective closure temperature of any sample depends on the cooling rate, grain size, mineralogy, and radiation damage (e.g., Farley, 2000; Flowers et al., 2009; Reiners et al., 2004, 2005). For this study, we collected 38 bedrock samples between ~18.5°S and 22.5°S along the fore arc of northern Chile, directly south of the Arica bend where a lack of thermochronometric data exist in previous literature (Figure 2a). The data tables are presented in the supporting information and a summary of all samples is in Table 1. Sampling in an area of long‐lasting arid climate, considerably uniform N‐S topography and subduction parameters, and a lack of accreted terrains (e.g., accreted microplates or island arcs) allows for near‐ideal conditions to explore potential exhumation gradients as a factor of spatial or temporal variations in tectonic factors (such as plate geometry or subduction angle) or climate.

Table 1

Summary Results of All Thermochronometric Data

Sample	Elevation (m)	Longitude (°)	Latitude (°)	Lithology	Unit namea	AHe age (Ma)	Error (2σ)	n	ZHe age (Ma)	Error (2σ)	n
Coastal Cordillera
15CL001	3	−70.20551	−19.58834	Granodiorite	JKg	49.24	10.95	5	‐	‐	‐
15CL029b	5	−70.13727	−20.18909	Granodiorite	JKg	57.61	5.13	5	137.87	19.27	3
15CL017	10	−70.32252	−18.52954	Andesite Porphory	J3i	43.02	4.16	3	57.39	94.83	3
15CL013b	89	−70.1888	−19.1555	Granodiorite	Jsg	24.53	3.92	5	115.58	7.32	3
15CL003	157	−70.19412	−19.55756	Granodiorite	JKg	12.62	5.16	5	‐	‐	‐
15CL004	265	−70.2013	−19.59009	Granodiorite	JKg	62.1	11.36	4	‐	‐	‐
15CL005	527	−70.18188	−19.56605	Granodiorite	JKg	55.29	30.2	5	‐	‐	‐
15CL030	737	−70.09938	−20.20582	Gabbro	JKg	4.87	1.53	5	‐	‐	‐
15CL048	954	−69.68369	−21.10682	Granite	CPg	109.51	29.29	5	‐	‐	‐
15CL052	979	−70.09074	−22.11668	Granodiorite	Jsg	90.39	5.55	4	‐	‐	‐
15CL007	1045	−70.08687	−19.55242	Granodiorite	JKg	61.58	5.84	5	‐	‐	‐
15CL051	1287	−69.71252	−22.28071	Granodiorite	Jig	57.91	14.28	5	‐	‐	‐
Precordillera
15CL037	1592	−69.44785	−19.86981	Metasandstone	Js1m	‐	‐	‐	‐	‐	‐
15CL038	1706	−69.38956	−19.86662	Granodiorite	KTg	46.73	11.51	5	‐	‐	‐
15CL068	1724	−69.17405	−20.96502	Granodiorite	Eg	32.87	30.08	5	‐	‐	‐
15CL039	1731	−69.37644	−19.86406	Granodiorite	KTg	46.29	8.292	4	‐	‐	‐
15CL060	1749	−69.39873	−20.11385	Diorite‐Gabbro	KTg	49.64	4.35	5	‐	‐	‐
15CL040	1758	−69.36977	−19.86345	Diorite‐Andesite	KTg	36.11	4.13	5	‐	‐	‐
15CL043	1924	−69.39466	−19.87144	Granodiorite	KTg	38.96	33.87	5	‐	‐	‐
15CL062	1955	−69.37035	−20.08307	Granite	KTg	54.09	40.98	4	60.52	6.96	3
15CL072b	1955	−69.18936	−22.28974	Granodiorite	Pag	41.87	3.21	5	60.33	3.75	3
15CL010	2006	−69.494	−19.34103	Metasandstone	J1 m	‐	‐	‐	62.19	3.79	3
15CL022b	2198	−69.79704	−18.46162	Granodiorite	KTg	28.36	19.02	4	51.09	3.06	3
15CL044	2224	−69.37527	−19.88301	Gabbro	KTg	32.6	14.78	3	‐	‐	‐
15CL066	2239	−69.06452	−20.92909	Metasandstone	JK1c	8.31	3.16	4	‐	‐	‐
15CL032	2271	−69.49072	−18.95063	Granodiorite	KTg	90.71	56.92	5	‐	‐	‐
15CL064	2273	−69.0533	−20.9265	Metasandstone	JK1c	12.89	9.72	5	‐	‐	‐
15CL023	2499	−69.76244	−18.44944	Granodiorite	KTg	34.43	26.83	5	‐	‐	‐
15CL073	2537	−69.0998	−22.30856	Granodiorite	Pag	36.67	3.3	5	‐	‐	‐
15CL024	2648	−69.75387	−18.4406	Granodiorite	KTg	23.83	8.45	5	‐	‐	‐
15CL076	2857	−68.96645	−22.32083	Granodiorite	Eg	16.89	7.16	5	‐	‐	‐
15CL056	2996	−69.18234	−20.08786	Granite	Og	11.96	2.52	5	‐	‐	‐
15CL074	3017	−68.97912	−22.31384	Granodiorite	Eg	22.47	5.62	5	‐	‐	‐

Note. n = number of grains measured.

As reported in SERNAGEOMIN (2003): CPg = Carboniferous‐Permian granite; Jig = Lower Jurassic granodiorite; Jsg = Middle to Upper Jurassic granodiorite; J3i = Jurassic volcanics; J1m = Jurassic sediments; Js1m = Middle to Upper Jurassic sediments; JK1c = Upper Jurassic to Lower Cretaceous sediments; JKg = Jurassic‐Cretaceous granodiorite; KTg = Upper Cretaceous to Lower Tertiary granodiorite; Pag = Paleocene granodiorite; Eg = Eocene granodiorite; Og = Oligocene granite.

Individual samples used in thermal‐kinematic modeling.

Methods: Sample Collection and Analysis

Samples were collected in deliberate N‐S near‐iso‐elevation (hereon referred to as iso‐elevation) profiles at sea level and 1,000, 1,500, 2,000, 2,500, and 3,000 m.a.s.l. (±250) in the CC and PC (Figure 4, see also supporting information). Sediment cover in the CD did not allow for bedrock (i.e., nondetrital/basin) samples located in between the two ranges. For age‐elevation relationships, three areas were additionally sampled as vertical transects in >200‐m intervals with minimal possible horizontal distances between samples. Due to limited exposure of suitable lithologies for exhumation studies, additional age‐elevation profiles, or profiles collected over large vertical distance, were not possible. Given this, the sampling strategy used in this study was designed to identify any regional‐scale variations in cooling ages with iso‐elevation profiles. Samples were collected from both endorheic and exorheic sections of the field area.

We follow standard mineral separation protocols and analytical procedures as those described in Stübner et al. (2016; detailed in the supporting information). Ages were calculated following the methods of Meesters and Dunai (2005), and error was propagated through our analysis to accurately estimate uncertainty in resulting ages. All grains were measured using spatially calibrated photos taken across and down the c axis to calculate alpha ejection correction factors (Ft; e.g., Farley et al., 1996; Farley & Stockli, 2002; Hourigan et al., 2005; Ketcham et al., 2011) using the geometric models of Ketcham et al. (2011). Throughout, we take a conservative error estimate of 5% on all Ft values derived from the maximum difference on select problematic grains between our Ft calculation and one modeled with a Matlab script, using a Monte Carlo alpha ejection simulation on 3‐D grain models built from our multiangle photos (Glotzbach et al., 2017).

Mean sample ages and associated error were derived following a modified version of the approach of Cooper et al. (2011) and Adams et al. (2015). Final error (reported at the 2σ level) on each sample was calculated after removing outliers based on qualitative exclusion (cooling ages > mapped crystallization age; SERNAGEOMIN, 2003, and references therein) and the Hampel identifier method (“median absolute deviation from the median” threshold of 4; Hampel, 1974; Pearson, 2011). Mean sample values were derived with a relative error weighted mean of the remaining ages (n > 3 for all but one sample) by propagating individual grain errors through the weighted mean equation (Bonamente, 2013). While it is common practice to weigh data by the inverse of their squared absolute error (the inverse of the variance), because analytical error correlates strongly with age, using a relative error weighting (%) prevents a strong and misleading bias toward younger samples.

Frequently, with (U‐Th)/He data, propagated analytical error is insufficient in accounting for the variation in a sample's single grain ages, which we assume part of one population. To address excess scatter from unknown human error or geologic processes, we calculate a mean square weighted deviation (MSWD) for each sample, estimating a standard deviation of the MSWD (σ MSWD) after the method of Wendt and Carl (1991), where σ MSWD = [2/(n − 1)]0.5. In cases where the MSWD >1.0 ± 2σ MSWD, we multiply the error by MSWD0.5 to expand it to an approximate 95% confidence interval matching the scatter in our data (Adams et al., 2015; Cooper et al., 2011; Ludwig, 2012; Ludwig & Titterington, 1994; Wendt & Carl, 1991). Previously published zircon and apatite (U‐Th)/He mean ages (Juez‐Larré et al., 2010) used in our analysis were recalculated from individual grain ages according to the same methods described and using individual grain Ft values of the original study.

Thermochronometry Results

Cooling Ages in the CC

The 12 CC samples (12 AHe and 3 ZHe ages) collected at low elevations (0–500 m.a.s.l.), hereafter referred to as “sea level” samples, show a large range in AHe ages of 12.6 to 62.1 Ma (Figure 5; see also supporting information). Most of these ages agree within error with the AHe sea level iso‐elevation profile of Juez‐Larré et al. (2010; 0–500‐m samples range from 37.6 to 83.8 Ma), which extends along the coast from latitudes ~20.9°S to ~27.0°S, overlapping considerably with our study area. Two of the three ZHe ages in our CC sea level transect had ages of 115.6 ± 7.3 and 137.9 ± 19.3 Ma, compared to a sea level zircon age of 170.2 ± 52.0 Ma from Juez‐Larré et al. (2010). Ages in the CC generally increase with elevation, though with significant local inconsistencies: the three ~1,000‐m iso‐elevation AHe ages in the CC yield a weighted mean of 77.3 Ma (Figure 5). All AHe samples are, however, younger than the AFT ages reported by Juez‐Larré et al. (2010; 118–145 Ma).

Cooling Ages in the PC

In the Precordillera, 21 samples (20 AHe and 4 ZHe ages) document a range of AHe ages from 8.3 to 90.7 Ma (Figure 6). Error‐weighted‐mean ages in the PC generally decrease with increasing elevation bins (48.5, 27.6, 33.1, and 14.0 Ma from the 1,500‐ to 3,000‐m bins), though individual sample cooling ages do not always locally follow this trend. In general, older ages also yield larger single grain scatter and resulting error. Furthermore, as the west‐to‐east increase in elevation is consistent throughout the fore arc, we note that ages decrease with easting longitude and with elevation, and both trends warrant further discussion (Figure 3). Four zircon ages from the PC were measured from the ~2,000‐m elevation bin and range from 51.1 to 62.2 Ma, more consistent than the AHe age range from this bin of 8.3 to 54.1 Ma. Three samples (15CL64, 15CL066, and 15CL056), two of which are sedimentary, yielded the youngest AHe ages (~10 Ma). AFT results (47.4–58.2 Ma) from Maksaev and Zentilli (1999), sampled at ~2,000 m.a.s.l. in the southernmost section of our field area, fall between our AHe and ZHe results. AHe results of McInnes et al. (1999) are omitted in our analyses due to their location in the Chuquicamata porphyry copper deposit, which has documented shallow (<4‐km) emplacement and is cross cut by the prominent West Fissure Fault, which may have had as much as 600‐m vertical and 25–40‐km sinistral displacement (Maksaev & Zentilli, 1999; McInnes et al., 1999). Samples follow a steady mean topography N‐S until a dip in elevation at ~22.5°S (Figure 6), and no clear latitudinal trend is discernable.

Age‐Elevation Relationships

To complement the iso‐elevation profiles, we analyzed two areas in the PC and one in the CC yielding sample ages defining vertical transects (VT, Figure 7) in addition to the Tocopilla VT (using recalculated ages for consistency) reported in Juez‐Larré et al. (2010). Assuming a constant AHe closure isotherm over the horizontal‐spatial and temporal extent of a given VT, we interpreted a simple age‐elevation relationship (AER) by fitting a straight line with elevation as the independent and age as the dependent variable, weighted according to the error (e.g., Glotzbach et al., 2011). While such an approach allows for some estimate of exhumation rates, this method is not ideal when applied to our samples because it does not quantify transient thermal conditions and we do not have a consistent suite of samples to robustly determine rates from them. Later, we improve upon this approach for calculating exhumation rates with a transient thermal model (see section 4).

Figure 7

Vertical transects (VTs) of apatite (U‐Th)/He data for four locations in the field area. Transect locations are shown in Figure 4. Age‐elevation relationships are plotted for samples from the Coastal Cordillera (Pisagua and Tocopilla) and the Precordillera (Putre and Calama). Age data and error for the Tocopilla VT are taken from Juez‐Larré et al. (2010) and recalculated for consistency with the analytical procedures used in this study. Mean (black‐solid) and maximum/minimum (black‐dotted) lines show linear regressed slopes weighted by their error for each data set, representing exhumation rate (Ø) estimates for that period (slope units are in km/Myr). Sample 15CL003 is excluded from the Pisagua VT as an outlier (see age‐elevation relationships, section 3.2.3). R 2 values are shown for each regression, for two subsections in the Tocopilla VT (red dashed lines) and for an age versus longitude correlation for the Calama VT (inset plot). Error bars represent 2σ uncertainties.

Linear regression of the ~2‐km Tocopilla VT (Figures 4 and 7c) yields a mean exhumation rate of 0.052 km/Myr (0.046 to 0.060 km/Myr at 2σ level) with an R 2 of 0.81 and an age range 37.6 to 83.6 Ma over a ~3‐km horizontal extent. Though Juez‐Larré et al. (2010) interpreted and modeled the section in two parts due to an apparent break in trend at ~50–60 Ma, we find that the lack of sufficient correlation of data in the subsections (R 2 = 0.15 for ~37–52 Ma; R 2 = 0.29 for ~56–84 Ma) suggests that a quantitative AER interpretation of anything besides the mean exhumation rate for the entire period is not warranted (Figure 7c).

The ~1,000 m Pisagua VT (Figures 4 and 7a), ~280 km north of Tocopilla, spans a greater horizontal extent of ~20 km and regresses to a mean exhumation rate of 0.12 km/Myr with an R 2 of 0.74 after the exclusion of a young outlier age, which we attribute to incision of the adjacent Tana‐Tiliviche canyon or the sample location along the Pisagua Fault (Coudurier‐Curveur et al., 2015). Ages ranging from 49.2 and 62.1 Ma are consistent with the Tocopilla VT, albeit shifted ~10–15 Myr toward older ages at the same elevation range. The minimum rate of 0.054 km/Myr (2σ) falls within range of the Tocopilla VT, whereas the maximum at this confidence interval yields a reversal in the regression trend.

In the Precordillera, both the Putre and Calama vertical transects (Figures 4, 7b, and 7d) yield AERs that are predominantly inverted, following the decrease in ages from west to east, and cannot be interpreted under consistent geothermal gradient conditions between these samples. In part, due to lower local relief, sampling at very close horizontal distances while maintaining meaningful elevation differences is not possible in the Calama area where the horizontal VT extent is ~23 km, compared to ~5 km in Putre. Schildgen et al. (2007, 2009) accounted for long‐wavelength lateral changes in the closure temperature isotherms in southern Peru by using depth below a paleosurface instead of elevation, but estimating an accurate paleosurface for samples outside of a clearly incised canyon is unrealistic in our case. With the available thermochronometer data, and the lack of a paleosurface, we are unable to conduct a meaningful 2‐D or 3‐D thermal model of the region to correct for lateral variations in the closure isotherm depth and therefore resort to calculating exhumation rates based on a 1‐D modeling approach (see below). Although PC AER results do not produce meaningful exhumation rates, a correlation between age and longitude is evident. The Calama transect, for example, demonstrates a cleaner age‐longitude correlation (R 2 = 0.97), than the inverted trend in the age‐elevation result (R 2 = 0.88; Figure 7d inset).

Thermal‐Kinematic Modeling

Methods: Model Setup

To evaluate the cooling histories of our multichronometer samples and relate them to exhumation rates, we employ a thermal‐kinematic model, which solves the transient advection‐diffusion heat equation. We follow the methods originally presented by Thiede and Ehlers (2013) in the Himalaya and subsequently applied by Adams et al. (2015) and Schultz et al. (2017) using a modified 1‐D version of the Pecube finite element scheme (Braun, 2003; Braun et al., 2012) with implicit time steps and a smart Monte Carlo search approach that identifies the range of exhumation histories that reproduce observed ages (Figures 9 and 10). Model free parameters are the exhumation rates, which are allowed to vary in discretely defined time steps. Between west (CC) and east (PC), we use spatially variable geothermal parameters but we do not vary them within a simulation for a given location (Table 2).

Figure 9

Pecube 1‐D thermal‐kinematic model results for exhumation in the Coastal Cordillera. Apparent exhumation histories since 200 Ma are depicted for three sea level samples and the inverse‐squared‐error‐weighted mean representative proxy run. Mean rate is represented by a solid black line. Cooling age ranges with their uncertainties that were used in modeling are shown in the sample run plots as red, cyan, and blue areas for ZHe, AFT, and AHe, respectively. 3 ranges on acceptable exhumation histories are depicted with dashed black lines and diagonal shading. The lower panel shows a synthesis of all mean apparent exhumation results and displays consistent histories across this section of the Coastal Cordillera.

Figure 10

Pecube 1‐D thermal‐kinematic model results for exhumation in the Precordillera. Apparent exhumation histories since 75 Ma are depicted for two of the 2,000‐m bin samples and the inverse‐squared‐error‐weighted mean representative proxy run. Mean rate is represented by a solid black line. Cooling age ranges with their uncertainties that were used in modeling are shown in the sample run plots as red, cyan, and blue areas for ZHe, AFT, and AHe, respectively. 3 ranges on acceptable exhumation histories are depicted with dashed black lines and diagonal shading. The lower panel shows a synthesis of all mean apparent exhumation results and displays consistent histories across the Precordillera since ~40–50 Ma.

Table 2

Parameters Used in Thermal‐Kinematic Modeling

Parameter	Value (Coastal Cordillera)	Value (Precordillera)	Source
Model thickness (km)	45	75	Springer (1999)
Thermal conductivity (W · m · K)	2.8	2.8	Schildgen et al. (2009) and Springer and Förster (1998)
Specific heat capacity (J · kg · K)	800	800	Schildgen et al. (2009)
Crustal density (kg/m3)	2750	2750	Schildgen et al. (2009) and Springer (1999)
Volumetric heat production of continental crust (uW/m3)	0.8	0.8	Schildgen et al. (2009) and Springer (1999)
Volumetric heat production of oceanic crust (uW/m3)	0.25	0.25	Springer (1999)
Temperature at base of model (°C)	250	500	Springer (1999)
Temperature at surface (°C)	15a	3b	aJuez‐Larré et al. (2010); bbased on 6.5 °C/km lapse rate, Schildgen et al. (2009)

Though as of late, 2‐D and 3‐D models have become fairly commonplace in evaluation of such data sets (e.g., Schildgen et al., 2009), our 1‐D approach is useful for isolating our parameter of interest, exhumation, without the need for simulating individual faults or assuming steady state topography. Our justification for using a 1‐D modeling approach comes from two factors. First, the regional sampling approach for data collection does not provide constraints on 2‐D (fault related) or 3‐D (topographic) processes. Considering faulting or topographic effects on the thermal field would be unconstrained in a 2‐D or 3‐D model and would result in over interpretating available data. Second, exhumation rates (demonstrated below) in the region are very low and below the threshold level of ~0.1–0.2 mm/yr where advective and 2‐D/3‐D thermal processes near faults become important (Ehlers, 2005; Ehlers et al., 2001; Mancktelow & Grasemann, 1997). In the Himalaya, 1‐D results have been shown to reproduce 2‐D and 3‐D models accurately in near‐surface thermochronometric systems above the ~350 °C isotherm (Thiede & Ehlers, 2013; Whipp et al., 2007). Though the slower erosion rates in the Andean fore arc may increase the relative effect of lateral variations in heat flow and produce slightly overestimated or underestimated exhumation rates, the added uncertainty is in our view compensated by the computational possibility of a much more complete range of derived exhumation histories (n = 140,000).

Thermal and spatial parameters for our model were obtained from Springer (1999) as previously done by Schildgen et al. (2009) in southern Peru and where applicable parameters were set to the same value given in the latter study (Table 2). We assigned constant thermophysical parameters and basal temperatures within each physiographic zone (Coastal Cordillera and Precordillera) and defined the depth and basal temperatures for each model space based on the subducting slab geometry of Springer (1999). Thus, although our modeling approach is 1‐D, it accounts for long‐wavelength lateral varitations in heat flow near a subduction zone by employing different basal temperatures in the model and hence spatial variability in upper plate thermal gradients as a function of distance from the subduction zone. Representative mean surface temperatures were used from the same previous studies for the two physiographic provinces and were held constant throughout the simulation as neither paleotemperatures nor reliable paleoelevation estimates are known.

To test exhumation histories between, as well as within, the CC and PC, one iso‐elevation bin in each domain was used as a representative zone for comparison of modeled cooling histories to observed ages in each given domain: we chose the sea level (CC) and 2,000‐m (PC) bins for our analysis. These bins had the highest amount of data and ZHe results allowing for a multichronometer approach. Both previously published (including AFT) and new ages were included. Samples from our study with high uncertainty (e.g., entirely overlapping 2σ range between thermochronometers) or which did not yield both AHe and ZHe ages were excluded. To incorporate the cooling ages from Juez‐Larré et al. (2010) into our model, we used values of 42.97 ± 2.36 Ma, 119.99 ± 11.66 Ma, and 170.9 ± 26.01 Ma for AHe, AFT, and ZHe, respectively, derived from inverse‐square‐error‐weighted means of sea level samples in the Tocopilla vertical transect of the study. Sample 15CL072 was sampled at the same location as AFT sample FT‐8 from Maksaev and Zentilli (1999; 47.4 ± 5.4 Ma) and modeled as one.

In each modeled elevation bin, we additionally calculated an error‐weighted mean age for both AHe and ZHe, using all sample ages for each chronometer within the bin. We use the AHe and ZHe mean bin ages as proxy values to derive a regional exhumation history model and compare it to local bedrock exhumation results (single samples). While the single‐multichronometer samples used are within 84 m vertically of each other in the CC and 243 m in the PC, mean proxy ages synthesize the full 500‐m elevation bin age range. Each sample was run individually for 20,000 iterations to maximize the number of unique solutions. For each model time step, model exhumation rates were selected randomly between 0 and 2 km/Myr assuming that mean exhumation rates did not exceed this range in the last ~200 Myr. A smart Monte Carlo inversion (genetic algorithm) was used to guide suitable exhumation histories with a built‐in randomization factor to break out of possible local minimums. Time steps were defined unique to each sample simulation, such that the number of time steps in a simulation matches the number of thermochronometers within that sample. The length of each time step is furthermore defined by the difference between thermochronometric ages, and each time step begins at a given chronometer age plus its calculated error. Thus, an individual exhumation rate represents the mean rate experienced by a sample between one closure isotherm and the next youngest closure isotherm (or the Earth's surface). Whereas choosing a longer time step size would average out signal in the data, a shorter time step size would result in an overinterpretation of the data.

In each Monte Carlo simulation, all predicted cooling ages were evaluated against the observed ages using a reduced chi‐squared () test between the observed sample age and the age predicted by the model against the uncertainty of each sample. In general, a higher uncertainty produces a larger number of acceptable fits and a more conservative estimate on possible exhumation histories. As a compromise, we define uncertainty in our model on the 1σ level for individual samples and on the 2σ level for the iso‐elevation bin proxy means because propagated error decreases with larger sample sizes. In each iteration, a cumulative chi‐squared ( ) was then calculated for a sample, where is the arithmetic mean of the values for each thermochronometer used in that sample. A threshold of 3, representing the 99.7% confidence interval, was considered to signify an acceptable erosion history, and samples within this range were saved to the model output.

Model Results

The 1‐D thermal model setup using the erosion rates based on our thermochronometric data and the boundary conditions and thermophysical parameters of Springer (1999) reproduces the 2‐D thermal structure of the original study well. A fixed temperature at depth in 1‐D defined by the angle of the subducting Nazca Plate in 2‐D results in a temperature inversion, slight in the CC and stronger in the PC, that mimics the thermal effect of heat being pulled down by a cooler subducting slab (Springer, 1999). Figure 8 demonstrates a range of thermal gradient solutions in the CC and PC based on 0.1 (“slow”) and 1.0 (“fast”) km/Myr exhumation rates in each region, indicating a minimum to maximum estimate of possible gradients in each region. Consistent with higher heat flow measurements recorded in the Precordillera compared to the Coastal Cordillera, the modeled PC geotherm is elevated, most importantly at depths of less than 20 km, above which all closure isotherms for the chronometers used in our study are found. Our results indicate near‐surface gradients of ~10–15 °C/km in the CC and ~15–25 °C/km in the PC. These ranges are consistent with previous thermal modeling in the Coastal Cordillera, which yielded gradients of under 20 °C/km since the Eocene for ~0.3 km/Myr erosion rates (Juez‐Larré et al., 2010). Modern borehole heat flow measurements have also calculated thermal gradients in our field area of ~10 °C/km in the CC and ~20 °C/km in the PC (Springer & Förster, 1998), consistent with the gradients associated with our slower exhumation rate simulations.

Figure 8

Geothermal gradient results from the Pecube 1‐D thermal‐kinematic model for exhumation rates of 0.1 km/Myr (representing “slow erosion” in blue) and 1.0 km/Myr (representing “fast erosion” in red) in the Coastal Cordillera (dashed) and Precordillera (solid). Depth to the subducting slab (gray dotted lines) is derived from the model results of Springer (1999). All four results display a temperature inversion at depth, and the Precordillera displays an elevated gradient at the surface and records higher present‐day heat flow measurements. Sample thermal gradients are indicated in the bottom left for reference.

Derived rates in the CC and PC are reasonably consistent within each physiographic zone between all accepted exhumation histories, including the iso‐elevation mean proxies (Figures 9 and 10). For each model simulation, apparent exhumation rates in the precursory time steps leading up to ZHe ages (and the AFT age for the Juez‐Larré et al., 2010, simulation) are consistently elevated (~0.8–1.5 km/Myr) with uncertainties always reaching the 2‐km/Myr upper limit allowed by our model setup. Beyond this first time step, the CC exhibits mean apparent exhumation rates between 0.04 and 0.16 km/Myr for all histories since ~120 Ma. The 10 m.a.s.l. sample from the Tocopilla VT, adopted from Juez‐Larré et al. (2010), has an older ZHe age of ~180 Ma and an AFT age of ~120 Ma matching the ZHe ages of our samples. Despite this, model results from this simulation are consistent within error to the others with a lower apparent exhumation rate in the AFT‐AHe time step (0.04 km/Myr) than that of the ZHe‐AHe time steps. In sample 15CL013 an increase in the apparent exhumation rate is evident since ~30–40 Ma but still holds under the 0.2‐km/Myr limit. In the PC, the three modeled exhumation rates are less consistent among themselves but largely maintain a mean exhumation rate between 0.06 and 0.21 km/Myr since ~50–60 Ma. Only a spike of higher rates in sample 15CL072 is apparent for a period between 44 and 53 Ma (AFT to AHe time step), though the exact value is dependent on the size of the time step over which this mean rate is given and is associated with significantly higher uncertainties (acceptable 3 range of ~0–0.8 km/Myr) . Nonetheless, this increase is not seen in other samples without AFT age constraints. Generally, in the PC, modeling points to a potential spatially variable decrease in mean exhumation rates since 50 Ma if all of the cooling signal is exhumation related (as opposed to cooling during crystallization).

Variation in apparent exhumation rates does not appear latitude dependent in either the PC or CC (compare north to south bedrock samples in Figures 9 and 10). Only the increase in exhumation rate of the northernmost CC sample, 15CL013, demonstrates the possibility of a spatial variation in exhumation. This trend, however, while hypothetically attributable to tectonic factors (e.g., location near plate syntaxis) may just as well be a result of proximity to an incised canyon, which could have formed in the last ~10 Myr (Jordan et al., 2010; Lease & Ehlers, 2013; Schildgen et al., 2009).

Discussion

Justifying Cooling‐Age‐Derived Exhumation Rates in the Andean Fore Arc

In the tectonically and magmatically active Andean fore arc, the interpretation of cooling ages is highly sensitive to the range of thermochronometers used, particularly in plutonic samples. Thermal effects of the magmatic arc and shallow pluton emplacement may limit the range of temperatures that track cooling through a rock's exhumation path, and care in interpretation must be exercised accordingly. We commonly find apatite U‐Th/He ages with large scatter at equal elevations even on small spatial scales due to the shorter‐wavelength topography‐induced thermal variability in the top few kilometers of the crust. ZHe and AFT ages, however, while representative of a deeper closure temperatures, which are less affected by topography, risk displaying the signal of a magmatic body's cooling unrelated to rock uplift and erosion. A pluton, emplaced at a shallower depth than the regional closure temperature of the investigated thermochronometer, will naturally display a signal of ambient conductive cooling despite never having been exhumed from the depth associated with that temperature. Attributing too much of the cooling signal to exhumation would, thus, result in modeled cooling histories overestimating exhumation rates. This also implies that barring any additional effects of magmatic intrusions elevating temperatures of our samples as they were being exhumed, model results are maximum estimates on exhumation rates.

A comparison of mineral cooling ages with crystallization ages for our data (SERNAGEOMIN, 2003, and references therein) shows an initial period of fast, followed by slow cooling, suggesting that magmatic bodies were emplaced near regional ZHe closure temperatures, cooling quickly to match the thermal surroundings. ZHe ages in both the CC and PC are generally within range or slightly younger than their mapped unit ages (SERNAGEOMIN, 2003), whereas lower temperature AHe ages are significantly younger. Thus, cooling ages more recent than ZHe dates should largely represent exhumational cooling through the top ~10 km of the crust (Figure 8). The poorly constrained apparent exhumation signal around the time of and especially prior to the ZHe age cannot be interpreted as exhumation cooling without deeper available thermochronometer data. In the first time step, most rates in the assigned 0–2 km/Myr erosion range yield an acceptable exhumation history in later time steps. This result is present even in the Juez‐Larré et al. (2010) sample simulation where a significantly older ZHe age with large error produces high (~1‐km/Myr mean) apparent exhumation rates albeit with a full 0–2‐km/Myr 3 acceptable erosion range (Figure 8).

Some lithologic units in the area, such as the Chuquicamata and El Abra copper porphyry deposits, have documented emplacement at ~2–3 km, depths currently above the AFT closure isotherm (Maksaev & Zentilli, 1999; McInnes et al., 1999). Though we have avoided sampling these higher elevation units for exhumation rates, the lower elevation Paleogene plutons of the Cerros de Monte Christo (samples 15CL072 and 15CL073) and Fortuna (samples 15CL074 and 15CL076) complexes run a similar risk. For sample 15CL072, a biotite 40Ar‐39Ar plateau age of 63.1 ± 0.3 Ma (Maksaev & Zentilli, 1999) derived from the same lithology and elevation is substantially older than the reported AFT age of 47.4 ± 5.4 Ma but concordant with our ZHe result of 60.3 Ma. At 2,800 m.a.s.l., however, Maksaev and Zentilli (1999) found concordant 40Ar‐39Ar and AFT ages indicating that while in the upper stretches the plutonic bodies may have been emplaced at AFT closure depths, the lower sections from which we derive exhumation rates were deposited below the AFT closure isotherm (though potentially near to the ZHe one). Higher elevation samples from the Fortuna complex demonstrate a similar problem but were not used for exhumation rate modeling (Maksaev & Zentilli, 1999).

Consistency between deeply emplaced plutonic cooling ages and buried then reexhumed volcanic or sedimentary deposits further suggests significant exhumation‐related cooling. In the PC, sample 15CL010, sampled from Jurassic sediments, yielded a ZHe cooling age of ~62 Ma, much younger than its deposition and consistent with the plutonic ZHe ages, suggesting that it was reset to ZHe closure temperatures. This trend is indicative of a period of burial, and subsequent exhumation through ZHe closure temperatures between the Jurassic and Cretaceous suggests that plutonic samples, despite their different provenance, experienced a similar exhumation history. Widespread deposition and burial in the Central Andes ending toward the Late Cretaceous have been recorded by others throughout this period in sedimentary sequences and thermal modeling (Michalak et al., 2016). In the CC, this ZHe trend is somewhat obscured by a poorly reproducing volcanic sample, 15CL017 (MSWD of 674). Scatter in the zircon data may imply significant time in the partial retention zone (e.g., Wolf et al., 1998) or maximal burial depths above the ZHe closure temperature. However, the AHe age in this sample (43 Ma) is consistent with other sea level samples and significantly younger than its Jurassic deposition age, suggesting burial and exhumation to at least the AHe closure temperatures.

It is also notable that the thermal effect of localized heat flow anomalies, generally located within ~10 km of active volcanoes, has been associated with shallow magma chambers at depths of ~4–6 km (Schildgen et al., 2009; Springer, 1999). In the PC, which is close to the present‐day arc, these depths correspond to depths between AHe and ZHe closure temperatures in our modeled thermal gradients. Although within the bounds of our study we have little constraint on where such anomalies may be, the reproducibility of ZHe ages across the PC and the agreement between the AHe ages used for modeling and mean iso‐elevation bin AHe ages suggest that the effect of such anomalies is likely inconsequential in the results. Long‐term net eastward migration of the arc has been estimated at ~1 km/Myr since the Jurassic (Andriessen & Reutter, 1994), which includes a period of westward retreat to its present location since around the early Miocene (Allmendinger et al., 1997); however, we see no clear imprint of thermal anomalies in the section of the fore arc sampled.

Interpretation of Low Exhumation Rates

Overall, bedrock thermochronometer‐derived exhumation rates in the Andean fore arc of northern Chile largely indicate background mean erosion rates (< 0.2 km/Myr) in both the Precordillera and Coastal Cordillera since at least ~50 Ma. While these rates may appear fast compared to other rates from cosmogenic isotopes reported in the west central Andes (e.g., Abbühl et al., 2010, 2011; Carretier et al., 2015; Kober et al., 2007, 2009; Placzek et al., 2010; Starke et al., 2017), we note that the rates reported here are averages over long timescales of 107 years and could contain within them punctuated periods of activity with long periods of little exhumation. Also, though 3 acceptable maximum values rise to ~0.8 km/Myr in the PC (Figure 10), these changes occur on shorter timescales and earlier in the mineral cooling histories where apparent exhumation is a less trustworthy indicator for mean regional exhumation or erosion. Moreover, because discrete time steps based on the cooling age data artificially define when the model allows erosion rates to vary, exact timing of changes is not easily constrained and direct comparison of mean rates integrated over time steps of different lengths should also be approached with care. With this said, uniformity of the exhumation history starts after closure of the ZHe system, such that the exhumation rates remain similar between closure of the two thermochronometer systems, and until present day. Rates prior to ZHe closure are not interpreted here, because they are more likely to influence the timing of pluton emplacement and postcrystallization cooling. In the CC this uniformity occurs between ~100 and 60 Ma, during the Cretaceous, and in the PC, between ~40 to 50 Ma, coeval with estimates of early Eocene Andes mountain building (Barnes & Ehlers, 2009; McQuarrie et al., 2005; Sempere et al., 1990). A possible decrease in exhumation in the PC could be indicated by the slightly lower rates in the AHe‐to‐present time steps as compared to the ZHe‐to‐AHe time steps. On the contrary, in the CC, sample 15CL013, sampled at the base of the Camarones canyon in the exorheic section of our field area shows the alternative possibility of increasing local rates in the north of the field area since the Miocene. Inhomogeneity in space and time is expected, however, as long‐term low erosion on a regional scale is driven locally by river incision while arid climates with little overland flow can yield negligible erosion elsewhere during these times (Jordan et al., 2010). Additional AFT data where it lacks would be necessary to further resolve the possible higher rates earlier in the sample cooling histories.

Results between AER regressions in the CC and thermal‐kinematic modeling also show general agreement. The magnitude of cooling events over AER time scales is ~0.1 km/Myr for a period between ~40 and 80 Ma (Figure 7). A punctuated exhumation history with temporal and spatial spurts of increased erosion, which is qualitatively suggested by the Juez‐Larré et al. (2010) data and expected in intermittently arid climates (Jordan et al., 2014), is indistinguishable at the resolution of our data from slow and steady erosion but not excluded because such spurts are subsumed in the modeled mean. First‐order estimates for such exhumation spurts can be assumed from maximum apparent exhumation rates in shorter time steps defined in our Pecube model but were still likely under 1 km/Myr (Figures 9 and 10). Using the linear regression results of the AERs is more complicated in this regard because the maximum slope within error yields an inverted trend: short‐term spurts on unrealistically short timescales cannot be decoupled from longer lower magnitude spurts. AER trends for lower temperature thermochronometers also generally overestimate exhumation rates because closure temperatures are affected by surface topography (Braun, 2002), although we have no evidence for significant relief change and do not expect this effect to be strong.

Longitudinal gradients are apparent in the raw ages but not in the derived exhumation rates. Thus, for proper interpretation, the documented trenchward increase in cooling ages must be decoupled from the regional thermal structure. Though exhumation in the CC and PC is resolved on different timescales, when accounting for a higher heat flow and geothermal gradient in the PC, apparent exhumation rates since the Eocene are consistent between both physiographic zones. As we compare only one iso‐elevation profile in each the CC and PC, only a bulk longitudinal trend is evident, though we do not expect results to vary at elevations between our profiles given the similarity between sea level and ~2,000 m.a.s.l. models. It is yet unclear, however, whether the longer‐standing slow exhumation rates since ~100 Ma in the CC are also slow during this time in the PC. While these rates suggest that aridity may have been present in the CC for much longer than many onset times have reported (Dunai et al., 2005; Hartley & Chong, 2002; Jordan et al., 2014; Rech et al., 2006), constraints on paleotopography are necessary to fully determine the origin of slow exhumation rates, which have topographic and geologic controls (slope and lithology) and climatic ones. Cosentino and Jordan (2017) used Sr isotope variations from coastal fog as a proxy for paleo‐elevations in the CC though an extensive and robust regional understanding of paleotopography is not yet established. Because the low closure temperature of AHe is highly susceptible to changes in topography, AHe signals have been used in the past as indicators of paleotopography and relief (e.g., Ehlers et al., 2006; Valla et al., 2010).

Apart from horizontal exhumational gradients and variable thermal structure, inverted age‐elevation trends on larger spatial scales (>10 km), as we find in our bulk data and PC AERs, can be a sign of decaying relief over time (Braun, 2002), that is, relatively higher denudation rates on ridge crests than valley bottoms. Though we do not have high spatial resolution constraints on exhumation rates across the entire fore arc, we have no other evidence to suggest a decrease in Andean fore‐arc relief and maintain that a substantial portion of the decreasing age signal comes from an elevated thermal gradient eastward, approaching the arc (e.g., Figure 8).

Latitudinal gradients in exhumation in the Andean fore arc may exist but are not extreme enough on the spatial and temporal scales of bedrock low‐temperature thermochronology and unlikely to be resolved with this technique. While exhumation rates are intrinsically linked with tectonically induced uplift rates, the low apparent exhumation rates that we document would not reflect changes in tectonic forcings if erosion rates cannot catch up to rock uplift rates on the relevant timescales. An erosional response to uplift is predicated in the presence of precipitation and depends on whether transport‐ or detachment‐limited conditions are assumed. In a long‐standing arid setting, tectonic forcings would not be well recorded in a slow‐to‐respond exhumational signal. Furthermore, documented knickpoints in the river profiles of the northern Chilean fore arc are an indicator of transient topography (e.g., Cooper et al., 2016; Hoke et al., 2007). Vertical fault motion could imprint local changes in exhumation and has been demonstrated to be the case in southern Peru (Noury et al., 2016). Our limited distribution of samples across individual faults precludes any meaningful calculation of fault kinematic history. However, we do not see consistent large‐scale thermochronometric differences in relation to mapped fault blocks in northern Chile, none of which are individually associated with enough displacement in the Miocene to be documented via low‐temperature thermochronometry. Brittle deformation in the Precordillera in the Paleogene was largely strike slip with a transpressive component that is unlikely to have resulted in the rates of exhumation that we observe, averaged over the scale that we address in this study (Andriessen & Reutter, 1994; Maksaev & Zentilli, 1999; Reutter et al., 1996; Tomlinson et al., 2001). Thus, spatial differences in regional cooling‐age‐derived exhumation rates in our field area are likely to be climate controlled. No consistent latitudinal trend is evident from our data, though the mean decrease of AHe and ZHe ages in the CC toward the north (Figure 5) and corresponding faster exhumation of sample 15CL013 is correlative to lower CC elevations and the presence of river canyons transporting Western Cordillera precipitation.

Comparison to Other Exhumation Rate Estimates

Previous studies in the Central Andean fore arc have quantified local exhumation rates in the Coastal Cordillera and Precordillera regions, though synthesis between ages across the physiographic areas in the fore arc has not been a principal focus. Vertical transect estimates using AHe data from the northern Chilean Coastal Cordillera (Juez‐Larré et al., 2010) and AFT from the Precordillera (Maksaev & Zentilli, 1999) reveal similar exhumation rates to ours of ~0.24–0.36 and ~0.1–0.2 km/Myr, respectively, for a period of ~10–20 Myr during the middle to late Eocene. Potentially higher exhumation in the Eocene is furthermore evident in our data for both the CC and PC (e.g., CC vertical transect data and PC Pecube sample run 15CL072, Figures 7a, 7c, and 9). Studies in southern Peru, north of the Arica bend, have yielded lower maximum estimates of ~0.08 km/Myr albeit averaged over a period spanning the Late Cretaceous to the middle‐late Miocene (Michalak et al., 2016). Schildgen et al. (2009), modeling AHe data sampled from deep southern Peruvian canyons, decoupled regional background exhumation rates from those in deeply incised canyons and concluded much lower background denudation rates of ~0.01 km/Myr from ca. 60 to 9 Ma (Schildgen et al., 2007, 2009). While Schildgen et al. (2007, 2009) suggested that background exhumation rates in the coastal regions of southern Peru were slower than the inland higher elevation rates based on a similar age distribution to what we find in northern Chile, when accounting for the variable thermal structure within the fore arc of northern Chile, we find that our data do not exclude a more uniform background exhumation across the entire present‐day fore arc over similar timescales. Using the same thermal model setup of Springer (1999) and the resulting 1‐D vertical thermal gradients that are possible in the CC and PC, we demonstrate the possibility of equally slow rates in the present‐day northern Chilean Precordillera since the Eocene. We should note, however, that whereas Schildgen et al. (2009) modeled background rates as those outside of and prior to areas of increased canyon incision, some of the increased incisional erosion signal is incorporated into our regional mean rates in the northern parts of our field area. General agreement between northern and southern field area samples implies that this difference should not be drastic. As compared to thermochronometric data in the eastern Andes where a distributed exhumation history has yielded exhumation rates ranging from 0.1 km/Myr in the eastern Cordillera to up to 0.4 and 0.6 km/Myr in the Interandean zone and Subandes (Barnes et al., 2008; Safran et al., 2006), the fore arc reveals a generally slower and older exhumation history driven largely by the arid climate that has developed on the leading edge of the Central Andean subduction zone.

More recent, shorter timescale denudation rates (<~106 years) based on catchment‐averaged 10Be, 26Al, and 21Ne cosmogenic nuclide concentrations and exposure ages have been reported, as low as <1 m/Myr, up to several orders of magnitude lower than bedrock thermochronometric exhumation rates in the southern Peruvian and northern Chilean fore arc but range up to 300 m/Myr (~0.3 km/Myr) at higher elevations and slopes, particularly in the modern‐day Precordillera (Abbühl et al., 2010, 2011; Carretier, Regard, Vassallo, Aguilar, et al., 2015; Carretier, Regard, Vassallo, Martinod, et al., 2015; Kober et al., 2007, 2009; McPhillips et al., 2013; Placzek et al., 2010; Starke et al., 2017). Generally, higher rates are reported in southern Peru than in northern Chile but correlations with precipitation or hillslope angles are inconsistent (Abbühl et al., 2011; Kober et al., 2009; Reber et al., 2017; Starke et al., 2017) and extreme climatic events (e.g., El Niño) may be important (Abbühl et al., 2010). While local erosion rates on hyperarid surfaces in the Central Depression of northern Chile may be negligible on millennial and even million‐year timescales (Dunai et al., 2005; Kober et al., 2007), intermittent hyperaridity (e.g., Jordan et al., 2014) and the wide range of local and catchment‐averaged erosion rates in the fore arc indicated by these methods is consistent with long‐term semiaridity to aridity suggested by mean exhumation rates on low‐temperature thermochronometry timescales. We further note that rates derived from samples in the Central Depression (Pampa de Tamurgal in our field area) represent values in a depositional basin where up to ~1,500 m of Miocene sediments (Jordan et al., 2010) would be stripped away at exhumation rates prescribed by our thermochronometric models. Thus, it is important to emphasize that our results cannot be extrapolated across the entire Atacama region and apply specifically to the two physographic areas where they are sampled (CC and PC). Some Miocene volcanic cover in the central section of our field area extends into the PC and records possible local areas of net deposition, though the original thickness of these units is unconstrained and exhumation cannot be directly inferred locally from our results to the north and south. Though thermochronometry can be a powerful tool in constraining changes in the top several kilometers of the crust, it is also imperfect in the slowly eroding environments such as the Central Andean fore arc, where rocks that potentially carry an AHe signal of a transition to even slower exhumation rates have yet to be exhumed. As such, even AHe data are not sensitive to more recent differences between aridity and hyperaridity. Nonetheless, if decreases to erosion rates as low as 0.2 m/Myr (e.g., Placzek et al., 2010 ) in the Atacama region since the Miocene were stable, regional, and long‐lived, higher rates than those modeled in our means during the Eocene and Oligocene would be necessary to produce the observed AHe ages.

Conclusions

To better understand the effect of long‐term tectonic and climatic controls on the evolution of the fore arc in the Central Andes of northern Chile, we interpreted apatite and zircon U‐Th/He ages and apatite fission track ages with a 1‐D thermal‐kinematic model and derived exhumation rates over 10‐Ma timescales. Interpretation of cooling ages in the Andes is heavily dependent on the geothermal gradient, which is fixed in neither time nor space and particularly sensitive to tectonic parameters in a subducting margin. Thermal‐kinematic modeling accounts for changes in these gradients based on boundary conditions derived from more complex 3‐D models. Higher heat flow and elevated thermal gradients in the Precordillera contribute significantly to the younger documented cooling ages east of the Coastal Cordillera. Any unconstrained thermal anomalies from the magmatic arc would be nonetheless difficult to decouple.

While much of the earlier cooling signal from plutonic samples is controlled by postcrystallization cooling of the magmatic body, later stages of apparent exhumation are likely representative of advection of material through rock uplift and erosion. Volcanic and sedimentary samples are in agreement with plutonic ones and were buried to AHe or ZHe closure temperatures by the Late Cretaceous. Since then, relatively consistent mean exhumation rates of <0.2 km/Myr have been pervasive throughout the fore arc of this study area. These rates nonetheless represent maximum estimates since some of the cooling signal is still potentially a result of ambient cooling of the plutonic bodies. Modern and local erosion rates in the region have been recorded orders of magnitude lower and at times at a virtual standstill. Areas of the Precordillera, much like the Central Depression, even record large sections of deposition through much of this time (approximately Miocene), and mean exhumation rates should be taken at a regional scale.

We find no evidence for significant gradients in exhumation in the field area either latitudinally or between physiographic zones (CC and PC). Such gradients would likely be tectonically controlled on these timescales if they existed. Differences that appear in the data can most easily be explained by increased local incision. As slow rates imply low erosion and aridity, it is unlikely that exhumation would respond significantly to tectonic uplift under these conditions. However, we have no evidence to suggest an exhumational or topographic steady state. As such, we see no significant effect of the Arica Bend syntaxis on exhumation rates in the northern Chilean fore arc and maintain that climate has been similar in the fore arc on long timescales (up to ~50 Ma).

Supporting Information S1

Click here for additional data file.