Sediment size and abrasion biases in detrital thermochronology

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.1 (1-12)

Earth and Planetary Science Letters ••• (••••) ••••••

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Sediment size and abrasion biases in detrital thermochronology Claire E. Lukens a,b,∗ , Clifford S. Riebe a , Leonard S. Sklar c,d , David L. Shuster e,f a

Department of Geology and Geophysics, University of Wyoming, Laramie, WY 82071, United States of America Department of Geography, Victoria University of Wellington, Wellington, 6140, New Zealand 1 c Department of Earth and Climate Sciences, San Francisco State University, CA 94132, United States of America d Department of Geography, Planning, and Environment, Concordia University, Montreal, Quebec, H3G 1M8, Canada 1 e Department of Earth and Planetary Sciences, University of California, Berkeley, CA 94720, United States of America f Berkeley Geochronology Center, Berkeley, CA 94709, United States of America b

a r t i c l e

i n f o

Article history: Received 12 June 2018 Received in revised form 22 October 2019 Accepted 27 October 2019 Available online xxxx Editor: A. Yin Keywords: geochronology thermochronometry U-Pb dating cosmogenic nuclides (U-Th)/He dating Inyo Creek

a b s t r a c t Detrital thermochronology has revolutionized the study of sediment provenance at orogen scales and spatial patterns in erosion at catchment scales. A strength of the method is that a handful of stream sand can be used to represent processes over landscape scales. However, it relies on both the widely recognized assumption that sand is supplied from every part of the catchment, and the implicit assumption that the fraction of eroded material that is sand-sized does not change across the landscape. These assumptions may be violated when the catchment contains spatial variations in the initial sizes of sediment produced by hillslope weathering or when abrasion during transport causes size reduction of sediment sourced from distal parts of the catchment. In either case, a detrital sample spanning a narrow range of sediment sizes (e.g., sand) may fail to represent the catchment as a whole, leading to bias in thermochronology of erosional and tectonic processes. We used forward modeling to quantify biases that can arise due to plausible abrasion rates and spatial variations in initial sediment size. Our results reveal significant reductions in the chance of detecting age populations originating from the highest, most distal parts of the landscape, leading to potentially erroneous interpretations in provenance studies. In tracer thermochronology, the biases distort detrital age distributions, leading to potentially profound misinterpretation of spatial patterns in erosion rates. Our analysis shows that the sediment size and abrasion biases increase with catchment area and relief but can be significant in small catchments (<10 km2 ) with moderate relief (>0.5 km). We show that the biases can be mitigated by analyzing a sufficient number of grains in multiple sediment size classes. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Detrital thermochronology has revolutionized studies of regional tectonics and mountain landscape evolution (Reiners and Brandon, 2006; Gehrels, 2014). Uranium-lead, fission-track, and (U-Th)/He ages in a handful of sand from a riverbed can reveal sediment provenance, providing insight into uplift and exhumation of catchments over the 104 – 108 yr timescales of tectonic processes (e.g., Fedo et al., 2003; Ruhl and Hodges, 2005; Lease et al., 2007; Naylor et al., 2015). Detrital ages can also reveal sediment source elevations and thus quantify spatial patterns in erosion rates over the 101 – 105 yr timescales of geomorphic

*

Corresponding author at: Department of Geography, Victoria University of Wellington, Wellington, 6140, New Zealand. E-mail address: [email protected] (C.E. Lukens). 1 Current address. https://doi.org/10.1016/j.epsl.2019.115929 0012-821X/© 2019 Elsevier B.V. All rights reserved.

processes (e.g., Stock et al., 2006; Vermeesch, 2007), providing insight into climatic control of sediment production (Riebe et al., 2015) and spatial variations in erosion by landsliding, glaciation, and rock fall (McPhillips and Brandon, 2010; Tranel et al., 2011; Ehlers et al., 2015; Fox et al., 2015; Lang et al., 2018; Whipp and Ehlers, 2019). An overarching assumption in many applications of detrital thermochronology is that a small sample of stream sediment can adequately represent the full population of ages in sediment that has been produced in the contributing area. However, this assumption may be violated in catchments with spatial variations in the production, transformation, and delivery of sediment from source to sampling point. Previous work on biases in detrital thermochronology has focused on issues specific to the datable minerals within the sediment. For example, the concentration and crystal size of the mineral of interest can vary in bedrock across a catchment (Moecher and Samson, 2006; Yang et al., 2012), leading to preferential sam-

JID:EPSL

2

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.2 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

pling of areas with high concentrations of analyzed sizes (Tranel et al., 2011; Avdeev et al., 2011). Likewise, preferential sampling can arise due to differences in mineral density and weathering susceptibility, which can reduce or enhance mineral survival during erosion from hillslopes and transport to the stream sampling point (e.g., Gleadow and Lovering, 1974; Lawrence et al., 2011; Sláma and Košler, 2012). Even a sample that is not affected by mineral-related biases can yield an incomplete perspective on sediment provenance if the number of analyzed grains is too small (Vermeesch, 2004; Avdeev et al., 2011). Numerous approaches have been proposed to minimize mineralrelated errors in detrital thermochronology. For example, to avoid biases due to spatial variations in mineral concentrations, studies select catchments with uniform lithology and either assume or confirm (by sampling) that the target mineral is uniformly distributed in bedrock (e.g., Stock et al., 2006). Biases due to differential mineral weathering are minimized by selecting highly resistant zircon or titanite as the target mineral (Gleadow and Lovering, 1974) or by sampling larger particles that protect less-resistant minerals such as apatite inside clast interiors (Reiners et al., 2007; Riebe et al., 2015). In addition, statistical guidelines have been developed for the minimum number of grains needed to quantify the full spectrum of ages in sediment provenance studies (Vermeesch, 2004). Even when study sites and sampling plans are carefully chosen to avoid mineral-related errors, applications of detrital thermochronology can be confounded by spatial variations in erosion rates over geomorphic timescales (Ruhl and Hodges, 2005). Areas with relatively high short-term erosion rates will be overrepresented in the sample, leading to misinterpretation of the catchment’s long-term exhumation history. The potential for this kind of bias can be evaluated by comparing the measured age distribution with catchment hypsometry. It has been proposed that a good match can be interpreted to reflect both steady exhumation of the landscape over millions of years and spatially uniform erosion over the shorter timescales of sediment production and delivery (Ruhl and Hodges, 2005). In landscapes where patterns in age distributions and catchment hypsometry do not match, detrital thermochronology can be used to quantify spatial variations in erosion rates (Stock et al., 2006; Vermeesch, 2007). Quantifying mismatches in detrital age and elevation distributions is the basis of powerful new applications of detrital thermochronology in geomorphology (McPhillips and Brandon, 2010; Tranel et al., 2011; Ehlers et al., 2015; Lang et al., 2018), including the simultaneous analysis of both temporal variations in exhumation rates and spatial variations in hillslope erosion rates using Bayesian inference (Avdeev et al., 2011). All of the applications outlined above rely on a previously underappreciated assumption—that the sampled sediment sizes are present in the same proportion in the size distribution of eroded material from all parts of the landscape. This assumption can be violated, for example, when finer sediment is preferentially produced at lower elevations (Fig. 1A-B). In such a scenario, a sample that only contains finer sediment (e.g., a grab sample of sand) would overrepresent the lower elevations in its age distribution (Fig. 1C). This is what we found at Inyo Creek in the Sierra Nevada, California, when we compared apatite ages from sand and fine gravel (Stock et al., 2006) with ages from coarse gravel (Riebe et al., 2015). The finer sediment had a markedly younger age distribution, suggesting that it is produced more readily at lower elevations (Riebe et al., 2015). This is one of many published examples of spatial variations in the initial size of sediment produced on hillslopes (e.g., Attal and Lavé, 2006; Whittaker et al., 2010; Marshall and Sklar, 2012; Attal et al., 2015; Genetti, 2017; Sklar et al., 2017; Roda-Boluda et al., 2018). Collectively, these results suggest that spatial variations in initial sedi-

