# Spatial correlation bias in late-Cenozoic erosion histories derived from thermochronology

Taylor F. Schildgen, Pieter A. van Der Beek, Hugh D. Sinclair, Rasmus C. Thiede

Thiede

► To cite this version:

Taylor F. Schildgen, Pieter A. van Der Beek, Hugh D. Sinclair, Rasmus C. Thiede. Spatial correlation bias in late-Cenozoic erosion histories derived from thermochronology. Nature, 2018, 559, pp.89-93. 10.1038/s41586-018-0260-6. insu-03595905

HAL Id: insu-03595905

https://insu.hal.science/insu-03595905v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spatial correlation bias in late-Cenozoic erosion histories derived from thermochronology

Taylor F. Schildgen<fnref n="1" />,²,⁶*, Pieter A. van der Beek³,⁶, Hugh D. Sinclair⁴ & Rasmus C. Thiede²,⁵

The potential link between erosion rates at the Earth's surface and changes in global climate has intrigued geoscientists for decades<fnref n="1" />,² because such a coupling has implications for the influence of silicate weathering³,⁴ and organic-carbon burial⁵ on climate and for the role of Quaternary glaciations in landscape evolution<fnref n="1" />,⁶. A global increase in late-Cenozoic erosion rates in response to a cooling, more variable climate has been proposed on the basis of worldwide sedimentation rates⁷. Other studies have indicated, however, that global erosion rates may have remained steady, suggesting that the reported increases in sediment-accumulation rates are due to preservation biases, depositional hiatuses and varying measurement intervals⁸⁻¹⁰. More recently, a global compilation of thermochronology data has been used to infer a nearly twofold increase in the erosion rate in mountainous landscapes over late-Cenozoic times⁶. It has been contended that this result is free of the biases that affect sedimentary records¹¹, although others have argued that it contains biases related to how thermochronological data are averaged¹² and to erosion hiatuses in glaciated landscapes¹³. Here we investigate the 30 locations with reported accelerated erosion during the late Cenozoic⁶. Our analysis shows that in 23 of these locations, the reported increases are a result of a spatial correlation bias—that is, combining data with disparate exhumation

histories, thereby converting spatial erosion-rate variations into temporal increases. In four locations, the increases can be explained by changes in tectonic boundary conditions. In three cases, climatically induced accelerations are recorded, driven by localized glacial valley incision. Our findings suggest that thermochronology data currently have insufficient resolution to assess whether late-Cenozoic climate change affected erosion rates on a global scale. We suggest that a synthesis of local findings that include location-specific information may help to further investigate drivers of global erosion rates.

Thermochronology allows the reconstruction of the cooling histo-
ries of minerals, which in turn can be used to infer when and how
quickly rocks were exhumed to the surface¹⁴. The averaging timescale
of a thermochronological system is equivalent to its cooling age. That
age is (to first order) the time since a mineral has cooled below its
closure temperature and has begun to retain radiogenic daughter
products¹⁴. The depth of the closure-temperature isotherm depends
on both the geothermal gradient and the exhumation rate. Therefore,
the different closure temperatures and resulting closure depths asso-
ciated with various thermochronological systems imply that erosion
rates have been averaged over different length scales and timescales¹⁴.
To assess whether changes in erosion rates have occurred, multiple data

<figure><img src="image_2.png" /><figcaption>**Fig. 1** | Common scenarios in which we find steep spatial gradients in thermochronological ages and exhumation rates. a–f, For each tectono-geomorphic setting, the crustal structure is sketched as a three-dimensional block and the expected thermochronological age and exhumation rate patterns are qualitatively indicated with coloured circles. Dashed lines represent the position of a crustal isochron. Half arrows</figcaption></figure>

show relative motion across faults. Arrows in e indicate particle motion pathways. The orange arrow in f represents upwelling middle-lower crust. The scale bar indicates the characteristic length scale of the different structures. Typical case studies discussed either in the main text or in Supplementary Information are indicated in parentheses for each tectonic setting.

<footnote n="1">¹GFZ German Research Centre for Geosciences, Telegrafenberg, Potsdam, Germany. ²Department of Earth and Environmental Sciences, University of Potsdam, Potsdam, Germany. ³Université Grenoble Alpes, Institut des Sciences de la Terre (ISTerre), Grenoble, France. ⁴School of Geosciences, The University of Edinburgh, Edinburgh, UK. ⁵Present address: Institut für Geowissenschaften, Christian Albrecht Universität zu Kiel, Kiel, Germany. ⁶These authors contributed equally: Taylor F. Schildgen, Pieter A. van der Beek. *e-mail: tschild@gfz-potsdam.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** | Constructing erosion histories from multiple samples or multiple thermochronometers. The lower figure illustrates thermochronological data from a tilted normal-fault-bounded block and its hanging wall. Coloured symbols correspond to thermochronological ages, as in Fig. 1. Half arrows show relative motion across the fault. **a, c, e,** Schematic age–depth (or age–temperature; see Methods) plots for the data. **b, d, f,** Erosion histories inferred by averaging data with disparate exhumation histories—a potential problem highlighted in ref. 12 (**a, b**); combining data with disparate exhumation histories, leading to the spatial correlation bias (and spurious increases) described here (**c, d**); or combining only data with common exhumation histories, resulting in correct exhumation histories (**e, f**). T~c~, closure temperature.</figcaption></figure>

points with common exhumation histories must be combined. These data points may be derived either from multiple thermochronometric systems within a single sample or from multiple samples collected from steep altitudinal transects¹⁵. However, if the samples are separated by some horizontal distance, they are more likely to have differing histories¹⁶. Several common scenarios illustrated in Fig. 1 show how thermochronological ages and exhumation rates can differ over distances of a few kilometres in regions that are tectonically deformed or deeply incised by rivers or glaciers.

The inverse relationship between erosion rates and thermochronological ages will lead to an apparent increase in erosion rates through time if erosion rates are simply averaged¹² (Fig. 2a, b). To avoid this bias in interpreting a global compilation of thermochronological data, ref. ⁶ used a linear inversion scheme¹⁷ in which proximal samples were combined using a prescribed correlation length scale to derive exhumation histories (see below). However, the magnitude of the potential bias that results from combining data points that do not share a common exhumation history (for example, across a tectonic boundary), which we refer to as a 'spatial correlation bias', can be substantial. Combining rapidly exhumed samples, which tend to be younger, with slowly exhumed samples, which tend to be older, can lead to spurious accelerations in erosion rate for the analysed region (Fig. 2c, d).

To investigate the impact of spatial correlation biases on the reported acceleration in late-Cenozoic erosion⁶, we first describe how a spatial correlation function, which is used to combine data points in the linear inversion model of ref. ¹⁷, is defined. Next, we conduct a review of all localities that have been reported to show an increase in late-Cenozoic erosion rates using this technique. In our review of each locality, we summarize previous interpretations and assess whether or not samples

are likely to share a common exhumation history. We also evaluate the performance of the linear inversion model with a series of synthetic tests. Finally, we consider cases in which accelerated late-Cenozoic erosion appears to be unrelated to a spatial correlation bias and discuss whether climate change was the likely driver.

To construct erosion histories, the inversion model described in ref. ¹⁷ combines data points according to an assumed spatial correlation function, which gives greater weight to samples in close proximity to the location where an erosion history is determined. The spatial correlation function is based on the construction of an empirical semi-variogram, which shows the variance (mean squared difference) in erosion rates between pairs of samples, derived from calculating a steady-state erosion rate for each sample, as a function of their separation distance¹⁸. The empirical semi-variogram constructed for the European Alps in ref. ¹⁸, and recreated in Extended Data Fig. 1, approximates a negative exponential, with a steep increase in variance with increasing distance up to about 80 km (or an e-folding length, defined as the 'spatial correlation length', of about 25 km), after which the variance remains approximately constant. Complex spatial variations in erosion rates are expected in tectonically active areas and probably explain the high variance for the Alps data (Extended Data Fig. 1). Anisotropy in the data (that is, the data do not vary in a similar way in all directions) also contributes to the scatter. In datasets that contain anisotropy, different correlation functions are typically applied in different directions.

