the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Optimising global landscape evolution models with ^{10}Be
Gregory A. Ruetenik
John D. Jansen
Pedro Val
Lotta YläMella
By simulating erosion and deposition, landscape evolution models (LEMs) offer powerful insights into Earth surface processes and dynamics. Streampowerbased LEMs are often constructed from parameters describing drainage area (m), slope (n), substrate erodibility (K), hillslope diffusion (D), and a critical drainage area (A_{c}) that signifies the downslope transition from hillslope diffusion to advective fluvial processes. In spite of the widespread success of such models, the parameter values are highly uncertain mainly because the advection and diffusion equations amalgamate physical processes and material properties that span widely differing spatial and temporal scales. Here, we use a global catalogue of catchmentaveraged cosmogenic ^{10}Bederived denudation rates with the aim to optimise a set of LEMs via a Monte Carlobased parameter search. We consider three model scenarios: advectiononly, diffusiononly, and an advection–diffusion hybrid. In each case, we search for a parameter set that best approximates denudation rates at the global scale, and we directly compare denudation rates from the modelled scenarios with those derived from ^{10}Be data. We find that optimised ranges can be defined for many LEM parameters at the global scale. In the absence of diffusion, n∼1.3, and with increasing diffusivity the optimal n increases linearly to a global maximum of n∼2.3. Meanwhile, we find that the diffusiononly model yields a slightly lower misfit when comparing model outputs with observed erosion rates than the advectiononly model and is optimised when the concavity parameter is raised to a power of 2. With these examples, we suggest that our approach provides baseline parameter estimates for largescale studies spanning long timescales and diverse landscape properties. Moreover, our direct comparison of modelpredicted versus observed denudation rates is preferable to methods that rely upon catchmentscale averaging or amalgamation of topographic metrics. We also seek to optimise the K and D parameters in LEMs with respect to precipitation and substrate lithology. Despite the potential bias due to factors such as lithology, these optimised models allow us to effectively control for topography and specifically target the relationship between denudation and precipitation. All models suggest a general increase in exponents with precipitation in line with previous studies. When isolating K under globally optimised models, we observe a positive correlation between K or D and precipitation > 1500 mm yr^{−1}, plus a local maximum at ∼300 mm yr^{−1}, which is compatible with the longstanding hypothesis that semiarid environments are among the most erodible.
 Article
(7029 KB)  Fulltext XML

Supplement
(1700 KB)  BibTeX
 EndNote