Fig. 1. Conceptual model for grain-size and abrasion bias. If the median sediment size (D50 ) produced on hillslopes increases with elevation across a catchment, then size distributions at low (blue), middle (red), and high (green) elevations (A-C) will contribute differing amounts of sediment to a sand-sized sample (gray bar in C) at the catchment outlet (star in A). In this scenario, lower elevations will be overrepresented in the sample. Sediment abrasion creates silt-sized sediment that is not sampled, reducing the contribution of all sediment sizes in proportion to their distance upstream from the sampling point. The example shown in (D) illustrates the combined effects of abrasion and the altitudinal increase in sediment size shown in B. Spatial variations in erosion rates can also influence the contributions of different elevations to the sample. C and D assume uniform erosion (gray horizontal line in E). If erosion rates instead increase with elevation (black line in E), then the higher elevations will contribute more sediment to the sample (F), but can still be underrepresented in sampled sand (gray bar in F) both with abrasion (dashed lines) and without it (solid lines). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

ment size may be the norm rather than the exception in mountain landscapes due to variations in climatic, lithologic, and geomorphic factors that control the conversion of bedrock to mobile sediment by weathering (Riebe et al., 2017; Sklar et al., 2017). Thus a handful of sand—or any narrow size range—from a stream may commonly fail to represent patterns in erosion from the catchment upstream, particularly in mountain settings where initial sediment size may range from silt to boulders. The use of a narrow range of sediment sizes is further complicated by abrasion during transport from hillslope sediment sources to the sampling point in the stream (Fig. 1D): Age distributions can be distorted because abrasion preferentially affects sediment eroded from more distal portions of the catchment (Sklar et al., 2006; Attal and Lavé, 2009; Lavarini et al., 2018). To determine whether sediment size and abrasion biases are large enough to affect interpretations in detrital thermochronology, we quantified how they vary across a range of plausible scenarios using a forward model of sediment production, mixing, and abrasion. Our results, which build on our previous analysis of similar

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.3 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

biases in cosmogenic studies of erosion rates (Lukens et al., 2016), show that sediment abrasion and initial size variations can lead to substantial bias in detrital thermochronology. We discuss the implications of these biases for tectonic studies of sediment provenance and geomorphologic studies of sediment production. Our model results suggest that the biases can be reduced by increasing the number of analyses and by sampling multiple sediment sizes at the catchment outlet. 2. Forward modeling Our goal was to quantify the magnitude of sediment size and abrasion biases in detrital thermochronology and identify the conditions under which they may arise in landscapes. We focused on catchment area and relief as morphometric factors that are likely to regulate these biases and their implications for detrital thermochronology. Catchments with larger areas should have more bias due to sediment abrasion, because they have more source locations with longer travel distances to the outlet and thus more potential for abrasion during transport. Similarly, catchments with a wider distribution of elevations (i.e., greater relief) should have larger ranges in sediment size if size varies with elevation, leading to greater bias. The effects of area and relief are determined by the distribution of travel distances to the outlet and the distribution of elevations within the catchment, known as the width function and hypsometry, respectively. To isolate the effects of catchment area and relief from potentially confounding variations in catchment shape and tributary structure, which can also influence the width function and hypsometry (Sklar et al., 2016), we generated joint distributions of elevation and travel distance for synthetic catchments derived from landscape scaling relationships following the approach of Sklar et al. (2016), which is briefly summarized in Appendix A. By modifying parameters in the relationships (Equations (A.1)-(A.4); Table A.1), we are able to independently change area and relief while keeping the shape of the hypsometry and width functions roughly constant. We followed the same approach in a companion study of sediment size and abrasion biases in detrital cosmogenic nuclide measurements (Lukens et al., 2016). To gauge bias in a sampled age distribution from a catchment, we applied a forward model of sediment production, mixing, and abrasion to the catchment’s joint distribution of elevation and travel distance. In our model, each point in the catchment (i.e., each combination of elevation and travel distance) has an initial sediment size distribution, an erosion rate, and a bedrock age. We can therefore impose different spatial patterns in sediment production across slopes. For a sampling point at the catchment outlet, we calculate the distribution of sediment sizes and associated ages of sediment delivered from bedrock upstream according to the distribution of erosion rates and account for size reduction along the travel path to the outlet (Lukens et al., 2016). We simulate ages measured at the outlet by sampling sand from the modeled sediment mixture, mimicking protocols often employed in the field. In the analyses that follow, we focus on the age distribution in sand. We then quantify bias in two ways: by comparing ages from the sampled sand with the distribution in bedrock upstream and with the “actual” age distribution of all sediment produced on catchment hillslopes. To identify conditions under which significant bias may arise in mountain catchments, we considered a range of scenarios in which erosion rates and sediment sizes vary across hillslopes where sediment is produced. Numerous studies have shown that erosion rates depend on climate and topography (e.g., Montgomery and Brandon, 2002; Reiners et al., 2003; Moon et al., 2011), but the factors that influence initial sediment sizes on hillslopes have rarely been investigated (Sklar et al., 2017). Recent studies have shown that

3

both erosion rate (Riebe et al., 2015) and sediment size (Genetti, 2017) can increase exponentially with altitude in steep catchments due to climatic and topographic effects (e.g., Riebe et al., 2004; Attal et al., 2015). We therefore used an exponential relationship to model an altitudinal increase in median sediment size on hillslopes (Equation (A.5)), modulating the exponential scaling factor to explore the sensitivity of bias to a wide but plausible range of altitudinal gradients in initial sediment size (gray lines, Fig. 1B). Erosion rates either remained constant or increased exponentially with elevation (Equation (A.6)) according to a tunable scaling factor that kept erosion rates within a plausible range (Fig. 1E-F). To simulate size reduction during transport, we employed the Sternberg exponential function of travel distance from sources on slopes to the catchment outlet (Equation (A.7)), using abrasion coefficients that are consistent with field observations (Morris and Williams, 1999; Attal and Lavé, 2009). We assumed that the abrasion coefficient does not vary with particle size, and that the mass lost by abrasion becomes silt (Kodama, 1994; Attal and Lavé, 2009), and thus is not represented in the measured detrital sample of sand. We considered two distinct spatial patterns of variation in bedrock age across hillslopes: one for studies of sediment provenance, and the other for tracer thermochronology. 3. Quantifying bias 3.1. Bias in provenance studies In sediment provenance studies, the goal is often to obtain a representative inventory of bedrock ages in the catchment upstream based on a sample of detrital mineral ages. According to one widely cited prescription for the minimum number of measured ages, at least 117 single-grain analyses are needed to have 95% confidence that the measured age distribution includes all age components representing at least 5% of the total distribution (Vermeesch, 2004). However, this will only yield a representative age population of the catchment upstream if the sampled sediment is itself representative of sediment produced across the catchment. In some rivers, this may be the case (Copeland et al., 2015). However, as shown in Fig. 1, collecting sediment from a narrow size range—a common sampling strategy—may not yield representative samples, if sediment size varies across the catchment or if abrasion leads to preferential losses of distally sourced sediment. We explored this possibility using our forward model to quantify how the chance of observing a representative inventory of ages is affected by the number of analyses when sediment size and abrasion biases are present in a catchment. For simplicity, we simulated sediment production and delivery when there are just two distinct age populations, representing a lithologic contact at a constant elevation (Fig. 2A). Sediment from the higher-elevation population, which is also farther from the outlet, is more strongly affected both by abrasion during transport and by altitudinal variations in the sediment sizes produced on slopes. We therefore focused on quantifying the probability of observing the higherelevation population in sediment collected from the outlet. Our simulations spanned ranges in abrasion rates, coverage by the higher-elevation lithology, sediment size variation, and the number of ages measured. For simplicity, we assumed that erosion rates are spatially uniform. We varied the abrasion coefficient (Equation (A.7)) from 0 to 0.2 km−1 , where 0.1 km−1 yields moderate size reduction during transport (e.g., Attal and Lavé, 2009). In successive simulations, we set the fractional coverage by the higher lithology to 5, 10, and 15 percent of the total catchment area. Finally, the exponential scaling factor on the altitudinal increase in sediment size (Equation (A.5)) ranged from 0 to 2.7 km−1 , producing reasonable ranges in particle sizes across the catchment in each of our simulations. In the most extreme case modeled here,

JID:EPSL

4

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.4 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

