Articles | Volume 11, issue 5
Research article
18 Sep 2023
Research article |  | 18 Sep 2023

Optimising global landscape evolution models with 10Be

Gregory A. Ruetenik, John D. Jansen, Pedro Val, and Lotta Ylä-Mella

By simulating erosion and deposition, landscape evolution models (LEMs) offer powerful insights into Earth surface processes and dynamics. Stream-power-based LEMs are often constructed from parameters describing drainage area (m), slope (n), substrate erodibility (K), hillslope diffusion (D), and a critical drainage area (Ac) 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 catchment-averaged cosmogenic 10Be-derived denudation rates with the aim to optimise a set of LEMs via a Monte Carlo-based parameter search. We consider three model scenarios: advection-only, diffusion-only, 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 10Be 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 diffusion-only model yields a slightly lower misfit when comparing model outputs with observed erosion rates than the advection-only 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 large-scale studies spanning long timescales and diverse landscape properties. Moreover, our direct comparison of model-predicted versus observed denudation rates is preferable to methods that rely upon catchment-scale 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 long-standing hypothesis that semi-arid environments are among the most erodible.

1 Introduction

To appreciate short-term changes in Earth surface processes, such as those induced by humans (Brown, 1981; Hooke, 2000), it is first necessary to understand long-term rates of denudation and deposition. Recognising this, some recent studies (e.g. Simoes et al., 2010) have derived erosion-transport rules from topography with an aim to predict macro-scale patterns of denudation and sediment flux. At more restricted spatial scales, denudation rates based on cosmogenic nuclides (e.g. 10Be) show a modest exponential correlation with catchment-averaged 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; Garcia-Castellanos 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 palaeo-topography, we settle for optimising LEM parameters with respect to catchment-averaged denudation rates estimated with cosmogenic 10Be. 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 10Be-derived 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 stream-power-plus-diffusion 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 trade-offs 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).

Catchment-averaged denudation rates from cosmogenic 10Be

Rates of catchment-scale denudation can be estimated by measuring the abundances of cosmogenic radionuclides, such as 10Be, in quartz-bearing 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 10Be 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 10Be-derived denudation rates and how we parameterise LEMs. Estimating catchment-averaged 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 long-term, steady bedrock erosion distributed uniformly across the catchment and that sediments have experienced continuous exposure to cosmic rays at/near the surface. In detail, long-term steady erosion refers to at least one attenuation length (∼0.6 m) of surface lowering integrated over a 103–105-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 cosmic-ray exposure. Other sources of potential discord relate to lithology, catchment size, and hypsometry, which are known to affect sediment transport dynamics and grain-size 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 10Be-derived denudation rates (e.g. Portenga and Bierman, 2011; Harel et al., 2016). Catchment-wide 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, 10Be-derived data currently offer the most widely distributed insight into long-term denudation on a global scale.

2 Methods

2.1 Stream power and hillslope diffusion

Stream power is represented by a non-linear 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, Ac) 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

(1) E predicted , advective = K A m S n ,

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 non-linearity 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 fastest-eroding 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 main-stem 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:

(2) E predicted , diffusive = D d 2 z d x 2 + d 2 z d y 2 p ,

where D is diffusivity, which is reported to range from 4.4×10-4 to 3.6×10-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 non-real components. Thus, we calculate an average catchment-wide 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 catchment-averaged 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:

(3) E predicted K = A m S n + D K d 2 z d x 2 + d 2 z d y 2 .

With our advection–diffusion model formulation (Eq. 3), we set out to solve simultaneously for globally optimised values of D/K, n, and Ac. 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 advection-dominant 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 Epredicted,diffusive/Etotal, where Etotal is the sum of Epredicted,diffusive and Epredicted,advective. The ratio Epredicted,diffusive/Etotal 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 Epredicted,diffusive/Etotal close to unity, diffusive processes will dominate, while values closer to zero represent advection dominance (Fig. 1).

