Low-temperature thermochronology in the Peruvian Central Andes: implications for long-term continental denudation, timing of plateau uplift, canyon incision and lithosphere dynamics

Abstract: In Peru, the western edge of the 4.5 km high Western Cordillera is cut by a >3 km deep canyon. To understand incision by the Cotahuasi–Ocoña River and the regional uplift history of this orogenic plateau capped by volcanic rocks, 26 crystalline rock samples were collected for low-temperature thermochronology from vertical profiles parallel and perpendicular to the canyon. Rock cooling histories confirm that most plateau denudation had occurred prior to 24 Ma but plateau incision peaked after c. 14–9 Ma in response to rapid surface uplift. The abrupt occurrence of a rock heating event is also detected during middle Miocene time. This was either a response to the emplacement of low-conductivity, regionally extensive ignimbritic caprock or a response to crustal-scale fluid circulation caused by wet melting of the overriding plate when magmatism resumed c. 24 Ma. The potential for thermochronology to provide information on past geothermal gradients is discussed, showing how it can be used as a proxy for understanding change in subducting slab dynamics, with oscillations in subduction angle having perhaps been the main on–off switch for magmatism in this Cordilleran setting. Supplementary material: Fission-track and apatite (U–Th)/He results and apatite chemistry indicators are available at http://www/geolsoc.org.uk/SUP18405.

Plate convergence zones are the focus of intense denudation in response to the gravitational potential afforded by crustal deformation, often resulting in high topographic relief. Despite being situated in an area of plate convergence for c. 200 Ma, the orocline of the Central Andes exhibits large expanses of lowrelief, elevated topography ( Fig. 1), collectively termed 'Altiplano', that would instead seem typical of more stable tectonic environments. In the Western Cordillera of southwestern Peru, the orogenic plateau edge is c. 4.5 km high and fluvial incision by the Cotahuasi-Ocoña River has exposed Cretaceous plutonic rocks in a c. 3 km deep canyon system. This results in exceptionally high topographic relief amidst the otherwise low-angle slope system of the Cordillera and Pre-Cordilleran forearc (Fig. 1). The scale of incision is large enough for low-temperature thermochronology to be used to constrain, from rock cooling histories, the timing and depths of erosion prior to and since canyon incision. Unlike previous thermochronological studies of the Central Andes, the data and methods used in this study also provide scope for discussing what proportion of the temperature signal is attributable to erosion rather than to any other non-erosional crustal process such as hydrothermal fluid circulation or change in crustal heat flow, thus also contributing insight into the dynamics and history of magmatism in the South American lithosphere.
Low-temperature thermochronology, which in this study integrates zircon fission-track (ZFT), apatite fission-track (AFT), and apatite (U-Th)/He (or AHe) analysis, is used to constrain the thermal history of rocks in the upper few kilometres of the Earth's crust, and is suited to monitoring the depth and timing of continental denudation over time scales of 1-100 Ma. Cooling histories derived from the observed data are typically converted to a rock exhumation rate, assumed to be driven by surface erosion, by assuming palaeotemperature gradients. Factors that can lead to uncertainty in this calculation include spatially nonuniform gradients caused by lateral variations in valley incision depths (Stüwe & Hintermüller 2000) or by variations in rock thermal conductivities, or unsteady gradients owing to local resetting by magmatic activity, to burial beneath low-conductivity cover rocks ), or to advective cooling or heating by crustal-scale fluid flow. The possibility that vertical temperature gradients may have changed over the course of orogeny is far from negligible, and is important because it affects the accuracy of computed depths of long-term denudation. This, in turn, can distort our understanding of palaeoclimatic and tectonic forcing factors in orogenic growth, and ultimately influence our interpretation of relief development and exhumation rates.
In this study we exploit modelling of single samples (Ketcham 2005), and joint inverse modelling of thermal histories from multiple samples (Gallagher et al. 2005), to assess the nature of palaeotemperature gradients. In total 26 samples were collected from vertical profiles for analysis, oriented parallel and perpendicular to the Cotahuasi-Ocoña canyon system and located at elevations between 0.5 and 3.8 km (Fig. 2). The interpretations obtained from our new thermochronological data and modelling approaches are compared with independent canyon incision chronologies based on a separate study of in-canyon (U-Th)/He ages (Schildgen et al. 2007) and on a large dataset of 40 Ar/ 39 Ardated volcanic materials that either drape the plateau surface or infill the canyon (Thouret et al. 2007). The results of this study emphasize the benefits of using coupled AFT and AHe approaches, and show that, as a result of the respective sensitivities of these two methods, important geological information can remain undetected if only one of them is used.