In the global analysis performed in ref. ⁶, the spatial correlation function used to combine and invert the data is an isotropic negative exponential with a uniform spatial correlation length of 30 km. In this model, data points close to a grid node are more strongly weighted when constructing an erosion history for the node's location; however, data points up to about 90 km away may still be used to derive an erosion history if no closer data points are available. Additional empirical semi-variograms were not reported to evaluate whether or not a negative exponential characterizes all parts of the global dataset, nor were the data evaluated for anisotropy.

A small subset (32 regions) of the global data compiled in ref. ⁶ produced erosion rates that were reported to be well resolved in the time bins from 6 to 4 million years (Myr) ago and from 2 to 0 Myr ago. The authors used these 32 regions (Extended Data Fig. 2) to assess global patterns in late-Cenozoic erosion-rate changes. Out of those regions, 27 produced erosion rates that are well resolved across each of four 2-Myr time bins from 8 to 0 Myr ago, which were in turn used to infer a nearly twofold late-Cenozoic increase in erosion rates across mountainous landscapes worldwide⁶. For each of the 32 regions with erosion histories that were reported to be well resolved (30 of which show accelerations), we conduct a literature review to compare the original data and interpretations to the erosion-rate histories reported in ref. ⁶. The regions include areas affected by normal faulting, transpression, thrust ramps, crustal pop-up structures and accretionary orogenic wedges. Other regions were affected by headward migration of river incision and glacial erosion. In each case, we explore whether the assumed spatial correlation function is justified, specifically considering whether the data points have a common exhumation history. Two of the regions that we reviewed are presented below; these were chosen because they represent an important part (21%) of the dataset that was used to infer a global increase in late-Cenozoic erosion rates⁶ and they are examples of common sources of the spatial correlation bias. Further discussion and reviews of all remaining regions can be found in Supplementary Information.

Our first case study comprises the western European Alps. The thermochronological database for this region defines two tectonic domains: the external zone, in which thermochronological ages are relatively young (<10 Myr to 15 Myr for apatite fission-track (AFT) ages; <20 Myr for zircon fission-track (ZFT) ages), and the internal zone, where AFT and ZFT ages are typically >20 Myr (Fig. 3a). The limit between these two zones corresponds to the Penninic frontal thrust (PFT), which coincides with a steep spatial gradient in thermochronological ages¹⁹⁻²¹. Interpretations of the data from the internal zone
<figure><img src="image_4.png" /><figcaption>**Fig. 3** | Thermochronology data and modelled erosion-rate changes for the western European Alps. **a**, Thermochronology data are derived from 52 different sources compiled by refs. ⁶,²⁵. Equivalent colours for thermochronometers correspond to equivalent one-dimensional steady-state erosion rates (see Methods for details). **b**, Normalized difference</figcaption></figure>

(see Methods) of erosion rates resolved from 2–0 Myr ago compared to those from 6–4 Myr ago, from the supplementary information of ref. ⁶. AHe, apatite (U–Th)/He system; MB, Mont Blanc massif; EP, Ecrins–Pelvoux massif; Rho, Rhône Valley; Rhi, Rhine Valley.

invoke steady, slow exhumation since early Miocene times, after a phase of rapid exhumation in the Oligocene²¹,²². By contrast, the external crystalline massifs in the external zone show rapid exhumation in late Miocene–early Pliocene times, followed by slower rates and finally localized exhumation of valley bottoms; the latter is attributed to glacial valley incision since the mid-Pleistocene²³.

In contrast to these interpretations, linear inversion models⁶,¹⁸,²⁴ show a broad band of erosion-rate increases in the western Alps during the Quaternary, which were attributed to a combination of slab breakoff and glacial erosion (Fig. 3b). An analogous earlier analysis reached similar conclusions²⁵. However, the band of increased erosion rates straddles the PFT (Fig. 3b), which marks the boundary between two blocks that have experienced contrasting exhumation histories. The inferred erosion-rate increases occur throughout the internal zone, east of the PFT, where no thermochronological ages younger than 4 Myr occur. Thus, the data cannot resolve any variations in erosion rates later than 4 Myr ago. Our synthetic tests show that the reported increases in erosion rates result from a combination of data across the PFT (Extended Data Fig. 3). When we impose different exhumation rates across this tectonic boundary but hold them constant through time, spurious increases are produced by the linear inversion in both the external and the internal Alps (see Methods, Extended Data Fig. 3).

In our second case study, the Mount Cook region of the Southern Alps in New Zealand (Fig. 4, Extended Data Fig. 4), oblique convergence between the Australian and Pacific plates is concentrated along the Alpine fault. This region has long been considered as a type example of a doubly vergent thrust wedge in which exhumation is concentrated in the retro-wedge owing to both the tectonic and the erosional asymmetry of the system²⁶,²⁷. Rocks are exhumed from progressively greater depths towards the Alpine fault. As a result, metamorphic grade increases and thermochronological ages decrease systematically towards the fault²⁸,²⁹. Each thermochronological system shows young (reset) ages close to the Alpine fault; the width of the reset age zone depends on the dip of the fault and the closure temperature of the

system (Fig. 4a–d). Thermo-mechanical²⁷ and thermo-kinematic²⁹,³⁰ modelling has shown that the age pattern can be explained by steady convergence at rates of around 10 mm yr⁻¹ since about 5–6 Myr ago²⁸, consistent with plate-tectonic reconstructions indicating a switch from pure strike-slip to oblique convergence along the Alpine fault at that time³¹.

In contrast to these earlier findings, the analysis reported in ref. ⁶ suggests a 20-fold increase in late-Cenozoic erosion rates along the crest of the Southern Alps. The increases are only resolved in a zone between about 20 km and 30 km southeast of the Alpine fault. Closer to the fault, all ages are considerably younger than 6 Myr, leading to loss of resolution for the older age bins, whereas farther away from the fault, all ages are substantially older than 6 Myr, leading to loss of resolution in the younger age bins. The inferred increases in erosion rates are due to the systematic spatial variations in ages, which reflect the laterally varying exhumation pathways, but are translated by the linear inversion model into temporal increases in erosion rates, as illustrated in another synthetic test (see Methods, Extended Data Fig. 5).

Quantifying how much of the proposed worldwide increase in late-Cenozoic mountain erosion rates results from a true increase recorded by samples that share a common exhumation history, rather than from a bias linked to the combination of samples with disparate exhumation histories, is a challenge. We have attempted to do so qualitatively by summarizing our conclusions from all of the areas for which increases in late-Cenozoic erosion rates were inferred⁶. We find that the reported increases can be explained by the combination of data with disparate exhumation histories in 23 cases. Increases in four cases can be explained by changes in tectonic boundary conditions, two are related to both changes in tectonics and glacial activity, and one appears to be associated with glacial valley incision alone (Extended Data Table 1, Extended Data Fig. 2). Hence, only in three regions do thermochronology data reflect an impact of late-Cenozoic cooling.

