Abstract

Deciphering erosion rates over geologic time is fundamental for understanding the interplay between climate, tectonic, and erosional processes. Existing techniques integrate erosion over different time scales, and direct comparison of such rates is routinely done in earth science. On the basis of a global compilation, we show that erosion rate estimates in glaciated landscapes may be affected by a systematic averaging bias that produces higher estimated erosion rates toward the present, which do not reflect straightforward changes in erosion rates through time. This trend can result from a heavy-tailed distribution of erosional hiatuses (that is, time periods where no or relatively slow erosion occurs). We argue that such a distribution can result from the intermittency of erosional processes in glaciated landscapes that are tightly coupled to climate variability from decadal to millennial time scales. In contrast, we find no evidence for a time scale bias in spatially averaged erosion rates of landscapes dominated by river incision. We discuss the implications of our findings in the context of the proposed coupling between climate and tectonics, and interpreting erosion rate estimates with different averaging time scales through geologic time.

Keywords

Erosion

climate variability

power-law

timescale dependence

glaciers

river incision

INTRODUCTION

Erosion drives landscape evolution and forms the link between climate and tectonics (1). Understanding changes in erosion rates through time is thus of paramount importance for deciphering Earth’s topographic evolution and evaluating the relationships between climate and tectonics (2–5). Estimated rates of landscape evolution (either erosion rates of upland catchments or sedimentation rates in depositional basins) integrate over a wide range of time scales, and controversy exists regarding the evidence for temporal variation in erosion rates and their interpretation (2, 3, 5–12). For example, the sedimentary record has commonly been interpreted as a measure of the erosional history of upland catchments, and it has been widely argued that global increases in sediment accumulation rates reveal the response of mountain erosion to climate cooling during Late Cenozoic time (3, 5). However, several confounding factors, such as stratigraphic incompleteness (8, 10, 12), sediment storage within a landscape, and internal dynamics of sedimentary basins (13–15), can introduce time-dependent biases that mask, alter, and mimic temporal changes in sediment accumulation rates. These time-dependent biases commonly manifest as observations of high sedimentation rates over short averaging time scales, which decay as a power law with increasing averaging time scale. Local measurements of river incision inferred from terraces also exhibit a time-dependent bias over millennial time scales (6), but it is unknown whether similar confounding factors exist at larger landscape scales.

Because climate and tectonics are the primary drivers of landscape-scale erosion, variations in erosion rates can give insights into the climatic and tectonic history of landscapes. However, it has been proposed that erosion is not a continuous process but occurs during discrete events (that is, “pulses of erosion”) (6, 16, 17) that are separated by intervals with little or no erosion, which we term “erosional hiatuses” (Fig. 1). Estimates of landscape-scale erosion (hereafter referred to as “estimated erosion rates”) typically average over both erosional pulses and hiatuses; therefore, it is important to evaluate how changes in estimated erosion rates reflect changes in the magnitude or temporal distribution of erosional pulses with age, measured as time before present.

Schematic of a time series of landscape-scale erosion highlighting the erosional pulses (blue bars) and hiatuses. In this framework, the functional dependence of estimated erosion rates on averaging time scale is determined by the probability distributions of the magnitudes of erosional pulses and hiatuses (fig. S1 and Supplementary Note 1).

Over geologic time, both the magnitudes of erosional pulses and durations of erosional hiatuses are expected to vary. In landscapes where erosion is dominated by river incision, erosion rates are often positively correlated with tectonic uplift rates (18–20) and, in some cases, precipitation rates (21–23), whereas the characteristic recurrence time of processes that drive landscape erosion varies between different climatic and tectonic settings. For example, in unglaciated bedrock landscapes like the San Gabriel Mountains, CA, some of the variability in erosion rates arises from wildfires, which have a recurrence interval of ~30 years (24, 25), and from changes in uplift rates due to tectonics over millions of years (26). In some landscapes, the bulk of the erosional flux has been attributed to landslides triggered by infrequent large earthquakes that occur on centennial to millennial time scales (27). Similarly, large-magnitude storm events may introduce intermittency in river erosion and landsliding over centennial time scales (28). In glaciated landscapes, previous workers have argued that the magnitudes of erosional pulses have increased over the Late Cenozoic in response to climate cooling, onset of glaciation, and large oscillations in climate (2, 3, 5, 7). The intermittency in glacial erosion is postulated to occur via changes in glacier extent and sliding velocity, which should reflect climate forcing over a wide range of time scales from months to millions of years (2, 29).

The probabilistic structure of the magnitude and recurrence of erosional pulses determines the nature of time scale dependence of the estimated erosion rates over long time scales—a relationship that can be compared to multiscale geochronological estimates of erosion rates from different landscapes (fig. S1 and Supplementary Note 1) (10). However, most existing techniques provide estimated erosion rates that integrate erosional pulses and hiatuses from some time in the past to the present, and paleo-erosion rates that integrate erosional pulses and hiatuses over time periods in the past are rare (30–33). Therefore, the age that corresponds to an estimated erosion rate is often necessarily correlated to the time span over which it averages. This correlation between the age and averaging time scale of estimated erosion rates is problematic because assessing changes in the frequency and magnitude of erosional pulses might require erosion rates measured over equal intervals of time (fig. S3) (10). Typically, sediment yield data average over decades, measurements of cosmogenic nuclides (CRNs) average over millennia, and thermochronology averages erosional pulses and hiatuses over millions of years. Thus, to compare erosion rates derived from different methods and to attribute changes in estimated erosion rates to secular changes in climate and tectonics, it is critical to quantify the time scales over which transience in landscape-scale erosion rates persists. This is particularly important for assessing the proposed feedbacks between climate and tectonic processes (2, 3, 5).