To appreciate shortterm changes in Earth surface processes, such as those induced by humans (Brown, 1981; Hooke, 2000), it is first necessary to understand longterm rates of denudation and deposition. Recognising this, some recent studies (e.g. Simoes et al., 2010) have derived erosiontransport rules from topography with an aim to predict macroscale patterns of denudation and sediment flux. At more restricted spatial scales, denudation rates based on cosmogenic nuclides (e.g. ^{10}Be) show a modest exponential correlation with catchmentaveraged slope, as does normalised steepness in stream profiles (Portenga and Bierman, 2011; Harel et al., 2016). Nevertheless, it is widely observed that steepness and stream power parameters are subject to considerable variation wherever climate and/or lithology differ (Harel et al., 2016; Gailleton et al., 2021; Marder and Gallen, 2023) and that the parameters also covary. A robust analysis must accommodate such interactions.
Earth's surface undergoes continuous modification through uplift and denudation over timescales too long to observe directly; hence landscape evolution models (LEMs) are crucial tools for building knowledge. Note that we use the term “denudation” to denote the mass loss leading to the lowering of Earth's surface. LEMs are often employed over expansive scales of space and time in order to study topographic response to changes in tectonics (e.g. Kooi and Beaumont, 1996; GarciaCastellanos et al., 2003), climate (e.g. Temme et al., 2009, Adams et al., 2020), and sea level (e.g. Pico et al., 2019; Ruetenik et al., 2019). And yet, large spatial and temporal scales require generalisation of model parameters that reliably accounts for processes of hillslope diffusion and advective fluvial erosion. Using LEMs to estimate denudation rates delivers the key advantage of bridging scales and defining an empirically derived mechanism at the local (grid cell) scale. This demands that denudation rates are integrated over scales matching the topographic changes they describe. At the local to regional scale, recent studies have focused on constraining LEM parameters via inversions that optimise rates of denudation, deposition, and topographic observations (e.g. Miller et al., 2013; Croissant and Braun, 2014; Pedersen et al., 2018; Barnhart et al., 2020). However, implementing many of these approaches at a global scale is challenging in terms of computational cost and because it often requires a compilation of a large set of observables (such as knickpoints, depositional patterns, and denudation rates). In the absence of computational power that can accurately simulate stratigraphy at the global scale, and without constraints on global palaeotopography, we settle for optimising LEM parameters with respect to catchmentaveraged denudation rates estimated with cosmogenic ^{10}Be. While LEM parameters can be estimated via direct topographic analysis (e.g. Wobus et al., 2006; Clubb et al., 2014; Mudd et al., 2018), this approach can lead to bias, as we discuss below.
Here, we determine LEM parameter values that minimise the variance among ^{10}Bederived apparent denudation rates in the OCTOPUS v.2 global catalogue (Codilean et al., 2022), and we analyse the capacity of LEMs to predict the denudation rate given those optimised parameters. Our LEM employs the common streampowerplusdiffusion formulation, which is subject to important limitations, such as the neglect of fluvial deposition and mass wasting processes (e.g. Whipple and Tucker, 1999). The tradeoffs involved in this simplified approach, we believe, are justified by the record of success with simulating landscape processes at large scales and across a wide range of lithologies, drainage areas, and steepnesses (e.g. Gallen et al., 2013; Miller et al., 2013; Fox et al., 2014).
Catchmentaveraged denudation rates from cosmogenic ^{10}Be
Rates of catchmentscale denudation can be estimated by measuring the abundances of cosmogenic radionuclides, such as ^{10}Be, in quartzbearing river sand (Granger et al., 1996; Von Blanckenburg, 2005). Such nuclides accumulate within minerals exposed to secondary cosmic rays in the upper few metres of the bedrock subsurface and are lost via erosion and radioactive decay (Lal, 1991). The attenuation of cosmic rays with depth causes the nuclide production rate to decrease exponentially (at 2 m depth the ^{10}Be production rate is <5 % that at the surface); hence, nuclide abundances measured in sediment are an inverse function of denudation rate.
The spatial variations observed in denudation rates across a range of climates and lithologies (Portenga and Bierman, 2011; Starke et al., 2020) suggest that the erosional processes driving the evolution of landscapes also vary. This has important implications for the interpretation of ^{10}Bederived denudation rates and how we parameterise LEMs. Estimating catchmentaveraged denudation rates from nuclide abundances in river sand depends on a multitude of assumptions (Von Blanckenburg, 2005; Mudd, 2016). For example, we assume that sediments were produced via longterm, steady bedrock erosion distributed uniformly across the catchment and that sediments have experienced continuous exposure to cosmic rays at/near the surface. In detail, longterm steady erosion refers to at least one attenuation length (∼0.6 m) of surface lowering integrated over a 10^{3}–10^{5}year timescale. Abrupt bedrock erosion events, for instance, caused by bedrock landsliding or glacial quarrying may bias denudation rate estimates. Similarly, long intervals of ice cover or intermittent deep sediment burial contradict the requirement for continuous cosmicray exposure. Other sources of potential discord relate to lithology, catchment size, and hypsometry, which are known to affect sediment transport dynamics and grainsize yields (Carretier et al., 2015; Riebe et al., 2015; Lukens et al., 2016; Zavala et al., 2020). The sources of deviation noted above are collectively responsible for the considerable variability observed in large compilations of ^{10}Bederived denudation rates (e.g. Portenga and Bierman, 2011; Harel et al., 2016). Catchmentwide denudation rates are commonly determined and published for settings that do not strictly comply with the method's premises; these estimates are best referred to as apparent denudation rates (Mudd, 2016). Nevertheless, ^{10}Bederived data currently offer the most widely distributed insight into longterm denudation on a global scale.
2.1 Stream power and hillslope diffusion
Stream power is represented by a nonlinear advection equation derived from observations of river morphology and generalised relationships for bed shear stress (e.g. Howard et al., 1994; Whipple and Tucker, 1999; Lague, 2014). It affords a description of channel incision as a function of upstream drainage area (A) and local slope (S) for the portion of the catchment (above the critical drainage area, A_{c}) where fluvial advection dominates over hillslope diffusion and debris flow processes (e.g. Lague and Davy, 2003; Fontana and Marchi, 2003; Whipple and Tucker, 1999). The stream power equation takes the form
where K is the advection coefficient or erosional efficiency that is (in our formulation) uniform across a catchment, and m and n determine the relative dependence of incision on drainage area and slope. The ratio $m/n$ defines how elevation scales with drainage area at a steady state relative to the longitudinal channel profile and typically varies from 0.3 to 0.6 (Wobus et al., 2006; Whipple and Tucker, 1999); n determines the erosional nonlinearity of the model, which is thought to be determined by regional discharge variability, as well as processes governing incision thresholds, and typically ranges from 1 to 4 (e.g. Lague, 2014). A global compilation of stream power parameters constrained by topographic metrics (Harel et al., 2016) reports an optimised n∼2.6, albeit with considerable scatter (also observed by Gailleton et al., 2021). The value of n may also vary depending on the location in the channel network; the steepest and fastesteroding locales such as knickpoints can have values closer to unity (Lague, 2014). In general, higher n results in larger erosional flux from steep terrain, while higher m results in larger erosional flux from big rivers. However, if m and n both increase (keeping their ratio and K constant), a larger fraction of erosional flux will be sourced from steeper mainstem channels.
The amalgamated outcomes of hillslope transport processes, such as rain splash, soil creep, and bioturbation, are primarily diffusive. Hence, to simulate hillslope processes, we include a diffusion equation:
where D is diffusivity, which is reported to range from $\sim \mathrm{4.4}\times {\mathrm{10}}^{\mathrm{4}}$ to $\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{2}}$ for linear diffusion (e.g. Fernandes and Dietrich, 1997). We use an exponentiated form of the diffusion equation in which concavity is raised to an exponent, p, in order to harmonise with Gabet et al. (2021), who posit that denudation rate scales approximately with hillslope concavity squared. For linear diffusion to satisfy mass conservation, the deposited and eroded sediment should be balanced. However, in the exponentiated formulation, negative concavities raised to a power of p can produce nonreal components. Thus, we calculate an average catchmentwide denudation rate in which we ignore deposited sediment (areas of negative concavity) and take an average based on eroded sediment only.
For our joint advection–diffusion model, we follow the common approach of combining stream power with linear diffusion (i.e. p=1 in Eq. 2). When Eqs. (1) and (2) are combined, D and K covary, and it has been noted that a higher $D/K$ ratio results in a lower sensitivity of erosion rate to catchmentaveraged slope (Forte et al., 2016a). Hence, we divide Eq. (3) by K, which allows the $D/K$ ratio to be optimised with respect to predicted erosion rate:
With our advection–diffusion model formulation (Eq. 3), we set out to solve simultaneously for globally optimised values of $D/K$, n, and A_{c}. $D/K$ and n determine the relative importance of advective versus diffusive processes in driving erosion; lower $D/K$ and higher n (which implies higher m given uniform $m/n$) will result in the dominance of advection, whereas higher $D/K$ and lower n will promote diffusion. While varying m and n, we keep their ratio constant at 0.45, a widely applied average channel concavity (e.g. Wobus et al., 2006; Harel et al., 2016) and in line with the global average of 0.42 reported by Gailleton et al. (2021) (despite considerable variability).
Based on previous modelling (e.g. Roering et al., 2007), one might expect advectiondominant landscapes to be rougher in outline relative to the smoothing effects of diffusion. However, Theodoratos et al. (2018) show that multiple sets of parameters can give rise to equifinality. In our modelling framework, the $D/K$ ratio covaries with n to determine the ratio of hillslope denudation and total (fluvial + hillslope) denudation, denoted here as ${E}_{\mathrm{predicted},\mathrm{diffusive}}/{E}_{\mathrm{total}}$, where E_{total} is the sum of E_{predicted,diffusive} and E_{predicted,advective}. The ratio ${E}_{\mathrm{predicted},\mathrm{diffusive}}/{E}_{\mathrm{total}}$ is therefore a function of both n and $D/K$. In principle, this metric is inversely proportional to an effective Péclet number for net denudation (Perron et al., 2008). For values of ${E}_{\mathrm{predicted},\mathrm{diffusive}}/{E}_{\mathrm{total}}$ close to unity, diffusive processes will dominate, while values closer to zero represent advection dominance (Fig. 1).
2.2 ^{10}Bederived apparent denudation rates
We conduct modelling experiments that employ randomly selected sets of LEM parameter values and then compare our model outputs with the global catalogue of ^{10}Bederived catchmentaveraged denudation rates, OCTOPUS v.2 (Codilean et al., 2022). In addition to apparent denudation rates (N=4631), OCTOPUS includes topographic data and catchment outlines. Using the catchment boundaries in OCTOPUS, we clipped rasters from the Hydrosheds Shuttle Radar Topography Mission (SRTM) dataset, a global digital elevation model (DEM) with 3 arcsec resolution (Lehner et al., 2008), plus the National Elevation Dataset (Gesch et al., 2002) for catchments in Alaska north of 60^{∘}. Local pits (i.e. lakes) within the catchments were filled using the priorityflood method of Barnes et al. (2014). In building the network of local upstream drainage areas for each cell, runoff is assumed to flow down the steepest descent in accordance with the D8 flowrouting algorithm. Slopes for every cell are computed along this steepestdescent flow path.
To determine the effect of DEM resolution, we test our models on 1 arcsec Copernicus DEM data, although for reasons of computational capacity we restrict our analyses to catchments with $<\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{6}}$ grid cells (N=3414) and to the diffusiononly and advectiononly scenarios. We find that DEM resolution is not decisively important to our results (Figs. S1 and S2 in the Supplement).
About 24 % of the OCTOPUS dataset is not exploitable for our purposes: 68 catchments are too small to be processed by the LEM (<3 DEM cells in any dimension), and we exclude 33 of the very largest catchments due to the extreme computational cost. Multiple denudation rate measurements included in OCTOPUS refer to samples from different locations within the same larger catchment. For such cases, a separate catchment is defined only where the drainage area differs by >5 %; otherwise, we amalgamate the data and derive a single average denudation rate. In total, 3640 catchmentwide apparent denudation rates (E_{apparent}) are used in our modelling experiments, ranging from 0.028 to 430 000 km^{2} (11 to 53 ×10^{6} grid cells). We do not separate ^{10}Be measurements conducted on different grainsize fractions.
Processing DEMs begins with computing the drainage network. To minimise the effects of disequilibrium perturbations (e.g. Schwanghart and Scherler, 2017), we smooth river profiles using a ∼1 km averaging window. Unsmoothed profiles are also used for comparison and yield similar results, where the largest discrepancies lie in the advection–diffusion model. After catchment slopes and drainage areas are computed, the diffusion and stream power equations can be solved. The LEM is run for exactly one time step on DEMs representing each of the 3640 catchments; the sediment flux to the catchment outlet is then averaged over the total drainage area to yield a LEMpredicted denudation rate (E_{predicted}). Because different values of input LEM parameters typically yield different denudation rates (with the exception of highly diffusive models, as described below), such models are then optimised by comparison with our catalogue of E_{apparent} data.
2.3 Monte Carlo simulations
We use a bruteforce Monte Carlo approach to investigate the parameter space by running randomly selected sets of parameters and testing the fit of modelled versus observed (^{10}Bederived) denudation rates. We adopt the philosophy of equifinality (e.g. Beven and Binley, 1992) to evaluate the model parameters applied in our LEMs; implicit in these assumptions is that multiple sets of parameters may give rise to a similar or equifinal result (e.g. Csilléry et al., 2010). Hence, we report both the range of optimal parameters and the bestfit model parameters.
In our framework, no more than three parameters are modified and compared at any one time (Table 1). This is possible thanks to several simplifying steps (detailed below) that require fewer modelling runs (e.g. Theodoratos et al., 2018). The performance of the model with a given set of parameters is evaluated based on the mismatch between ${E}_{\mathrm{predicted}}^{*}$ and E_{apparent} with respect to the likelihood function (Beven and Binley, 2014). Modelled and observed rates are compared directly, and no regression is involved. In so doing, one or more local maxima representing an optimised parameter set may be identified in the space defined by parameter values versus the likelihood function. A range is then defined within 1 % of the peak (for example, if the bestfit model has a score of 0.500, we report the range of parameters from models with scores > 0.495).
For each randomly selected set of parameter values (Table 1), the LEM computes a single time step, and the erosion in each grid cell is integrated. ${E}_{\mathrm{predicted}}^{*}$ is then scaled by employing a log transformation on all modelled catchments:
where N is the number of observations, and K^{*} is K prior to log transformation (see Eq. 7 below). The performance of $\mathrm{log}\left({E}_{\mathrm{predicted}}^{*}\right)$ against log (E_{apparent}) is then calculated for each parameter set. This setup offers the advantage of limiting the number of parameters varied; it is not our aim to determine the absolute values of coefficients because these covary; instead, we focus on the ratio $D/K$ (see Sect. 2.1).
After optimal values of A_{c}, $D/K$, and/or n are found, the associated value of K can be corrected for log transformation using an unbiased estimator (after Ferguson, 1986).
where s^{2} is an estimate of the variance:
An equivalent form of Eqs. (4)–(7) is used for the diffusiononly model, simply by replacing K with D. Although we believe this log transformation is justified, we also provide coefficients without Eqs. (6) and (7) (Fig. S7c).
Given that our denudation rate data span several orders of magnitude, we compare log (E_{apparent}) and $\mathrm{log}\left({E}_{\mathrm{predicted}}^{*}\right)$ using the Nash–Sutcliffe coefficient of efficiency (NSE):
Optimised values are defined as the maximum NSE value only where they are surrounded by local maxima in the likelihood, which is the case in all the experiments below. We also tested a version of NSE that incorporates measurement uncertainty, based on Harmel and Smith (2007), but it had little effect on the optimised parameters (Figs. S3 and S4). Additionally, we calculate the mean absolute error (MAE) to gauge the sensitivity to the likelihood function (Figs. S5 and S6).
2.4 The influence of lithology and precipitation on denudation
We subdivide the OCTOPUS catchments according to (1) areally dominant lithology based on the GLiM global geologic map (Hartmann and Moosdorf, 2012), which gives a vectorised description of lithology compiled from a number of regional highresolution geologic maps at a target resolution of 1:1 000 000, and (2) spatially averaged mean annual precipitation (MAP) using the CHELSA dataset (Karger et al., 2017).
Precipitation differences between lithologic subgroups can be significant; for instance, averages range from 851 to 2077 mm yr^{−1} for unconsolidated sediment and metamorphic rocks, respectively. To address lithological variations in the presence of climatic differences between lithologic subgroups in the advectiononly model, we attempt to isolate substrate effects with the form of the stream power equation given by Kooi and Beaumont (1996), which explicitly includes precipitation variations that are normally folded into K under the assumption that precipitation (P) scales linearly with discharge:
Although many factors influence K besides precipitation, we use K_{lith} in this case to denote the variable we are attempting to isolate. We do not attempt to correct for variable precipitation in calculating D, for instance, by devising an equivalent D_{lith} from Eq. (9), and it is only applied when calculating K.
3.1 Advectiononly model
We apply a streampowerbased advectiononly model (Eq. 1) (excluding hillslope diffusion), with two free parameters: a slope exponent (n) and critical drainage area (A_{c}). Variations in m are fixed to n such that concavity is held constant at $m/n=\mathrm{0.45}$. We report the optimised values in terms of maximum value and an optimised range of values (Q_{0.01}) that are within 1 % of the maximum. The advectiononly model (Fig. 2) is globally optimised at n∼1.28 (Q_{0.01}=1.23–1.43; Fig. 2a) and at A_{c}∼0.05 km^{2} (Q_{0.01}=0.03–0.07 km^{2}; Fig. 2b). We note that n changes by ∼6 % (n∼1.36) using the higherresolution 1 arcsec DEMs (Fig. S1). The slight differences in n are likely the result of the sensitivity to higher catchmentaveraged slopes, which naturally arises from higherresolution topographic data (Table S1).
3.2 Diffusiononly model
The diffusiononly model (Fig. 3) is globally optimised with the hillslope diffusion exponent, p∼2.0 (NSE = 0.51; Fig. 3b), and with negligible dependence on DEM resolution, p∼2.0 for the 1 arcsec models (Fig. S2). We also find that the 1 arcsec Copernicus DEM has a lower overall performance than the 3 arcsec SRTM DEM (see discussion in Fig. S2).
3.3 Advection–diffusion model
While the optimisation of our diffusiononly model with p=2 (Fig. 3) is an intriguing result that invites further investigation (cf. Gabet et al., 2021), we retain linear diffusion (p=1) in our advection–diffusion experiment because for p≠1 the model is (1) numerically unstable when implemented in LEMs, and (2) it fails to accommodate hillslope deposition and hence does not conserve mass.
The advection–diffusion model (Fig. 4) is globally optimised at n∼2.3 (Q_{0.01}=1.59–2.70; Fig. 4c) and $D/K\sim \mathrm{1.79}\times {\mathrm{10}}^{\mathrm{6}}$ m^{0.9n+1} (${Q}_{\mathrm{0.01}}=\mathrm{2.97}\times {\mathrm{10}}^{\mathrm{4}}$ to 3.44×10^{7}; Fig. 4d). For n and $D/K$, the Q_{0.01} ranges are quite broad in part because both parameters are codependent (Fig. 4a). Optimum A_{c}∼0.03 km^{2} (Q_{0.01}=0.01–0.17; Fig. 4e) is similar to that from the advectiononly models. The models are more diffusive when n is low, and $D/K$ is high, and E_{predicted} is dependent mainly on catchment slope. This results in values clustering at NSE ∼ 0.34 (Fig. 4c–e) for models with high diffusion. Because $D/K$ covaries with n (and m; Fig. 4a), we find that sediment transport derived from diffusional processes is maximised when ${E}_{\mathrm{predicted},\mathrm{diffusive}}/{E}_{\mathrm{total}}$ is ∼0.43 (Fig. 4f).
In their benchmark study, Portenga and Bierman (2011) employ stepwise regression to relate their compilation of ^{10}Bederived denudation rates to a range of factors embracing topography, climate, lithology, and seismicity. That study, along with the later inclusion of normalised steepness (e.g. Harel et al., 2016; Marder and Gallen, 2023), added substantially to our knowledge of how and why denudation rate varies. Our alternative approach here focuses on the erosional processes at play in terms of advective and diffusive mass flux rather than attempting to interpret the machinations of landscape response to internal and external agents. A significant advantage is that explicit relations between m, n, topography, denudation, and LEM parameters are derived at the scale of the DEM grid cell within each catchment, and success is gauged from the absolute difference between modelled (${E}_{\mathrm{predicted}}^{*}$) and ^{10}Bederived (E_{apparent}) denudation rates. In other words, we evaluate LEM parameters as they are commonly implemented in the models.
4.1 Optimised parameters for landscape evolution models
We first consider some comparisons with previous work regarding advectiononly approaches. Our optimised A_{c}∼0.05 km^{2} (Fig. 2b) for the 3 arcsec resolution models falls near the minimum of the range applied in previous studies, such as Whipple and Tucker (1999), who suggest 0.059–0.14 km^{2}. Our optimised n∼1.3 (n∼1.4 for the 1 arcsec models) is much lower than the 2.6 reported by Harel et al. (2016), which is derived from regression of denudation rate and normalised steepness (k_{s,ref}). Harel et al. (2016) then use the product of k_{s,ref} and a scaling drainage area to calculate M_{χ}. In principle, M_{χ} and n should be similar to our K and n values; however, n is derived from a regression of E_{apparent} against M_{χ}, and it is integrated based on each pixel within the catchment. The large discrepancy between our globally optimised values of n and those of Harel et al. (2016) may stem from the inability of the latter method to accommodate inherent nonlinearities at the subcatchment or subreach scale when n≠1, whereas our approach is designed to capture some of these nonlinear effects. This is particularly important in transient catchments, as spatial heterogeneity in denudation rate is often controlled by steep areas in the catchment, such as knickpoints, and higher values of n amplify the proportion of denudation derived from steep areas relative to the rest of the catchment (Fig. 1).
While linear diffusion (p=1) is commonly applied in landscape evolution studies (e.g. Forte et al., 2016b), our optimised p∼2 for the diffusiononly model is consistent with Gabet al. (2021), in which denudation rate correlates best with the square of hillslope convexity. In response to Gabet et al. (2021), Struble and Roering (2021) point to a systematic underestimation of curvature in natural landscapes that may be an artefact of the numerical methods used for estimating curvature from DEMs. Gabet et al. (2021) employ highresolution (∼1 m) lidar data, but the broader point made by Struble and Roering (2021) poses a serious limitation for largescale LEM analyses that are typically restricted to lowerresolution DEMs. In such cases, the need for mass conservation and numerical stability is an important consideration. And yet, a diffusion equation with the exponent p≠1 is numerically unstable and physically unexplained and does not accommodate deposition (the result of negative curvature). What does it say about the utility of running LEMs on natural landscapes if the optimised parameter value (p∼2) cannot be implemented? Struble and Roering (2021) suggest that p∼2 enhances the influence of steep, rapidly eroding areas on average curvature, which is commonly underestimated by many methods. Below, we discuss how the influence of these steep areas may be approximated by the stream power equation coupled with linear diffusion.
Our advection–diffusion model allows us to explore aspects of how hillslope and river processes govern sediment flux in river catchments. Theoretically, for a catchment in a perfect state of massflux equilibrium (or steady state), hillslopes and rivers erode at the same rate, so either one should be equally useful as a proxy for denudation rate. There would be no apparent advantage to combining advection and diffusion in the same model since they would both yield the same average denudation rate. Landscapes are, however, more often not at a steady state (at least over timescales integrated by cosmogenic ^{10}Be), and the slight dominance of advective denudation in our optimised model (${E}_{\mathrm{predicted},\mathrm{diffusive}}/{E}_{\mathrm{total}}=\mathrm{0.43}$; Fig. 4f) suggests that transient signals disproportionately affect catchmentaveraged denudation rates.
The positive relationship we observe between optimised n and relative diffusivity gives rise to a compelling possibility: as diffusivity increases, advective denudation becomes less important as a proxy for the total average denudation rate within the catchment and more of a proxy for transience focused on the most rapidly eroding zones. However, in the absence of diffusion, river incision must account for all sediment that would otherwise be eroded diffusively from hillslopes. The increase in optimised n with diffusivity therefore represents the expanding role of steep transient zones in dictating the catchmentscale denudation rate. Our optimised result for the advection–diffusion model, n∼2.3 (Fig. 4c), is compatible with previous work suggesting that, typically, n>1 (e.g. Lague, 2014; Harel et al., 2016).
The best correlation between predicted and apparent denudation rates occurs when $D/K\sim \mathrm{1.79}\times {\mathrm{10}}^{\mathrm{6}}$ (Fig. 4d). This outcome broadly agrees with other studies that use K values in the range of $\sim {\mathrm{10}}^{\mathrm{8}}$ to 10^{−5} m^{(n−1)} yr^{−1} and D values in the expected range noted by Fernandes and Dietrich (1997) of $\mathrm{4.4}\times {\mathrm{10}}^{\mathrm{4}}$ to $\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{2}}$ m^{2} yr^{−1}. Whipple et al. (2017) report an optimal $D/K$ ratio of 5×10^{2} from Himalayan catchments, although fixing n=1 in their models is a limiting assumption because $D/K$ covaries with n, as we show. The diffusion model employed here assumes that the longterm flux of hillslope material is similar to the amount transported in one time step. In reality, catchments may not be at a steady state, and the hillslope denudation rate may change notably so as to change the rate of hillslope flux within individual catchments.
4.2 Erosion and precipitation
Correlating topographically derived metrics with mean annual precipitation (MAP) on a global scale has been a longstanding goal (e.g. Ahnert, 1970). Harel et al. (2016) examine correlations between stream power variables and climate as defined by the Köppen–Geiger scheme. They find that aridity yields the lowest n and that, in general, k_{s,ref} and temperature covary inversely: warm deserts yield the highest k_{s,ref} and polar regions the lowest, on average, albeit with large uncertainties. Because our approach compares model results to denudation rates, rather than using regression, we can more directly correlate K and n as they are implemented in LEMs under differing climate. This also means that our method is potentially better suited for transient catchments.
We demonstrate a general increase in n and p with precipitation for both the advection model and the diffusion model (Fig. 5b and c), with the exception of the highest MAP bin (∼2300–6500 mm yr^{−1}). Our advection–diffusion model shows a similar increase in n, aside from the lowest MAP bin (Fig. 5d). In the streampowerbased models, the rising n with respect to precipitation suggests that wetter environments favour a nonlinear erosional response perhaps tied to heavytailed floodfrequency distributions. Extremely variable flow regimes are characteristic of drylands (Zaman et al., 2012), and yet, in two of the three models (advectiononly and diffusiononly), the exponent (n or p) is low. By contrast, only the advection–diffusion model yields a relatively high n (∼2.8) for drylands. This may indicate that the advection–diffusion model, which we suggest is more sensitive to transience, does a better job at simulating dryland catchments. However, we emphasise the relatively low performance of all models (NSE = 0.2–0.3) for the driest settings (Fig. 5b–d), which may reflect the challenge of capturing the erosional dynamics of drylands with such simple models.
It is difficult to differentiate changes in exponents from changes in model coefficients due to their inherent coupling in powerlaw functions (e.g. Syvitski et al., 2000). In response to this issue, we calculate the coefficient D (for diffusiononly) or K (advection) for each of the optimised models (i.e. with constant exponents) using Eq. (5) in catchments represented by 20 MAP bins distributed regularly ($N=\mathrm{182}\pm \mathrm{1}$; Fig. 5e). By looking at how the coefficients vary within each bin, we effectively reduce the influence of topography and isolate the relationship between erosion rate and precipitation. We emphasise that, due to the wide range of data used, optimised models can have large uncertainties (dry regions having the highest uncertainty relative to the mean), with residuals that are generally lognormally distributed (Fig. S8). We also trialled different drainagearea bin sizes, where a similar (albeit noisier) pattern remains with higher bin density (Fig. S7).
Optimised coefficients show a local peak in denudation rates centred around 300 mm yr^{−1}, subsequently dipping overall from around 1100 to 1600 mm yr^{−1} before increasing again for extremely wet regions (Fig. 5e). These results agree well with the classic work of Langbein and Schumm (1958), which suggests that the fastesteroding environments are semiarid (MAP ∼ 250 mm yr^{−1}). This relationship is thought to be a product of the interplay between denudation, vegetation, total precipitation, and storm frequency in semiarid regions – an outcome reproduced by Istanbulluoglu and Bras (2006), who show a positive relationship between sediment transport and the effects of reduced sediment cover and increased runoff during drought.
While Langbein and Schumm (1958) had scant access to data from wetter settings, our results reveal an upward trend in coefficient values for MAP > 1500 mm yr^{−1} (Fig. 5e). This is in line with the global study of sediment yield and climate of Walling and Kleo (1979), which also shows a further major peak at ∼800 mm yr^{−1} and may be attributed to the most expansive agricultural production globally (e.g. Hyman et al., 2016). In contrast to Walling and Kleo (1979), who do not isolate the effects of variable land use, topography, and geology on sediment yields, our use of denudation rates based on ^{10}Be means that we can largely ignore the effects of land use. This may explain the subdued peak at ∼ 800 mm yr^{−1} in our data.
Marder and Gallen (2023) find a nonlinear relationship between ^{10}Bederived denudation rate and k_{sn} via regression; they also find that the exponent relating k_{sn} to erosion rate (here related to n) increases with precipitation. They took steps to exclude transient catchments, as transience may influence n due to the nonlinear effects of stream power, as noted above. Likewise, we find that n generally increases with precipitation, giving credence to the idea that our approach may account to some degree for transience. In particular, we envisage cases in which steep zones, such as adjusting areas downstream of knickpoints, are responsible for a large proportion of the total catchment erosional flux, even when they may represent a small fraction of the drainage area (Willenbring et al., 2013). Where denudation is a nonlinear function of slope and drainage area, an average of these denudation rates is unlikely to be proportional to the integrated k_{sn}.
One acknowledged shortcoming of our approach (noted in Sect. 2.4) is that some lithologies may be overrepresented in areas with higher or lower MAP. As we discuss in the next section, it is a challenge to discriminate biases due to lithology from those linked to precipitation in bins that span heterogeneous lithologies.
4.3 Erosion and lithology
We assessed the variability in LEM parameters within each lithological bin following a similar approach to Sect. 4.2, with the additional step of employing K_{lith} (Eq. 9) to isolate the effects of lithology on K in the advectiononly model. The threefold range in erodibility (Fig. 6e) is much lower than that reported elsewhere, in some cases by several orders of magnitude (Sklar and Dietrich, 2001; GarciaCastellanos and O'Connor, 2018). This may be due to the greater focus on the differential erodibility within individual sites and to the spatial scale of the analysis. Despite our efforts to account for some of the covariation with MAP, our analysis inevitably smooths out some variability owing to the diversity of catchments incorporated within each lithological bin.
We find that for the advectiononly model (Fig. 6b), n is higher in metamorphic and intermediate plutonic rocks, which also tend to be more resistant (e.g. Moosdorf et al., 2018). This is not surprising given that n is thought to be influenced by higher thresholds for rock detachment, which in turn can lead to more nonlinear behaviour. More complex relationships involving hillslope diffusion may result in the general lack of a comprehensible pattern emerging in the advection–diffusion models (Fig. 6c and d). The higher number of parameters used in the advection–diffusion model leads to broader distributions (Fig. 6d).
When looking at the variability in K, our overall results are somewhat expected: all three models agree on the general tendency of sedimentary rocks being most erodible and plutonic and volcanic being least, although the relative magnitudes of these differences vary between models. Unconsolidated sedimentary rock is the first or secondmost erodible of all according to the diffusion and advectiononly models (Fig. 6e). Models disagree about the least erodible subcategories: basic plutonic, intermediate volcanic, or pyroclastic. The pyroclastic, intrusive plutonic, and unconsolidatedsediment rock types show the most variance. For pyroclastic rock types, the advection models suggest relatively low erodibility, whereas the diffusiononly models suggest moderate erodibility. Most models agree on the moderate erodibility of plutonic–intrusive rocks, but due to the low sample size, the uncertainty is very high. Less expected is the relatively high erodibility of carbonate sedimentary rocks shown for all our models. For example, other studies argue that carbonates should be less erodible (Ott, 2020) and that volcanic rocks may be relatively more erodible than we show here (Moosdorf et al., 2018). Some of these variable findings possibly arise due to factors such as fracture density and weathering condition, which are difficult to accommodate in largescale analyses (Neely et al., 2019), and these factors also contribute to the large uncertainty (Fig. S9).
4.4 Erosion and drainage area
Are LEM parameters sensitive to catchment drainage area? Violin plots (Fig. 7) show exponent (n and p) values binned according to drainage area. The advection–diffusion model shows no dependence on area, but for the advectiononly and diffusiononly models, the exponent values generally increase with drainage area, n ranging from ∼1.1 to ∼1.5, and p ranging from ∼1.8 to 2.2. Additionally, we note the extreme range in n (∼1.3 to ∼3) for small catchments in the advection–diffusion models, potentially indicating the influence of landslides in headwaters (e.g. Yanites et al., 2010). This dependence of n and p on drainage area may be difficult to discern from other effects, but one possibility may be the inherent link with precipitation. Larger catchments are more likely to cross large rainfall gradients, which as we show, gives rise to higher n or p.
4.5 Limitations and future considerations
An important limitation of our approach is that it fails to unravel the $m/n$ ratio, which is known to vary widely (Gailleton et al., 2021). We ran several trials (using optimised A_{c}=0.05 km^{2} and 20 values of n linearly spaced from 0 to 4) and found that $m/n=\mathrm{0.3}$, 0.45, and 0.6 yielded negligible differences in fit (NSE = 0.48, 0.47, and 0.45, respectively) despite the similar values of optimised n of 1.4, 1.2, and 1.2, respectively (Fig. S10). One option for improvement is to determine $m/n$ from the river profile concavity using recent integralbased techniques (e.g. Harel et al., 2016; Gailleton et al., 2021) before running the suite of models. This may be especially useful for tuning the $m/n$ ratio within specific regions.
Our approach cannot account for several factors that are known to bias ^{10}Bederived denudation rates (e.g. Von Blanckenburg, 2005; Dingle et al., 2018; Hippe et al., 2019; Struck et al., 2018a). Such factors violate two key underpinning premises of the method: steady and uniform erosion across the catchment and continuous exposure of sediment at/near the surface (noted in the section “Catchmentaveraged denudation rates from cosmogenic ^{10}Be”). For example, sudden pulses of sediment from sources of deepseated mass wasting can shift apparent denudation rates downward by diluting the nuclide abundances in the river sediment sample. Given that our LEM is deterministic, it is difficult to model the impact of landslides on nuclide inventories, particularly in small headwater catchments where their influence will be greater (e.g. Yanites et al., 2010). A second pertinent issue is that our modelling assumes the erosional flux is transported instantly to the catchment outlet without intermediate storage. And yet, river sediment is likely to experience multiple episodes of erosion and deposition, especially in large lowland catchments where the volume of sediment storage expands greatly. Intermediate sediment storage can push ^{10}Bederived denudation rates either upwards or downwards depending on the duration and depth of sediment burial (Wittmann and von Blanckenburg, 2016; Struck et al., 2018b). The extent of bias in the OCTOPUS dataset (Codilean et al., 2022) is essentially unknown, but it no doubt contributes to the data scatter we report here.
In addition to the complications with streampower and diffusionbased models (noted in the Introduction), the inherent assumption that channel width increases monotonically with stream discharge (or its proxy, drainage area) is one that is violated widely. Channel width narrows at knickpoints (Whitbread et al., 2015; Yanites, 2018) and across diverse substrate erodibilities (Jansen et al., 2010; Croissant et al., 2017). The value of n extracted from transient landscapes will be influenced by this process, and a goal of future global analyses should be to incorporate more sophisticated rules for channel width evolution (e.g. Yanites, 2018).
The approach described here is sufficiently robust to incorporate many different types of models. We would like to see, for instance, an exploration of the way in which drainage is routed through catchments. Different flowrouting schemes can give rise to notably different drainage networks and corresponding drainage areas (e.g. Endreny and Wood, 2003). Here we have used unidirectional “D8” flow, but alternatives such as multipleflow routing can produce alternative results and perspectives (e.g. Pelletier, 2004). Future efforts may test a range of flowrouting methods against denudation rate in order to test their efficacy.
We have examined the most widely used parameters applied to a set of three landscape evolution model setups: (1) a streampowerbased advectiononly model, (2) a diffusiononly model, and (3) an advection–diffusion hybrid model. We optimised the parameter values by directly comparing the catchmentaveraged denudation rates predicted by our three models with a global catalogue of ^{10}Bederived apparent denudation rates (Codilean et al., 2022).
The diffusiononly model outperformed the advectiononly model when applying p∼2. However, the physical implications and numerical limitations of this result make it impractical for implementation in LEMs. Instead, we propose that linear diffusion coupled with fluvial denudation (advection–diffusion) captures a high proportion of sediment derived from rapidly eroding, steep areas in a similar sense to a diffusion model with the exponent p∼2. In the advection–diffusion hybrid model, the best agreement between the predicted and apparent denudation rates is observed with n∼2.3 (assuming a fixed concavity, $m/n=\mathrm{0.45}$), while the ratio of diffusivity $/$ advection coefficient ($D/K$) is optimised at $\sim \mathrm{1.79}\times {\mathrm{10}}^{\mathrm{6}}$.
The Monte Carlo method employed here is a simple and powerful means of identifying ideal parameter sets over large spatial scales and is especially useful for dealing with sparse datasets. We applied the same approach to elucidate differences in optimal LEM parameters when considering lithology and precipitation. By looking at the LEM coefficients, we were able to better account for the influence of topography when isolating the relationship between denudation rate and precipitation/lithology. Of particular interest was a general upward trend in the coefficients (K and D) with respect to precipitation and a local maximum centred at ∼300 mm yr^{−1}. This local maximum may represent the higher erodibility of semiarid environments identified by Langbein and Schumm (1958).
Nevertheless, many other influences on denudation are yet to be explored in a robust way. At the local and regional scale, optimised values will differ from what we have inferred here, but future studies may use these parameter ranges as a baseline to inform largescale landscape evolution studies. Moreover, our methodology could be extended to incorporate more complexity into the canonicaladvectionbased and diffusionbased equations applied here.
Code related to this article is available online at https://doi.org/10.5281/zenodo.8317033 (Ruetenik et al., 2023). Maps were produced using PyGMT (https://doi.org/10.5281/zenodo.8303186, Tian et al., 2023).
The supplement related to this article is available online at: https://doi.org/10.5194/esurf118652023supplement.
GR conceived the study, performed data analysis, and devised the code. JDJ and PV assisted with framing the study, and LYM performed code analysis. All coauthors contributed to manuscript production.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Many thanks to Bruce Wilkinson for his guidance in framing the study and to Daniel GarciaCastellanos and Kim Huppert for their constructive advice. We are grateful for the reviews from Richard Ott and Boris Gailleton, who encouraged us to explore new datasets and methods and greatly improved the manuscript. Several calculations were performed on the CSDMS Blanca HPC at the University of Colorado, Boulder.
This paper was edited by Tom Coulthard and reviewed by Richard Ott and Boris Gailleton.
Adams, B. A., Whipple, K. X., Forte, A. M., Heimsath, A. M., and Hodges, K. V.: Climate controls on erosion in tectonically active landscapes, Sci. Adv., 6, 3166, https://doi.org/10.1126/sciadv.aaz3166, 2020.
Ahnert, F.: Functional relationships between denudation, relief, and uplift in large midlatitude drainage basins, Am. J. Sci., 268, 243–263, 1970.
Barnes, R., Lehman, C., and Mulla, D.: Priorityflood: An optimal depressionfilling and watershedlabeling algorithm for digital elevation models, Comput. Geosci., 62, 117–127, 2014.
Barnhart, K. R., Tucker, G. E., Doty, S. G., Shobe, C. M., Glade, R. C., Rossi, M. W, and Hill, M. C.: Inverting topography for landscape evolution model process representation: 1. Conceptualization and sensitivity analysis, J. Geophys. Res.Earth, 125, 1–31, https://doi.org/10.1029/2018JF004961, 2020.
Beven, K. and Binley, A.: The future of distributed models: model calibration and uncertainty prediction, Hydrol. Process., 6, 279–298, 1992.
Beven, K. and Binley, A.: GLUE: 20 years on, Hydrol. Process., 28, 5897–5918, 2014.
Brown, L. R.: World population growth, soil erosion, and food security, Science, 214, 995–1002, https://doi.org/10.1126/science.7302578, 1981.
Carretier, S., Regard, V., Vassallo, R., Aguilar, G., Martinod, J., Riquelme, R., Christophoul, F., Charrier, R., Gayer, E., Farías, M., and Audin, L.: Differences in 10Be concentrations between river sand, gravel and pebbles along the western side of the central Andes, Quater. Geochronol., 27, 33–51, 2015.
Clubb, F. J., Mudd, S. M., Milodowski, D. T., Hurst, M. D., and Slater, L. J.: Objective extraction of channel heads from highresolution topographic data, Water Resour. Res., 50, 4283–4304, 2014.
Codilean, A. T., Munack, H., Saktura, W. M., Cohen, T. J., Jacobs, Z., Ulm, S., Hesse, P. P., Heyman, J., Peters, K. J., Williams, A. N., Saktura, R. B. K., Rui, X., ChishiroDennelly, K., and Panta, A.: OCTOPUS database (v.2), Earth Syst. Sci. Data, 14, 3695–3713, https://doi.org/10.5194/essd1436952022, 2022.
Croissant, T. and Braun, J.: Constraining the stream power law: a novel approach combining a landscape evolution model and an inversion method, Earth Surf. Dynam., 2, 155–166, https://doi.org/10.5194/esurf21552014, 2014.
Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid postseismic landslide evacuation boosted by dynamic river width, Nat. Geosci., 10, 680–684, 2017.
Csilléry, K., Blum, M. G., Gaggiotti, O. E., and François, O.: Approximate Bayesian computation (ABC) in practice, Trends Ecol. Evol., 25, 410–418, 2010.
Dingle, E. H., Sinclair, H. D., Attal, M., Rodés, Á., and Singh, V.: Temporal variability in detrital ^{10}Be concentrations in a large Himalayan catchment, Earth Surf. Dynam., 6, 611–635, https://doi.org/10.5194/esurf66112018, 2018.
Endreny, T. A. and Wood, E. F.: Maximizing spatial congruence of observed and DEMdelineated overland flow networks, Int. J. Geogr. Inform. Sci., 17, 699–713, 2003.
Ferguson, R. I.: River loads underestimated by rating curves, Water Resour. Res., 22, 74–76, 1986.
Fernandes, N. F. and Dietrich, W. E.: Hillslope evolution by diffusive processes: The timescale for equilibrium adjustments, Water Resour. Res., 33, 1307–1318, 1997.
Fontana, G. D. and Marchi, L.: Slope–area relationships and sediment dynamics in two alpine streams, Hydrol. Process., 17, 73–87, 2003.
Forte, A. M., Whipple, K. X., Bookhagen, B., and Rossi, M. W.: Decoupling of modern shortening rates, climate, and topography in the Caucasus, Earth Planet. Sc. Lett., 449, 282–294, 2016a.
Forte, A. M., Yanites, B. J., and Whipple, K. X.: Complexities of landscape evolution during incision through layered stratigraphy with contrasts in rock strength, Earth Surf. Proc. Land., 41, 1736–1757, 2016b.
Fox, M., Goren, L., May, D. A., and Willett, S. D.: Inversion of fluvial channels for paleorock uplift rates in Taiwan, J. Geophys. Res.Earth, 119, 1853–1875, 2014.
Gabet, E. J., Mudd, S. M., Wood, R. W., Grieve, S. W. D., Binnie, S. A., and Dunai, T. J.: Hilltop curvature increases with the square root of erosion rate, J. Geophys. Res.Earth, 126, 1–16, https://doi.org/10.1029/2020JF005858, 2021.
Gailleton, B., Mudd, S. M., Clubb, F. J., Grieve, S. W., and Hurst, M. D.: Impact of changing concavity indices on channel steepness and divide migration metrics, J. Geophys. Res.Earth, 126, 1–31, https://doi.org/10.1029/2020JF006060, 2021.
Gallen, S. F., Wegmann, K. W., and Bohnenstiehl, D. R.: Miocene rejuvenation of topographic relief in the southern Appalachians, GSA Today, 23, 4–10, 2013.
GarciaCastellanos, D. and O'Connor, J. E.: Outburst floods provide erodability estimates consistent with longterm landscape evolution, Sci. Rep., 8, 1–9, 2018.
GarciaCastellanos, D., Vergés, J., Gaspar Escribano, J., and Cloetingh, S.: Interplay between tectonics, climate, and fluvial transport during the Cenozoic evolution of the Ebro Basin (NE Iberia), J. Geophys. Res.Solid, 108, 2347, https://doi.org/10.1029/2002JB002073, 2003.
Gesch, D., Oimoen, M., Greenlee, S., Nelson, C., Steuck, M., and Tyler, D.: The national elevation dataset, Photogram. Eng. Remote Sens., 68, 5–32, 2002.
Granger, D. E., Kirchner, J. W., and Finkel, R.: Spatially averaged longterm erosion rates measured from insitu produced cosmogenic nuclides in alluvial sediment, J. Geol., 104, 249–257, 1996.
Harel, M. A., Mudd, S. M., and Attal, M.: Global analysis of the stream power law parameters based on worldwide 10Be denudation rates, Geomorphology, 268, 184–196, 2016.
Harmel, R. D. and Smith, P. K.: Consideration of measurement uncertainty in the evaluation of goodnessoffit in hydrologic and water quality modeling, J. Hydrol., 337, 326–336, 2007.
Hartmann, J. and Moosdorf, N.: The new global lithological map database GLiM: A representation of rock properties at the Earth surface, Geochem. Geophy. Geosy., 13, 1–37, https://doi.org/10.1029/2012GC004370, 2012.
Hippe, K., Gordijn, T., Picotti, V., Hajdas, I., Jansen, J. D., Christl, M., Vockenhuber, C., Maden, C., Akçar, N., and IvyOchs, S.: Fluvial dynamics and ^{14}C^{10}Be disequilibrium on the Bolivian Altiplano, Earth Surf. Proc. Land., 44, 766–780, 2019.
Hooke, R. L.: On the history of humans as geomorphic agents, Geology, 28, 843–846, 2000.
Howard, A. D., Dietrich, W. E., and Seidl, M. A.: Modeling fluvial erosion on regional to continental scales, J. Geophys. Res.Solid, 99, 13971–13986, 1994.
Hyman, G., Barona, E., Biradar, C., Guevara, E., Dixon, J., Beebe, S., Castano, S. E., Alabi, T., Gumma, M. K., Sivasankar, S., and Rivera, O.: Priority regions for research on dryland cereals and legumes, F1000 Research, 5 pp., https://doi.org/10.12688/f1000research.8657.2, 2016.
Jansen, J. D., Codilean, A. T., Bishop, P., and Hoey, T. B.: Scaledependence of lithological control on topography; bedrock channel geometry and catchment morphometry in western Scotland, J. Geol., 118, 223–246, 2010.
Istanbulluoglu, E. and Bras, R. L.: On the dynamics of soil moisture, vegetation, and erosion: Implications of climate variability and change, Water Resour. Res., 42, W06418, https://doi.org/10.1029/2005WR004113, 2006.
Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., SoriaAuza, R. W., Zimmermann, N. E., Linder, H. P., and Kessler, M.: Climatologies at high resolution for the earth's land surface areas, Sci. Data, 4, 1–20, 2017.
Kooi, H. and Beaumont, C.: Largescale geomorphology: Classical concepts reconciled and integrated with contemporary ideas via a surface processes model, J. Geophys. Res.Solid, 101, 3361–3386, 1996.
Lague, D.:. The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61, 2014.
Lague, D. and Davy, P.: Constraints on the longterm colluvial erosion law by analyzing slopearea relationships at various tectonic uplift rates in the Siwaliks Hills (Nepal), J. Geophys. Res.Solid, 108, 2129, https://doi.org/10.1029/2002JB001893, 2003.
Lal, D.: Cosmic ray labeling of erosion surfaces: in situ nuclide production rates and erosion models, Earth Planet. Sc. Lett., 104, 424–439, 1991.
Langbein, W. B. and Schumm, S. A.: Yield of sediment in relation to mean annual precipitation, Eos Trans. Am. Geophys. Union, 39, 1076–1084, 1958.
Lehner, B., Verdin, K., and Jarvis, A.: New global hydrography derived from spaceborne elevation data, Eos Trans. Am. Geophys. Union, 89, 93–94, 2008.
Lukens, C. E., Riebe, C. S., Sklar, L. S., and Shuster, D. L.: Grain size bias in cosmogenic nuclide studies of stream sediment in steep terrain, J. Geophys. Res.Earth, 121, 978–999, 2016.
Marder, E. and Gallen, S. F.: Climate control on the relationship between erosion rate and fluvial topography, Geology, 51, 424–427, 2023.
Miller, S. R., Sak, P. B., Kirby, E., and Bierman, P. R.: Neogene rejuvenation of central Appalachian topography: Evidence for differential rock uplift from stream profiles and erosion rates, Earth Planet. Sc. Lett., 369, 1–12, 2013.
Moosdorf, N., Cohen, S., and von Hagke, C.,: A global erodibility index to represent sediment production potential of different rock types, Appl. Geogr., 101, 36–44, 2018.
Mudd, S. M.: Detection of transience in eroding landscapes, Earth Surf. Proc. Land., 42, 24–41, 2016.
Mudd, S. M., Clubb, F. J., Gailleton, B., and Hurst, M. D.: How concave are river channels?, Earth Surf. Dynam., 6, 505–523, https://doi.org/10.5194/esurf65052018, 2018.
Neely, A. B., DiBiase, R. A., Corbett, L. B., Bierman, P. R., and Caffee, M. W.: Bedrock fracture density controls on hillslope erodibility in steep, rocky landscapes with patchy soil cover, southern California, USA, Earth Planet. Sc. Lett., 522, 186–197, 2019.
Ott, R. F.: How lithology impacts global topography, vegetation, and animal biodiversity: A globalscale analysis of mountainous regions, Geophys. Res. Lett., 47, 1–11, https://doi.org/10.1029/2020GL088649, 2020.
Pedersen, V. K., Braun, J., and Huismans, R. S.: Eocene to midPliocene landscape evolution in Scandinavia inferred from offshore sediment volumes and preglacial topography using inverse modelling, Geomorphology, 303, 467–485, 2018.
Pelletier, J. D.: Persistent drainage migration in a numerical landscape evolution model, Geophys. Res. Lett., 31, L20501, https://doi.org/10.1029/2004GL020802, 2004.
Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of firstorder valleys, J. Geophys. Res.Earth, 113, F04016, https://doi.org/10.1029/2007JF000977, 2008.
Pico, T., Mitrovica, J. X., Perron, J. T., Ferrier, K. L., and Braun, J.: Influence of glacial isostatic adjustment on river evolution along the US midAtlantic coast, Earth Planet. Sc. Lett., 522, 176–185, 2019.
Portenga, E. W. and Bierman, P. R.: Understanding Earth's eroding surface with ^{10}Be, GSA Today, 21, 4–10, 2011.
Riebe, C. S., Sklar, L. S., Lukens, C. E., and Shuster, D. L.: Climate and topography control the size and flux of sediment produced on steep mountain slopes, P. Natl. Acad. Sci. USA, 112, 15574–15579, 2015.
Roering, J. J., Perron, J. T., and Kirchner, J. W.: Functional relationships between denudation and hillslope form and relief, Earth Planet. Sc. Lett., 264, 245–258, 2007.
Ruetenik, G. A., Jansen, J. D., Val, P., and YläMella, L.: Code and data for Ruetenik et al., (2023): Optimising global landscape evolution models with 10Be (v0.13), Zenodo [data set and code], https://doi.org/10.5281/zenodo.8317033, 2023.
Ruetenik, G., Moucha, R., and de Boer, B.: Deformation in response to landscape evolution during glacial cycles on the US Atlantic passive margin, Earth Planet. Sc. Lett., 526, 115759, https://doi.org/10.1016/j.epsl.2019.115759, 2019.
Schwanghart, W. and Scherler, D.: Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques, Earth Surf. Dynam., 5, 821–839, https://doi.org/10.5194/esurf58212017, 2017.
Simoes, M., Braun, J., and Bonnet, S.: Continentalscale erosion and transport laws: A new approach to quantitatively investigate macroscale landscapes and associated sediment fluxes over the geological past, Geochem. Geophy. Geosy., 11, Q09001, https://doi.org/10.1029/2010GC003121, 2010.
Sklar, L. and Dietrich, W. E.: Sediment supply, grain size and rock strength controls on rates of river incision into bedrock, Geology, 29, 1087–1090, 2001.
Starke, J., Ehlers, T. A., and Schaller, M.: Latitudinal effect of vegetation on erosion rates identified along western South America, Science, 367, 1358–1361, 2020.
Struble, W. T. and Roering, J. J.: Hilltop curvature as a proxy for erosion rate: wavelets enable rapid computation and reveal systematic underestimation, Earth Surf. Dynam., 9, 1279–1300, https://doi.org/10.5194/esurf912792021, 2021.
Struck, M., Jansen, J. D., Fujioka, T., Codilean, A. T., Fink, D., Egholm, D. L., Fülöp, R. H., Wilcken, K. M., Price, D. M., and Kotevski, S.: Soil production and transport on postorogenic desert hillslopes quantified with ^{10}Be and ^{26}Al, Geol. Soc. Am. Bull., 130, 1017–1040, 2018a.
Struck, M., Jansen, J. D., Fujioka, T., Codilean, A. T., Fink, D., Fülöp, R.H., Wilcken, K. M., Price, D. M., Kotevski, S., Fifield, L. K., and Chappell, J.: Tracking the ^{10}Be–^{26}Al sourcearea signal in sedimentrouting systems of arid central Australia, Earth Surf. Dynam., 6, 329–349, https://doi.org/10.5194/esurf63292018, 2018b.
Syvitski, J. P., Morehead, M. D., Bahr, D. B., and Mulder, T.: Estimating fluvial sediment transport: the rating parameters, Water Resour. Res., 36, 2747–2760, 2000.
Temme, A. J. A. M., Baartman, J. E. M., and Schoorl, J. M.: Can uncertain landscape evolution models discriminate between landscape responses to stable and changing future climate? A millennialscale test, Global Planet. Change, 69, 48–58, 2009.
Theodoratos, N., Seybold, H., and Kirchner, J. W.: Scaling and similarity of a streampower incision and linear diffusion landscape evolution model, Earth Surf. Dynam., 6, 779–808, https://doi.org/10.5194/esurf67792018, 2018.
Tian, D., Uieda, L., Leong, W. J., Schlitzer, W., Fröhlich, Y., Grund, M., Jones, M., Toney, L., Yao, J., Magen, Y., JingHui, T., Materna, K., Belem, A., Newton, T., Anant, A., Ziebarth, M., Quinn, J., and Wessel, P.: PyGMT: A Python Interface for the Generic Mapping Tools, Zenodo [data set], https://doi.org/10.5281/zenodo.8303186, 2023.
Von Blanckenburg, F.: The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment, Earth Planet. Sc. Lett., 237, 462–479, 2005.
Walling, D. E. and Kleo, A. H. A.: Sediment yields of rivers in areas of low precipitation: a global view, The Hydrology of areas of low precipitation, International Association of Hydrological Sciences (IAHS) Publication, Vol. 128, 479–493, 1979.
Whipple, K. X. and Tucker, G. E.: Dynamics of the streampower river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs, J. Geophys. Res.Solid, 104, 17661–17674, 1999.
Whipple, K. X., Forte, A. M., DiBiase, R. A., Gasparini, N. M., and Ouimet, W. B.: Timescales of landscape response to divide migration and drainage capture: Implications for the role of divide mobility in landscape evolution, J. Geophys. Res.Earth, 122, 248–273, 2017.
Whitbread, K., Jansen, J. D., Bishop, P., and Attal, M.: Substrate, sediment, and slope controls on bedrock channel geometry in postglacial streams, J. Geophys. Res.Earth, 120, 779–798, 2015.
Willenbring, J. K., Gasparini, N. M., Crosby, B. T., and Brocard, G.: What does a mean mean? The temporal evolution of detrital cosmogenic denudation rates in a transient landscape, Geology, 41, 1215–1218, 2013.
Wittmann, H. and von Blanckenburg, F.: The geological significance of cosmogenic nuclides in large lowland river basins, EarthSci. Rev., 159, 118–141, 2016.
Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., Sheehan, D., and Willett, S. D.: Tectonics from topography: Procedures, promise, and pitfalls, Special papers 398, Geological Society of America, 55 pp., https://doi.org/10.1130/2006.2398(04), 2006.
Yanites, B. J.: The dynamics of channel slope, width, and sediment in actively eroding bedrock river systems, J. Geophys. Res.Earth, 123, 1504–1527, 2018.
Yanites, B. J., Tucker, G. E., Mueller, K. J., and Chen, Y. G.: How rivers react to large earthquakes: Evidence from central Taiwan, Geology, 38, 639–642, 2010.
Zaman, M. A., Rahman, A., and Haddad, K.: Regional flood frequency analysis in arid regions: A case study for Australia, J. Hydrol., 475, 74–83, 2012.
Zavala, V., Carretier, S., and Bonnet, S.: Influence of orographic precipitation on the topographic and erosional evolution of mountain ranges, Basin Res., 32, 1574–1599, 2020.