Theoretically, regional changes in exhumation rates could be correctly resolved using linear modelling approaches in regions that have
<figure><img src="image_5.png" /><figcaption>**Fig. 4** | Thermochronology data and modelled erosion-rate changes for the Mount Cook region in New Zealand. **a–d**, Thermochronology data from refs 28,33 (AFT, ZFT), ref. 29 (ZHe) and ref. 30 (AHe, AFT, ZHe, ZFT). Equivalent colours for thermochronometers correspond to equivalent one-dimensional steady-state erosion rates. Coloured lines in</figcaption></figure>

a–c and e mark the approximate boundaries between reset and un-reset thermochronometer ages (red, ZFT; orange, ZHe; yellow, AFT). AHe, apatite (U–Th)/He system; ZHe, zircon (U–Th)/He system. e, Normalized difference of erosion rates resolved from 2–0 Myr ago compared to the those from 6–4 Myr ago, from the supplementary information of ref. 6.

not experienced changes in short-wavelength relief, with data available from either (1) multiple thermochronometers analysed from the same sample in an area of mostly vertical exhumation or (2) steep elevation transects of samples that do not cross spatial gradients in exhumation¹⁷. Reducing the spatial correlation length in the linear inversion model results in fewer cases of unwarranted data combination, as illustrated in our synthetic tests (see Methods and Extended Data Figs. 3, 5, 6). However, as the correlation length is decreased, fewer areas yield well-resolved erosion histories. More problematically, our synthetic tests show that the spatial correlation bias can create spurious increases even in areas of high data density, multiple thermochronometers available from individual samples and steep sample transects (Extended Data Fig. 3).

On the basis of our literature review and synthetic tests, we find that thermochronology data do not resolve a worldwide late-Cenozoic increase in mountain erosion rates. Both the magnitude of the inferred increase and its purported occurrence throughout mountainous landscapes appear to be largely artefacts related to combining data with disparate exhumation histories. These results, combined with previous analyses of bias in the sedimentary record⁸⁻¹⁰, call into question the evidence presented to date for a worldwide late-Cenozoic increase in erosion rates.

Considering the pitfalls of a global approach to data analysis that neglects the local context of the data, an approach that synthesizes detailed local work and includes location-specific aspects in data analysis is probably our best way forward. Such an approach, as summarized here, shows only a few locations in which resolved increases in late-Cenozoic erosion rates are unambiguously linked to late-Cenozoic climate change, and in those areas, such increases are mostly restricted to localized glacial valley incision²³,³². Hence, even in cases where data do resolve late-Cenozoic changes in landscape relief, the importance of such changes with respect to global sediment and carbon budgets must be carefully considered.

## Online content

Any Methods, including any statements of data availability and Nature Research reporting summaries, along with any additional references and Source Data files, are available in the online version of the paper at https://doi.org/10.1038/s41586-018-0260-6.

Received: 15 October 2017; Accepted: 8 May 2018; Published online 4 July 2018.

1. Molnar, P. & England, P. Late Cenozoic uplift of mountain ranges and global climate change: chicken or egg? *Nature* **346**, 29–34 (1990).
2. Raymo, M. E. & Ruddiman, W. F. Tectonic forcing of late Cenozoic climate. *Nature* **359**, 117–122 (1992).

3. Berner, R. A., Lasaga, A. C. & Garrels, R. M. The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years. *Am. J. Sci.* **283**, 641–683 (1983).

4. Kump, L. R., Brantley, S. L. & Arthur, M. A. Chemical weathering, atmospheric CO₂, and climate. *Annu. Rev. Earth Planet. Sci.* **28**, 611–667 (2000).

5. France-Lanord, C. & Derry, L. A. Organic carbon burial forcing of the carbon cycle from Himalayan erosion. *Nature* **390**, 65–67 (1997).

6. Herman, F. et al. Worldwide acceleration of mountain erosion under a cooling climate. *Nature* **504**, 423–426 (2013).

7. Zhang, P., Molnar, P. & Downs, W. R. Increased sedimentation rates and grain sizes 2–4 Myr ago due to the influence of climate change on erosion rates. *Nature* **410**, 891–897 (2001).

8. Sadler, P. M. The influence of hiatuses on sediment accumulation rates. *Geores. Forum* **5**, 15–40 (1999).

9. Schumer, R. & Jerolmack, D. J. Real and apparent changes in sediment deposition rates through time. *J. Geophys. Res.* **114**, F00A06 (2009).

10. Willenbring, J. K. & von Blanckenburg, F. Long-term stability of global erosion rates and weathering during late-Cenozoic cooling. *Nature* **465**, 211–214 (2010).

11. Herman, F. & Champagnac, J.-D. Plio-Pleistocene increase of erosion rates in mountain belts in response to climate change. *Terra Nova* **28**, 2–10 (2016).

12. Willenbring, J. K. & Jerolmack, D. J. The null hypothesis: globally steady rates of erosion, weathering fluxes and shelf sediment accumulation during Late Cenozoic mountain uplift and glaciation. *Terra Nova* **28**, 11–18 (2016).

13. Ganti, V. et al. Time scale bias in erosion rates of glaciated landscapes. *Sci. Adv.* **2**, e1600204 (2016).

14. Reiners, P. W. & Brandon, M. T. Using thermochronology to understand orogenic erosion. *Annu. Rev. Earth Planet. Sci.* **34**, 419–466 (2006).

15. Hurford, A. J. Uplift and cooling pathways derived from fission track analysis and mica dating: a review. *Geol. Rundsch.* **80**, 349–368 (1991).

16. Stephenson, J., Gallagher, K. & Holmes, C. C. Low temperature thermochronology and strategies for multiple samples 2: partition modelling for 2D/3D distributions with discontinuities. *Earth Planet. Sci. Lett.* **241**, 557–570 (2006).

17. Fox, M., Herman, F., Willett, S. D. & May, D. A. A linear inversion method to infer exhumation rates in space and time from thermochronometric data. *Earth Surf. Dynam.* **2**, 47–65 (2014).

18. Fox, M., Herman, F., Willett, S. D. & Schmid, S. M. The exhumation history of the European Alps inferred from linear inversion of thermochronometric data. *Am. J. Sci.* **316**, 505–541 (2016).

19. Seward, D. & Mancktelow, N. S. Neogene kinematics of the central and western Alps: evidence from fission-track dating. *Geology* **22**, 803–806 (1994).

20. Fügenschuh, B. & Schmid, S. M. Late stages of deformation and exhumation of an orogen constrained by fission-track data: a case study in the Western Alps. *Geol. Soc. Am. Bull.* **115**, 1425–1440 (2003).

21. Tricart, P., van der Beek, P., Schwartz, S. & Labrin, E. Diachronous late-stage exhumation across the western Alpine arc: constraints from apatite fission-track thermochronology between the Pelvoux and Dora-Maira Massifs. *J. Geol. Soc. Lond.* **164**, 163–174 (2007).

22. Beucher, R., van der Beek, P., Braun, J. & Batt, G. E. Exhumation and relief development in the Pelvoux and Dora-Maira massifs (western Alps) assessed by spectral analysis and inversion of thermochronological age transects. *J. Geophys. Res.* **117**, F03030 (2012).

23. Valla, P. G., Shuster, D. L. & van der Beek, P. A. Significant increase in relief of the European Alps during mid-Pleistocene glaciations. *Nat. Geosci.* **4**, 688–692 (2011).

24. Fox, M., Herman, F., Kissling, E. & Willett, S. D. Rapid exhumation in the Western Alps driven by slab detachment and glacial erosion. *Geology* **43**, 379–382 (2015).

25. Vernon, A. J., van der Beek, P. A., Sinclair, H. D. & Rahn, M. K. Increase in late Neogene denudation of the European Alps confirmed by analysis of a fission-track thermochronology database. *Earth Planet. Sci. Lett.* **270**, 316–329 (2008).