To assess the impact of the averaging time scale on estimated erosion rates, we compiled estimated erosion rates from both fluvial and glaciated landscapes across the globe (table S2). First, we show that there is not a systematic trend in estimated erosion rates with averaging time scale in landscapes dominated by river incision. Second, we show that estimated erosion rates in glaciated landscapes exhibit a systematic scaling relation with averaging time scale (time scale bias) from decades to millions of years—a pattern that is consistent with theoretical expectations arising from heavy-tailed erosional hiatuses. Third, using numerical simulations, we assess whether changes in magnitudes of erosional pulses can be masked by changes in estimated rates because of the averaging time scale bias. Finally, we discuss the possible origin and implications of heavy-tailed erosional hiatuses in interpreting estimated erosion rates with age.

RESULTS

We compiled a global data set of estimated erosion rates from landscapes in diverse climate and tectonic settings that covers time scales from 101 to 108 years (Materials and Methods and table S2). We used previously published data of sediment yield to assess decadal-scale estimated erosion rates; volumetric rates of lake sediment fill, fjord sedimentation, and CRNs to assess millennial-scale estimated erosion rates; and mineral cooling ages from low-temperature thermochronology to assess million-year time scale estimated erosion rates (Materials and Methods). Measurements of bed load and suspended sediment fluxes provide an estimate of basin-averaged erosion rate over the time scale of measurement, which is typically on the order of decades. Because of their relatively short measurement time, these rates may be affected by inadequate sampling of large events and also by anthropogenic interference (17). To compare these decadal-scale with millennial-scale estimates of erosion rates, we used previously published basin-averaged erosion rates derived from CRNs in alluvial sediments, which average over the time required to erode roughly ~60 cm of rock (34–39). To assess millennial-scale erosion rates in glaciated landscapes, we used previously published measurements of lake sediment fill, fjord sedimentation, and CRN data that were corrected for shielding from snow and ice (40–42). Finally, we compared the decadal- and millennial-scale erosion rate estimates with million-year erosion rate estimates by using previously published bedrock thermochronology data. Thermochronological data are point measurements that record mineral cooling histories, which provide information of land surface lowering (43, 44). Because common thermochronometers [such as fission tracks or (U-Th)/He in apatite or zircon] integrate cooling histories over several kilometers depth, these point measurements reflect landscape-scale erosion. To make the comparison with the other two methods, we averaged these point measurements across the area of interest for each thermochronometer with a fixed closure temperature. Although a scaling exponent of −1 between estimated erosion rate and averaging time scale is expected because of spurious correlation for data derived from a single chronometer (45), we compare here the rates across multiple chronometers with a wide range of averaging time scales (Supplementary Note 2 and fig. S3).

Results show that the time scale bias in estimated erosion rates in glaciated landscapes (including landscapes that were predominantly shaped by glacial processes but are dominated by fluvial processes at present, such as the Rhône Valley in the Swiss Alps) significantly differs from locations where glaciers were not the dominant erosional agent (Fig. 2). Across most of the glaciated landscapes, we find that estimated erosion rates show a pronounced (often about two orders of magnitude) increase toward the present, which is characterized by an inverse power-law trend on the averaging time scale (Fig. 2A and Materials and Methods). Despite the diversity of tectonic settings (for example, active convergence in the St. Elias Range, AK versus relative tectonic quiescence in Svalbard, Norway), the scaling exponents of the inverse power-law trends fall within a range of −0.71 to −0.11 (Fig. 3 and table S1).

Fig. 2Worldwide compilation of erosion rate estimates over a range of averaging time scales.

(A) Estimated erosion rates in glaciated landscapes show marked increase with a decrease in averaging time scale, which correlates with age and is characterized by an inverse power-law trend (Materials and Methods). My, million years; NW, Northwest. (B) Estimated erosion rates in landscapes where glacial processes are not dominant. The colored boxes around the square markers indicate the SE in averaging time scale and erosion rate estimates (table S2). The colored dashed lines indicate the best power-law fit to the data of estimated erosion rates and the averaging time scale (Materials and Methods). The dashed black lines denote a slope of −1, which could be expected because of plotting a rate against time (45).

Plot showing the estimated scaling exponents (including SE) of the best-fitting power laws, which describe the time scale bias in erosion rates of glaciated landscapes. Results indicate that the scaling exponents across a wide range of glaciated landscapes lie within a relatively narrow range of −0.71 to −0.11.