Figure 1(a) Catchment example (Swakop River, Namibia) clipped from a Hydrosheds DEM based on the shapefile provided in OCTOPUS v.2 (Codilean et al., 2022). The lower panels (b–e) show corresponding relative denudation rates (colour ramp spans 0 %–98 % of the range) for differing parameter values. No diffusion is included in (b) and (d); hence erosion is focused in the channels. In (c) and (e), a moderately high (107) diffusivity is used relative to advection, which causes erosion to be more focused on hillslopes.


2.210Be-derived 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 10Be-derived catchment-averaged 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 priority-flood 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 flow-routing algorithm. Slopes for every cell are computed along this steepest-descent 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 <1.3×106 grid cells (N=3414) and to the diffusion-only and advection-only 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 catchment-wide apparent denudation rates (Eapparent) are used in our modelling experiments, ranging from 0.028 to 430 000 km2 (11 to 53 ×106 grid cells). We do not separate 10Be measurements conducted on different grain-size fractions.

Table 1Parameter values and ranges for the three model set-ups.

Download Print Version | Download XLSX

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 LEM-predicted denudation rate (Epredicted). 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 Eapparent data.

2.3 Monte Carlo simulations

We use a brute-force Monte Carlo approach to investigate the parameter space by running randomly selected sets of parameters and testing the fit of modelled versus observed (10Be-derived) 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 best-fit 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 Epredicted* and Eapparent 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 best-fit 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. Epredicted* 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 log(Epredicted*) against log (Eapparent) is then calculated for each parameter set. This set-up 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 Ac, 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).

(6) K = K * e s 2 2 ,

where s2 is an estimate of the variance:

(7) s 2 = 1 N - 1 log E apparent - log E predicted * 2 .

An equivalent form of Eqs. (4)–(7) is used for the diffusion-only 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 (Eapparent) and log(Epredicted*) using the Nash–Sutcliffe coefficient of efficiency (NSE):

(8) NSE = 1 - E apparent - E predicted * 2 E apparent - mean E apparent 2 .

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 high-resolution 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 advection-only 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:

(9) E predicted , advective = K lith ( P A ) m S n .

Although many factors influence K besides precipitation, we use Klith 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 Dlith from Eq. (9), and it is only applied when calculating K.

3 Results

3.1 Advection-only model

We apply a stream-power-based advection-only model (Eq. 1) (excluding hillslope diffusion), with two free parameters: a slope exponent (n) and critical drainage area (Ac). Variations in m are fixed to n such that concavity is held constant at m/n=0.45. We report the optimised values in terms of maximum value and an optimised range of values (Q0.01) that are within 1 % of the maximum. The advection-only model (Fig. 2) is globally optimised at n∼1.28 (Q0.01=1.23–1.43; Fig. 2a) and at Ac∼0.05 km2 (Q0.01=0.03–0.07 km2; Fig. 2b). We note that n changes by ∼6 % (n∼1.36) using the higher-resolution 1 arcsec DEMs (Fig. S1). The slight differences in n are likely the result of the sensitivity to higher catchment-averaged slopes, which naturally arises from higher-resolution topographic data (Table S1).

Figure 2The advection-only model. (a) Optimised n=1.28 (Q0.01=1.23–1.43), (b) optimised Ac=0.05 km2 (Q0.01=0.03–0.07 km2), (c) apparent versus predicted erosion rate (NSE = 0.48; no regression is performed, and the black line indicates a 1:1 fit).


3.2 Diffusion-only model

The diffusion-only 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).

Figure 3The diffusion-only model. (a) The sole free parameter (p) is optimised at p=2.00. (b) Apparent versus predicted erosion rate (NSE = 0.51; the black line represents a perfect 1:1 fit).


3.3 Advection–diffusion model