26. Beaumont, C., Kamp, P. J. J., Hamilton, J. & Fullsack, P. The continental collision zone, South Island, New Zealand: comparison of geodynamical models and observations. *J. Geophys. Res.* **101**, 3333–3359 (1996).

27. Batt, G. E. & Braun, J. The tectonic evolution of the Southern Alps, New Zealand: insights from fully thermally coupled dynamical modelling. *Geophys. J. Int.* **136**, 403–420 (1999).

28. Batt, G. E., Braun, J., Kohn, B. P. & McDougall, I. Thermochronological analysis of the dynamics of the Southern Alps, New Zealand. *Geol. Soc. Am. Bull.* **112**, 250–266 (2000).

29. Herman, F., Braun, J. & Dunlap, W. J. Tectonomorphic scenarios in the Southern Alps of New Zealand. *J. Geophys. Res.* **112**, B04201 (2007).

30. Herman, F., Cox, S. C. & Kamp, P. J. J. Low-temperature thermochronology and thermokinematic modeling of deformation, exhumation, and development of topography in the central Southern Alps, New Zealand. *Tectonics* **28**, TC5011 (2009).

31. Walcott, R. I. Modes of oblique compression: Late Cenozoic tectonics of the South Island of New Zealand. *Rev. Geophys.* **36**, 1–26 (1998).

32. Ehlers, T. A., Farley, K. A., Rusmore, M. E. & Woodsworth, G. J. Apatite (U-Th)/He signal of large-magnitude accelerated glacial erosion, southwest British Columbia. *Geology* **34**, 765–768 (2006).

33. Tippett, J. M. & Kamp, P. J. J. Fission track analysis of the Late Cenozoic vertical kinematics of continental Pacific crust, South Island, New Zealand. *J. Geophys. Res.* **98**, 16119–16148 (1993).

**Acknowledgements** We thank P. Valla for discussions on the topic. D. Jerolmack and A. Ault provided feedback on an early version of the manuscript. T.F.S. acknowledges support from the Emmy Noether Program of the Deutsche Forschungsgemeinschaft (DFG), grant number SCHI 1241/1-1, and the Helmholtz Association's Initiative and Networking Fund. P.A.v.d.B acknowledges support from the Institut Universitaire de France (IUF). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). R.C.T. was funded by the DFG grant number TH 1371/5-1.

**Reviewer information** *Nature* thanks R. Flowers, E. Kirby and J. Spotila for their contribution to the peer review of this work.

**Author contributions** T.F.S., P.A.v.d.B. and H.D.S. discussed early versions of the ideas presented here, which were refined following discussion among all the authors. T.F.S., P.A.v.d.B. and R.C.T. compiled data for the figures. T.F.S. and P.A.v.d.B. created all figures and constructed the *edot2age.m* Matlab script. P.A.v.d.B. carried out the synthetic model runs with the help of T.F.S. T.F.S., P.A.v.d.B. and H.D.S. analysed individual case studies. All authors discussed the results and contributed to the writing of the text.

**Competing interests** The authors declare no competing interests.

**Additional information**

*Extended data* is available for this paper at https://doi.org/10.1038/s41586-018-0260-6.

*Supplementary information* is available for this paper at https://doi.org/10.1038/s41586-018-0260-6.

*Reprints and permissions information* is available at http://www.nature.com/reprints.

*Correspondence and requests for materials* should be addressed to T.F.S.

*Publisher's note:* Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
# METHODS

## Obtaining normalized-difference maps of erosion rates from the global analysis.

The full results of ref. ⁶ are provided in a 207-page data table in the supplementary information of that paper. That table provides the longitude, latitude, erosion rate and resolution for each of the grid points in the model⁶, divided into time bins (that is, in 2-Myr increments between 8 Myr ago and the present). The analysis of ref. ⁶ employed ratios of erosion rates from 2–0 Myr ago and compared to 6–4 Myr ago to assess global patterns of erosion-rate changes. A disadvantage of ratios (final rate divided by initial rate) is that the values increase substantially for very slow initial erosion rates. To avoid this problem, we instead calculate normalized differences (ND) between the final erosion rates ($e_f$ from 2–0 Myr ago) and initial erosion rates ($e_i$, from 6–4 Myr ago):

$$ND = \frac{e_f - e_i}{\max(e_f, e_i)} \qquad (1)$$

The difference between the final and initial erosion rate is divided by the maximum of the two rates. Thus, increases are scaled between 0 and 1 and decreases are scaled between 0 and -1. These values track fractional changes in erosion rates, so an increase in erosion rates from 0.5 mm yr⁻¹ to 1 mm yr⁻¹ yields the same normalized difference as an increase from 0.05 mm yr⁻¹ to 0.1 mm yr⁻¹. Increases and decreases are also symmetric, so an increase from 0.05 mm yr⁻¹ to 0.1 mm yr⁻¹ has a normalized difference of 0.5, whereas a decrease from 0.1 mm yr⁻¹ to 0.05 mm yr⁻¹ has a normalized difference of -0.5. Calculating these values for each 'resolved' location (following the definition in ref. ⁶) was automated using a Matlab script. The data were then exported to ArcMap GIS (Esri) to be plotted in map view. Such maps were not shown in ref. ⁶, but they are critical for assessing the quality of the model results; all the data needed to make those maps can be found in the supplementary information of ref. ⁶.

### Finding the steady-state thermochronological age corresponding to a steady erosion rate.
When constructing maps for each case study, we coloured each thermochronological data point according to its corresponding one-dimensional steady-state erosion rate. We emphasize that these rates are calculated only to facilitate visual comparison of data from different thermochronometers; the rates do not reflect an attempt to reconstruct detailed exhumation histories. If the data points of multiple thermochronometers from a small region have the same colour, the data are consistent with minimal changes in exhumation rates through time. Hence, this approach provides a simple way to visually assess mean, time-averaged exhumation rates from multiple thermochronological systems. However, because our predicted erosion rates do not consider sample elevation, any comparison of samples should be among those at a similar elevation.

In detail, we determined what age corresponds to a given exhumation rate for each thermochronometer (Extended Data Table 2) and coloured data in each map according to a range of specified rates. Changing the parameters in this model will modify the absolute rates, but not the relative rates, which are of interest here. In the following, we describe the steps used to make the conversion from erosion rate to age on the basis of the methodology developed in refs ³⁴,³⁵, which are encoded in the Matlab script edot2age.m (available in Supplementary Information).

The age of a thermochronological system is, to the first order, equal to the time since that system crossed its closure-temperature isotherm as it was exhumed to the surface. The closure temperature ($T_c$) of a monotonously cooling thermochronological system can be modelled as³⁶:

$$T_c = \frac{E_a}{R \ln(A \tau D_0 / a^2)} \qquad (2)$$

where R is the universal gas constant, A is a geometry factor, $E_a$, $D_0$ and a are experimentally determined diffusion parameters and $\tau$ is a characteristic time that depends on the cooling rate, $\delta T/\delta t$, where t is the time:

$$\tau = - \frac{R T_c^2}{E_a (\delta T / \delta t)} \qquad (3)$$

We calculate the closure temperatures for the apatite and zircon (U–Th)/He systems (AHe and ZHe, respectively) and for the AFT and ZFT systems as a function of cooling rate using diffusion parameters given in ref. ¹⁴ (see Extended Data Table 3). We note that although AFT and ZFT annealing cannot strictly be treated as linear Arrhenius diffusion³⁷, the above equations provide a good first-order approximation for this process as well, as discussed in ref. ¹⁴.