In contrast, in landscapes where glaciers are subordinate agents of erosion, we do not observe a negative power-law relationship between estimated erosion rates and the age or averaging time scale (Fig. 2B). In Taiwan, the central Himalayas, and the San Gabriel Mountains, estimated erosion rates from decades to millions of years are comparable (Fig. 2B), suggesting that erosion rates may not have changed significantly in these landscapes (18, 46, 47). This concordance between estimated erosion rates and age was also reported in fluvial landscapes in tectonically inactive regions like the Appalachians, where previous work has documented relatively constant rates of estimated erosion from decades to several million years (48); however, we note that spatial variability of estimated erosion rates within the orogen has also been documented (49–51). Although the San Gabriel Mountains harbor approximately an order of magnitude of spatial variability in erosion rates (18), our data do not show a sustained increase toward the present. Moreover, although variability in erosion rates within both the San Gabriel and Himalayan mountains provides evidence of transient topographic adjustment (52–55), it does not apparently manifest itself in a time scale dependence of estimated erosion rates in either landscape. Despite known variability in magnitudes of erosional pulses in fluvial landscapes due to intermittent floods or landslides (6, 17, 28), our data show no evidence for a time scale bias in estimated erosion rates from decades to millions of years.

Age and time scale bias in estimated erosion rates of glaciated landscapes

The most striking finding from our analysis is the observed power-law decrease of estimated erosion rates with averaging time scale or age in glaciated landscapes (Fig. 2A). Although previous workers observed a similar decrease in estimated erosion rates with age (41, 56), it is unknown whether this represents changes in the magnitude and frequency of erosional pulses with age or it emerges because of the systematic increase in averaging time scale that correlates with age. We can consider the causes of this apparent time scale dependence of estimated erosion rates using two end-member frameworks, which we term as “deterministic” and “probabilistic.” The first and possibly most direct interpretation of this trend can be obtained within a deterministic framework, where the temporal trends in estimated erosion rates result directly from an increase in magnitudes of erosional pulses toward the present (Fig. 4A). Within this deterministic framework, the magnitude of erosional pulses should increase as a power-law function in time toward the present (Fig. 4A), such that the estimated erosion rates show a power-law decrease with age (Fig. 4C). The implications and requirements of this deterministic model are as follows: (i) the magnitude of erosional pulses has to be high during the present, (ii) the magnitude of erosional pulses should decrease monotonically as a power-law function of age, and (iii) these relationships must be consistent globally across glaciated landscapes because the range of scaling exponents that describe the relationship between estimated erosion rates and the averaging time scale is relatively narrow (Figs. 2A and 3). A scaling exponent of −0.5, for example, implies that erosion rates double in time toward the present every 0.25T, where T is the age. For example, this model predicts that the erosion rates estimated 6 years ago are double those estimated 25 years ago, which in turn are double those estimated by sediment yield data ~100 years ago, and so on.

(A) Representative time series of the magnitudes of erosion with age for the deterministic model, where the magnitudes of erosion decline as a power-law function with age. (B) Representative time series of the magnitudes of erosion with age for the probabilistic model. The magnitudes of erosional pulses in this simulation were held constant with age and were separated by erosional hiatuses drawn from a truncated Pareto distribution (tail index, 0.5; upper bound, 200 ky). (C) Estimated erosion rates for both the time series shown in (A) (deterministic; gray markers) and (B) (probabilistic; red markers) as a function of increasing averaging time scale. In these simulations, the estimated erosion rates average from some time in the past to the present, such that the averaging time scale is equal to the age. Both models result in an inverse power-law relationship of estimated erosion rates with averaging time scale (or age), but for different reasons. For the deterministic model, the power-law trend occurs because the magnitudes of erosional pulses are imposed to decrease as a power-law function of age, and averaging time scale correlates with age. For the probabilistic model, the power-law trend occurs because of the intrinsic dependence of estimated erosion rates on averaging time scale for heavy-tailed erosional hiatuses. We generated 1000 independent realizations of the probabilistic model and estimated the slope and the intercept of the best-fitting power law between estimated erosion rates and averaging time scales less than the upper bound. The black and gray lines show the best-fitting power law function with a mean power-law slope (−0.50) along with the SE (0.004) evaluated by averaging over these 1000 different realizations. The estimated erosion rates saturate for time scales greater than the upper bound on erosional hiatuses (shaded gray area).

Alternatively, previous work with regard to explaining similar observations of power-law decrease in sedimentation rates with age has shown that inverse power-law trends in estimated rates can arise from averaging multiscale variability in hiatuses, which are characterized by heavy-tailed distributions (fig. S1 and Supplementary Note 1) (10, 14). In this probabilistic framework, erosional hiatuses with a probability density function that decays as a power law with a tail index α < 1 will result in a power-law dependence of estimated erosion rates on the averaging time scale (scaling exponent of α − 1; Supplementary Note 1 and fig. S1) (6, 10). Unlike the case of thin-tailed erosional hiatuses (for example, exponentially decaying), heavy-tailed erosional hiatuses result in systematic divergence of estimated erosion rates over averaging time scales that are characterized by the aforementioned power law (6, 10). Figure 4B shows a simple end-member schematic representation of pulses of erosion with age that have a constant magnitude but are interspersed with erosional hiatuses, which are described by a truncated heavy-tailed distribution [truncated Pareto distribution with tail index of 0.5 and an upper bound of 200 ky (thousand years)]. The estimated erosion rates within this probabilistic framework also show a negative power-law dependence on averaging time scale (Fig. 4C), for averaging time scales that are less than the longest possible hiatus (given by the upper bound on the Pareto distribution), consistent with observations from our global compilation of estimated erosion rates in glaciated landscapes (Fig. 2A). For averaging times greater than the longest hiatus, the estimated erosion rates are no longer a function of the averaging time scale (Fig. 4C). Because most measurement techniques average erosion from some time in the past to the present, estimated erosion rates often average over a time interval that correlates with age; thus, the power-law trend of estimated erosion rate with averaging time scale in Fig. 4C could be mistaken as a relation between estimated erosion rate and age (fig. S2). The implications and requirements of this probabilistic model are as follows: (i) erosional hiatuses exhibit variability over a wide range of time scales and can be described by a heavy-tailed distribution (Supplementary Note 1) (6, 10); (ii) the magnitudes of erosional pulses can, but need not, change with age; and (iii) any relationship that describes the magnitudes of erosional pulses need not be globally consistent because the relationship between estimated erosion rate and the averaging time scale is dominated by the heavy-tailed nature of the erosional hiatuses rather than the magnitude of erosional pulses (10). We note that in this simulation the magnitudes of erosional pulses were selected to be constant for ease of demonstration. However, the time scale bias in estimated erosion rates persists even if the magnitudes of erosional pulses are variable and thin-tailed (fig. S4 and Supplementary Note 3) (10); only the erosional hiatuses need to be heavy-tailed to reconcile the observed averaging time scale–dependent power-law behavior.