Fig. 2. Forward modeling results for a catchment that has two distinct, normally distributed peaks in bedrock age representing chronologically distinct lithologies (A) that occur in upper and lower portions of the catchment (inset). Black line shows actual bedrock age distribution. Red line shows age distribution in sand at catchment mouth when median grain size increases from 1 to 100 mm across the catchment (with no sediment abrasion). Gray shaded area in A shows 95% confidence interval when n=117. Solid lines in B show how the probability of observing one or more ages from the upper population (left axis) depends on altitudinal increase in median sediment size (bottom axis) and the percentage of area underlain by the upper age population (contour labels). Dotted lines in B show combined effects of sediment abrasion and altitudinal increases in median sediment size for the same scenarios. Squares and circles denote examples discussed in the text. Top axis in B shows median sediment size at highest elevation of a 2.3 km-relief catchment when the initial median sediment size is 1 mm and the exponential scaling factor on the altitudinal increase in median sediment size is given by bottom axis. Input values are: relief = 2.3 km, area = 18 km2 , lower age = 25±2.5 My (mean±std. dev.), and upper age = 50±2.5 My.

a gradient of 2.7 km−1 over our modeled catchment (relief = 2.3 km) corresponds to an increase in median grain size (D 50 ) from 1 mm (sand) at the catchment outlet to ∼500 mm (boulders) at the divide, similar to our observations at Inyo Creek on the steep eastern flanks of the High Sierra (e.g., Riebe et al., 2015). When both the abrasion coefficient and the altitudinal increase in sediment size equaled 0 km−1 , we simulated the commonly assumed case of negligible abrasion and spatially uniform initial sediment size. We varied the number of randomly sampled ages (n) from 30 to 300 and replicated this random sampling process 10,000 times for each combination of parameter values in our model. The fraction of outcomes that produced at least one age from the higher-elevation population is the probability of observing that population in a random sample of n ages in sand. Across the range of conditions outlined above, the probability of observing one or more ages from the upper age population ranges from 0 to 1 (Fig. 2B). For example, the probability is 0.84 for n = 117 ages when abrasion is negligible, the higher-elevation population covers 5% of the catchment, and the exponential scaling factor on median sediment size is 1 km−1 , corresponding to an increase from 1 mm at the outlet to 10 mm at the highest ridge (solid square, Fig. 2B). In comparison, the probability is

Fig. 3. Chances of detecting a higher-elevation (upper) age population with different abrasion coefficients and different sample sizes for measured ages. The catchment modeled here has the same area and relief as in Fig. 2. Erosion rates are uniform (0.25 mm/yr everywhere). The upper age population occurs in the top 10% of the catchment. (A) Probability of observing the upper population using a range of abrasion coefficients (contour labels, with units of km−1 ). (B) Chance of detecting the upper age population for a range in the number of single grain analyses (n) from a sample of sand when the abrasion coefficient is 0.1 km−1 .

only 0.1 when median sediment size increases from 1 mm to 100 mm (solid circle, Fig. 2B). Increasing the fractional coverage of the higher-elevation population increases the probability of observing it in the sample. Adding sediment abrasion has the opposite effect, dropping the probability from 0.85 to 0.22 and from 0.1 to less than 0.01 in the examples described above (open symbols, Fig. 2B). In addition to variability in sediment production, we considered a range of abrasion coefficients and different numbers of analyses (Fig. 3). Increasing the altitudinal gradient in sediment size substantially reduces the chances of detecting the upper age population, even if abrasion is negligible and sample sizes are large (Fig. 3). The probability of observing a higher-elevation age is also sensitive to abrasion rate (Fig. 3A), decreasing to nearly 0 for all sediment production scenarios when the abrasion coefficient is 0.2 km−1 or greater. When sediment abrades readily during transport, the probability of detecting the upper age will be very low, regardless of whether sediment size varies. Conversely, when clasts are abrasion resistant, sediment from distal locations will survive the journey to the outlet, and the probability of detecting the upper age will be much higher. The chance of detecting the upper age also increases with n (Pullen et al., 2014), growing especially quickly when n is low (Fig. 3B). 3.2. Bias in tracer thermochronology In the relatively new field of tracer thermochronology, which focuses on quantifying spatial variations in erosion rates over geomorphic timescales, the goal is to obtain a distribution of ages that represents where sediment comes from (Stock et al., 2006). If erosion rates are spatially uniform, all points on the land-

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.5 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

5

scape are equally likely to contribute sediment to the stream, and the predicted age distribution can be obtained by combining catchment hypsometry with the bedrock age-elevation relationship (Brewer et al., 2003; Ruhl and Hodges, 2005). According to established protocols of tracer thermochronology, a statistically significant difference between the measured and predicted age distributions can be interpreted to reflect spatial variations in erosion rates (Stock et al., 2006). However, as we show next, spatial variations in initial sediment size and particle abrasion during transport—both of which may be common in mountain landscapes—may distort or obscure the true spatial pattern of erosion rates within a catchment (Riebe et al., 2015; Lukens et al., 2016; Lavarini et al., 2018). To explore these complications and identify the conditions under which they may arise, we simulated them in our forward model across a range of plausible scenarios. Bedrock ages increase linearly with elevation in all scenarios (Fig. 4A). Thus, we simulated typical conditions in tracer thermochronology: The age-elevation relationship in bedrock is known and can be used to detect spatial variations in erosion rates. Our analysis differs, however, in that we relaxed implicit assumptions about negligible abrasion and negligible variations in the initial sediment sizes produced on slopes. Our results show that when the altitudinal variation in initial sediment size is present without spatial variations in erosion and abrasion during transport, the age distribution in stream sand is shifted towards younger ages—i.e., derived from lower elevations— relative to the hypsometry (Fig. 4B). In the absence of a priori information about sediment size variations across the catchment, such a shift might be erroneously interpreted as evidence of faster erosion from low elevations, even though erosion rates are spatially uniform. Adding the complication of an altitudinal increase in erosion rates shifts the observed age distribution (“Sand” in Fig. 4C) towards older ages (i.e., higher elevations), but it still fails to capture the true spatial pattern in erosion rates (“Actual” in Fig. 4C). The observed pattern in ages will be misleading even when the altitudinal increase in erosion rates is very steep; under the circumstances illustrated in Fig. 4D, for example, the measured pattern of ages in sand might lead to the false conclusion that erosion rates are spatially uniform, even though they actually span an order of magnitude across the modeled catchment’s 2.0 km of relief. Adding abrasion shifts the curves in the opposite direction, towards younger ages relative to patterns shown in Figs. 4B-D, because of greater abrasion over long travel distances of sediment from more distal locations (which have higher elevations and thus older ages). 3.2.1. Quantifying the magnitude of bias In all simulations, an altitudinal increase in sediment size produces an offset between the measured age distributions in stream sand and actual age distributions in sediment produced on hillslopes. The size of the shift should depend on many factors. In general, the biases are bigger when the altitudinal increase in erosion rate is steeper and the abrasion coefficient is higher. The size of the bias should also depend on catchment area, because of its influence on the distribution of travel distances to the outlet. Bias should also vary with relief, because of its influence on the range in particle size for a fixed altitudinal gradient in initial sediment size (Equation (A.5)). To quantify the effects of area and relief on tracer thermochronology, we varied each factor separately while holding the other constant. We limited catchment area to 5–90 km2 , to evaluate the significance of the abrasion bias in catchments with relatively short travel distances to the outlet. In larger catchments, even the strongest rocks suffer substantial size reduction due to long travel distances (Sklar et al., 2006; Attal and Lavé, 2009). Likewise, we limited relief to 0.5–3 km, to focus on catchments where altitudinal gradients in initial sediment

Fig. 4. Influence of spatial patterns in erosion rate on age distributions measured in sand. (A) Model inputs in these examples: bedrock ages (color gradient on catchment map) increase linearly with elevation; median sediment size increases from 1 to 24 mm over 2.3 km of catchment relief (black line in plot); and sediment abrasion is omitted for simplicity. If erosion rates are spatially uniform (inset in B), then the age distribution in sand (red line) will reflect the young ages present at lower elevations in the catchment (B), while the age distribution produced on hillslopes (blue line) will lie directly on top of the distribution predicted by hypsometry (dashed black line). If erosion rates increase modestly with elevation (i.e., with Z in inset of C), then the age distribution produced on hillslopes will reflect higher erosion rates at high elevations (C), and neither the distribution in sand (red line) nor the distribution produced on hillslopes (blue line) will mimic the distribution predicted by hypsometry (dashed black line). If erosion rates increase steeply with elevation (i.e., with Z in inset of D), the age distribution in sand may mimic the distribution predicted by hypsometry (D), while the actual distribution of ages in eroded sediment reflects increasing erosion rates at higher elevations.

size could produce a substantial but realistic range of variations in sediment sizes produced on hillslopes. In addition to avoiding both trivial and extreme examples of sediment size and abrasion bias,