To transform the closure temperature into a closure depth, the geothermal gradient needs to be known. The equation describing the temperature, $T$, (in kelvin) as a function of depth, $z$, (in kilometres) in a one-dimensional system including vertical rock advection (but no heat production) is³⁵,³⁸.

$$T(z) = T_0 + (T_L - T_0) \frac{1 - e^{-z e / \kappa}}{1 - e^{-L e / \kappa}} \qquad (4)$$

where L (km) is the 'characteristic length' of the system (in our case, the crustal thickness), $\epsilon$ (km yr⁻¹) is the vertical exhumation rate, $\kappa$ (km² yr⁻¹) is the thermal diffusivity, and $T_0$ and $T_L$ (K) are the temperatures at the top and the base of the system (at depths 0 and L), respectively.

We can rearrange this equation to give the depth as a function of temperature (that is, to derive a closure depth from a closure temperature):

$$z(T) = \frac{\kappa}{\epsilon} \ln \left[ 1 - \frac{T - T_0}{T_L - T_0} \left( 1 - e^{-L e / \kappa} \right) \right] \qquad (5)$$

We take the derivative of equation (4) with respect to depth to obtain the geothermal gradient:

$$\frac{dT}{dz} = \frac{\epsilon(T_L - T_0)}{\kappa(1 - e^{-L e / \kappa})} e^{-z e / \kappa} \qquad (6)$$

Equations (2)–(6) are solved iteratively. The geothermal gradient is used together with the input exhumation rate to calculate the cooling rate at the closure depth and thereby the closure temperature. The iteration quickly converges to a steady-state geothermal gradient, closure temperature and depth, and thermochronologic age as a function of exhumation rate.

### Synthetic tests of the linear inversion model.
We conducted synthetic tests to evaluate whether the spatial correlation bias may exist in various tectonic settings and to test the impact of different values for the spatial correlation length and the a priori erosion rate used in the model inversion. To create synthetic datasets, we extracted elevation data and sample locations from three different areas: the western European Alps, the Mount Cook area of New Zealand and the Wasatch Mountains. We predicted the ages at the sample locations by running forward models with the three-dimensional thermal-kinematic code Pecube³⁹, assuming spatially variable but temporally constant exhumation rates. In these forward models, we used exhumation rates that produce a range of ages similar to what is found in the real data in each area. We added a random 10% error to these ages to obtain a more natural spread and then applied the linear inversion model of ref. ¹⁷ to the predicted ages. For both the thermal–kinematic forward model and the linear inversion, we used the same thermal and thermochronometric parameters as in ref. ⁶; these are summarized in Extended Data Table 3. For each synthetic test, we show maps (Extended Data Figs. 3, 5, 6) of the input erosion rates (constant in time) and the erosion rates predicted by the inversion model from 6–4 Myr ago and from 2–0 Myr ago. Model-predicted maps include contour lines showing the resolution. We also present maps showing the normalized difference between the erosion rates of the two time windows. Importantly, because the input exhumation rates are constant in time for all our model runs, any predicted change between the erosion rates of different time bins is an artefact of the linear inversion model.

For the synthetic test of the western European Alps, we extracted the topography and sample locations from a portion of the range spanning the northeast–southwest striking, sub-linear segment of the Penninic frontal thrust (PFT) between the Ecrins–Pelvoux massif and the Rhône valley (Fig. 3). We modelled the PFT as a vertical tectonic boundary separating two regions of contrasting but temporally constant exhumation rates, and we set the exhumation rates for the external and internal zones to be 1.0 mm yr⁻¹ and 0.25 mm yr⁻¹, respectively (Extended Data Fig. 3a). We ran the model at these rates for 20 Myr before predicting ages for various thermochronometers (AHe, AFT and ZFT) at the sample locations for which real data are available. Our choice to use real data locations reflects our aim to perform the linear inversion with a realistic distribution of data points. Predicted AHe ages (18 in total) range from 2.0 Myr to 3.3 Myr in the external zone; a single age in the internal zone is 8.6 Myr. Predicted AFT ages (176 in total) range from 2.5 Myr to 5.6 Myr in the external zone and from 11.3 Myr to 22.5 Myr in the internal zone. Predicted ZFT ages (114 in total) range from 6.4 Myr to 8.8 Myr in the external zone; in the internal zone, the ZFT ages are unreset (that is, they have not been exhumed deep enough to yield ages representative of the imposed exhumation rate) and they range from 32.5 Myr to 34.5 Myr (Extended Data Fig. 3b). Using a spatial correlation length of 30 km (as used in ref. ⁶), the linear inversion of the data predicts 'resolved' increases in exhumation rates through time across an approximately 110-km-wide zone (using the resolution cut-off of 0.25, as in ref. ⁶) that spans both sides of the PFT; that is, both the internal and the external zones (Extended Data Fig. 3e). Using a correlation length of only 10 km, the inversion results in a narrower (about 30-km-wide) zone of resolved increases, but those increases still span both the internal and external zones (Extended Data Fig. 3h). The decreased area of resolved changes for the 10 km correlation length occurs because fewer samples are combined when determining exhumation histories. This synthetic test illustrates the strength of the spatial correlation bias, even in regions
with a high data density, multiple thermochronometers available for individual samples, and steep age-elevation transects.

For the New Zealand synthetic test, we extracted topographic data and sample locations from the Mount Cook area (Fig. 4). We predicted exhumation rates from a model that assumes overthrusting for 10 Myr at a constant rate of 5 mm yr⁻¹ along a fault with a geometry taken from ref. ³⁰; dipping at 45° from the surface to a depth of 15 km and at 9.5° below that depth, before soling out in a horizontal detachment at a depth of 20 km. Predicted exhumation rates associated with this setup vary stepwise, from 2.2 mm yr⁻¹ adjacent to the fault to zero about 47 km away from it (Extended Data Fig. 5a). From this model, we predict 13 AHe ages between 1.0 Myr and 5.4 Myr, 50 reset AFT ages between 0.8 Myr and 5.4 Myr, 19 reset and 3 partially reset ZHe ages between 1.7 Myr and 7.7 Myr and between 10 Myr and 34 Myr, respectively, and 30 reset ZFT ages between 1.8 Myr and 5.8 Myr. In total, 7 AFT, 5 ZHe and 33 ZFT ages are predicted to be unreset and vary from 80 Myr to 100 Myr (Extended Data Fig. 5b). The linear inversion of the data using a 30 km spatial correlation length predicts very low rates (close to the a priori erosion rate of 0.35 mm yr⁻¹) from 6-4 Myr ago throughout the model domain, with a predicted resolution greater than 0.25 in a large part of it (Extended Data Fig. 5c). From 2-0 Myr ago, the model predicts rates that increase towards the fault (Extended Data Fig. 5d); as a result, large and resolved increases (normalized differences ranging between 0.7 and 0.9) are predicted over most of the model domain (Extended Data Fig. 5e). Choosing a shorter correlation length of 10 km leads to similar results, but with lower resolution, in particular from 6-4 Myr ago, leading to predicted resolved increases over a smaller part of the model domain (Extended Data Fig. 5h).