While the optimisation of our diffusion-only 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 (Q0.01=1.59–2.70; Fig. 4c) and D/K1.79×106 m0.9n+1 (Q0.01=2.97×104 to 3.44×107; Fig. 4d). For n and D/K, the Q0.01 ranges are quite broad in part because both parameters are co-dependent (Fig. 4a). Optimum Ac∼0.03 km2 (Q0.01=0.01–0.17; Fig. 4e) is similar to that from the advection-only models. The models are more diffusive when n is low, and D/K is high, and Epredicted 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 Epredicted,diffusive/Etotal is ∼0.43 (Fig. 4f).

Figure 4Model parameters representing variations in the relative dominance of advection versus diffusion. (a) Covariance of D/K with n; when D/K is low (no diffusion), optimal n approaches ∼1.3 (y intercept). (b) The best correspondence between Epredicted* and Eapparent is achieved with NSE = 0.52, where (c) n∼2.3 (Q0.01=1.59–2.7), (d) D/K1.79×106 (Q0.01=2.97×104 to 3.44×107), and (e) Ac∼0.03 km2 (Q0.01=0.01–0.17 km2). Clustering at NSE  0.34 in panels (c)(e) represents parameter sets where diffusion dominates over advection. (f) Sediment transport derived from diffusional processes is maximised when Epredicted,diffusive/Etotal is ∼0.43.


4 Discussion

In their benchmark study, Portenga and Bierman (2011) employ stepwise regression to relate their compilation of 10Be-derived 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 (Epredicted*) and 10Be-derived (Eapparent) 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 advection-only approaches. Our optimised Ac∼0.05 km2 (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 km2. 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 (ks,ref). Harel et al. (2016) then use the product of ks,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 Eapparent 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 non-linearities at the sub-catchment or sub-reach scale when n≠1, whereas our approach is designed to capture some of these non-linear 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 diffusion-only 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 high-resolution (∼1 m) lidar data, but the broader point made by Struble and Roering (2021) poses a serious limitation for large-scale LEM analyses that are typically restricted to lower-resolution 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 mass-flux 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 10Be), and the slight dominance of advective denudation in our optimised model (Epredicted,diffusive/Etotal=0.43; Fig. 4f) suggests that transient signals disproportionately affect catchment-averaged 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 catchment-scale 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/K1.79×106 (Fig. 4d). This outcome broadly agrees with other studies that use K values in the range of 10-8 to 10−5 m(n−1) yr−1 and D values in the expected range noted by Fernandes and Dietrich (1997) of 4.4×10-4 to 3.6×10-2 m2 yr−1. Whipple et al. (2017) report an optimal D/K ratio of 5×102 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 long-term 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 long-standing 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, ks,ref and temperature covary inversely: warm deserts yield the highest ks,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 stream-power-based models, the rising n with respect to precipitation suggests that wetter environments favour a non-linear erosional response perhaps tied to heavy-tailed flood-frequency distributions. Extremely variable flow regimes are characteristic of drylands (Zaman et al., 2012), and yet, in two of the three models (advection-only and diffusion-only), 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.

Figure 5(a) Global distribution of catchments in the OCTOPUS v.2 catalogue of 10Be-derived apparent denudation rates (Codilean et al., 2022) coloured by mean annual precipitation (MAP). (b–d) Violin plots (symmetric kernel density plots) of n or p values within the top 1 % of model runs for the advection (b), diffusion (c), and advection–diffusion (d) models. Blue stars correspond to the NSE value for each bin. (e) Coefficients for diffusion (circle), advection (square), and advection–diffusion (cross) calculated for each globally optimised model per MAP bin. Diffusion model plots represent only the global maximum within each bin due to insufficient model runs to form a distribution. All panels use the same colour ramp, which corresponds to the MAP bin; blue shading represents the average bootstrapped 90 % confidence interval spanned by the three models, and the thick blue line is the mean. The grey line represents the range spanned by the 90 % confidence interval averaged among all three models relative to the mean (see Fig. S8).