JID:EPSL

6

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.6 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

abrasion and altitudinal increases in erosion rate yield the largest values of the D statistic. D is less sensitive to relief when abrasion is included, because abrasion affects all catchments equally when area is constant (orange lines in Fig. 5A). We also explored how D varies with catchment area while holding relief constant (Fig. 5B). D is insensitive to area when the abrasion coefficient is zero, and the magnitude of D depends on the altitudinal gradients in both initial sediment size and erosion rate (cyan lines in Fig. 5B). When the abrasion coefficient is nonzero, D increases with area because longer travel distances lead to more size reduction of sediment from distal slopes (black and orange lines in Fig. 5B). When initial sediment size increases with elevation, the increase in D with area is less steep than when sediment size is spatially uniform (cf. orange and black lines in Fig. 5B), because catchments with small and large area share the same range in initial sediment size when relief is constant. The combined effects of sediment abrasion and altitudinal increases in both initial sediment size and erosion rate yield the largest values of D at the largest catchment areas (solid orange line in Fig. 5B). To identify the conditions needed for D to imply a statistically significant difference, we calculated p=0.05 significance thresholds of the K-S test for a range of sample sizes (gray lines in Figs. 5A and B). With a sample size of 30, a D statistic of 0.22 or more suggests a statistically significant difference between two populations. Nearly all the modeled scenarios exceed this criterion (Fig. 5), suggesting that the sediment size and abrasion biases may often be a confounding factor in tracer thermochronology of erosion patterns in mountain catchments.

Fig. 5. Magnitude of bias resulting from abrasion and altitudinal increases in sediment size and erosion rate across varying catchment relief and area. The D statistic from the K-S test is the vertical distance between the cumulative age distributions of ages in sand at the catchment mouth and ages produced on hillslopes, which reflects the magnitude of bias. In A, relief varies, and area is held constant at 26 km2 . In B, area varies, and relief is held constant at 2 km. Both panels show three scenarios: abrasion bias only (black); sediment size bias only (cyan); and both abrasion and sediment size biases (orange). Solid lines show amplifying effect of spatially varying erosion rates (relative to dashed lines that show uniform erosion case). Gray horizontal lines are critical D values for different numbers of ages in each distribution (n, contour labels); errors larger than these D values are statistically significant with 95% confidence.

our chosen ranges in area and relief are also consistent with the geometry of catchments studied using geomorphic applications of tracer thermochronology (e.g., Stock et al., 2006; Vermeesch, 2007; McPhillips and Brandon, 2010; Tranel et al., 2011; Ehlers et al., 2015). To quantify differences between the age distributions of stream sand and actual sediment produced on hillslopes, we used D, the Kolmogorov-Smirnov (K-S) test statistic, which is calculated as the maximum difference between the cumulative age distributions (“CAD”, after Vermeesch, 2007). Values of the D statistic range from 0, reflecting perfect agreement between two CADs, to 1, indicating no overlap in the range of the CADs. The D statistic thus quantifies the magnitude of error introduced by either (or both) of the abrasion and sediment size biases in our modeling scenarios. We first held area constant and varied relief (Fig. 5A). When sediment size is spatially uniform, the D statistic is insensitive to relief and depends only on the magnitude of the abrasion coefficient (black lines in Fig. 5A). When sediment size increases with elevation, D increases with increasing relief at a rate that depends on the abrasion coefficient and the spatial pattern in erosion rates (cyan and orange lines in Fig. 5A). The increase in D with relief is steeper, and D is larger (implying a more statistically significant difference), when erosion rates also increase with elevation (solid cyan and orange lines in Fig. 5A). The combined effects of sediment

3.2.2. Identifying the location of bias To correctly interpret measured age distributions and design studies that minimize bias, it is crucial to not only quantify the magnitude of bias but also identify which parts of the catchment are misrepresented. However, because the D statistic is based on cumulative distributions, it does not reveal which parts of the catchment are over- or underrepresented in a sample (i.e., the spatial distribution of bias). Instead it only shows whether the difference is large or small. We are unaware of a statistical metric that would quantify the spatial distribution of bias and how it varies with different sediment size and abrasion effects. Therefore, we developed a new metric that incorporates the measured, actual, and hypsometric age distributions, into a single value called the “focusing error” (F E).



FE =

i∈I

( S i − H i ) z∗i − R



i∈ J

( A i − H i ) z∗i (1)

Here, S i , H i , and A i and are the heights of the normalized age distributions of sand, hypsometry, and the actual distribution of ages in sediment produced on catchment slopes, respectively (Fig. 6). The subscript i refers to the number of a specific elevation bin (an integer) that varies from 1 to the total number of elevation bins; I is the index set of integers corresponding to positive values of S i − H i ; and J is the index set of integers corresponding to positive values of A i − H i . To calculate F E, we first convert the normalized empirical distribution function (EDF) of ages measured in sand at the outlet (S) and the EDF of ages actually produced on hillslopes ( A) into EDFs of elevation using the age-elevation relationship of the catchment bedrock. We then subtract from each of these EDFs the EDF of hypsometry (H ). The set notation in Equation (1) indicates that we only evaluate S i − H i and A i − H i at elevations where the hypsometry falls below the curve in question (Fig. 6). In Equation (1), z∗i = zi − zmedian , where zi is the midpoint of each elevation bin. Thus z∗i is negative when zi is below the median elevation (zmedian ) and positive when zi is above it. R is to-

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.7 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

7

Fig. 6. Conceptual model for focusing error (F E). Red lines show the age distribution in sand, blue lines show the actual age distribution of sediment produced on hillslopes, and black dashed lines show the age distributions predicted by hypsometry when erosion rates are assumed to be spatially uniform. Focusing error is calculated by summing positive differences between curves S (sand) and H (hypsometry) after weighting the differences by z*, the elevation relative to catchment midpoint, and subtracting the sum of all z*-weighted positive differences between curves A (actual) and H . Text at top right shows example calculation for scenario at top left. Panels A-F illustrate different patterns in the distributions S and A to show how differences in age data translate into differences in F E. Insets in panels A-F show the modeled patterns of erosion rate and sediment size variations that yield each result. While hypothetical, these results are all based on forward modeling using realistic distributions of elevation and travel distance and a wide range of reasonable spatial variations in erosion rates and sediment sizes.

tal catchment relief—a normalizing factor that keeps values of F E between −1 and 1 and therefore makes F E comparable between catchments regardless of differences in relief. The sign of F E reflects the spatial distribution of bias. The quantities (S i − H i ) and ( A i − H i ) are always positive, so the sign of each term in the sums of Equation (1) is determined by the sign of the local z∗i . If sampled sand originates mostly at eleva-

tions below the median, the first sum in the numerator will be negative (Fig. 6). Likewise, if the actual sediment produced on hillslopes comes mostly from high elevations, the second sum in the numerator will be negative, and the overall result will be negative (Fig. 6). In general, F E will be negative when upper elevations are underrepresented in sand (Fig. 6A-C) and positive when upper elevations are over-represented in sand (Fig. 6E-F). When the mea-

JID:EPSL

AID:115929 /SCO

8

[m5G; v1.261; Prn:14/11/2019; 9:24] P.8 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