For the Wasatch synthetic test, we extracted the topography and sample locations from the portion of the mountain range for which thermochronological data are available (Supplementary Fig. 1a). We set the exhumation rate on the left side of the model (along the Wasatch fault) to 1.2 mm yr⁻¹ and on the right side of the model to zero, creating a smooth spatial variation in exhumation rates across the model domain. We ran Pecube at these rates for 12 Myr before predicting AHe and AFT ages at the sample locations where real data are available (Extended Data Fig. 6). Twenty-seven predicted AHe ages range from 1.9 Myr to 7.3 Myr, 33 reset AFT ages range from 2.7 Myr to 10 Myr, whereas 4 unreset AFT ages in the east are between 25 Myr and 30 Myr. Linear inversions of these synthetic data produce a spatial gradient in exhumation rates that is damped compared to the input gradient. The linear inversion of the data, using a 30-km spatial correlation length and an a priori erosion rate of 0.35 ± 0.1 mm yr⁻¹, predicts resolved moderate increases in exhumation rates throughout the model domain, with normalized differences varying from 0.08 in the east to 0.46 in the west (Extended Data Fig. 6e). As in the previous synthetic tests, reducing the correlation length to 10 km does not strongly change the pattern but reduces the resolution, in particular from 6-4 Myr ago (Extended Data Fig. 6f-h). When the a priori erosion rate is reduced to 0.1 ± 0.1 mm yr⁻¹, the model predicts strong decreases (normalized difference up to -0.42) in the east and strong increases (up to 0.61) in the west (Extended Data Fig. 6k). These synthetic tests illustrate again how samples with multiple thermochronometers do not counteract the spatial correlation bias. Moreover, this example shows the dependence of the linear inversion results on the assumed a priori erosion rate, as the model tends towards these rates where the resolution is relatively low.

The linear inversion model found spurious erosion-rate increases for all of the synthetic tests, which were assigned spatially variable but temporally constant exhumation rates. Overall, shorter spatial correlation lengths restrict the area for

which spurious increases were interpreted, but they do not eliminate spurious increases, even in the synthetic test of the western Alps, where data are dense, multiple thermochronometer ages are available for several individual samples, and numerous steep elevation transects characterize the dataset.

Our synthetic tests furthermore illustrate that choosing a higher resolution cut-off (higher than 0.25) reduces the regions for which spurious increases are resolved; however, even at resolutions of 0.5-0.6, spurious increases occur. On a global scale, increasing the resolution cut-off to 0.5 would result in only seven locations in the world with resolved late-Cenozoic increases in exhumation rates (the western European Alps, the Apennines, Taiwan, Fiordland, a single point in the Wasatch range and a single point in the central Himalaya), whereas a cut-off of 0.6 reduces the number of locations to three (the western Alps, the Apennines and Taiwan). Nevertheless, as we argued earlier and in Supplementary Information, most of those increases are spurious.

We also find that the magnitude of the acceleration is dependent on the chosen a priori erosion rate. In their discussion of the sensitivity of the inversion results to the a priori erosion rate, the authors of ref. ⁶ note that “choosing a prior erosion rate that is substantially different from the actual solution will lead to a wrong solution.” However, they used a constant and relatively low (0.35 ± 0.1 mm yr⁻¹) a priori erosion rate for all ‘tectonically active’ regions, potentially further exacerbating the inferred recent increase in erosion rates.

**Data availability.** Data for the European Alps and New Zealand are based on the global thermochronology data compilation published in ref. ⁶. All other thermochronology data shown were compiled by the authors, but can also be obtained from the supplementary information of ref. ⁶. All shaded relief maps (for example, in Figs. 3 and 4 and in Extended Data Figs. 3, 5 and 6) were created in ArcGIS 10.3.1 from 3-arcsecond Shuttle Radar Topography Mission (SRTM) digital elevation data available from the US Geological Survey. Input files for Pecube that were used to predict data for our synthetic tests and input files for the linear inversion model are available at https://github.com/TaylorSchildgen/SpatialCorrelationBias. The linear inversion model is available at https://zenodo.org/record/1215685.

34. Brandon, M. T., Roden-Tice, M. K. & Garver, J. I. Late Cenozoic exhumation of the Cascadia accretionary wedge in the Olympic Mountains, northwest Washington State. *Geol. Soc. Am. Bull.* **110**, 985–1009 (1998).

35. Braun, J., van der Beek, P. & Batt, G. *Quantitative Thermochronology: Numerical Methods for the Interpretation of Thermochronological Data* (Cambridge Univ. Press, Cambridge, 2006).

36. Dodson, M. H. Closure temperature in cooling geochronological and petrological systems. *Contrib. Mineral. Petrol.* **40**, 259–274 (1973).

37. Laslett, G. M., Green, P. F. & Duddy, I. R. Thermal annealing of fission tracks in apatite 2. A quantitative analysis. *Chem. Geol. Isot. Geosci. Sect.* **65**, 1–13 (1987).

38. Mancktelow, N. S. & Grasemann, B. Time-dependent effects of heat advection and topography on cooling histories during erosion. *Tectonophysics* **270**, 167–195 (1997).

39. Braun, J. et al. Quantifying rates of landscape evolution and tectonic processes by thermochronology and numerical modeling of crustal heat transport using PECUBE. *Tectonophysics* **524–525**, 1–28 (2012).

40. Cressie, N. *Statistics for Spatial Data* (Wiley Interscience, Hoboken, 1993).

41. Farley, K. A. Helium diffusion from apatite: general behavior as illustrated by Durango fluorapatite. *J. Geophys. Res.* **105**, 2903–2914 (2000).

42. Reiners, P. W., Spell, T. L., Nicolescu, S. & Zanetti, K. A. Zircon (U-Th)/He thermochronometry: He diffusion and comparisons with <sup>40</sup>Ar/<sup>39</sup>Ar dating. *Geochim. Cosmochim. Acta* **68**, 1857–1887 (2004).
<figure><img src="image_6.png" /></figure>

**Extended Data Fig. 1** | Empirical semi-variogram of erosion rates derived from AFT data from the European Alps. a, Gamma plot, showing the mean of the squared differences (variance) in erosion rates for each distance bin. The plot shows a negative exponential pattern, as described in ref. ¹⁸. The variance at a distance of zero, defined as the ‘nugget’, is non-negligible. The nugget is considered to represent the measurement error plus the variance at distances smaller than the smallest measurable distance⁴⁰. The ‘sill’ is the relatively steady variance value

reached at distances beyond the initial increase⁴⁰. The nugget-to-sill ratio, which in this case is about 0.35, indicates that approximately 35% of the variance is unaccounted for with the spatial correlation function⁴⁰.

b, Box plot of the range of squared differences in erosion rates, with boxes outlining quartiles of data, and error bars corresponding to 1.5 times the inner quartile distance. Horizontal lines within boxes show medians of data within each distance bin.
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 2 | Global coverage and inferred cause of resolved changes in late-Cenozoic erosion rates.** The data are based on the inversion results of ref. ⁶. Locations where erosion rates are resolved both from 2–0 Myr ago and from 6–4 Myr ago are shown with thin black circles,

whereas locations where erosion rates are resolved in all four 2-Myr time bins between 8 Myr ago and the present are shown with thick black circles. Circles filled with multiple colours indicate contributions from multiple factors.
<figure><img src="image_8.png" /></figure>

**Extended Data Fig. 3** | Synthetic test for the western Alps. **a**, Input exhumation rates are 1 mm yr⁻¹ in the external zone and 0.25 mm yr⁻¹ in the internal zone; they are held constant for 20 Myr. **b**, Predicted thermochronological ages; different symbols for thermochronological systems and colour scale for ages are as in Fig. 3a. **c-e**, Predicted erosion rates for a correlation length (λ) of 30 km from 6-4 Myr ago (**c**), from

2-0 Myr ago (**d**), and normalized difference in erosion rates from the comparison of **c** to **d** (as defined in equation (1)), shown only where the resolution in each time bin is greater than 0.25 (e). Contours in **c-e** show the predicted resolution. **f-h**, Same as **c-e**, but for a correlation length of 10 km. Ma, million years ago.
<figure><img src="image_9.png" /><figcaption>**Extended Data Fig. 4** | Overview map of the South Island, New Zealand.</figcaption></figure>