It is difficult to differentiate changes in exponents from changes in model coefficients due to their inherent coupling in power-law functions (e.g. Syvitski et al., 2000). In response to this issue, we calculate the coefficient D (for diffusion-only) 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=182±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 log-normally distributed (Fig. S8). We also trialled different drainage-area 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 fastest-eroding environments are semi-arid (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 semi-arid 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 10Be 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 non-linear relationship between 10Be-derived denudation rate and ksn via regression; they also find that the exponent relating ksn 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 non-linear 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 non-linear function of slope and drainage area, an average of these denudation rates is unlikely to be proportional to the integrated ksn.

One acknowledged shortcoming of our approach (noted in Sect. 2.4) is that some lithologies may be over-represented 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 Klith (Eq. 9) to isolate the effects of lithology on K in the advection-only 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; Garcia-Castellanos 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.

Figure 6(a) Global distribution of catchments in the OCTOPUS v.2 catalogue of 10Be-derived apparent denudation rates (Codilean et al., 2022) coloured by dominant lithology. (b–d) Violin plots of n or p values within the top 1 % of model runs for the advection (b), diffusion (c), and advection–diffusion (d) models. Blue stars correspond to the NSE for each bin. (e) Coefficients (normalised by their maximum values) for different best-fit models within 12 lithologic subsets (from Hartmann and Moosdorf, 2012). Coefficients for diffusion (circle), advection (square), and advection–diffusion (cross) calculated for each globally optimised model per lithologic bin; blue shading represents average bootstrapped 90 % confidence interval spanned by the three models, and the thick blue line is the mean. All panels use the same colour ramp, which corresponds to the lithologic bin. The grey line represents the range spanned by the 90 % confidence averaged among all three models relative to the mean (see Fig. S9).


We find that for the advection-only 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 non-linear 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 second-most erodible of all according to the diffusion- and advection-only models (Fig. 6e). Models disagree about the least erodible subcategories: basic plutonic, intermediate volcanic, or pyroclastic. The pyroclastic, intrusive plutonic, and unconsolidated-sediment rock types show the most variance. For pyroclastic rock types, the advection models suggest relatively low erodibility, whereas the diffusion-only 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 large-scale 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 advection-only and diffusion-only 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.

Figure 7Distributions of n or p within each area bin for (a) advection-only, (b) diffusion-only, and (c) advection–diffusion models. Blue stars represent the corresponding NSE for each bin. The top 1 % of model runs are used.


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 Ac=0.05 km2 and 20 values of n linearly spaced from 0 to 4) and found that m/n=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 integral-based 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 10Be-derived 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 “Catchment-averaged denudation rates from cosmogenic 10Be”). For example, sudden pulses of sediment from sources of deep-seated 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 10Be-derived 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 stream-power- and diffusion-based 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 flow-routing 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 multiple-flow routing can produce alternative results and perspectives (e.g. Pelletier, 2004). Future efforts may test a range of flow-routing methods against denudation rate in order to test their efficacy.

5 Conclusions

We have examined the most widely used parameters applied to a set of three landscape evolution model set-ups: (1) a stream-power-based advection-only model, (2) a diffusion-only model, and (3) an advection–diffusion hybrid model. We optimised the parameter values by directly comparing the catchment-averaged denudation rates predicted by our three models with a global catalogue of 10Be-derived apparent denudation rates (Codilean et al., 2022).

The diffusion-only model outperformed the advection-only 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=0.45), while the ratio of diffusivity / advection coefficient (D/K) is optimised at 1.79×106.

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 semi-arid 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 large-scale landscape evolution studies. Moreover, our methodology could be extended to incorporate more complexity into the canonical-advection-based and diffusion-based equations applied here.

Code and data availability

Code related to this article is available online at (Ruetenik et al., 2023). Maps were produced using PyGMT (, Tian et al., 2023).