sured and actual age distributions match, F E will be 0, even if they differ from hypsometry (e.g., Fig. 6D), indicating that stream sand accurately reflects the origins of the sediment across the catchment, with no associated bias in tracer thermochronology. Larger values of F E, whether more negative or more positive, indicate greater overall deviation between the measured and actual age distributions. The magnitude of F E also reflects the spatial distribution of bias. Because z∗i is effectively a moment arm in Equation (1); it weights departures from the hypsometry more when they are farther from the midpoint of the catchment. Thus F E will be larger when measured and actual age distributions are more strongly mismatched and when departures from hypsometry are farther apart, on opposite sides of the catchment midpoint. In Fig. 6, the largest F E occurs when both the measured and actual distributions differ from the hypsometry and have large differences in elevation between large departures (Fig. 6A). In contrast, F E is relatively small when departures are at similar locations (Fig. 6C-D). Two scenarios shown in Fig. 6 have similar F E but differing spatial patterns in erosion rates and initial sediment size (Fig. 6E-F). This shows that F E does not uniquely identify specific sediment production patterns that are responsible for observed spatial distributions in bias. To explore how the magnitude and spatial distribution of bias varies with relief and area, we calculated F E across the same range of conditions used to explore the magnitude of error using the D statistic (see Fig. 5). In all cases, F E was negative (Fig. 7), because of abrasion and the assumed pattern of increasing sediment size with elevation; both factors make upper elevations underrepresented in sampled stream sand. In scenarios with abrasion but without variation in initial sediment size, F E is insensitive to relief (black lines in Fig. 7A) but becomes increasingly negative with increasing area (black lines in Fig. 7B) for both uniform and increasing erosion rates (dotted and solid lines). Conversely, when sediment size increases with altitude but abrasion is negligible, F E becomes increasingly negative with increasing relief (cyan lines in Fig. 7A) and is insensitive to area both when erosion is uniform and when it increases with altitude (cyan lines in Fig. 7B). F E is largest (most negative) for a given relief and area when all of the factors (size reduction by abrasion and altitudinal gradients in both erosion rate and initial sediment size) are present (solid orange lines in Fig. 7A-B). Unlike the D statistic from the K-S test, our analysis of F E is limited by the lack a method for judging the statistical significance of F E values. However, F E is most likely to be significant when its absolute magnitude is large. Therefore, by showing how F E varies with relief and area (Fig. 7), our simulations illustrate the conditions under which thermochronology of detrital sand is at the greatest risk of bias due to the sediment size and abrasion effects. 4. Discussion Our analyses show that detrital thermochronology is susceptible to two previously underappreciated sources of bias—spatial variations in initial sediment size and size reduction by abrasion during transport. Where these effects occur, they can have implications for sediment provenance and tracing studies, even when erosion rates are spatially uniform across the contributing area. Here we consider whether these effects are common in landscapes and how sensitive predicted bias is to model assumptions. We also outline sampling strategies to reduce the risk of sediment size and abrasion biases in detrital thermochronology. 4.1. Prevalence of bias due to spatial variations in initial sediment size Although quantifying initial hillslope sediment size is a new area of research, available data suggest spatial variations may

Fig. 7. Focusing error (F E) resulting from abrasion and altitudinal increases in sediment size and erosion rate across varying catchment relief and area. In A, relief varies, and area is held constant at 26 km2 . In B, area varies, and relief is held constant at 2 km. Both panels show three scenarios: abrasion bias only (black); sediment size bias only (cyan); and both abrasion and sediment size biases (orange). Solid lines show amplifying effect of spatially varying erosion rates (relative to dashed lines that show uniform erosion case). F E is always negative, indicating an overrepresentation of the lower elevations and an underrepresentation of the higher elevations in measured ages from detrital sand, because erosion rate and grain size are either constant or increase with elevation in these scenarios. The most negative F E values occur in the largest, highest-relief catchments in which erosion rates increase with elevation and sediment sizes are reduced by abrasion during transport.

be common in mountain landscapes due to variations in lithology, erosion, and climate (Whittaker et al., 2010; Marshall and Sklar, 2012; Attal et al., 2015; Riebe et al., 2015; Roda-Boluda et al., 2018). For example, lithologies can differ in fracture spacing, mechanical strength, and chemical weathering susceptibility, leading to differences in initial sediment size (Sklar et al., 2017; Roda-Boluda et al., 2018). In addition, crustal deformation can drive spatial gradients in erosion rates, which cause differences in topographic slope, soil residence time, and sediment transport processes that correlate with observed variations in sediment sizes produced on hillslopes (Whittaker et al., 2010; Attal et al., 2015; Riebe et al., 2015; Sklar et al., 2017; Roda-Boluda et al., 2018). Spatial gradients in precipitation, temperature, and vegetation due to differences in elevation, for example, affect the relative importance of chemical, physical, and biological weathering that help explain altitudinal gradients in initial sediment size (Marshall and Sklar, 2012; Riebe et al., 2015; Sklar et al., 2017). Lithologic, erosional, and climatic gradients can be pronounced in mountain landscapes, suggesting that spatial variations in initial sediment size may be the norm rather than the exception. Although spatial variations in initial sediment size are likely common, there is no general theory for modeling sediment size

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.9 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

variations across landscapes (Sklar et al., 2017). In the absence of theory, we employed an exponential trend in sediment size with elevation (Equation (A.5)) and used scaling factors that produce sediment sizes consistent with field observations from Inyo Creek, California (Lukens, 2016; Genetti, 2017). Many other relationships are possible in mountain landscapes: talus production might result in a bimodal size distribution at high elevations, erosion rates might change dramatically across a fault, and changes in lithology might create patterns in sediment size that do not correlate with elevation. Although patterns in bias could differ from what we have presented here, the modeled scenarios are plausible and perhaps even conservative in our analysis of potential bias in detrital age distributions. Our results therefore suggest that a high risk of sediment size biases may be the norm rather than the exception in detrital thermochronology. Our analyses indicate that sediment size bias is most likely to be large in catchments with high relief and large area. These effects can be pronounced when both initial sediment size and erosion rate increase with elevation (e.g., Fig. 6A). Conversely, the two effects can be offsetting (e.g., Fig. 6C), for example, when higher elevations produce coarser sediment due to altitudinal variations in climate (Riebe et al., 2015) and lower elevations erode faster due to a transient wave of incision (Willenbring et al., 2013). However, such offsetting scenarios may be rare due to positive correlations between erosion rates and sediment size (Sklar et al., 2017). For example, faster erosion can lead to landsliding, which produces coarser sediment (Whittaker et al., 2010), while slower erosion leads to longer hillslope residence times and finer sediment sizes (Callahan et al., 2019). These interactions illustrate the importance of evaluating how catchment-scale variations in lithology, tectonics, and climate jointly influence the potential for bias due to spatial variations in initial sediment size. 4.2. Sensitivity of abrasion bias to modeling assumptions Particle size reduction during transport is ubiquitous in mountain landscapes (Attal and Lavé, 2009; Miller et al., 2014), and thus important to consider in detrital thermochronology. However, more work is needed to improve models of particle size reduction in rivers (Szabó et al., 2015). There are many sources of uncertainty in quantifying particle size reduction during transport, including: whether abrasion produces silt or sand; the relative importance of abrasion versus fracturing; and how rock properties influence size reduction by abrasion and fracturing. Most models assume that size reduction occurs due to abrasion resulting from particle collisions during transport (e.g., Lukens et al., 2016; Dingle et al., 2017), with the lost mass converted to silt or sand. We assumed that all mass is lost due to abrasion that produces silt, as is commonly observed in laboratory experiments (Kodama, 1994; Attal and Lavé, 2009), and that size reduction is an exponential function of distance (Morris and Williams, 1999). We neglected fracturing by high-energy collisions (Arabnia and Sklar, 2016), which is rarely modeled in studies of downstream fining in rivers due to a lack of data to parameterize it (see Le Bouteiller et al., 2011 for an exception). When abrasion produces silt rather than sand, abrasion products are not part of sand-sized samples collected from the channel. However, the sampled sand does include contributions from coarser clasts that are reduced to sand as mass is lost by abrasion during transport. Where lithologic and transport conditions favor production of sand rather than silt, the chance of sampling ages from more distal, higher elevations should increase (Lavarini et al., 2018), shifting the location of bias compared to when abrasion produces silt. In addition, because mass loss by abrasion scales with particle mass, bigger particles would contribute more sand to the sample, leading to overrepresentation of areas that produce coarser initial size distributions.

9