To adequately explain the data presented in Fig. 2A, we argue that the deterministic model is far more restrictive than the probabilistic model because it requires the magnitudes of erosional pulses to be globally high during the present and also to decrease as a power-law function with age, with unrealistically accelerated erosion rates predicted in the immediate future. In contrast, the probabilistic model requires only that the erosional hiatuses be heavy-tailed. Thus, we favor the interpretation that the observed power-law decrease of estimated erosion rates with age (Fig. 2A) is a consequence of sampling heavy-tailed erosional hiatuses over averaging time scales that increase with age. Nevertheless, it is likely that the magnitudes of erosional pulses have changed over time, in contrast to the inputs imposed in our end-member model. If this is correct, then this raises the question of whether the increase in estimated erosion rates due to heavy-tailed erosional hiatuses masks changes in the magnitudes of erosional pulses through time. Moreover, we implicitly implied in our end-member model in Fig. 4 that erosional hiatuses are spatially coherent across the landscape, which may not be the case in natural landscapes where erosional pulses can be discrete in both space and time. This raises the question of whether spatial averaging affects the dependence of estimated erosion rates on averaging time scale.

Sensitivity analysis of age and time scale bias in estimated erosion rates

Using numerical simulations, we determined whether heavy-tailed distributions of erosional hiatuses can overshadow the signal of a systematic increase in the magnitudes of erosional pulses with age. These numerical simulations were designed to evaluate end-member cases that highlight the role of heavy-tailed erosional hiatuses and how their probabilistic structure affects estimated erosion rates over different averaging time scales. We generated a time series of magnitudes of erosional pulses, which were separated by heavy-tailed erosional hiatuses (truncated Pareto distribution; tail index, 0.5; upper bound, 200 ky). Similar to most geochronological techniques, we calculated the estimated erosion rate by averaging erosional pulses and hiatuses from some time in the past to the present; thus, the age of the estimated erosion rate also corresponds to the averaging time scale. As a basis for evaluating the influence of increasing erosion rates over time, we considered a control case in which the magnitudes of erosional pulses were held constant through time (10 mm), and the resulting decrease in estimated erosion rates with age is solely a consequence of time scale bias, equivalent to results shown in Fig. 4.

We explored two additional scenarios, which we compare to the control case, where the magnitudes of erosional pulses during 0 to 50 ka (thousand years ago) were by 2-fold [approximately the median of reported worldwide increase in the study by Herman et al. (2); Fig. 5A] and 10-fold of their value before 50 ka. Note that the age of the increase at 50 ka was arbitrarily chosen for demonstration purposes and that the model behavior is insensitive to this choice. The estimated erosion rates for the first 50 ky are identical for the control case and the cases where there is an imposed increase in magnitudes of erosional pulses (Fig. 5B). This is expected because no erosion rate proxy would sample the increase in magnitudes of erosional pulses over this time frame. Although there is an increase (2- and 10-fold) in magnitudes of erosional pulses at 50 ky, our results indicate that the estimated erosion rates are indistinguishable from the control case up to an averaging time scale of 200 ky (Fig. 5B). The change in magnitude of erosional pulses by 2-fold or even 10-fold has little effect on the estimated erosion rates because the change in averaging time scale with age spans five orders of magnitude, which in turn results in an estimated erosion rate that changes by more than a factor of 100 (fig. S2). For averaging time scales greater than the upper bound of the erosional hiatuses, the change in magnitudes of erosional pulses is evident when compared against each other. Thus, when erosional hiatuses have a heavy-tailed distribution, temporal trends in estimated erosion rates can confidently be interpreted to reflect changes in magnitudes of erosional pulses only when the averaging time scale is held constant with age (fig. S2) or when the averaging time scales are longer than the longest possible erosional hiatus.

Fig. 5Effect of heavy-tailed erosional hiatuses on masking changes in magnitudes of erosional pulses with age.