Circles are locations where erosion rates are resolved (resolution >0.25) from both 2–0 Myr ago and from 6–4 Myr ago from the extended data

of ref. ⁶. Black boxes outline the regions described in the case studies of Mount Cook (Fig. 4), Fiordland (Supplementary Fig. 8) and Marlborough fault system (Supplementary Fig. 16).
<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 5** | Synthetic test for the Mount Cook region in New Zealand. **a**, Input exhumation rates are predicted from a kinematic model of overthrusting along a fault with a ramp-flat geometry³⁰; rates vary from 2.2 mm yr⁻¹ adjacent to the fault to zero far from it, and they are held constant for 10 Myr. **b**, Predicted thermochronological ages; different symbols for thermochronological systems and the colour scale for ages are

as in Fig. 4a. **c-e**, Predicted erosion rates for a correlation length of 30 km from 6-4 Myr ago (**c**), from 2-0 Myr ago (**d**), and normalized difference in erosion rates from the comparison of **c** to **d** (as defined in equation (1)), shown only where the resolution in each time bin is greater than 0.25 (e). Contours in **c-e** show the predicted resolution. **f-h**, Same as plots **c-e**, but for a correlation length of 10 km.
<figure><img src="image_11.png" /><figcaption>**Extended Data Fig. 6** | Synthetic test for the Wasatch Mountains.</figcaption></figure>

a, Input exhumation rates increase smoothly from zero at the eastern boundary to 1.2 mm yr⁻¹ at the western boundary and are held constant for 12 Myr. b, Predicted thermochronological ages. Symbols for thermochronological systems and colour scale for ages are as in Supplementary Fig. 1a. c–e, Predicted erosion rates for a correlation length

of 30 km from 6–4 Myr ago (c), from 2–0 Myr ago (d), and normalized
difference in erosion rates from the comparison of c to d (as defined in
equation (1)), shown only where the resolution in each time bin is >0.25
(e). Contours in c–e show the predicted resolution. f–h, Same as c–e, but
for a correlation length of 10 km. i–k, Same as f–h, but for a prior erosion
rate (e_p) of 0.1 mm yr⁻¹ instead of 0.35 mm yr⁻¹.
<table><caption>Extended Data Table 1 | Global inversion results from ref. ⁶ and likely causes for erosion-rate increases</caption><thead><tr><th rowspan="2">Location</th><th colspan="2">Bins 6-4 Ma and 2-0 Ma</th><th colspan="2">All four bins</th><th colspan="4">Normalized difference 2-0 Ma minus 6-4 Ma</th><th rowspan="2">Tectonic setting</th><th rowspan="2">Cause of increase*</th></tr><tr><th>Points</th><th>% total</th><th>Points</th><th>% total</th><th>min</th><th>max</th><th>average</th><th>median</th></tr></thead><tbody><tr><td>Alps</td><td>360</td><td>21.8</td><td>328</td><td>30.0</td><td>-0.44</td><td>0.93</td><td>0.31</td><td>0.30</td><td></td><td></td></tr><tr><td>Eastern Alps</td><td>70</td><td>4.2</td><td>69</td><td>6.3</td><td>0.00</td><td>0.13</td><td>0.07</td><td>0.08</td><td>Post-orogenic</td><td>No significant increase</td></tr><tr><td>Western Alps</td><td>290</td><td>17.6</td><td>259</td><td>23.7</td><td>-0.44</td><td>0.93</td><td>0.37</td><td>0.38</td><td>Normal/thrust faulting</td><td>Spurious/Glacial</td></tr><tr><td>New Zealand</td><td>221</td><td>13.4</td><td>130</td><td>11.9</td><td>-0.67</td><td>0.99</td><td>0.51</td><td>0.57</td><td></td><td></td></tr><tr><td>Marlborough Fault System</td><td>28</td><td>1.7</td><td>14</td><td>1.3</td><td>0.15</td><td>0.99</td><td>0.53</td><td>0.49</td><td>Transpression</td><td>Spurious</td></tr><tr><td>Mt. Cook region</td><td>56</td><td>3.4</td><td>11</td><td>1.0</td><td>-0.67</td><td>0.98</td><td>0.70</td><td>0.83</td><td>Accretionary wedge</td><td>Spurious</td></tr><tr><td>Fiordland</td><td>137</td><td>8.3</td><td>105</td><td>9.6</td><td>-0.14</td><td>0.90</td><td>0.43</td><td>0.47</td><td>Accretionary wedge</td><td>Spurious/Tectonic</td></tr><tr><td>W. Himalaya</td><td>200</td><td>12.1</td><td>122</td><td>11.2</td><td>-0.43</td><td>0.81</td><td>0.21</td><td>0.20</td><td>Thrust ramp</td><td>Spurious</td></tr><tr><td>Patagonia</td><td>146</td><td>8.8</td><td>60</td><td>5.5</td><td>-0.06</td><td>0.85</td><td>0.27</td><td>0.22</td><td>Transpression</td><td>Spurious/Tectonic</td></tr><tr><td>Northern Apennines</td><td>117</td><td>7.1</td><td>91</td><td>8.3</td><td>-0.14</td><td>0.99</td><td>0.43</td><td>0.46</td><td>Accretionary wedge</td><td>Spurious/Tectonic</td></tr><tr><td>Taiwan</td><td>87</td><td>5.3</td><td>27</td><td>2.5</td><td>-0.94</td><td>0.99</td><td>0.61</td><td>0.74</td><td>Accretionary wedge</td><td>Spurious/Tectonic</td></tr><tr><td>Central Himalaya</td><td>75</td><td>4.5</td><td>53</td><td>4.8</td><td>-0.83</td><td>0.91</td><td>0.20</td><td>0.10</td><td>Thrust ramp</td><td>Spurious</td></tr><tr><td>Namche Barwa</td><td>66</td><td>4.0</td><td>49</td><td>4.5</td><td>-0.95</td><td>0.89</td><td>0.29</td><td>0.35</td><td>Crustal-scale pop-up</td><td>Spurious</td></tr><tr><td>Bhutan</td><td>53</td><td>3.2</td><td>38</td><td>3.5</td><td>-0.29</td><td>0.48</td><td>0.08</td><td>0.09</td><td>Thrust ramp</td><td>Spurious</td></tr><tr><td>Gongga Shan</td><td>49</td><td>3.0</td><td>42</td><td>3.8</td><td>-0.09</td><td>0.98</td><td>0.41</td><td>0.28</td><td>Transpression</td><td>Spurious</td></tr><tr><td>Olympic Mts.</td><td>39</td><td>2.4</td><td>37</td><td>3.4</td><td>0.13</td><td>0.56</td><td>0.31</td><td>0.28</td><td>Accretionary wedge</td><td>Spurious</td></tr><tr><td>Merida Andes</td><td>36</td><td>2.2</td><td>6</td><td>0.5</td><td>0.32</td><td>0.79</td><td>0.61</td><td>0.63</td><td>Transpression</td><td>Spurious</td></tr><tr><td>British Columbia</td><td>33</td><td>2.0</td><td>21</td><td>1.9</td><td>0.09</td><td>0.61</td><td>0.34</td><td>0.29</td><td>Valley incision</td><td>Glacial</td></tr><tr><td>Nanga Parbat</td><td>24</td><td>1.5</td><td>13</td><td>1.2</td><td>0.52</td><td>0.91</td><td>0.74</td><td>0.75</td><td>Crustal-scale pop-up</td><td>Spurious</td></tr><tr><td>Longmen Shan</td><td>21</td><td>1.3</td><td>16</td><td>1.5</td><td>0.20</td><td>0.65</td><td>0.43</td><td>0.43</td><td>Thrust ramp</td><td>Spurious</td></tr><tr><td>San Gabriel Mts.</td><td>18</td><td>1.1</td><td>17</td><td>1.6</td><td>0.36</td><td>0.78</td><td>0.56</td><td>0.55</td><td>Transpression</td><td>Spurious</td></tr><tr><td>Eritrea</td><td>18</td><td>1.1</td><td>6</td><td>0.5</td><td>0.00</td><td>0.98</td><td>0.47</td><td>0.38</td><td>Normal-fault block</td><td>Spurious</td></tr><tr><td>Alborz</td><td>16</td><td>1.0</td><td>9</td><td>0.8</td><td>0.23</td><td>0.81</td><td>0.46</td><td>0.42</td><td>Accretionary wedge</td><td>Tectonic</td></tr><tr><td>Wasatch</td><td>14</td><td>0.8</td><td>9</td><td>0.8</td><td>0.05</td><td>0.29</td><td>0.20</td><td>0.21</td><td>Normal-fault block</td><td>Spurious</td></tr><tr><td>St. Elias</td><td>11</td><td>0.7</td><td>0</td><td>0.0</td><td>0.51</td><td>0.62</td><td>0.59</td><td>0.60</td><td>Accretionary wedge</td><td>Tectonic/Glacial</td></tr><tr><td>Papua New Guinea</td><td>10</td><td>0.6</td><td>3</td><td>0.3</td><td>0.00</td><td>0.80</td><td>0.15</td><td>0.03</td><td>Accretionary wedge</td><td>Tectonic</td></tr><tr><td>San Juan Mts.</td><td>7</td><td>0.4</td><td>0</td><td>0.0</td><td>0.59</td><td>0.99</td><td>0.83</td><td>0.85</td><td>Post-orogenic</td><td>Spurious</td></tr><tr><td>Pamir</td><td>7</td><td>0.4</td><td>6</td><td>0.5</td><td>0.02</td><td>0.08</td><td>0.05</td><td>0.05</td><td>Normal-fault block</td><td>No significant increase</td></tr><tr><td>Southern Peru</td><td>6</td><td>0.4</td><td>4</td><td>0.4</td><td>0.27</td><td>0.50</td><td>0.43</td><td>0.47</td><td>Valley incision</td><td>Spurious</td></tr><tr><td>Bolivia</td><td>5</td><td>0.3</td><td>5</td><td>0.5</td><td>0.12</td><td>0.28</td><td>0.18</td><td>0.15</td><td>Valley incision</td><td>Spurious</td></tr><tr><td>Southeast Tibet</td><td>3</td><td>0.2</td><td>2</td><td>0.2</td><td>0.05</td><td>0.12</td><td>0.07</td><td>0.05</td><td>Thrust ramp</td><td>Tectonic</td></tr><tr><td>Aconquija</td><td>3</td><td>0.2</td><td>0</td><td>0.0</td><td>0.48</td><td>0.52</td><td>0.50</td><td>0.50</td><td>Thrust block</td><td>Spurious</td></tr><tr><td>Greater Caucasus</td><td>2</td><td>0.1</td><td>0</td><td>0.0</td><td>0.56</td><td>0.65</td><td>0.60</td><td>0.60</td><td>Thrust ramp</td><td>Tectonic</td></tr><tr><td>Tien Shan</td><td>1</td><td>0.1</td><td>0</td><td>0.0</td><td>0.35</td><td>0.35</td><td>0.35</td><td>0.35</td><td>Thrust block</td><td>Tectonic/Glacial</td></tr></tbody></table>