Regional geology and topographic setting
Important geological and geophysical variations occur at all scales along the strike of the Andes. The study area (Fig. 1) lies astride the Coastal Cordillera, the Western Cordillera and the westernmost edge of the Altiplano between 15 and 178S, at the northern end of the modern Central Andean steep subduction segment. We begin with an overview of the geological and topographic setting of the Peruvian canyon territory along the Western Cordillera (13-158S and 71-73830'W), in which we demonstrate that this area is not just an extension of better studied regions of Andean subduction located immediately to the south (e.g. Chile; Oncken et al. 2006).
Subduction associated with Andean mountain building began at least 185 Ma ago. After an extended period of Mariana-type subduction dominated by back-arc extensional subsidence between c. 200 Ma and c. 90 Ma, the construction of the present- day Andean Cordillera began with a switch to an Andean-style subduction regime, involving intense magmatic episodes interspersed with compressive tectonics (Benavides-Cáceres 1999;Gregory-Wodzicki 2000). Parallel to the strike of the Pacific trench, the marine volcanic arc changed to a continental magmatic arc, and a forearc basin began to form during the late Mesozoic on the western flank of the current Western Cordillera. Uplift began in the Western and Coastal Cordilleras c. 50 Ma ago, perhaps as early as 70 Ma (McQuarrie et al. 2005), and continued later and more slowly in the Eastern Cordillera (40-10 Ma) while undergoing a marked interruption in the Western Cordillera. Uplift resumed in both Cordilleras and in the Altiplano c. 30-20 Ma in phase with periods of higher convergence rates between the Nazca and South American plates (Sébrier et al. 1988), but accelerated significantly after 10 Ma (Schildgen et al. 2007(Schildgen et al. , 2009Thouret et al. 2007;Garzione et al. 2008) in a context of abruptly diminishing convergence rates .
In this Peruvian canyon territory, the 1.5-2.0 km high Coastal Cordillera consists of a Precambrian to Palaeozoic metamorphic basement, the Arequipa Massif (Ramos 2008), covered by a sequence of Cenozoic forearc sediments. Precambrian gneisses ( Fig. 1) crop out wherever they are not overlain by Lower Palaeozoic metamorphic conglomerates and sandstones. James & Sacks (1999) have emphasized the general coincidence between the cratonic Arequipa terrane, which extends eastward, and the modern Altiplano, suggesting that its mechanical strength has been sufficient to withstand successive Mesozoic and Cenozoic orogenic cycles without undergoing the deformation observed elsewhere in the Andes.
The Cenozoic inner-forearc sediments that separate the risen Western Cordillera from the Pacific coastline are known as the Moquegua Formation. They were deposited in an extensional setting west of the magmatic arc and comprise 0.5 km of Eocene to Oligocene continental sandstone and siltstone (Huamán 1985;Roperch et al. 2006). The 0.7 km thick upper sequence of the Moquegua Formation lies unconformably on the lower sequence of the same formation, and forms an upward coarsening unit of nonmarine conglomerate, sandstone and siltstone with an increasing frequency of tuff intercalations. This accumulation of alluvial debris and lahar deposits was sourced by the Western Cordillera between 34 and 16 Ma during an early phase of progressive uplift. Clasts in the lower Moquegua 'A' sequence (Roperch et al. 2006) include coastal batholith rocks (Huamán 1985), demonstrating that exposure of the plutons sampled in this study had already occurred by Oligocene time.
Finally, the volcanic arc itself and its suite of igneous rocks spans three geological eras of intermittent magmatism (Figs 1  and 3). The oldest arc rocks belong to the Chocolate  Ma median ages given by Mamani et al. 2009) and Toquepala (78-50 Ma median ages) Formations. The age range of those arcs is, however, not precisely constrained because of the effects of rock weathering on K-Ar age determinations. In the study area, where rocks belonging to the Chocolate Formation are absent (Fig. 1), granodiorite and tonalite units, collectively known as the Andean batholith, have been ascribed to middle Cretaceous and late Cretaceous-Palaeogene magmatism (e.g. Pecho Gutierrez 1983;Beckinsale et al. 1985;Mukasa 1986). These crystalline rocks crop out where the studied canyons and their tributaries have cut deepest into the crust (Fig.  1). After a c. 20 Ma volcanic lull between c. 50 and 30 Ma (Sandeman et al. 1995;Benavides-Cáceres 1999), the Tacaza Arc (30-24 Ma) produced volcanic rocks that are preserved within the fluvial sequences of the forearc sedimentary sequence. The Alpabamba and Huaylillas ignimbrites (24-13 Ma) of the Huaylillas Arc (24-10 Ma: Mamani et al. 2009;Figs 1 and 3) are stratigraphic equivalents to the Oxaya ignimbrites of northern Chile, and predate the major phase of deep canyon incision (Thouret et al. 2007). These large volumes of welded rhyodacitic ignimbrite are preserved as undeformed, plateau-like cover up to 1 km thick north of the town of Cotahuasi (Fig. 1), and they thin out towards the west over the forearc prism, mostly as a result of inherited depositional thickness. Their extensive preservation supports the idea of low erosion rates since 13 Ma (Thouret et al. 2007). The Lower Barroso Arc (10-3 Ma; Mamani et al. 2009) corresponds to eroded stratovolcanoes whereas the Upper Barroso Arc (3-1 Ma) corresponds to better preserved stratovolcanoes, domes and cones (Fig. 3). Another magmatic phase during the Pliocene produced the less voluminous Sencca ignimbrites, which are inset in valleys and therefore post-date the onset of canyon incision (Thouret et al. 2007). Finally, the modern Frontal Arc is situated c. 120 km above the Wadati-Benioff zone (Fig. 1). It formed less than 1 Ma ago on the west side of the Western Cordillera, and the study area includes several coeval composite cones shown in Figure 3 (Thouret et al. 2007;Mamani et al. 2009).

Topographic setting
Instead of an abrupt topographic step typical of many mountain fronts, the Western Cordillera forms a large monocline (Fig. 4) that extends from Nazca to the Chilean border (Dollfus 1973). Thus, in sharp contrast to the much more rugged topography of the eastern flank of the Andes, here the orogenic plateau connects to the Pacific coastline via a 60 km long topographic ramp, which descends from c. 4.4 km above sea level (a.s.l.) to sea level, and where mean gradients vary between c. 88 and 28 (Cuno Cuno; Fig. 2). Where rivers have incised the monocline the canyon drainage tends to expand laterally along the hinge, developing a WNW-ESE-trending anticlinal valley and exposing  Figure 1. Location of mines exploiting hydrothermal ore deposits is according to Pecho Gutierrez (1983) and Olchauski & Dávila (1994).  Gutierrez (1983); Olchauski & Dávila (1994); updates according to Thouret et al. (2007). the strike-parallel batholith (Fig. 1). Although it mimics the geometry of mountain-front pediments and fans (Tosdal et al. 1981(Tosdal et al. , 1984Laharie 1985;Clark et al. 1990), the slope of such a large, continuous monocline is primarily an effect of crustal deformation of the pre-Huaylillas land surface in response to crustal mechanisms that will be discussed below. Flexural models of the overriding plate can explain this feature by invoking a steep coast-to-Cordillera lithospheric rigidity gradient (Tassara 2005;Fig. 4).

Analytical FT and AHe methods
The outcrops that were sampled occur along the Ocoña valley system in the middle to late Cretaceous plutonic units termed Tiabaya (Pecho Gutierrez 1983) and Incahuasi (Olchauski & Dávila 1994). Fission-track analysis on zircon and apatite (ZFT and AFT, respectively) exploits information contained in the thermally activated shortening (or annealing) process of metastable fission tracks, which are physical damage zones generated in the crystal lattice by spontaneous fission of 238 U atoms (Dumitru 2000). A sample's ZFT age captures cooling through the zircon closure temperature, which varies as a function of radiation damage in crystal grains and ranges between c. 270-300 8C according to laboratory experiments (Yamada et al. 1995) and c. 187-223 8C according to field-calibrated studies (Bernet 2009). A sample's measured AFT age and spontaneous track length distribution together record upper crustal residence history below total annealing temperatures (,110-135 8C for exposure times .10 6 -10 7 years for common chemical compositions of apatite). In most cases, AFT thermal histories record sample cooling driven by rock uplift and exhumation processes. The AHe thermochronometer is based on the accumulation in apatite crystals of radiogenic 4 He from the decay of 235 U, 238 U and 232 Th series nuclides. Helium retentivity in the crystal allows a temperature-dependent closure age, normally younger than the AFT age, to be calculated.

Modelling of thermal histories and inference of palaeotemperature gradients
AFT partial track annealing zones and AHe partial helium retention zones (PAZ and PRZ, respectively) overlap. For lowchlorine apatites, this property can help to reconstruct rock thermal histories between c. 120 8C and c. 40 8C. AFT and AHe ages can be interpreted as cooling ages if cooling is monotonic and fairly rapid because annealing and diffusive losses remain limited. However, when cooling is slower or locally complicated, 4 He is partially retained at temperatures between c. 80 and 40 8C (Wolf et al. 1998), fission track length distributions reflect higher levels of track annealing, and rock time-temperature (t-T) paths are less straightforward to interpret. Furthermore, because AFT annealing kinetics are sensitive to apatite chemical composition,  Tassara (2005), reaches a regional minimum in the study area.
we measured the chlorine content (wt %) and the kinetic parameter D par on a range of samples. Fission tracks are optically revealed by chemical etching of the crystals, and D par (in ìm) is the mean fission-track etch pit diameter parallel to the crystal c-axis. It serves as a proxy for determining chemistry-dependent AFT annealing behaviour at crystal level (Donelick et al. 1999).
To recover information about thermal histories from measured AFT and AHe data, searches of randomly generated cooling histories were performed by two methods. First, thermal histories of single samples were generated with the HeFTy modelling package (Ketcham 2005), which includes the multi-compositional AFT annealing algorithm of Ketcham et al. (2007). This software performs inverse modelling of the AFT and AHe data with a Monte Carlo search for better fitting cooling paths through a set of user-defined t-T boxes. The outputs are thermal histories that predict AFT and AHe parameters that most closely match the measured data. Goodness-of-fit between predicted and observed parameter distributions is based on the exceedance of a statistical merit value defined by Kuiper's Statistic (Ketcham 2005). Because no thermal history ever uniquely defines a set of AFT or AHe data, the better fitting thermal histories must always be considered in the context of local geological constraints.
The second approach is based on the method presented by Gallagher et al. (2005). This method incorporates essentially the same annealing and diffusion algorithms as that of Ketcham (2005). However, the main difference is that a suite of samples from a vertical profile are modelled jointly, using both the AFT and AHe data as available. The key assumptions in this approach are that all samples have the same form of thermal history, and that the highest elevation sample has always been colder than the lowest. In this approach, we use the highest elevation sample as the reference point for inferring the thermal history and we introduce additional parameters into the modelling procedure; that is, the temperature difference (or offset) between the highest and lowest elevation samples. The thermal histories for all other samples are determined by linear interpolation. The method is formulated in a Bayesian framework, and uses a trans-dimensional approach (e.g. Gallagher et al. 2009) so that we do not explicitly choose the number of time-temperature points in the thermal history. Rather the data themselves determine the complexity of the solution.

New constraints on continental denudation in the Central Andes of SW Peru
The ZFT ages (Fig. 5) range from 91 AE 4 Ma to 126 AE 7 Ma. These are consistent with other ZFT data for this Mesozoic batholith (Wipf 2006), and with Rb/Sr, K/Ar and U/Pb granite crystallization ages (which peaked c. 109-100 Ma; Pecho Gutierrez 1983;Beckinsale et al. 1985;Mukasa 1986;Clark et al. 1990). This also suggests regionally uniform cooling conditions. If cooling through the c. 230 8C zircon closure temperature had been caused by significant erosional unroofing, a geothermal gradient of, for example, 25-30 8C km À1 would position these rocks at a palaeodepth of 8-10 km during late Cretaceous time. However, petrological studies have shown that the pluton was emplaced at a depth of c. 3 km, spreading laterally to form tabular bodies (Bussell & Pitcher 1985;Haederle & Atherton 2002). The ZFT ages, therefore, mostly reflect the terminal cooling phases of magma crystallization during Cretaceous time. This interpretation is supported by the fact that cratonic rocks exposed near the coast (Figs 1 and 3 (2ó); Wipf 2006), also confirming that zircon crystals at those locations were not reset by heating associated with emplacement of the Andean batholith.
The AFT ages range between c. 105.0 AE 3.9 and 11.9 AE 1.8 Ma, whereas AHe ages range between 54.4 AE 3.1 Ma and 1.6 AE 0.1 Ma (Fig. 5). Unlike the ZFT data, these age ranges are similar to AFT age ranges in the coastal Precambrian rocks obtained by Wipf (2006), and unambiguously record regionalscale post-magmatic unroofing and/or changes in crustal temperature gradient. The ages fall into two elevation-dependent groups with distinct cooling signatures. The more elevated samples, hereafter termed canyon-rim samples (OCO06, OCO07, OCO08, OCO09, OCO23, OCO04-01; Fig. 2), all exhibit long mean track lengths (MTL . 14 ìm). Such signatures are typical of rocks that underwent rapid cooling followed by a long residence time at near-surface temperatures (Fig. 6). With distinctly shorter tracks (typically 12-14 ìm) but an age range poorly distinguishable from that of the canyon-rim samples (i.e. Age-elevation plots of AFT, ZFT and AHe data. AFT age distribution suggests a thermally stratified crustal section. The crustal section exhibits an uplifted former AFT partial annealing zone (PAZ), where higher temperatures in the past caused greater fission-track shortening, and hence smaller MTL values, whereas the shallower and cooler samples exhibit larger MTL values. AHe ages show significant dispersion, interpreted here as having been caused by an unevenly distributed crustal heating event. This occurred at relatively low temperatures, and thus affected the helium thermochronometer more than the fission-track system. 55-105 Ma), the lower elevation, in-canyon samples (Figs 2 and 5) represent exhumation from deeper crustal levels. Dispersion of AFT ages as a function of crustal depth is typical of a fissiontrack PAZ. Along with the clearly stratified MTL distribution (Fig. 5), the pattern reveals that the Ocoña canyon has incised a thermally stratified crust in which apatites have resided at a range of partial annealing temperatures over the last c. 100 Ma.

Age of canyon incision brackets the timing of major Andean plateau uplift
In-canyon bedrock cooling signatures imply that uplift-driven valley incision in the Western Cordillera began after 14 Ma (Schildgen et al. 2007;Thouret et al. 2007;Picard et al. 2008). AHe valley-floor ages also decrease upstream from c. 13 Ma (OCO04-04, OCO11, OCO12) to c. 2-4 Ma (OCO17, OCO18) (Fig. 5). Furthermore, upstream from the confluence between the Marán and Cotahuasi rivers (OCO17, OCO18), the canyons form V-shaped bedrock gorges, whereas downstream the valley becomes a c. 1 km wide braided channel (Fig. 2). These features indicate asynchronous canyon incision: as relief increased, steepened river reaches, or knickzones, migrated upstream in successive waves, causing delayed deepening of the V-shaped canyons, while also eroding successive generations of pyroclastic and other unwelded canyon-fill deposits (Thouret et al. 2007). Using the Huaylillas ignimbrite deposits as a marker horizon (Thouret et al. 2007), it follows that depths of fluvial incision since 13 Ma are approximately equal to canyon depth.

Constraints on pre-canyon denudation rates
In contrast to the in-canyon samples, the canyon-rim samples record the geomorphological history of the low-relief plateau surface. After a period of rapid cooling from ZFT to AFT critical Fig. 6. Typical single-sample cooling histories. (a) Cooling history for a sample with no AHe data. The AFT thermometer appears insensitive to the magnitude of reheating that occurred after 20-15 Ma, highlighting the information lost when AFT and AHe data are not or cannot be used jointly. (b) Cooling history typical of samples from the Ocoña west bank, where coupled AFT and AHe data do not detect any occurrence of heating after 20-14 Ma. This explains the older AHe ages for samples from that entire area, suggesting that the rocks have been close to the surface since c. 50 Ma and that Cenozoic denudation has been limited. (c) Example of an in-canyon sample that was subjected to reheating c. 14 Ma. The heating event explains the difference in AHe age between OCO11 and samples such as OCO24. In all panels, darker and lighter grey shading define good and acceptable model thermal histories, respectively. Dashed and dotted horizontal lines mark the approximate temperature ranges of the AFT partial annealing and AHe partial retention zones, respectively. Black curve adjusted to AFT length data (histograms) illustrates the best model fit obtained for each sample. The two peaks of canyon incision (black vertical bars) determined independently from 40 Ar/ 39 Ar dating of canyon-fill volcanic deposits (Thouret et al. 2007) reveal the insensitivity to incision of the most elevated canyonmargin rocks such as OCO07, but show a good match with AHe cooling ages within the deeply incised area (OCO11).
temperatures, which we can infer to have occurred between c. 120 and c. 90 Ma based on the gap between ZFT and AFT age clusters shown in Figure 5, the decline in cooling and, by inference, in erosion rates after 80-70 Ma (see Fig. 6) indicates either a change in surface conditions (i.e. a more arid climate, or a reduction in topographic slope gradients; e.g. Gunnell et al. 2009) or an abrupt reduction in heat flow. James & Sacks (1999), for example, have argued that the Palaeogene period in SW Peru coincided with a period of rapid crustal cooling correlating with the onset of low-angle subduction. We discuss in the next section how a reduction in heat flow almost certainly contributed to trends in the cooling signal, but here we suggest that aridity and low-gradient topography also conspired to limit continental denudation. There is also independent evidence for long-term geomorphological quiescence during Cenozoic time in the Western Cordillera (Hoke et al. 2007). For example, near Arequipa ( Fig. 1) 40 Ar/ 39 Ar ages from alunite-group minerals preserved within Neogene weathering profiles suggest negligible depths of mechanical erosion since 36-39 Ma at least in the forearc (Quang et al. 2003). Although such evidence of low geomorphological activity mostly documents the period during the 50-25 Ma volcanic lull, the AFT and AHe data suggest that the period of low denudation may have started earlier, and that canyon-rim samples remained at low temperatures (and hence shallow depths beneath the topographic surface), between c. 60 and c. 20 Ma (Fig. 7).
Early exhumation of samples to shallow depths is supported by the relatively old AHe ages collected on the west bank of the Ocoña canyon, irrespective of elevation (OCO24, OCO25, OCO21, OCO04-07, OCO14; see Fig. 2). The difference between these older AHe ages (.35 Ma), which are predominant to the west, and the younger AHe ages (,25 Ma), which prevail to the east of the canyon's axis, is potentially problematic. The effect of variable grain size is taken into account in the modelling, so we do not consider this to be an explanation of the older ages. Furthermore, there is no obvious global correlation between the concentration of helium's parent elements (U, Th) and the ages, so we do not consider the older ages to reflect preferential retention of helium as a result of radiation damage (e.g. Schuster et al. 2006). For example, the oldest ages in sample OCO04-07 (at 1690 m) have a range of effective U (eU) of 25-60 ppm, whereas the much younger ages in sample OCO05 (2200 m) are associated with a range in eU of 15-30 ppm. Also, sample OCO04-09 (1200 m) has the highest range of eU (160-280 ppm), but the youngest of the three ages has the highest value.
Overall, the consistency in spatial distribution between these two populations seems to override any other consideration. Unless a so far undetected structural divide exists in the study area, the pattern suggests a role played after 30-25 Ma by spatially heterogeneous reheating of the apatites between areas situated east and west of the canyon, respectively. In this region, where ore mineralizations are widespread and mined (Fig. 2), anisotropic fluid circulation in the plutonic rocks could be responsible for the younger AHe ages observed in canyon-floor and east-bank samples, with only some older AHe ages having been rejuvenated by fluid circulation (e.g. OCO11, Fig. 6). The spatial heterogeneity of AHe age resetting by hot fluid circulation may well have occurred on a finer scale than this simple eastwest divide corresponding roughly to the canyon axis. For example, sample OCO26, which occurs on the Ocoña right bank amidst a cluster of c. 55-35 Ma sample ages, has yielded two 20-22 Ma aliquot ages as well as one aliquot age of 43 Ma (Fig.  5). Clearly, at current sampling density levels the resolution of the data cannot provide further insight into this hypothesis.
A c. 20-14 Ma crustal heating signal Use of the combined sample approach to thermal modelling (Gallagher et al. 2005) was motivated by the high relief of the study area, which allowed a selection of samples to be treated as a profile that is nearly vertical and has not been structurally disrupted by major post-cooling fault displacements (Fig. 3). We selected 10 samples distributed up the canyon between 0.5 and 3.6 km a.s.l. In this approach, temperature offsets between the highest and lowest samples are inferred along with a variable number of time-temperature points that define the thermal history. In this approach the data determine how many timetemperature points are required (i.e. the model complexity), and the temperature offsets can then be interpreted as palaeogeothermal gradients. Results based on three model scenarios combining AFT and AHe data are presented in Figure 7.
The prior information built into all models assumed that the temperature history can pass through any point between 0 and 200 8C from 200 Ma to the present day, and the temperature offset between the uppermost and lowermost samples was allowed to vary between 0 and 100 8C. These broad bounds encourage the solutions to have fairly simple forms unless the data require the thermal histories to incorporate more complexity. The three models we ran had slightly different parameters. For model 1 the present-day temperatures were fixed (using a value of 5 8C for the shallowest with an atmospheric temperature lapse rate of about 6.5 8C km À1 ), whereas in model 2 the present-day temperatures were free to vary between 0 and 40 8C. In both of these models the palaeotemperature gradient was assumed to be unknown but constant over time. In model 3 we fixed the present-day temperatures as in model 1, but allowed the palaeotemperature gradient to vary over time.
The results for all models (Fig. 7) imply rapid cooling around 100-80 Ma, as starting the thermal history at low temperatures (e.g. model 3) requires rapid cooling from temperatures in excess of 120-150 8C immediately prior to this time. These results are consistent with the interpretation of the ZFT data, and broadly correspond to the range in AFT ages. The solutions shown in Figure 8 either have progressive cooling up to about 15 Ma followed by a short-lived (,1 Ma) reheating event, to temperatures of 80-100 8C, or slow heating to temperatures around 70-80 8C followed by cooling to surface temperatures starting about 20 Ma. In terms of the fit to the observed data, it is clear that although we fit the general elevation trend for the AFT age (and the MTL data are scattered around what is essentially the mean value), the older AHe ages are not well predicted by any of these joint models. In all cases, the inferred temperature offsets are around 10-15 8C, which implies temperature gradients of 3-5 8C km À1 . Such gradients are clearly very low, even allowing for the uncertainties of AE2-4 8C km À1 .
The results from the combined data models all reveal two key characteristics: the gradients have remained consistently very low during the Cenozoic, and all models imply an apparent increase in gradient c. 20-14 Ma. This increase for models 1 and 3 (Fig.  7) could be the result of fixing the present-day temperatures to values that require a higher gradient than that inferred for the past. However, the results from model 2, in which we allowed the present-day temperatures to vary, imply a similar increase also. Given this, we also note that the resolution of the increase is close to the limit of the uncertainties, except perhaps for model 1. The low geothermal gradients produced by all models are nevertheless consistent with data from other forearc settings, where heat flow is not only low but also at least 60% less in the overriding plate than in the subducting plate (Stein 2003). Palaeogradients of 6-10 8C km À1 , for example, have been computed in the Sierra Nevada, California (Dumitru 1990(Dumitru , 1991. The lowest heat flow values in the Andes are currently recorded in the low-angle subduction areas between 4 and 158S and 27 and 338S, respectively. These values of 20-40 mW m À2 (Henry & Pollack 1988) coincide with the currently non-volcanic segments of the Andes, suggesting that volcanic lulls should correlate with low-angle subduction and with anomalously low geothermal gradients. Given the volcanic lull between 50 and 25 Ma in the study area (James & Sacks 1999;Mamani et al. 2009), it is reasonable to assume that low-angle subduction prevailed at the time. The abnormally low heat flow values resulted perhaps from advective cooling caused by slab dewatering and hydration of the overriding plate (James & Sacks 1999). The model results suggest that after 20 Ma temperature gradients increased somewhat, although they remained low enough to be still compatible with subduction settings (Stein 2003). The reheating event after 20 Ma remains undetected in thermal models obtained from samples that yielded either only AFT results (Fig. 6) or only AHe results (Schildgen et al. 2007(Schildgen et al. , 2009, respectively. Abrupt reheating by 40-60 8C during the middle Miocene seems, however, to be required when thermal histories are constrained by both thermochronometers.
In methodological terms, we find that results produced by each of the two modelling approaches are broadly similar. However, the multiple-sample modelling approach will not provide as a good a fit to the observed data as modelling single samples. It tries instead to find a compromise solution that provides an acceptable fit to all the available data for all samples simultaneously. This means that we can use samples with no track length data, samples with only AHe or only AFT ages, as well as all the samples with a complete set of data. Importantly, we also obtain from the offset parameters a direct estimate of the temperature gradient over time, which is required for converting the temperature histories to an equivalent denudation depth.

Causes of middle Miocene crustal heating
Both the magnitude and timing of the heating event detected by the data point to two possible causes: the rapid response to instant burial beneath the freshly erupted Alpabamba (22-18 Ma) and Huaylillas (14.3-12.7 Ma) ignimbrites around that time (Thouret et al. 2007), or the delayed response to a steepened subduction angle. In the first scenario, the relatively low thermal conductivity of the ignimbrite caprock (which can be as low as 0.35 W m À1 K À1 presumably as a consequence of air-filled pore-space; e.g. Keating 2005), with preserved thicknesses of 350 m on the canyon rim (Fig. 3), could explain why the temperature gradient in the upper crust steepened. In the second scenario, the termination of the Cenozoic volcanic lull c. 25 Ma ushered in the well-documented Huaylillas to present era of volcanic activity (Schildgen et al. 2007;Thouret et al. 2007;Mamani et al. 2009). This resumption of magmatism, which might be attributed to a return to steeper subduction, would have been promoted by wet melting of the overlying lithosphere. Magmatism was itself fuelled by an influx of melting asthenosphere driven by downward rotation of the subducting Nazca slab (James & Sacks 1999). Fluid circulation could have been responsible for heating the upper crust c. 20-15 Ma. The observed asymmetry in AHe age resetting on either side of the canyon could, furthermore, support the idea of anisotropic fluid circulation in the upper crust. This would explain the inability of the joint data models to simultaneously fit the younger (OCO05, OCO04-04, OCO11) and the older (OCO04-07, OCO04-09) AHe ages (Fig. 7c).
To summarize, given that thermochronology is sensitive only to temperature, the data cannot easily discriminate between the two most likely causes of crustal heating after 20 Ma: either the root cause of regional volcanism (a shift from lower-angle to higher-angle subduction c. 25 Ma) or the proximate cause attributable to rock burial beneath low-conductivity ignimbrites, or a combination of both. However, rapid rock cooling caused by canyon incision after 10 Ma soon competed with, and overwhelmed, the heating signal (Fig. 6).

Rise of the Western Cordillera and crustal response by flexural deformation
The thermochronological data clearly detect rapid canyon incison after 10 Ma, a process directly ascribable to rapid crustal uplift in the Central Andes, but cannot by themselves constrain palaeoaltimetric change in the study area. Thin, c. 24 Ma marine layers (Camaná Formation) preserved c. 2 km a.s.l. south of Caravelí (see Fig. 1) and interfingering with the upper Moquegua Formation (Huamán 1985;Kennan 2000;Cruzado & Rojas 2007) provide an indication of the magnitude of rock uplift since early Miocene time. The current elevation of those Oligocene to early Miocene marine rocks in the piedmont (Tosdal et al. 1984) and near Caravelí suggests that this low-energy environment was close to sea level between at least 40 and 24 Ma. It follows that low plateau erosion during Cenozoic time, which we interpreted as being due to low topographic gradients, can also be explained by the low elevation of the Cordillera above oceanic base level until at least 24 Ma. The area was thus a volcanically active, low-elevation land surface covered by onlapping forearc sediments, all of these components forming a low-relief, composite landscape: partly erosional, partly depositional, and partly eruptive after c. 25 Ma.
The cause of uplift after 24 Ma remains debated. Although the crust is 65-70 km thick (James 1971;Kono et al. 1989;Schmitz 1994;Beck et al. 1996;Masson et al. 2000) and despite prolonged tectonic convergence, crustal foreshortening in the Western Cordillera has been negligible. Instead, N808E strike-slip and normal faults prevail in the forearc, arc and western Altiplano of SW Peru (Pecho Gutierrez 1983;Olchauski & Dávila 1994;Sempere et al. 2008;Schildgen et al. 2009). Reports of compressional deformation in the Western Andes exist but invariably come from further south (usually 19-218S, in northern Chile), where interpretations concerning the nature, intensity, sense and surface expression of deformation in Precordilleran rocks have also differed significantly (compare Muñoz & Charrier 1996;Wörner et al. 2002;Victor et al. 2004;Elger et al. 2005;Farias et al. 2005;Zeilinger et al. 2005). Temporal patterns of deformation in the Central Andes between 46 Ma and the present indicate that deformation in the Western Cordillera after 7 Ma involved crustal shortening as well as oblique extensional to transtensional kinematics, both occurring synchronously in different locations (Oncken et al. 2006). This broadly agrees with both observed and modelled World Stress Map maximum horizontal compressional stresses (S H ): in the study area, these have led since 10 Ma to a normal faulting regime with S H oriented parallel to the strike of the Andes (Heidbach et al. 2008). Here we further suggest, in agreement with James & Sacks (1999), that limited compressional deformation of the crust in SW Peru can be explained by the strength of the Precambrian Arequipa block, which has survived magmatic infusion over several cycles of Andean subduction. This view is supported by the observation that, at this latitude, virtually all of the Andean deformational foreshortening has been confined to the eastern Altiplano and to the hotter and more ductile Eastern Cordillera and Sub-Andean fold-belt, beneath all of which the Arequipa block is absent.
An alternative explanation for crustal buoyancy in the 4-4.5 km high Western Cordillera is magmatic thickening (James & Sacks 1999;Mamani et al. 2009), whereby prolonged magmatic activity will have facilitated the topographic growth of the Precambrian basement (Arequipa craton). This interpretation is supported by balanced cross-sections across the Central Andes, based on which Schmitz (1994) modelled 300 km of shortening between the Western Cordillera and the Brazilian foreland. Comparing those results with seismic refraction data of Wigger et al. (1993), Schmitz concluded that whereas crustal thickening in the Altiplano and Eastern Cordillera could be explained by shortening and thrust stacking, the presence of 30 km of excess crustal material beneath the Western Cordillera could not be accounted for by this mechanism. Consistent with inferences originally made by James (1971), it was concluded that the 30 km of excess crust beneath the Western Cordillera could be attributed to magmatic addition. Detailed summaries of the magmatic succession of the Peruvian Andes (Mamani et al. 2009) clearly indicate that the Western Cordillera is an aggregation of seven Jurassic to Quaternary magmatic arcs, which can be approximately mapped on the basis of existing whole-rock and mineral-separate radiometric ages ( Fig. 1; see, e.g. Tosdal et al. 1981Tosdal et al. , 1984Beckinsale et al. 1985;Klinck et al. 1986;Mukasa 1986;Clark et al. 1990;Thouret et al. 2007). An indirect indication of occluded magma bodies at depth (as a consequence of magma ponding) is provided by the highly fractionated character of calc-alkaline lavas (relative to typical island arc subduction volcanism) throughout the history of Central Andean arc magmatism. Fractional crystallization thus appears to have been a recurring phenomenon in the study area, the surface area of exposed Andean batholith already indicating that the volumes of intrusive rocks are much greater than the volumes of their associated extrusive products (Francis & Rundle 1976).
Contributions to crustal thickness other than magmatic addition are plausible, and one of these is crustal thickening by lateral crustal flow. This has been previously proposed as a causal mechanism to explain the large monocline of eastern Tibet (Clark & Royden 2000), another vast, asymmetric orogenic plateau similar to the Altiplano. In the context of SW Peru, underthrusting of the Brazilian shield is thought to have caused shortening of the lower crust, with westward ingress of ductile crustal material resulting in excess mass on the western side of the Central Andes (Isacks 1988;Husson & Sempere 2003). The colder, rigid forearc would have acted as a passive indenter, thus allowing the monoclinal bulge portrayed in Figure 4 to form under the effect of hydraulic ramming by the inflow of ductile crust from the east. Crustal flow, however, would imply either hotter than normal crust (a situation contradicted by the data) or wetter than normal crust. Wet crust is thus the only condition compatible with our inference of crustal-scale fluid flow in the overriding South American plate. Such a process would have nevertheless contributed to only a fraction of the 30 km of excess crust detected by Schmitz (1994) in the Western Cordillera. This, therefore, still affords an overriding advantage to the magmatic addition proposition detailed above.
Another process conducive to monocline formation would be a strengthening of coupling between the subducting Nazca and the overriding South American plate (Sacks et al. 1978;. Two possible conditions could have increased the magnitude of resistive forces to shear stress at the plate interface. One is a deficit in trench lubrication along the sediment-starved continental margin (Lamb 2006), an interpretation partly supported by the limited Cenozoic denudation rates recorded in this study. The other is the excess gravitational potential energy of the 4 km high orogenic plateau, which rose rapidly after 10 Ma. Studies of the Nazca-South American plate convergence history indicate that not only did subduction rates peak between 25 and 20 Ma and steadily decline from 20 Ma to the present, but that a tangential strike-slip component prevalent during the earlier Cenozoic also decreased from 25 Ma onwards (Somoza 1998). As a result, the convergence angle became considerably more orthogonal to the South American plate margin. Finer detail shows that convergence rates decreased from 10.3 AE 0.2 cm a À1 c. 10 Ma ago to 6.7 AE 0.2 cm a À1 at the present day, and this 30% reduction in plate convergence coincided with the acceleration of topographic rise of the Andean orogen (Heidbach et al. 2008). This is also deemed sufficient to initiate the normal faulting regime reported in the Altiplano (Heidbach et al. 2008) and the Western Cordillera (Schildgen et al. 2009). Another consequence of this post-10 Ma stress distribution could be the formation, or at least the accentuation, of the monocline, which had perhaps already begun to form c. 25-20 Ma in response to the increase in orthogonality of convergence.

Conclusion
The thermal histories obtained from this multiple chronometer study of the Cotahuasi-Ocoña canyon provide an integrated history of crustal denudation since the igneous rocks, soon after crystallization, cooled through the zircon and apatite FT critical temperatures in late Cretaceous time (Figs 5 and 6). We infer that the Andean orogenic plateau in southwestern Peru began to rise at first after 24 Ma and again in late Miocene time, with uplift accelerating after 13-9 Ma. This is in agreement with increasing evidence obtained elsewhere that the uplift of the Central Andes increased both abruptly and recently (Garzione et al. 2008). Cooling histories also indicate low rates of background denudation during the earlier Cenozoic; that is, before the strong post-10 Ma cooling signal occurred as a consequence of canyon incision into the flank of a large monocline that is remarkably exempt of thrust tectonics. From the presence of granitic clasts in the Moquegua 'A' Formation in the forearc basin, we know that the shallow Andean batholith had been unroofed by Oligocene time at the latest.
The timing of Cordilleran plateau-edge deformation was probably triggered by steepening of the topographic slope facing the Pacific Ocean as a consequence of magmatic addition and crustal buckling. Crustal buckling was caused either by lower crustal flow, by frictional drag at the trench, or by the new balance between driving and resistive forces generated after 10 Ma in the forearc by the large, rising orogenic plateau. Given that canyon incision began c. 13 Ma and intensified after 9 Ma, we suggest that the giant Cordilleran monocline formed after 14 Ma. Flexural deformation may have occurred in several pulses, incrementally steepening the mountain front and driving successive waves of knickzone migration into the plateau (Thouret et al. 2007).
The data also capture some variation in the Cenozoic temperature gradient. One key aspect is that the geothermal gradient has remained extremely low since c. 80 Ma. Although this could reflect a lack of sensitivity in the annealing models, the values (3-8 8C km À1 ) are still within a range compatible with the lower bracket of cool forearc settings. Another important finding is that the temperature gradient steepened after 20-15 Ma. Given the magnitude of temperature change (heating by 40-60 8C) and its abrupt occurrence during the middle Miocene, we intepret this event either as a response to crustal heating after emplacement of the low-conductivity, regionally extensive Huaylillas ignimbritic caprock, or as a response to crustal-scale fluid circulation caused by wet melting of the overriding plate. Although the physical evidence remains indirect, we speculate, based on the arguments of James & Sacks (1999), that the observed change in heat flow would have been a response to change in subduction angle c. 25 Ma; that is, when the leading edge of the Nazca plate resumed a steeper gradient after a 25 Ma period of low-angle, non-volcanic subduction. Subduction of a seamount could also have switched off the magmatic system between c. 50 and c. 25 Ma without altering the subduction angle. However, as ridges and seamounts are more buoyant than the surrounding ocean floor, low-angle subduction of less than 458 is likely to be promoted by the gain in slab buoyancy (Wipf et al. 2008). Overall, either of the interpretations provided above is plausible, neither is geologically incompatible, and slab angle appears to be a key parameter.
To summarize, the data imply that, in subduction settings, (1) lower than normal geothermal gradients correlate with volcanic activity, with orthogonal plate convergence and, by association with existing theoretical models but also with present-day slab geometries (Fig. 1), with relatively steep subduction angles; and that (2) even lower geothermal gradients correlate with volcanic lulls and with greater tangential components to plate convergence. These latter conditions were prevalent prior to c. 25 Ma (Somoza 1998). The inferred chronology of change in geothermal gradients also indirectly supports the hypothesis that the subduction angle in Cordilleran settings can change over time in a given region, apparently abruptly. None of these correlations directly explain the driving mechanisms of change in magmatic activity. However, given that slab steepness controls opportunities for asthenospheric inflow into the region situated between the subducting slab and the overriding plate, the convergence angle may ultimately exert important controls over magmatic activity through controls on the slab angle.