For simplicity, we assume that a single abrasion coefficient applies to all sizes. Available data and theory suggest that this is reasonable in mountain landscapes, except for fine sand at the lower limit of sizes susceptible to abrasion (Kendall, 1978; Wright and Smith, 1993; Wright et al., 1998). Published abrasion coefficients span three orders of magnitude (Sklar et al., 2006; Attal and Lavé, 2009). Our results show that just one order of magnitude of variation can introduce marked differences in bias (Fig. 3A). Lower abrasion coefficients might be appropriate for mono-mineralic grains of resistant minerals like zircon and titanite, or mica traveling as washload, which are likely to survive for long distances and be less susceptible to abrasion bias. We used a middle value for the coefficient (α = 0.1 km−1 ) in most of our simulations because most sediment on mountain hillslopes is multi-mineralic. Thus, our finding that abrasion biases can have a significant effect on age distributions in sediment should be robust and broadly representative of the scale of errors that are likely in mountain landscapes. 4.3. Strategies for reducing the sediment size and abrasion biases Our analysis shows that the abrasion and sediment size biases can confound interpretations of detrital thermochronology in mountain landscapes. However, our results also point to a number of strategies for reducing the potential for misinterpretation. In provenance studies, biases can be reduced by measuring more ages and thus increasing the probability of sampling the full range of ages present in the catchment (Fig. 3B). In tracer studies of geomorphic processes, additional measurements will make measured age distributions easier to distinguish from age distributions predicted from catchment hypsometry, because higher n values correspond to lower critical values of the D statistic in K-S tests. However, more measurements alone will not be enough to eliminate either the abrasion or sediment size bias, because they arise from distortion in the shape of the measured age distribution relative to the true age distribution of sediment produced on slopes. One strategy for identifying whether bias may be a problem at a field site is to analyze ages from multiple size classes and determine whether they are different. However, this will not be sufficient to account for bias, because it would be impossible to know a priori which (if any) distribution of ages represents the true spatial pattern of erosion in the catchment. It may be possible to account for bias by quantifying the effects of abrasion and sediment size variability across landscapes. To accomplish this goal, future studies will need to focus efforts on (i) quantifying sediment size reduction by both abrasion and fracturing and (ii) predicting spatial variability in sediment size across a wide range of landscapes (Sklar et al., 2017). 5. Concluding remarks We showed that detrital thermochronology is prone to two previously underappreciated sources of bias in studies that use a narrow range of sediment sizes sampled from streams: (i) spatial variations in the initial size distribution of sediment produced on hillslopes; and (ii) sediment abrasion during transport. When present, these biases can lead to misinterpretation in sediment provenance studies and in tracer thermochronology. In provenance studies, the age distribution measured in sand may not provide a complete inventory of bedrock ages. In tracer thermochronology, sediment abrasion and spatial variations in initial sediment size may obscure or distort inferred spatial patterns in erosion rates. In the examples presented here, bias was larger in higher-relief catchments where sediment size increases quickly with elevation and in larger-area catchments where abrasion reduces sediment size during transport.

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.10 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

10