<footnote>*See Supplementary Information for detailed descriptions. Ma, million years ago.</footnote>
<table><caption>Extended Data Table 2 | Thermochronological ages corresponding to steady erosion rates</caption><thead><tr><th>Erosion rate (mm yr<sup>-1</sup>)</th><th>AHe age (Ma)</th><th>AFT age (Ma)</th><th>ZHe age (Ma)</th><th>ZFT age (Ma)</th></tr></thead><tbody><tr><td>0.07</td><td>31.3</td><td>50.5</td><td>85.7</td><td>104.0</td></tr><tr><td>0.10</td><td>22.4</td><td>35.8</td><td>60.4</td><td>73.1</td></tr><tr><td>0.15</td><td>15.2</td><td>24.1</td><td>40.4</td><td>48.7</td></tr><tr><td>0.20</td><td>11.5</td><td>18.1</td><td>30.2</td><td>36.3</td></tr><tr><td>0.30</td><td>7.6</td><td>11.9</td><td>19.8</td><td>23.8</td></tr><tr><td>0.40</td><td>5.7</td><td>8.7</td><td>14.6</td><td>17.4</td></tr><tr><td>0.50</td><td>4.4</td><td>6.8</td><td>11.4</td><td>13.6</td></tr><tr><td>0.60</td><td>3.6</td><td>5.6</td><td>9.2</td><td>11.0</td></tr><tr><td>0.70</td><td>3.0</td><td>4.6</td><td>7.7</td><td>9.2</td></tr><tr><td>0.80</td><td>2.6</td><td>3.9</td><td>6.5</td><td>7.7</td></tr><tr><td>0.90</td><td>2.2</td><td>3.4</td><td>5.6</td><td>6.7</td></tr><tr><td>1.20</td><td>1.5</td><td>2.3</td><td>3.9</td><td>4.6</td></tr><tr><td>1.50</td><td>1.1</td><td>1.7</td><td>2.8</td><td>3.4</td></tr><tr><td>2.00</td><td>0.7</td><td>1.1</td><td>1.8</td><td>2.2</td></tr><tr><td>3.00</td><td>0.4</td><td>0.6</td><td>1.0</td><td>1.1</td></tr></tbody></table>
<table><caption>**Extended Data Table 3** | Diffusion parameters used for closure-temperature calculations</caption><thead><tr><th>Thermochronometer<br>(reference)</th><th>E<sub>a</sub><br>(kJ mol<sup>-1</sup>)</th><th>D<sub>0</sub><br>(cm<sup>2</sup> s<sup>-1</sup>)</th><th>a<br>(µm)</th><th>T<sub>c</sub>; ṇ=0.1<br>mm yr<sup>-1</sup><br>(°C)†</th><th>T<sub>c</sub>; ṇ=1.0<br>mm yr<sup>-1</sup><br>(°C)</th></tr></thead><tbody><tr><td>AHe <sup>41</sup></td><td>138</td><td>50</td><td>100</td><td>69</td><td>88</td></tr><tr><td>ZHe <sup>42</sup></td><td>169</td><td>0.46</td><td>100</td><td>189</td><td>219</td></tr><tr><td>AFT <sup>14‡</sup></td><td>147</td><td>3.73 x 10<sup>4</sup></td><td></td><td>110</td><td>134</td></tr><tr><td>ZFT <sup>14§</sup></td><td>208</td><td>1.82 x 10<sup>6</sup></td><td></td><td>230</td><td>258</td></tr></tbody></table>

In all calculations, R = 8.3145 J K⁻¹ mol⁻¹, A = 27, L = 30 km, κ = 30 km² Myr⁻¹, and the initial geothermal gradient is 30°C km⁻¹ (T₀ = 0°C; T_L = 900°C).

†Predicted closure temperature (T_c) for slow (0.1 mm yr⁻¹) and rapid (1 mm yr⁻¹) exhumation rates (e).

‡Values corresponding to the 'average apatite composition' of ref. ¹⁴, with a value of D₀/a² (s⁻¹).

§Values corresponding to the 'natural, radiation damaged' zircon of ref. ¹⁴, with a value of D₀/a² (s⁻¹).