(A) Representative time series of magnitudes of erosional pulses showing a twofold increase at 50 ky. (B) Estimated erosion rates as a function of the averaging time scale for the cases when the upper bound on erosional hiatuses was 200 ky and the erosional pulses had a 2-fold (blue markers) and 10-fold (red markers) increase at 50 ka. The estimated erosion rates average from some time in the past to the present, such that the averaging time scale is equal to the age. The estimated erosion rates for the cases with an increase in the magnitudes of erosional pulses are indistinguishable from the control case for averaging time scales that are shorter than the longest possible hiatus (red markers plot on top of the gray and blue markers). The black and gray lines show the best-fitting power-law function with mean power-law slope (−0.50) along with the SE (0.004) evaluated by averaging over these 1000 different realizations. (C) Plot showing the results of time scale dependence of estimated erosion rates where time series of erosion with equal magnitudes of erosional pulses (10 mm) and the same truncated Pareto distribution of erosional hiatuses (tail index, 0.5; upper bound, 200 ky) were spatially averaged (Materials and Methods). The estimated erosion rates average from some time in the past to the present, such that the averaging time scale is equal to the age. Results show that the time scale bias in estimated erosion rates is reduced with increasing number of series that are spatially averaged, that is, large values of n. The best-fitting power-law slope, along with SE resulting from averaging across multiple realizations, is indicated in the legend, where the power-law fits to estimated erosion rates were made for averaging time scales between 10 and 200 ky (upper bound on erosional hiatuses).

Our analysis of the probabilistic model suggests that the temporal trends in data from glaciated landscapes (Fig. 2A) can arise without changes in magnitudes of erosional pulses over time because of erosional hiatuses that have a heavy-tailed distribution. In the simulations in Fig. 5B, we assumed that erosional pulses and hiatuses are spatially correlated across the area of interest. This might be typical of glacial landscapes, as discussed below; however, estimated erosion rates might also average uncorrelated and spatially discrete erosional events across a landscape, which could affect the trend of estimated erosion rate with averaging time scale. To evaluate the effect of spatial averaging, we generated several (n = 1 to 100) independent time series of erosion pulses that each represent uncorrelated erosional events acting on spatially discrete portions of the landscape, and averaged these series to calculate spatially averaged erosion across the landscape (Materials and Methods). Numerical simulations show that when erosional hiatuses are drawn from the same probability distribution (truncated Pareto distribution; tail index, 0.5; upper bound, 200 ky), the time scale bias in estimated erosion rates is a function of the degree of spatial averaging across a landscape (that is, the value of n in Fig. 5C). The power-law slope flattens with increasing degree of spatial averaging until, beyond a certain large n, the time scale bias in estimated erosion rates is eliminated (Fig. 5C). The time scale bias in estimated erosion rates is eliminated when the product of n and Δt, where Δt is the averaging time scale, is larger than the upper bound of erosional hiatuses (Fig. 5C). For example, in Fig. 5C, n = 100 eliminates the power-law scaling for averaging time scales greater than Δt = 2 ky (that is, 200 ky/n).

Within the probabilistic framework, the lack of a time scale bias in estimated erosion rates for fluvial landscapes may reflect a high degree of spatial averaging of uncorrelated erosional events, even if local erosional events are characterized by heavy-tailed hiatuses, as has been inferred from the record of river terraces (6). It might also reflect spatial averaging of erosional pulses that are characterized by thin-tailed hiatuses (fig. S1). Moreover, different power-law scaling exponents in glaciated landscapes (Fig. 3) may be attributed to different degrees of spatial averaging of heavy-tailed erosional hiatuses that are coherent across a landscape. Nonetheless, for the probabilistic model to explain the time scale bias in estimated erosion rates of glaciated landscapes (Fig. 2A), either spatial averaging must be minimal or erosional variability must be reasonably correlated across landscapes over a wide range of time scales approaching millions of years—the possibility of which is discussed below.

DISCUSSION

We hypothesized that in glaciated landscapes, relatively high-magnitude erosional pulses interspersed with a heavy-tailed distribution of spatially correlated erosional hiatuses could arise from the influence of natural climate variability on erosional processes. The principal processes of glacial erosion are abrasion and quarrying, which both scale with the sliding velocity of glaciers (29, 57, 58). Although glacial sliding is complex and no generally accepted sliding law exists, in most existing models, the sliding velocity of glaciers scales with the ice thickness, surface slope, and effective pressure and is therefore tightly coupled to climate (29). Hence, it is reasonable to expect that glacial erosion is influenced by secular trends in climate. Moreover, recent work by Herman et al. (57) suggests that glacial erosion is nonlinearly related to sliding velocity of glaciers with an exponent greater than unity, and thus, the variability in glacial erosion can be amplified when compared to the variability in glacial sliding. For example, over relatively short time scales (seasonal to centennial time scales), fluctuations in glacial sliding velocities (59), including surges (60, 61), can introduce variability in subglacial erosion (62). Over millennial and orbital time scales, the waxing and waning of ice masses introduce variability in sliding velocities of glaciers, which affect erosion within a catchment. Finally, over millions of years, the absence or presence of ice sheets or the changing conditions at the bed, that is, frozen-based versus wet-based (63, 64), can itself result in long-term variability in erosion that affects erosion rate estimates.