Our analyses focused on underappreciated complications in traditional applications of detrital thermochronology, leading to the somewhat negative view that misinterpretation may be common in studies of mountain catchments. We conclude with a more optimistic view: that these complications can open new avenues in surface processes research. The factors that lead to the biases described here also reflect and influence a variety of processes that shape Earth’s surface. For example, sediment abrasion is a major factor in downstream fining in rivers (Miller et al., 2014), and the size distribution of sediment produced on hillslopes is a key regulator of channel geometry (Wohl, 2004), incision into bedrock (Sklar and Dietrich, 2004), and aquatic habitat (Riebe et al., 2014). With enough data from multiple sediment sizes, the sensitivity of detrital thermochronology to abrasion and spatial variations in initial sediment size could ultimately be used to quantify how these phenomena vary among catchments. The complications explored here therefore present opportunities for expanding knowledge of sediment production, transport, and size reduction processes, which are central to the interplay of tectonics, climate, and life in Earth surface dynamics and the sedimentary record. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work was supported by the National Science Foundation (EPS 1208909, EAR 1325033, and EAR 1331939); the David and Doris Dawdy Fund for Hydrologic Research; a Wyoming NASA Space Grant (NASA grant #NNX10AO95H); an American Geophysical Union Horton Research Grant; and the Ann and Gordon Getty Foundation. We thank Doug Burbank, Jeffrey Rahl, and Peter Copeland for helpful comments and suggestions on previous versions of the manuscript.

Equation (A.1) defines the lowest elevation for each travel distance in the stream channel. At travel distances on the long profile that fall above the fluvial part of the landscape, this relationship may not hold (e.g., Stock and Dietrich, 2003). We used a constant S h to approximate the channel slope in these reaches.

ks S h = θ ( L max − x)−θ H kA

(A.3)

This slope extends from the channel head, at L ch , to the catchment divide, at L max . We calculate the elevation of the unchanneled valley above the channel head according to Equation (A.4).





1−θ H 1−θ H + S h (x − xch ) zc = kc L max − L ch

(A.4)

Here, xch is the distance measured along the channel between the catchment divide and the channel head. Together, Equations (A.1) and (A.4) define the lowest elevation for each travel distance in the catchment. To define the highest elevation, we calculated a ridge profile using a power-law scaling relationship between travel distance and elevation. We then used a best-fit beta distribution to spread the area at each travel distance across the range of elevations between the channel and the ridge (Sklar et al., 2016). When applied across the entire catchment, this approach generates a modeled distribution of elevations at each travel distance. Integrated across all travel distances, these distributions yield the catchment hypsometry. Likewise, the distribution of travel distances can be integrated across all elevations to yield the catchment’s travel distance distribution. See Sklar et al. (2016) for details of this approach, and Lukens et al. (2016) for a recent application. To produce catchments with a range in area and relief, we vary only two parameters: k s , which determines the maximum elevation in the channel, and L max , the longest travel distance. The explored range of parameters (see Table A.1) produced realistic distributions of elevation and travel distance for mountain catchments ranging in size from 4 to 90 km2 , with relief ranging from 0.5 to 3 km.

Appendix A. Detailed methods A.2. Modeling changes in sediment size, erosion rate, and abrasion A.1. Generating distributions of elevation and travel distance To generate realistic distributions of elevation and travel distance for catchments spanning a range in area and relief, we followed the approach of Sklar et al. (2016), using landscape scaling relationships from the literature to calculate a channel long profile (Equation (A.1)), which defines the elevation of the channel (zc ) as a function of travel distance upstream of the outlet.



1−θ H zc = kc L max − ( L max − x)1−θ H



(A.1)

Here L max is the longest travel distance in the catchment, q and H are empirical constants, and the coefficient kc is calculated in Equation (A.2) using additional empirical constants k s and k A .

kc =

k s k−θ A

(A.2)

1−θH

Many of the details of our forward model are elaborated in Lukens et al. (2016). Here we describe how we varied three components of the model in our analysis of effects of bias: (i) altitudinal variations in the mean of the initial size-distribution of sediment produced on slopes; (ii) altitudinal variations in erosion rates of points on catchment slopes; and (iii) size reduction during transport from sediment sources on hillslopes to the sampling point at the catchment outlet. We set the sediment size distribution produced on hillslopes to be log-normal with a standard deviation of 2 φ units (consistent with Marshall and Sklar, 2012 and Genetti, 2017). For simplicity, this shape and spread in the modeled sediment size distribution does not change across the catchment. However, the median sediment size does change as a function of elevation according to Equation (A.5).

Table A.1 Parameters used to generate joint distributions of elevation and travel distance. Parameter

Value used or range explored

Description

ks kA L max q H L ch

3-25 1.28 5-30 km 0.31 1.75 600 m

Intercept coefficient for channel profile Hack’s Law coefficient Longest travel distance Channel curvature coefficient Hack’s Law exponent Distance from catchment divide to channel head

JID:EPSL

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.11 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

D i = D min ed( Z i − Z 0 )

(A.5)

Here, D i is the median grain size at point I , zi is the elevation at i, and D min is the median grain size at the lowest elevation (z0 ) in the catchment. We varied the exponential scaling factor, d, from 0 to 0.0028 m−1 , producing a realistic range of sediment sizes across different conditions (see main text). Note also that this parameter appears as the bottom horizontal axis in Figs. 2B and 3. This simple altitudinal increase in sediment size is supported by the limited empirical data that do exist (Riebe et al., 2015; Genetti, 2017), and the model allows for more complex variability that could be incorporated if additional empirical support arises in the future. Erosion rates are also modeled as an exponential function of elevation.

E i = aeb(zi −zavg )

(A.6)

Here E i is the erosion rate at each point, zavg is the median elevation (in km) in the catchment, a is the erosion rate at the median elevation (equal to 0.25 mm yr−1 in all models), and b is the scaling factor that modulates erosion rate across the catchment (0, 0.65 km−1 , or 1.3 km−1 ). The range of b used in our model produces plausible ranges in erosion rates across the range in relief that we considered in our analysis (e.g., Portenga and Bierman, 2011; Attal et al., 2015; Olen et al., 2015; Riebe et al., 2015; Lukens et al., 2016). Sediment abrasion during transport is modeled as a function of travel distance from the point of origin on hillslopes to the stream outlet (e.g., Kodama, 1994; Morris and Williams, 1999), according to Equation (A.7).

D p = D 0 e −α L

(A.7)

Here, D p is the diameter of the particle (mm) at the stream outlet, D 0 is the initial diameter, L is the travel distance (km), and α is the abrasion coefficient (ranging from 0 to 0.2 km−1 in models presented here). We assume that particle size is reduced via abrasion, and that the material produced by abrasion is less than 62 μm in diameter (and thus not part of samples collected from the channel). We also assume that all sediment sizes break down at the same rate. The abrasion coefficient used in Figs. 2, 4, 5, and 7 (0.1 km−1 ) is a moderate value for the abrasion of gravel (see main text). References Arabnia, O., Sklar, L.S., 2016. Experimental study of particle size reduction in geophysical granular glows. Int. J. Eros. Control Eng. 9, 122–129. https://doi.org/10. 13101/ijece.9.122. Attal, M., Lavé, J., 2006. Changes of bedload characteristics along the Marsyandi River (central Nepal): implications for understanding hillslope sediment supply, sediment load evolution along fluvial networks, and denudation in active orogenic belts. In: Geological Society of America Special Papers, vol. 398, pp. 143–171. Attal, M., Lavé, J., 2009. Pebble abrasion during fluvial transport: experimental results and implications for the evolution of the sediment load along rivers. J. Geophys. Res. 114, F04023. https://doi.org/10.1029/2009JF001328. Attal, M., Mudd, S.M., Hurst, M.D., Weinman, B., Yoo, K., Naylor, M., 2015. Impact of change in erosion rate and landscape steepness on hillslope and fluvial sediments grain size in the Feather River basin (Sierra Nevada, California). Earth Surf. Dyn. 3 (1), 201–222. https://doi.org/10.5194/esurf-3-201-2015. Avdeev, B., Niemi, N.A., Clark, M.K., 2011. Doing more with less: Bayesian estimation of erosion models with detrital thermochronometric data. Earth Planet. Sci. Lett. 305 (3–4), 385–395. https://doi.org/10.1016/j.epsl.2011.03.020. Brewer, I.D., Burbank, D.W., Hodges, K.V., 2003. Modeling detrital cooling-age populations: insights from two Himalayan catchments. Basin Res. 15 (3), 305–320. https://doi.org/10.1046/j.1365-2117.2003.00211.x. Callahan, R.P., et al., 2019. Arrested development: erosional equilibrium in the southern Sierra Nevada, California, maintained by feedbacks between channel incision and hillslope sediment production. Geol. Soc. Am. Bull. 131, 1179–1202. https://doi.org/10.1130/B35006.1.

11

Copeland, P., Bertrand, G., France-Lanord, C., Sundell, K., 2015. 40Ar/39Ar ages of muscovites from modern Himalayan rivers: Himalayan evolution and the relative contribution of tectonics and climate. Geosphere 11, 1837–1859. https:// doi.org/10.1130/GES01154.1. Dingle, E.H., Attal, M., Sinclair, H.D., 2017. Abrasion-set limits on Himalayan gravel flux. Nature 544, 471–474. https://doi.org/10.1038/nature22039. Ehlers, T.A., Szameitat, A., Enkelmann, E., Yanites, B.J., Woodsworth, G.J., 2015. Identifying spatial variations in glacial catchment erosion with detrital thermochronology. J. Geophys. Res., Earth Surf. 120 (6), 2014JF003432. https:// doi.org/10.1002/2014JF003432. Fedo, C.M., Sircombe, K.N., Rainbird, R.H., 2003. Detrital zircon analysis of the sedimentary record. Rev. Mineral. Geochem. 53 (1), 277–303. https://doi.org/10. 2113/0530277. Fox, M., Leith, K., Bodin, T., Balco, G., Shuster, D.L., 2015. Rate of fluvial incision in the Central Alps constrained through joint inversion of detrital 10 Be and thermochronometric data. Earth Planet. Sci. Lett. 411, 27–36. https://doi.org/10. 1016/j.epsl.2014.11.038. Gehrels, G., 2014. Detrital Zircon U-Pb geochronology applied to tectonics. Annu. Rev. Earth Planet. Sci. 42 (1), 127–149. https://doi.org/10.1146/annurev-earth050212-124012. Genetti, J.R., 2017. Evidence for Variations of Hillslope Sediment Size with Elevation at Inyo Creek. M.S. thesis. San Francisco State University. 144 pp. Gleadow, A.J.W., Lovering, J.F., 1974. The effect of weathering on fission track dating. Earth Planet. Sci. Lett. 22 (2), 163–168. https://doi.org/10.1016/0012-821X(74) 90077-6. Kendall, K., 1978. The impossibility of comminuting small particles by compression. Nature 272, 710–711. https://doi.org/10.1038/272710a0. Kodama, Y., 1994. Experimental study of abrasion and its role in producing downstream fining in gravel-bed rivers. J. Sediment. Res. 64 (1), 76–85. Lang, K.A., Ehlers, T.A., Kamp, P.J.J., Ring, U., 2018. Sediment storage in the Southern Alps of New Zealand: new observations from tracer thermochronology. Earth Planet. Sci. Lett. 493, 140–149. https://doi.org/10.1016/j.epsl.2018.04.016. Lavarini, C., Attal, M., da Costa Filho, C.A., Kirstein, L.A., 2018. Does pebble abrasion influence detrital age population statistics? A numerical investigation of natural data sets. J. Geophys. Res., Earth Surf. 123, 2577–2601. https://doi.org/10.1029/ 2018JF004610. Lawrence, R.L., Cox, R., Mapes, R.W., Coleman, D.S., 2011. Hydrodynamic fractionation of zircon age populations. Geol. Soc. Am. Bull. 123 (1–2), 295–305. https:// doi.org/10.1130/B30151.1. Le Bouteiller, C., Naaim-Bouvet, F., Mathys, N., Lavé, J., 2011. A new framework for modeling sediment fining during transport with fragmentation and abrasion. J. Geophys. Res., Earth Surf. 116 (F3), F03002. https://doi.org/10.1029/ 2010JF001926. Lease, R.O., Burbank, D.W., Gehrels, G.E., Wang, Z., Yuan, D., 2007. Signatures of mountain building: detrital zircon U/Pb ages from northeastern Tibet. Geology 35 (3), 239–242. https://doi.org/10.1130/G23057A.1. Lukens, C.E., 2016. Spatial Variation in Sediment Production Across Mountain Catchments: Insights from Cosmogenic Nuclides and Tracer Thermochronometry. Ph.D. Dissertation. University of Wyoming. 189 pp. Lukens, C.E., Riebe, C.S., Sklar, L.S., Shuster, D.L., 2016. Grain size bias in cosmogenic nuclide studies of stream sediment in steep terrain. J. Geophys. Res., Earth Surf. 121 (5), 2016JF003859. https://doi.org/10.1002/2016JF003859. Marshall, J.A., Sklar, L.S., 2012. Mining soil databases for landscape-scale patterns in the abundance and size distribution of hillslope rock fragments. Earth Surf. Process. Landf. 37, 287–300. https://doi.org/10.1002/esp.2241. McPhillips, D., Brandon, M.T., 2010. Using tracer thermochronology to measure modern relief change in the Sierra Nevada, California. Earth Planet. Sci. Lett. 296 (3–4), 373–383. https://doi.org/10.1016/j.epsl.2010.05.022. Miller, K.L., Szabó, T., Jerolmack, D.J., Domokos, G., 2014. Quantifying the significance of abrasion and selective transport for downstream fluvial grain size evolution. J. Geophys. Res., Earth Surf. 119, 2412–2429. https://doi.org/10.1002/ 2014JF003156. Moecher, D.P., Samson, S.D., 2006. Differential zircon fertility of source terranes and natural bias in the detrital zircon record: implications for sedimentary provenance analysis. Earth Planet. Sci. Lett. 247 (3–4), 252–266. https://doi.org/10. 1016/j.epsl.2006.04.035. Montgomery, D.R., Brandon, M.T., 2002. Topographic controls on erosion rates in tectonically active mountain ranges. Earth Planet. Sci. Lett. 201, 481–489. Moon, S., Page Chamberlain, C., Blisniuk, K., Levine, N., Rood, D.H., Hilley, G.E., 2011. Climatic control of denudation in the deglaciated landscape of the Washington Cascades. Nat. Geosci. 4, 469–473. https://doi.org/10.1038/ngeo1159. Morris, P.H., Williams, D.J., 1999. A worldwide correlation for exponential bed particle size variation in subaerial aqueous flows. Earth Surf. Process. Landf. 24 (9), 835–847. https://doi.org/10.1002/(SICI)1096-9837(199908)24:9<835::AIDESP15>3.0.CO;2-G. Naylor, M., Sinclair, H.D., Bernet, M., van der Beek, P., Kirstein, L.A., 2015. Bias in detrital fission track grain-age populations: implications for reconstructing changing erosion rates. Earth Planet. Sci. Lett. 422, 94–104. https://doi.org/10. 1016/j.epsl.2015.04.020. Olen, S.M., Bookhagen, B., Hoffmann, B., Sachse, D., Adhikari, D.P., Strecker, M.R., 2015. Understanding erosion rates in the Himalayan orogen: a case study from

JID:EPSL

12

AID:115929 /SCO

[m5G; v1.261; Prn:14/11/2019; 9:24] P.12 (1-12)

C.E. Lukens et al. / Earth and Planetary Science Letters ••• (••••) ••••••

the Arun Valley. J. Geophys. Res., Earth Surf. 120 (10), 2014JF003410. https:// doi.org/10.1002/2014JF003410. Portenga, E.W., Bierman, P.R., 2011. Understanding Earth’s eroding surface with 10 Be. GSA Today 21, 4–10. https://doi.org/10.1130/G111A.1. Pullen, A., Ibáñez-Mejía, M., Gehrels, G.E., Ibáñez-Mejía, J.C., Pecha, M., 2014. What happens when n = 1000? Creating large-n geochronological datasets with LAICP-MS for geologic investigations. J. Anal. At. Spectrom. 29 (6), 971. https:// doi.org/10.1039/c4ja00024b. Reiners, P.W., Ehlers, T.A., Mitchell, S.G., Montgomery, D.R., 2003. Coupled spatial variations in precipitation and long-term erosion rates across the Washington Cascades. Nature 426, 645–647. https://doi.org/10.1038/nature02111. Reiners, P.W., Brandon, M.T., 2006. Using thermochronology to understand orogenic erosion. Annu. Rev. Earth Planet. Sci. 34 (1), 419–466. https://doi.org/10.1146/ annurev.earth.34.031405.125202. Reiners, P.W., Thomson, S.N., McPhillips, D., Donelick, R.A., Roering, J.J., 2007. Wildfire thermochronology and the fate and transport of apatite in hillslope and fluvial environments. J. Geophys. Res., Earth Surf. 112 (F4), F04001. https:// doi.org/10.1029/2007JF000759. Riebe, C.S., Kirchner, J.W., Finkel, R.C., 2004. Sharp decrease in long-term chemical weathering rates along an altitudinal transect. Earth Planet. Sci. Lett. 218, 421–434. https://doi.org/10.1016/S0012-821X(03)00673-3. Riebe, C.S., Sklar, L.S., Overstreet, B.T., Wooster, J.K., 2014. Optimal reproduction in salmon spawning substrates linked to grain size and fish length. Water Resour. Res. 50, 898–918. https://doi.org/10.1002/2013WR014231. Riebe, C.S., Sklar, L.S., Lukens, C.E., Shuster, D.L., 2015. Climate and topography control the size and flux of sediment produced on steep mountain slopes. Proc. Natl. Acad. Sci. 112 (51), 15574–15579. https://doi.org/10.1073/pnas. 1503567112. Riebe, C.S., Hahm, W.J., Brantley, S.L., 2017. Controls on deep critical zone architecture: a historical review and four testable hypotheses. Earth Surf. Process. Landf. 42, 128–156. http://doi.wiley.com/10.1002/esp.4052. Roda-Boluda, Duna C., D’Arcy, Mitch, McDonald, Jordan, Whittaker, Alexander C., 2018. Lithological controls on hillslope sediment supply: insights from landslide activity and grain size distributions. Earth Surf. Process. Landf. 43, 956–977. https://doi.org/10.1002/esp.4281. Ruhl, K.W., Hodges, K.V., 2005. The use of detrital mineral cooling ages to evaluate steady state assumptions in active orogens: an example from the central Nepalese Himalaya. Tectonics 24 (4), TC4015. https://doi.org/10.1029/ 2004TC001712. Sklar, L.S., Dietrich, W.E., 2004. A mechanistic model for river incision into bedrock by saltating bed load. Water Resour. Res. 40, W06301. https://doi.org/10.1029/ 2003WR002496. Sklar, L.S., Dietrich, W.E., Foufoula-Georgiou, E., Lashermes, B., Bellugi, D., 2006. Do gravel bed river size distributions record channel network structure? Water Resour. Res. 42 (6), W06D18. https://doi.org/10.1029/2006WR005035. Sklar, L.S., Riebe, C.S., Lukens, C.E., Bellugi, D., 2016. Catchment power and the joint distribution of elevation and travel distance to the outlet. Earth Surf. Dyn. 4, 799–818. https://doi.org/10.5194/esurf-4-799-2016.

Sklar, L.S., Riebe, C.S., Marshall, J.A., Genetti, J., Leclere, S., Lukens, C.L., Merces, V., 2017. The problem of predicting the size distribution of sediment supplied by hillslopes to rivers. Geomorphology 277, 31–49. https://doi.org/10.1016/j. geomorph.2016.05.005. Sláma, J., Košler, J., 2012. Effects of sampling and mineral separation on accuracy of detrital zircon studies. Geochem. Geophys. Geosyst. 13 (5), Q05007. https:// doi.org/10.1029/2012GC004106. Stock, G.M., Ehlers, T.A., Farley, K.A., 2006. Where does sediment come from? Quantifying catchment erosion with detrital apatite (U-Th)/He thermochronometry. Geology 34 (9), 725–728. https://doi.org/10.1130/G22592.1. Stock, J., Dietrich, W.E., 2003. Valley incision by debris flows: evidence of a topographic signature. Water Resour. Res. 39, WR001057. https://doi.org/10.1029/ 2001WR001057. Szabó, T., Domokos, G., Grotzinger, J.P., Jerolmack, D.J., 2015. Reconstructing the transport history of pebbles on Mars. Nat. Commun. 6, 1–7. https://doi.org/10. 1038/ncomms9366. Tranel, L.M., Spotila, J.A., Kowalewski, M.J., Waller, C.M., 2011. Spatial variation of erosion in a small, glaciated basin in the Teton Range, Wyoming, based on detrital apatite (U-Th)/He thermochronology. Basin Res. 23 (5), 571–590. https:// doi.org/10.1111/j.1365-2117.2011.00502.x. Vermeesch, P., 2004. How many grains are needed for a provenance study? Earth Planet. Sci. Lett. 224 (3–4), 441–451. https://doi.org/10.1016/j.epsl.2004.05.037. Vermeesch, P., 2007. Quantitative geomorphology of the White Mountains (California) using detrital apatite fission track thermochronology. J. Geophys. Res., Earth Surf. 112, F03004. https://doi.org/10.1029/2006JF000671. Whipp, D.M., Ehlers, T.A., 2019. Quantifying landslide frequency and sediment residence time in the Nepal Himalaya. Sci. Adv. 5, eaav3482. https://doi.org/10. 1126/sciadv.aav3482. Whittaker, A.C., Attal, M., Allen, P.A., 2010. Characterising the origin, nature and fate of sediment exported from catchments perturbed by active tectonics. Basin Res. 22, 809–828. https://doi.org/10.1111/j.1365-2117.2009.00447.x. Willenbring, J.K., Gasparini, N.M., Crosby, B.T., Brocard, G., 2013. What does a mean mean? The temporal evolution of detrital cosmogenic denudation rates in a transient landscape. Geology 41, 1215–1218. https://doi.org/10.1130/G34746.1. Wohl, E., 2004. Limits of downstream hydraulic geometry. Geology 32 (10), 897–900. https://doi.org/10.1130/G20738.1. Wright, J., Smith, B., 1993. Fluvial comminution and the production of loess-sized quartz silt: a simulation study. Geogr. Ann., Ser. A, Phys. Geogr. 75A (1–2), 25–34. https://doi.org/10.2307/521050. Wright, J., Smith, B., Whalley, B., 1998. Mechanisms of loess-sized quartz silt production and their relative effectiveness: laboratory simulations. Geomorphology 23, 15–34. https://doi.org/10.1016/S0169-555X(97)00084-6. Yang, S., Zhang, F., Wang, Z., 2012. Grain size distribution and age population of detrital zircons from the Changjiang (Yangtze) River system, China. Chem. Geol. 296, 26–38.