The supplement related to this article is available online at:

Author contributions

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 co-authors contributed to manuscript production.

Competing interests

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 Garcia-Castellanos 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.

Review statement

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,, 2020. 

Ahnert, F.: Functional relationships between denudation, relief, and uplift in large mid-latitude drainage basins, Am. J. Sci., 268, 243–263, 1970. 

Barnes, R., Lehman, C., and Mulla, D.: Priority-flood: An optimal depression-filling and watershed-labeling 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,, 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,, 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 high-resolution 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., Chishiro-Dennelly, K., and Panta, A.: OCTOPUS database (v.2), Earth Syst. Sci. Data, 14, 3695–3713,, 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,, 2014. 

Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid post-seismic 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 10Be concentrations in a large Himalayan catchment, Earth Surf. Dynam., 6, 611–635,, 2018. 

Endreny, T. A. and Wood, E. F.: Maximizing spatial congruence of observed and DEM-delineated 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,, 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,, 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. 

Garcia-Castellanos, D. and O'Connor, J. E.: Outburst floods provide erodability estimates consistent with long-term landscape evolution, Sci. Rep., 8, 1–9, 2018. 

Garcia-Castellanos, 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,, 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 long-term erosion rates measured from in-situ 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 goodness-of-fit 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,, 2012. 

Hippe, K., Gordijn, T., Picotti, V., Hajdas, I., Jansen, J. D., Christl, M., Vockenhuber, C., Maden, C., Akçar, N., and Ivy-Ochs, S.: Fluvial dynamics and 14C-10Be 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.,, 2016. 

Jansen, J. D., Codilean, A. T., Bishop, P., and Hoey, T. B.: Scale-dependence 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,, 2006. 

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, 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.: Large-scale 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 long-term colluvial erosion law by analyzing slope-area relationships at various tectonic uplift rates in the Siwaliks Hills (Nepal), J. Geophys. Res.-Solid, 108, 2129,, 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,, 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 global-scale analysis of mountainous regions, Geophys. Res. Lett., 47, 1–11,, 2020. 

Pedersen, V. K., Braun, J., and Huismans, R. S.: Eocene to mid-Pliocene landscape evolution in Scandinavia inferred from offshore sediment volumes and pre-glacial 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,, 2004. 

Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of first-order valleys, J. Geophys. Res.-Earth, 113, F04016,, 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 mid-Atlantic coast, Earth Planet. Sc. Lett., 522, 176–185, 2019. 

Portenga, E. W. and Bierman, P. R.: Understanding Earth's eroding surface with 10Be, 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],, 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,, 2019. 

Schwanghart, W. and Scherler, D.: Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques, Earth Surf. Dynam., 5, 821–839,, 2017. 

Simoes, M., Braun, J., and Bonnet, S.: Continental-scale 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,, 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,, 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 10Be and 26Al, 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 10Be–26Al source-area signal in sediment-routing systems of arid central Australia, Earth Surf. Dynam., 6, 329–349,, 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 millennial-scale test, Global Planet. Change, 69, 48–58, 2009. 

Theodoratos, N., Seybold, H., and Kirchner, J. W.: Scaling and similarity of a stream-power incision and linear diffusion landscape evolution model, Earth Surf. Dynam., 6, 779–808,, 2018.  

Tian, D., Uieda, L., Leong, W. J., Schlitzer, W., Fröhlich, Y., Grund, M., Jones, M., Toney, L., Yao, J., Magen, Y., Jing-Hui, 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],, 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 stream-power 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, Earth-Sci. 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.,, 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. 

Short summary
We compare models of erosion against a global compilation of long-term erosion rates in order to find and interpret best-fit parameters using an iterative search. We find global signals among exponents which control the relationship between erosion rate and slope, as well as other parameters which are common in long-term erosion modelling. Finally, we analyse the global variability in parameters and find a correlation between precipitation and coefficients for optimised models.