Climate history and the influence of glacial processes on landscapes can be informed by the δ18O values of marine calcite (derived from well-preserved benthic foraminifera), which provide a proxy for continental ice volume and global temperature (Fig. 6) (65, 66). Spectral analysis of the global δ18O record (Materials and Methods) (65, 66) reveals major periodicities at 100 and 40 ky, which are widely attributed to orbital variations in solar forcing. However, substantial variability persists down to the smallest resolvable scale, which is characterized by the power-law decay in spectral density (Fig. 6) (66). Nested phenomena, such as Heinrich events, Dansgaard-Oeschger events, Bond events, and the glacial-interglacial transitions, contribute to observed climate variability over multiple time scales. In glaciated landscapes, erosional variability is therefore expected to occur over a wide range of time scales that overlap with the long characteristic time scales of tectonic processes—dynamics that may lead to a heavy-tailed distribution of erosional hiatuses in glaciated landscapes.

The δ18O carbonate record (blue curve; inset) from marine basins (65) shows global cooling associated with the onset of late Plio-Pleistocene glaciations. Increasing δ18O values reflect cooler temperatures and greater continental ice volume and vice versa. The power spectral density (red curve; Materials and Methods) of this climate proxy indicates that apart from the major climatic periodicities at 100 and 40 ky (paced by variation in solar insolation), several nested characteristic time scales exist down to the smallest scale as shown by the power-law decay of the spectral density (66). We hypothesize that this climate variability leads to heavy-tailed erosional hiatuses and thus a time scale bias in estimated erosion rates of glaciated landscapes.

Although glacial processes may not operate during interglacial periods, landscapes that are currently dominated by fluvial incision (for example, the Rhône Valley) still carry the signature of past glaciations (67). Moreover, integrating over several cycles of climate variability may introduce additional complexity in the erosional processes in glaciated landscapes. In particular, topographic adjustment can lead to a period of increased erosion rates when glaciers begin sculpting landscapes that are not in equilibrium with glacial processes (68) or with the size of glaciers (69). Furthermore, topographic adjustment may reverse during short interglacial periods, when ice retreat results in debuttressing of steep hillslopes and frequent slope failures (70, 71). In addition, erosion rates from glaciated landscapes can be influenced by lag times between sediment production and transport (31); much of the sediment that is produced in glacially carved valleys can be transiently stored during the interglacial period because of relatively low topographic slopes (for example, in overdeepenings) that characterize many glacial landscapes (67, 72, 73). Moreover, recent work by Willett et al. (74) highlights the long time scales involved in the transient adjustment of landscapes because of drainage divide migration and river capture. This can be pronounced in glaciated landscapes owing to the significant topographic adjustment during glaciated periods (67, 73), indicating that the erosional effects of glaciation outlast the intervals of glaciation itself. Unlike fluvial-dominated landscapes (with no oscillations between ice and ice-free conditions), glaciated landscapes are subject to topographic adjustment over a wide range of time scales leading to glacially conditioned erosional processes (67, 73), which can lead to long time scales of transient adjustment.

Our global analysis supports the notion that multiscale climate variability in glaciated landscapes can result in a significant increase in estimated erosion rates toward the present because of a heavy-tailed distribution of erosional hiatuses and measurement techniques that average over a time interval that correlates with age. Our numerical simulations suggest that the longest erosional hiatus (that is, the upper bound on the heavy-tailed distribution), which can be inferred from the saturation of the time scale dependence of estimated erosion rates (Figs. 4C and 5B and Supplementary Note 1), determines the minimum averaging time scale over which estimated erosion rates are insensitive to the averaging time scale. In our data compilation, apart from Patagonia, a negative power-law dependence of estimated erosion rates on averaging time scale continues to persist over million-year time scales (Fig. 2A). Moreover, our simulations reveal that systematic changes in magnitudes of erosional pulses with age may not be measurable because the estimated erosion rate can shift by several orders of magnitude in response to different averaging time scales (Fig. 5B). Thus, our analysis suggests that inferring changes in the rates of glacial landscape evolution needs measurements to average erosion over equal time intervals, regardless of age (fig. S2), or average over large temporal and spatial scales relative to the characteristic size and hiatus durations of erosional events so as to remove the time scale bias in estimated erosion rates.

Although estimated erosion rates in fluvial landscapes also likely bear the influence of climate variability, fluvial erosion rates are thought to increase nearly linearly with water discharge (75), which is a proxy of precipitation and, in turn, for climate; thus, the sensitivity of erosion to climate for these landscapes may be lower when compared to glaciated landscapes. Moreover, the impact of large floods may be lessened in fluvial systems because rivers can adjust their geometries to changing flows (76), leading to more uniform magnitudes of erosional pulses through time. Our numerical simulations suggest that spatial averaging of uncorrelated heavy-tailed erosional hiatuses over large spatial domains may reduce or eliminate the time scale bias in estimated erosion rates (Fig. 5C) (9). For example, the landscape-averaged erosion rates in Fig. 2B are in contrast to incision rates from dated fluvial terraces that show a time scale bias (6), which is consistent with estimated erosion rates from terraces being point measurements with minimal spatial averaging in contrast to the catchment-averaged erosion rates in Fig. 2B that appear to be insensitive to averaging time scale. Our results in Fig. 2 therefore support the possible interpretation of catchment-wide erosion rates in fluvial-dominated landscapes in terms of secular changes in climate and/or tectonic forcing (18, 20).

The power-law exponent that describes the time scale dependence of estimated erosion rates in different regions potentially yields additional insight into the intermittency of erosion (Fig. 2A). Power-law exponents closer to 0, for example, indicate less variability in erosional hiatuses (less intermittent) or a greater degree of spatial averaging of uncorrelated pulses of erosion. Power-law exponents closer to −1 indicate higher variability in erosional hiatuses (more intermittent) or spatially correlated pulses of erosion. In our current synthesis, there are no strong relationships between the power-law exponent and the expected thermal regime of glaciers (Fig. 3), suggesting that although the magnitude of erosional pulses may be higher in temperate glaciers when compared to subpolar glaciers (41, 64, 77), the intermittency in erosion may be comparable. However, the power-law exponent in landscapes shaped by tidewater glaciers is higher than that of other glaciers (Fig. 3), suggesting that erosion is more intermittent or spatially correlated in these regions, possibly because of additional marine controls, such as sea level and ocean temperatures. Moreover, landscapes with high uplift and erosion rates may respond more quickly to external forcing, which could reduce erosion variability over long time scales. Although this reasoning may hold for fluvial-dominated landscapes, glaciated landscapes with high uplift rates (for example, New Zealand) continue to show a time scale bias (albeit to a lesser degree; Fig. 3 and table S1) in estimated erosion rates over millions of years.

An additional source of complexity in interpreting temporal trends in estimated erosion rates from different thermochronometers is related to the nature of the thermal structure of Earth. Deeper isotherms are less sensitive to topography, whereas shallow isotherms are sensitive to spatial differences associated with topography or drainage divide migration (78). Thus, older thermochronological ages may inherently average over larger spatial domains. Moreover, thermochronology cannot resolve young, low-magnitude erosional pulses that do not exhume a sufficient depth of rock, potentially biasing global compilations of estimated erosion rates to be high for younger ages (11).

Our findings have implications for the ongoing debate on potential feedbacks between global climatic cooling during the Pleistocene and increases in erosion rates worldwide (2, 3, 5, 7, 11, 12). Previous workers have postulated that climate cooling and the onset of glaciation led to a global increase in erosion rates in glaciated landscapes (2, 3, 5). However, our results suggest that an increase in estimated erosion rates toward the present can be reasonably expected if the averaging time scale correlates with age and is shorter than the longest possible erosional hiatus. Although glaciated landscapes carry an unmistakable imprint of Pleistocene glaciations, our numerical simulations reveal that uneven sampling of heavy-tailed erosional hiatuses has the potential to dominate the relation between estimated erosion rate and age. Furthermore, the time scales of climate variability may overlap with tectonic time scales, which complicates untangling climatic versus tectonic effects on erosion rates in glaciated landscapes. Ultimately, our results suggest that it is important to characterize the stochastic nature of erosional events, their spatial extents, and the longest erosional hiatus to better interpret changes in estimated erosion rates with age.

MATERIALS AND METHODS

Estimated erosion rates from decades to millions of years

To assess the erosion rates on million-year time scales, we compiled thermochronological data. Thermochronology interprets isotopic ages in terms of thermal history (43), and accordingly provides a cooling age of a sample below a certain temperature. This temperature depends on the thermochronometric system applied and ranges from >900°C [(U-Th)/Pb (44)] to c. 60°C [apatite (U-Th)/He (43)]. We compiled available data from landscapes spanning both fluvial and glaciated landscapes (table S2), mostly Ar-Ar muscovite, apatite and zircon fission track, and (U-Th)/He dating, and computed the erosion rates using the depth to closure and the time since closure of each thermochronometer. Where data of surface elevation and temperature were available, we used the method proposed by Willett and Brandon (79). We then averaged the erosion rates estimated from the cooling ages of a particular thermochronometer to access million-year time scale erosion rates.

To assess millennial-scale estimated erosion rates, we used data of in situ produced terrestrial CRNs from river sediments. We compiled previously published estimates of basin-wide erosion rates (table S2). In situ produced terrestrial CRNs, which are rare isotopes that are exclusively produced by cosmic radiation in the uppermost meters of Earth’s surface, allow determining erosion rate estimates at time scales of 102 to 105 years (34–37). This approach is applied by dating landforms or deposits that record progressive incision (80, 81) or by estimating catchment-averaged erosion rates, derived from CRN abundances in river sediment (38, 39). Where available, we also used erosion rate estimates derived from lake and fjord sedimentation rates. The error on these erosion rates depended on the precision of mapping of the stratigraphic units in the lakes and fjords and was given in the original publications as c. 30% (41). These derived erosion rates average over entire catchments. Note that storage within the catchment, which might affect erosion rate estimates on decadal time scales, may be present in CRN dating, although sediment from large rivers like the Amazon had also been shown to preserve the erosion signal from the exhuming mountain range (82). To assess decadal-scale erosion rates, we used the reported values of erosion rates estimated from sediment yield data published in previous studies (see the Supplementary Materials).

We compiled data from landscapes for which estimates of erosion rates were available across varied time scales of integration, that is, sediment yield data, CRN-derived erosion rates or fjord sedimentation rates, and erosion rates derived from low-temperature thermochronometry. Erosion rates derived from sediment yield measurements and CRN concentrations in river sands were inherently spatially averaged over the catchment of interest; however, erosion rates derived from low-temperature thermochronometry were point measurements. To perform the comparison between estimated erosion rates across decadal, millennial, and million-year time scales, we averaged the erosion rates derived from thermochronology within the same area (figs. S5 and S6). When data within these locations were not available, we used the measurements from the nearest locations that overlapped with the erosion rate measurements on shorter time scales (see Supplementary Note 4 for additional notes on data compilation specific to each landscape). While averaging and presenting errors on the estimated erosion rates, we made sure that erosion rates derived from a single thermochronometer were used, that is, we averaged erosion rates derived from apatite fission track independently of erosion rates derived from zircon fission track, for example.

Power-law fitting

We performed linear regression analysis on the compiled data with their errors to determine the best-fitting power laws, where the data considered were the logarithms of estimated erosion rates and the averaging time scale. Further, we performed Monte Carlo error analysis for each data set where we computed 5000 linear fits between the logarithms of estimated erosion rates and the averaging time scales (83). These methods were used to determine the scaling exponent of the power-law fits along with their SE (reported in table S1 and Fig. 3). The observed negative power-law trend of estimated erosion rates on averaging time scale persisted over million-year time scales in all compiled glaciated landscapes except for Patagonia (Fig. 2A)—here, it appears that the millennial-scale erosion rate is comparable to the million-year time scale erosion rates in this landscape, suggesting that the upper bound of erosional hiatuses is on the order of several millennia. Our results indicated that the scaling exponent was statistically not different from zero for the San Gabriel Mountains, Taiwan, and the central Himalayas (table S1).

Numerical simulations of time series of erosion

We generated independent time series of erosion by drawing truncated Pareto random variables to represent erosional hiatuses and punctuated these hiatuses with a pulse of erosion of a given magnitude. Each erosional hiatus was attributed units of 1 year in time, which represents the minimum resolution of the generated time series. A Pareto random variable (xp) with a tail index α can be generated by raising a uniform random variable, u, to a power −α, that is, xp= u−α. A truncated Pareto random variable (xtp) with a tail index α and upper bound ν can be iteratively generated by retaining the Pareto random variable less than ν and discarding values greater than the upper bound, that is, xtp = {xp≤ν}. We generated a sufficiently long time series of erosion (>5 million years) and then excluded the first and last 500 ky of data (greater than twice of the upper bound of simulations shown here, that is, 200 ky). This not only avoids any initial or boundary condition effects but also allows for spatial points to be in the middle of a hiatus as opposed to always starting with an erosional pulse or a hiatus. The power-law fits to averaging time scale and the estimated erosion rate were then made provided that positive measurements of erosion rate were made at all averaging time scales, similar to data shown in Fig. 2. Further, to assess how different degrees of spatial averaging changed the averaging time scale dependence of estimated erosion rates, we fit power-law functions to estimated erosion rates for averaging time scales less than the upper bound on erosional hiatuses, that is, on the same time scales as for the time series without any spatial averaging.

Spectral analysis of the climate proxy using wavelets

We used wavelet analysis to characterize the multiscale variability of time-series δ18O values from well-preserved benthic foraminifera [a stack compiled in the study by Lisiecki and Raymo (65)], a commonly used high-resolution global climate proxy. Using Fourier analysis, previous workers have shown that the climate variability is coupled over diverse time scales ranging from seasonal to millennial time scales (66). Wavelet analysis is a localized multiscale analysis of a time series and can be more useful for time series that exhibit many scales of variability (84–86). We performed wavelet analysis on the climate proxy time series data with increasing order of derivatives of the Gaussian function (N = 1, 2, 3, …) until the computed wavelet power spectral density was the same for different orders of wavelets to ensure that our analysis was performed on stationary increments of the climate proxy data.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

,
Exhumation, crustal deformation, and thermal structure of the Nepal Himalaya derived from the inversion of thermochronological and thermobarometric data and modeling of the topography. J. Geophys. Res.115,
B06407 (2010).

,
Uplift and exhumation along the Arun River (eastern Nepal): Implications for the mechanism of uplift of the high Himalaya and the coupling between erosion and tectonics. Am. Geophys. Union Fall Meet.,
T31E-07 (2009).

Acknowledgments: We thank R. DiBiase, I. Larsen, and L. Malatesta for fruitful discussions and K. Whipple, F. Pazzaglia, R. Schumer, J. Braun, and anonymous reviewers for critical comments on an earlier version of this manuscript. Funding: We are grateful for financial support from Caltech (California Institute of Technology). Terrestrial Hazard Observation and Reporting Center (made possible by Foster and Coco Stanback) (M.P.L. and W.W.F.) and Tectonics Observatory (to C.v.H. and J.-P.A.) at Caltech provided additional financial support. National Center for Earth-Surface Dynamics 2 synthesis postdoctoral program, Imperial College London Junior Research Fellowship (V.G.), and Alexander von Humboldt Foundation (D.S.) are acknowledged for additional support. Author contributions: V.G., C.v.H., and D.S. compiled the data. M.P.L., V.G., and C.v.H. conceived this study. V.G. performed all the numerical simulations. V.G. and C.v.H. analyzed the compiled data and wrote the manuscript with input from other authors. All authors contributed to interpretation and discussion of the results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.