Quantifying erosion rates and processes remains a central focus of studying the Earth's surface. Measurement of in situ–produced cosmogenic radionuclides (CRNs) enables a level of quantification that would otherwise be impossible or fraught with uncertainty and expense. Remarkable success stories punctuate the field over the last decade as CRN-based methodologies are pushed to new limits. Inherent to all is an assumption of steady-state rates and processes. This paper focuses on the use of cosmogenic 10Be and 26Al, extracted from quartz in bedrock, saprolite, and detrital material to quantify sediment production or erosion rates and processes. Previous results from two very different field areas are reviewed to highlight the potential for non-steady-state processes in shaping soil-mantled landscapes. With this potential in mind, a numerical model is presented, following a review of the CRN conceptual framework, to test the effects of non-steady-state erosion rates and processes on CRN concentrations. Results from this model focus on 10Be concentrations accumulated under modeled variations in erosion rates with different ranges, frequencies, and styles of variability. In general, the higher the maximum erosion rate, the higher the impact on the CRN concentration and, therefore, the more likely that point measurements will capture the variable signal. Conversely, the higher the frequency of erosional variation, the less likely point measures are to accurately determine rates, but the closer the inferred rate is to the mean of the long-term erosion rate. Modeling results are applicable for point-specific erosion rates, but endorse the catchment-averaged approach for determining average rates. Potentially large uncertainties emphasize the need for careful sample selection, with adequate numbers of samples collected for quantifying the processes eroding the land. The two field examples show how analyzing enough samples can define a clear soil production function despite the potential for non-steady-state processes. The model presented here is ready for application to catchment-averaged processes, as well as modeling the role of muons in variable erosion rate scenarios.

Here I briefly review the state of using CRN measurements (focusing specifically on in situ–produced 10Be and 26Al) to infer soil production and erosion rates. Using this context, I present a numerical model to test the effect of stochastic sediment production and transport processes on the accumulation of CRNs in eroding bedrock that can be either exposed or beneath a mobile soil mantle. This work was motivated by field observations from two very different field sites. First, a steep, soil-mantled hilly landscape in the Oregon Coast Range (Fig. 1) where geomorphic processes of erosion and soil production were identified to be nonuniform and potentially catastrophic (Heimsath et al., 2001b; Montgomery et al., 1998; Roering and Gerber, 2005; Schmidt, 1999). Measurements of the in situ–produced CRNs, 10Be and 26Al, from weathered bedrock beneath an actively eroding soil mantle led to determining an apparent soil production function, where soil production decreased exponentially with increasing soil thickness (Heimsath et al., 2001b). The second site, a gentle upland landscape eroding almost an order of magnitude more slowly (Fig. 2), was likely to have undergone dramatic changes in the dominant geomorphic processes due to Pleistocene climate changes (Heimsath et al., 2001a). Conclusions reached for both sites depended in part on assuming steady-state conditions for both the overlying soil mantle and the erosion rate (equivalent to the soil production rate) of the underlying bedrock.

Figure 1. Photograph of a freshly clear-cut slope in the Oregon Coast Range near the field site used by Heimsath et al. (2001b). The clear-cut reveals exposed bedrock along the sharply convex-up ridge crests and the ridge-and-valley topography characteristic of the region. Light-colored patch at the uppermost extent of the valley to the far right is bedrock exposed by recent landsliding. Scale of the clear-cut at the ridgeline is ∼200 m. Average slopes are ∼45°.

Figure 1. Photograph of a freshly clear-cut slope in the Oregon Coast Range near the field site used by Heimsath et al. (2001b). The clear-cut reveals exposed bedrock along the sharply convex-up ridge crests and the ridge-and-valley topography characteristic of the region. Light-colored patch at the uppermost extent of the valley to the far right is bedrock exposed by recent landsliding. Scale of the clear-cut at the ridgeline is ∼200 m. Average slopes are ∼45°.

Figure 2. Photograph of Frogs Hollow field area, looking north from the south side of the Bredbo River. Unchanneled valley to the left of the mid-photo ridge was sampled to yield average erosion rate shown by C on Figure 4; average and bedrock incision rates are of the Bredbo, flowing right to left in the midst of the trees in the foreground. Note outcropping tors especially prevalent upon the convex-up noses. The “tor profile” used for the data shown in Figure 6 is the highest visible tor on the ridge crest of the leftmost ridge. Results from the nuclide concentrations suggested “stripping” of the landscape at ca. 150 ka (Heimsath et al., 2001a). Relief of the unchanneled valley is ∼17 m; the tor in the middle of the photograph is about three meters high and two and a half wide.

Figure 2. Photograph of Frogs Hollow field area, looking north from the south side of the Bredbo River. Unchanneled valley to the left of the mid-photo ridge was sampled to yield average erosion rate shown by C on Figure 4; average and bedrock incision rates are of the Bredbo, flowing right to left in the midst of the trees in the foreground. Note outcropping tors especially prevalent upon the convex-up noses. The “tor profile” used for the data shown in Figure 6 is the highest visible tor on the ridge crest of the leftmost ridge. Results from the nuclide concentrations suggested “stripping” of the landscape at ca. 150 ka (Heimsath et al., 2001a). Relief of the unchanneled valley is ∼17 m; the tor in the middle of the photograph is about three meters high and two and a half wide.

Observations of how both the overlying soil thickness and the local production rate of soil, or erosion rate of exposed bedrock, may vary in a stochastic (or even predictable) ways across the landscape led to the development of the numerical model presented here to calculate CRN concentrations under non-steady-state conditions. Importantly, the model here differs from those presented by others (Bierman and Steig, 1996; Lal, 1991; Small et al., 1997) in that it integrates the governing differential equation for the CRN concentrations under changing conditions of erosion and overlying soil thickness (i.e., shielding from cosmicray penetration), using numerical techniques rather than solving the equations analytically. I present results specifically with the deviations from steady state for the Oregon and southeastern Australia cases in mind, using extensive in situ measurements of 10Be and 26Al for comparison with the modeling results. Furthermore, the model is freely available1, and has broad application as a test for the potential use of CRN methodology for tackling the more complex geomorphic problems that are currently being pursued.

To develop this model, I first review the geomorphic conceptual framework and then explicitly derive the differential equations used in CRN applications. Brief field descriptions with a summary of results from the two field areas motivating this study are then presented, with the observations and data from the respective studies driving the modeling. Specific discussion into how soil production functions and erosion rates can vary across landscapes serves not only as a review of results from different field areas, but also as a launching point for future work. I conclude, therefore, by placing the modeling results into a broader context by suggesting what the next steps are for application of such modeling efforts. Specifically, the model as developed here is evaluating how variable erosion rates affect CRN concentrations for point-specific samples without accounting for muogenic production. The next steps thus involve determining how variable and stochastic erosion rates affect CRN concentrations in sediments and also how the muogenic signature might be affected.

CONCEPTUAL FRAMEWORK AND MODEL

The Geomorphology

The focus here is on hilly and mountainous soil-mantled landscapes where the bedrock is actively converted to a continuous soil mantle (also referred to as regolith—the key is that the soil is the mobile layer). Importantly, saprolite, or weathered bedrock, is conceptualized as bedrock with the criterion that geomorphic processes do not physically mobilize it: It retains relict rock structure. When conditions of local steady state are assumed, the soil production rate equals the erosion rate, and the local soil thickness remains temporally constant (Heimsath et al., 1999). The bedrock-soil interface lowers spatially at the soil production rate, and the soil acts as a continuously moving layer removing sediment produced locally and transported from upslope (Heimsath et al., 1997) such that the lowering rate of the soil-bedrock interface is equivalent to the landscape lowering rate. Importantly, this rate can vary across the landscape, as discussed below, such that the landscape is not lowering at the same rate everywhere and is, therefore, out of dynamic equilibrium (Ahnert, 1967, 1987).

Cosmogenic nuclides are produced within the mineral grains present in the soil column as well as in a profile with depth in the underlying bedrock. Here we are specifically interested in the in situ–produced CRN concentrations in the bedrock at the soil-bedrock boundary, not in the soil. The sloping soil mantle is treated like an additional buffer against cosmic-ray penetration, and CRN production rates for a specific sample location are corrected to normalize against a flat, unburied surface (Dunne et al., 1999). Concentrations of both 10Be and 26Al measured from the bedrock or saprolite (chemically weathered, inplace bedrock) beneath the soil layer determine soil production rates (e.g., Heimsath et al., 1997, 1999; Small et al., 1999). It is especially important to recognize that under steady-state conditions the application of CRN measurements to determining soil production rates is identical to the procedure more commonly used to determine exposed bedrock erosion rates except for the explicit correction of CRN production rates for the overlying soil mantle.

Sediment produced from exposed bedrock as well as that eroded from the mobile soil layer is transported to channels and eventually removed from the landscape. Physical processes mix this sediment during transport such that a grab sample of sand from a sandbar in a channel is likely to contain sediment grains from eroding parts of the basin draining to that point. Following the well-constrained studies of Granger et al. (1996) and Bierman and Steig (1996), which showed that the CRN concentration of such a grab sample reflects a spatial average of the erosion rates acting in the basin, extensive application of this technique suggests it is robust across a wide range of geomorphic conditions (Clapp et al., 2000, 2002; Hewawasam et al., 2003; Matmon et al., 2003a, 2003b; Nichols et al., 2005; Schaller et al., 2001; Wobus et al., 2005). There are several critical assumptions, however, that were well articulated by Granger et al. (1996) and Bierman and Steig (1996). Perhaps most important are the assumptions that the collected sediments integrate the erosional processes operating across the landscape and that these processes are occurring close to steady state. The role of landslides in contributing pulses of sediment that are not spatially integrated as well as not being close to steady state is a potentially serious problem for this methodology, but testing this is beyond the scope of this paper and is being modeled with a different approach (Niemi et al., 2005). In the discussion below, I will address specifically how modeling results presented here may be used in an approach similar to Niemi et al. (2005).

The Cosmogenic Nuclide Approach

Here I review the derivation of the equations commonly used in geomorphic applications, to clarify the integration of the governing differential equation and its implications. This may be redundant to most experienced users, but the growing number of students of this technique warrants an explicit review. Similarly, while all users of CRNs apply the steady-state solution of Lal (1991), the derivation of the governing equations remains poorly understood despite the thorough review of Gosse and Phillips (2001). This derivation also builds the foundation for the numerical model presented below. Muogenic production of nuclides will be assumed to be small and will be left out of the derivation for the sake of simplicity and also ease of comparison with previous studies where muogenic production was not accounted for. It will follow that an extra term for nuclide production by muons (Gosse and Phillips, 2001; Granger and Smith, 2000; Stone et al., 1998b) is simple to add to the numerical model. This will be a critical addition for further applications of this model, as the impacts of variable erosion rates on the muogenic component of CRN concentrations is currently unknown. Following Lal (1991), for a flat bedrock target with no soil cover, the production rate of a radionuclide, P(z, t) (atoms/g/yr), declines from the surface nuclide production rate, PS, with depth, z (cm), such that

where ρ is the density of the target material (g cm−3), and Λ is the absorption mean free path for the nuclear interacting particles in the target (g cm−2). An absorption coefficient (cm−1) is defined for bedrock as µ = ρ/Λ to simplify all further equations. An additional coefficient is defined for the overlying soil mantle further dampening the production of nuclides as µS = ρS/Λ, where ρS is the soil bulk density and is assumed constant based on field measurements in the well-mixed soils of the landscapes studied here. It is assumed that the nuclides are produced only by cosmic-ray nucleons and that the nuclide production rate is constant over time. The current debate over production rate uncertainties (Clark et al., 1995; Dunai, 2000; Nishiizumi et al., 1996; Stone et al., 1998a) and the likely contribution of muons to nuclide concentrations under moderate and high erosion rates (Brown et al., 1995a; Granger and Smith, 2000; Stone et al., 1998b) are not relevant to the derivation of this model because the application is to determine relative errors resulting from non-steady-state conditions of erosion.

All nuclides considered here decay radioactively at a rate proportional to the concentration, with a constant of proportionality, λ, which is specific to the nuclide. Thus, the concentration of nuclides in the rock horizon beneath a soil mantle with thickness h can be modeled as

For most applications of geomorphology, the concentration of nuclides is measured from the surface of the target material in an actively eroding environment. To evaluate equation 2 over the exposure history of the sample, the position with depth, z(t), at time t of the sample must therefore be defined as a function of the erosion rate. The simplest case is that exposed rock has eroded at a steady rate over its exposure history (Lal, 1991; Nishiizumi et al., 1989). If the original position, generally assumed to be where there is no penetration of cosmic rays, of the sample under the ground is some depth, z0, and the constant erosion rate is ϵ, then the sample position with time is

Substituting equation 3 into equation 2 provides a governing differential equation for a bedrock surface, either below soil cover or not (in which case h = 0), that is eroding at a constant rate over time, and considering N as a function of time only,

Grouping similar terms and simplifying for now by letting h = 0, yields

Substituting this quantity into the above equation and multiplying both sides by eλt yields

Grouping the exponentials on the right side of the equation into terms with and without the time variable,

and integrating both sides with respect to time, t, yields, with N0 as a constant,

Solving for N(t) thus results in an expression that only requires solving a simple integral,

Under the steady-state assumption (i.e., that the erosion rate, ϵ, is constant) the integral can be solved analytically. If ϵ is a function of time, however, the integral must be solved numerically, as described below.

Assuming, for now, that the steady-state assumption holds, the concentration at a specific time, T, is determined by evaluating equation 9 at t = T such that

At time T, a particle starting at depth z0 will have risen to a depth z =z0– ϵT. Substituting z0 = Z + ϵT into equation 10 leads to

The integral can now be solved analytically as follows:

Equation 12 provides a closed-form solution for concentration under the steady-state assumption (i.e., constant erosion, ϵ, over the sample exposure history). Finally, this reduces to the widely used solution as T >> (λ + µϵ)−1 such that (e.g., Lal, 1991)

which can be solved for the steady-state erosion rate by rearrangement:

Equation 14 is mostly used to interpret nuclide concentrations in the context of a surface eroding by processes assumed to be operating relatively constantly over time, such as grain-grain spallation of sandstone, thin exfoliation of exposed granite, and biogenic soil production from saprolite beneath an active soil mantle (in which case z is very close to zero, but the depth term for soil, h, will need to be accounted for in setting the nuclide production rate). This relationship is also used to interpret concentrations extracted from detrital sediments collected from channels to determine erosion rates averaged across the entire drainage basin. The model presented below does not, however, spatially integrate point-specific samples to determine basin-wide concentrations under variable erosion rate conditions. Recent modeling work suggests that the processes eroding the basin need not be in steady state as long as the basin is large enough to integrate signals from both the low-magnitude/high-frequency events (e.g., processes approaching steady state) and the low-frequency/high-magnitude events (e.g., landslides) (Niemi et al., 2005). While this model is not yet fully tested, the conclusions are intuitive and are discussed below.

The Numerical Model

I now extend the model equation 4, such that the effects of variable conditions in erosion and overlying soil depth can be evaluated. It will become clear that a temporally variable soil depth results in the same scenario as a variable erosion rate scenario, as the soil production rate, which depends on the soil depth, is equivalent to the local erosion rate. The model is then applied to hypothetical cases of temporally variable erosion rates, based on the field observations of Heimsath et al. (2001a, 2001b), as well as extended to potentially extreme scenarios to compare with the earlier models testing the effects of extreme events on nuclide concentrations (Bierman and Steig, 1996; Lal, 1991; Small et al., 1997).

Under conditions of non-steady-state erosion (i.e., the erosion rate, now denoted as ϵ(t), is a function of time and can be modeled with different approximations), the initial differential equation for nuclide concentration can be written as

Following the same derivation as done above leads to a similar expression as in equation 9:

Note that, unlike the steady-state condition, the above integral does not have a closed-form solution for nontrivial erosion rates, ϵ(t). This simply is because the erosion rate is a function of the integrating variable, t. To solve for the non-steady-state condition of equation 15, a numerical approximation of the differential equation is, therefore, used. Specifically, by summing both sides of the differential equation, the concentration at time T is approximated by the following:

where N(0) = N0 (the initial concentration, typically equal to zero), and the depth, z(t), at time is t0

With this numerical solution, the nuclide concentration can be estimated for any arbitrary erosion function, ϵ(t), that might approximate potential temporal variations in erosion observed in the field. It is important to note that as constructed here, the model is evaluating point-specific concentrations and inferred erosion rates. The obvious next step is to incorporate this point-specific model into a model that integrates the concentrations derived from all points upslope of any given catchment outlet. The approach could be similar to the one used by Niemi et al. (2005), except that every pixel in the catchment digital elevation model (DEM) would be assigned a CRN concentration based on its variable erosion rate history. Such a model would have to track the production and transport of sediment from all points and is, therefore, beyond the scope of this paper.

Equation 17 can also be modified easily to account for nuclide production by muons, as suggested above. In this case, an additional production term would be added to the right side and the modeling code adapted with the simple addition of the muogenic production term. For the purposes of this paper, this extension is not needed, which also makes results presented here easier to compare to other studies that have not accounted for muogenic production. Similarly, the dampening of nuclide production due to the overlying soil thickness, included in equation 1, is trivial in comparison to the oscillations in erosion rates modeled here and is therefore left out for simplicity. Specifically, if overlying soil thickness is varying stochastically or regularly, then the soil production rate will vary as a function of the change in thickness. Since soil production rates equal erosion rates, for the purposes of this paper it is sufficient to model variation in erosion rates. Specific forms for the variable erosion rate scenarios will be described following the field descriptions and summary of previous results.

FIELD AREAS AND PREVIOUS RESULTS

In this section, I briefly review the pertinent aspects of the field areas that motivated the model and present the soil production and erosion rates determined for both sites using both 10Be and 26Al concentrations from exposed and soil-mantled, weathered bedrock as well as stream sediments and bedrock. Both field sites differ in critical ways from those examined by similar studies across completely soil-mantled landscapes (e.g., Heimsath et al., 1997, 2000). Some of these differences will be reviewed here and a comparison between the sites will be drawn after application of the model, with specific discussion of the form and magnitude of the respective soil production functions.

Oregon Coast Range

The Oregon Coast Range is a well-studied region of the Pacific Northwest that exemplifies humid-temperate, soil-mantled landscapes where stochastic as well as steady-state erosional processes drive landscape development (e.g., Dietrich et al., 2003; Heimsath et al., 2001b; Montgomery et al., 1997, 1998, 2000; Roering and Gerber, 2005; Roering et al., 1999, 2001). The actively eroding landscape is characterized by ridge-and-valley topography (Fig. 1), and while there is some hint of long-term equilibrium in erosion rates, there is also clear indication of local disequilibrium across the landscape (Heimsath et al., 2001b; Reneau and Dietrich, 1991; Roering and Gerber, 2005). Bedrock is a relatively undeformed Eocene turbidite sandstone and siltstone called the Tyee Formation (Heller et al., 1985; Lovell, 1969; Snavely et al., 1964) and outcrops in unweathered units along sharply convex ridge crests. The soil-mantled ridges (a few meters wide at the crest and over one hundred meters down the crest to the channel) have high spatial variability of soil thicknesses such that there is little systematic variation between soil depth and topographic position (Heimsath et al., 2001b). Soil depths on slopes vary between zero and one to two meters in the upland areas with relatively unweathered bedrock often exposed in recently evacuated debris-flow or landslide scars in the hollows. Creeks and rivers draining the region are mixed alluvial-bedrock and are thought to be incising at roughly the same rate as the long-term average uplift rate of 100–300 m/m.y. documented by marine terrace records (Kelsey and Bockheim, 1994; Kelsey et al., 1994). Rainfall is high, averaging about two meters a year, and the dense forest, dominated by Douglas fir that can grow to be over forty meters tall, is actively cleared for timber harvesting.

Results relevant to this paper include a cosmogenic nuclide (both 26Al and 10Be) determined soil production function (Fig. 3) (Heimsath et al., 2001b). Apparent soil production is shown to be an inverse exponential function of overlying soil thickness, with a maximum inferred soil production rate of 268 m/m.y. occurring under soil depth that approaches zero. This soil production function is determined from saprolite samples (27 samples and eight duplicates) collected under soils thicker than ∼20 cm, while erosion rate data from exposed bedrock (three samples and one duplicate) as well as from under soils thinner than 20 cm (three samples and two duplicates) average 160 m/m.y. Two samples and one duplicate of stream sediments determined a catchment-averaged erosion rate of 117 m/m.y. for the field area, which agreed well with other estimates of average erosion rates for the region (Reneau and Dietrich, 1991), but is much lower than the postfire rates determined from recently reanalyzed data (Roering and Gerber, 2005). Stochastic processes of tree throw and shallow landsliding can change dramatically local soil thicknesses over short timescales and, therefore, potentially affect the inferred rates of soil production, as well as the inferred long-term erosion rates (Heimsath et al., 2001b). These rates, as well as the apparent soil production function, enable modeling of landslide susceptibility for the region, which is especially important given the potentially fatal implications of debris flows on local homeowners. Similarly, long-term erosion rates from cosmogenic nuclides are significantly lower than short-term rates from sediment trap studies and suggest an increase in sediment removal from the landscape associated with timber harvesting. It is therefore both scientifically interesting and relevant to land management to assess the robustness of these data.

Figure 3. Apparent soil production rates (m/m.y.), calculated from the in situ–produced radionuclides 10Be and 26Al extracted from bedrock samples versus the observed normal soil depth, H, in cm. The soil production function, ϵ(H) = (268 ± 25) e−(0.03 ± 0.02)H, is the variance-weighted least squares fit to the soil production rates under soil depths greater than 15 cm (plotted with black filled circles). Soil production rates for the shallower samples and the exposed bedrock samples (plotted with gray filled squares) were not used to determine this function because of lithologic differences, as discussed by Heimsath et al. (2001b). The average erosion rates for the catchment from nuclide concentrations from three stream sediment samples are plotted with open squares on the far right of the plot, as labeled. Rates are calculated from the concentrations of both nuclides with a few exceptions. Error bars are 1σ propagated from accelerator mass spectrometry (AMS), atomic absorption (AA), bulk density, absorption mean free path, and soil depth uncertainties. Modified from Heimsath et al. (2001b).

Figure 3. Apparent soil production rates (m/m.y.), calculated from the in situ–produced radionuclides 10Be and 26Al extracted from bedrock samples versus the observed normal soil depth, H, in cm. The soil production function, ϵ(H) = (268 ± 25) e−(0.03 ± 0.02)H, is the variance-weighted least squares fit to the soil production rates under soil depths greater than 15 cm (plotted with black filled circles). Soil production rates for the shallower samples and the exposed bedrock samples (plotted with gray filled squares) were not used to determine this function because of lithologic differences, as discussed by Heimsath et al. (2001b). The average erosion rates for the catchment from nuclide concentrations from three stream sediment samples are plotted with open squares on the far right of the plot, as labeled. Rates are calculated from the concentrations of both nuclides with a few exceptions. Error bars are 1σ propagated from accelerator mass spectrometry (AMS), atomic absorption (AA), bulk density, absorption mean free path, and soil depth uncertainties. Modified from Heimsath et al. (2001b).

Southeastern Australian Highlands

Rolling soil-mantled hills dotted with outcropping tors are characteristic of the southeastern Australian highlands (Fig. 2). Land use results often in clearing of the dominant sclerophyll forests and subsequent gullying across the stony hills. Where undisturbed, shallow soils across the ridges are typically free from rilling and grade into thick, undissected colluvial fills in the hollows. Bedrock varies across the region between Ordovician metasediments and Devonian granites. Only the granites contain enough quartz for cosmogenic nuclide analyses, though abundant quartz veins in the metasediments can also be used. Gradients are gentle, averaging less than 20° in comparison with the nearly 45° slopes common in the Oregon Coast Range forests. Rainfall is a fraction of the Oregon Coast Range site and falls throughout the year to average between half to three-quarters of a meter a year. Temperatures average ∼18–22 °C in the summer and 5–8 °C in the winter, with short periods of freezing temperatures and minor snowfall events. Soil production and erosion are primarily due to biogenic processes and creep, with some evidence of overland flow (Heimsath et al., 2001a).

Results relevant here are the soil production and erosion rates determined from cosmogenic nuclide (both 26Al and 10Be) concentrations (Fig. 4). Five samples of the granitic saprolite from under different soil depths led to the inference of an apparent soil production function for the Frogs Hollow site. Noticeably, the slope of the function is significantly steeper than the function determined in Oregon, as well as those reported for more completely soil-mantled landscapes in both California (Heimsath et al., 1997, 1999) and the coastal lowlands of southeastern Australia (Heimsath et al., 2000). This steep function was subsequently expanded with the collection of significantly more data from another highland site, which also led to the confirmation of the soil production function first reported from the lowland site (Fig. 5) (Heimsath et al., 2006). Importantly, the addition of new data from a greater range of soil depths and additional landscapes showed the robustness of the soil production function and its applicability across the soil-mantled landscapes of southeastern Australia. In any case, the initial quantification of soil production rates from the highland site enabled numeric modeling of landscape evolution and soil development for the region that was supported by field observations (Heimsath et al., 2001a). Soil production rates inferred from nuclide concentrations ranged from ∼50 m/m.y. under 25 cm of soil to ∼11 m/m.y. under 65 cm of soil. Average erosion rates for the landscape, determined from both 26Al and 10Be concentrations in sediments, were ∼15 m/m.y. (denoted by C and Riv on Fig. 4), roughly an order of magnitude slower than the Oregon Coast Range. Incision rates for the river draining the landscape were determined from three samples (both 26Al and 10Be) from the fluvially sculpted and polished bedrock of the active channel bed and averaged 9 m/m.y., suggesting that the soil-mantled tributary catchment is eroding slightly more rapidly. The slightly higher rates from the tributary were explained by suggesting the landscape is responding to a pulse of incision, inferred from knickpoints in the channel's long profile, potentially driven by some period of increased erosion.

Figure 4. Apparent soil production rates versus observed normal soil depths, H, for Frogs Hollow. Filled circles show the soil production rates versus soil depth from both 26Al and 10Be concentrations. The variance-weighted best fit to the filled circles is ϵ(H) = (143 ± 20) e−(0.042 ± 0.003)H, where soil production is in m/m.y. and soil depth is in cm. Gray filled squares plotted at zero depth are exposed bedrock samples at the ground surface, inferred to be emerging core stones. Error bars show 1σ error propagated from all sources of error (uncertainty in AA, AMS, bulk density, and soil depth measurements, and the attenuation length of the cosmic rays) except the uncertainty in nuclide production rates. Rates plotted on the far right of the plot, as labeled, show average erosion rates from (1) sediments from the Frogs Hollow catchment area (half-filled black square, labeled C), (2) sediments from the Bred-bo River upstream of Frogs Hollow (gray diamond, labeled Riv), and (3) the bedrock incision rates of the Bredbo (four samples, gray halffilled squares, labeled Brk). Modified from Heimsath et al. (2001a).

Figure 4. Apparent soil production rates versus observed normal soil depths, H, for Frogs Hollow. Filled circles show the soil production rates versus soil depth from both 26Al and 10Be concentrations. The variance-weighted best fit to the filled circles is ϵ(H) = (143 ± 20) e−(0.042 ± 0.003)H, where soil production is in m/m.y. and soil depth is in cm. Gray filled squares plotted at zero depth are exposed bedrock samples at the ground surface, inferred to be emerging core stones. Error bars show 1σ error propagated from all sources of error (uncertainty in AA, AMS, bulk density, and soil depth measurements, and the attenuation length of the cosmic rays) except the uncertainty in nuclide production rates. Rates plotted on the far right of the plot, as labeled, show average erosion rates from (1) sediments from the Frogs Hollow catchment area (half-filled black square, labeled C), (2) sediments from the Bred-bo River upstream of Frogs Hollow (gray diamond, labeled Riv), and (3) the bedrock incision rates of the Bredbo (four samples, gray halffilled squares, labeled Brk). Modified from Heimsath et al. (2001a).

Figure 5. Apparent soil production rates plotted for all four southeastern Australia field sites discussed in Heimsath et al. (2006). Note the remarkable overlap between the two highland sites, Brown Mountain and Frogs Hollow (FH), which are ∼50 km apart, but at roughly the same elevation and in the same physiographic region. Also note that when soil production rates from these highland sites are combined with the data from Heimsath et al. (2000) at the base of the escarpment (NR), as well as from the coastal lowlands (Snug), the soil production function appears to be quite robust. Erosion rates of bedrock (BRK) exposed at the soil surface are plotted for all sites with different symbols. Variance-weighted best fit to all the soil production data plotted here is ϵ(H) = (53 ± 2) e−(0.022 ± 0.001)H. Error bars show the same sources of uncertainty as in Figures 3 and 4. Modified from Heimsath et al. (2006).

Figure 5. Apparent soil production rates plotted for all four southeastern Australia field sites discussed in Heimsath et al. (2006). Note the remarkable overlap between the two highland sites, Brown Mountain and Frogs Hollow (FH), which are ∼50 km apart, but at roughly the same elevation and in the same physiographic region. Also note that when soil production rates from these highland sites are combined with the data from Heimsath et al. (2000) at the base of the escarpment (NR), as well as from the coastal lowlands (Snug), the soil production function appears to be quite robust. Erosion rates of bedrock (BRK) exposed at the soil surface are plotted for all sites with different symbols. Variance-weighted best fit to all the soil production data plotted here is ϵ(H) = (53 ± 2) e−(0.022 ± 0.001)H. Error bars show the same sources of uncertainty as in Figures 3 and 4. Modified from Heimsath et al. (2006).

What makes this highland site so interesting from a non-steady-state erosional perspective are the results from a profile of bedrock samples collected from the side of an outcropping granite tor. Our interest in long-term landscape development led to the development of a “Kuniometer,” named for Kuni Nishiizumi, who came up with the idea, which is a profile sampled for nuclide concentrations from ground level to the highest point on a bedrock tor exposed in a soil-mantled landscape (Heimsath et al., 2000, 2001a). The idea is simple. Concentrations of nuclides from the profile sampled will depend on the relative rates of lowering between the soil-mantled landscape and the eroding bedrock that is outcropping. Steady-state lowering will lead to an increase in nuclide concentration with height above ground level, while some dramatic change in erosion that might have left the bedrock abruptly exposed will lead to a constant nuclide concentration with height. While the lowland site led to support for steady state (Heimsath et al., 2000), results from Frogs Hollow led to the suggestion of a dramatic stripping event in the late Pleistocene (Heimsath et al., 2001a). Specifically, the observed nuclide concentration profile suggests that a period of increased erosion ca. 150,000 yr ago led to the rapid exhumation of the tor by the removal of over two meters of easily erodible saprolite (Fig. 6). This evidence for dramatic changes in erosional regime leads to obvious concern for the potential effects on nuclide concentrations in the saprolite and, therefore, the inferred soil production rates.

Figure 6. Observed concentrations of 10Be (note, results from 26Al yielded the same story) for tor profile sampled for exposed bedrock. Nuclide concentrations are plotted as closed black circles against height above the present ground surface, with measured nuclide concentrations normalized to sea level and the error bars showing 1σ propagated from AMS, AA, bulk density, and absorption mean free path uncertainties. The black dashed line plots the best-fit model prediction for the steady-state scenario from the upper Bega Valley site of Nunnock River of Heimsath et al. (2000), scaled to fit the Frogs Hollow data more closely and showing the marked deviation from observed concentrations. These observations of concentrations against height above ground are best explained by a model that posits complete emergence of the tor early in its exposure history, ca. 150 ka. Modified from Heimsath et al. (2001a).

Figure 6. Observed concentrations of 10Be (note, results from 26Al yielded the same story) for tor profile sampled for exposed bedrock. Nuclide concentrations are plotted as closed black circles against height above the present ground surface, with measured nuclide concentrations normalized to sea level and the error bars showing 1σ propagated from AMS, AA, bulk density, and absorption mean free path uncertainties. The black dashed line plots the best-fit model prediction for the steady-state scenario from the upper Bega Valley site of Nunnock River of Heimsath et al. (2000), scaled to fit the Frogs Hollow data more closely and showing the marked deviation from observed concentrations. These observations of concentrations against height above ground are best explained by a model that posits complete emergence of the tor early in its exposure history, ca. 150 ka. Modified from Heimsath et al. (2001a).

RUNNING THE MODEL

Given these two dramatically different field examples (Oregon Coast Range and southeastern Australian highlands) for how erosion rates might vary over time, there are reasonable constraints for how to model the effects of potential non-steady-state scenarios on point samples used to determine soil production or erosion rates. Several forms of oscillating input erosion rates can be chosen from and are used to model (which is encoded in Matlab® and freely available upon request2) conditions of variable erosion. Additionally, there are two parameters that govern the resolution and run time of the model: the time step and the maximum time. The maximum time should be set such that the end of the model run reaches a steady-state nuclide concentration, or, in the case of an oscillating concentration, a steady-state oscillation with the imposed variable erosion rate. Conceptually, for steady-state erosion scenarios, the faster the erosion rate, the shorter the time to steady-state nuclide concentrations. The time step is the interval over which the model evaluates the numerical integration. Setting the time step governs, therefore, the efficiency of the numerical computations as well as the resolution of the integration. Naturally, the shorter the time step, the longer the computations take to complete, but the greater the resolution of the model output. For higher rates of erosion, the time step needs to be shortened to insure analytical accuracy. Setting both the time step and the maximum time involves editing the source code for the model, and should be done prior to choosing the form of the model and the parameters governing the input for the model.

Setting the model parameters and the modeled erosion scenario is done by the command line in Matlab®. Each variable erosion rate functional form is based on field observations as well as a range of studies that have either speculated or documented the ways in which erosion rates and processes can change across a landscape. To check the mathematics of the numerical integration, the model also enables input of either zero erosion or steady-state erosion rates, functions that result in steadily decreasing steady-state nuclide concentrations with increasing erosion rate (Fig. 7) (Lal, 1991; Cerling and Craig, 1994). The variable erosion models are mathematically represented as the following: (1) step function, (2) rectified full wave, (3) square wave, (4) exponential, and (5) sawtooth. Examples of each of these functions are plotted in Figure 8 for erosion oscillations between 100 m/m.y. and 1000 m/m.y., showing two full cycles. Parameters governing the frequency and magnitude of each of these functions can be adjusted according to predicted, or hypothesized, variations in erosion rates.

Figure 7. Modeled 10Be accumulation in rock surface samples as a function of time and steady-state erosion rate. Concentration curves are labeled with the input erosion rate written above the resulting curve, varying from 0 to 1000 m/m.y. Note that the faster the erosion rate, the shorter the time needed for the nuclide concentration to reach steady state, and that this concentration is lower. Concentrations here and for the modeling results below are calculated for a sea-level 10Be production rate of 6 atoms/g/yr, and would scale for production rates depending on sample location and shielding. Note that this figure is similar to ones commonly shown in textbooks describing the cosmogenic nuclide methodology and is implicit in the plots of Lal (1991), then shown plotted in this form by Cerling and Craig (1994).

Figure 7. Modeled 10Be accumulation in rock surface samples as a function of time and steady-state erosion rate. Concentration curves are labeled with the input erosion rate written above the resulting curve, varying from 0 to 1000 m/m.y. Note that the faster the erosion rate, the shorter the time needed for the nuclide concentration to reach steady state, and that this concentration is lower. Concentrations here and for the modeling results below are calculated for a sea-level 10Be production rate of 6 atoms/g/yr, and would scale for production rates depending on sample location and shielding. Note that this figure is similar to ones commonly shown in textbooks describing the cosmogenic nuclide methodology and is implicit in the plots of Lal (1991), then shown plotted in this form by Cerling and Craig (1994).

Figure 8. Input erosion rate functions plotted against time to show representative examples. Here the maximum erosion of the oscillations is set for 1000 m/m.y., while the “pedestal” is set at 100 m/m.y. Period is set for 15 k.y. in this example. A: Full-wave rectification. Erosion increases gradually to a peak, then drops back down to the background rate before rising again. B: Square wave. Similar to A, except that rise and fall of erosion is abrupt, best modeled with a square wave. Also, the time spent at the higher erosion rate is longer. C: Exponential. Erosion rate increases exponentially with time to a peak and then abruptly drops back to the base rate. D: Sawtooth. Erosion rate increases linearly with time to a peak and then drops abruptly back to the base rate. Both of the latter two functions can also be reversed such that the increase is abrupt and the return to base conditions is more gradual. For the purposes of the modeling examples here, the conclusions do not change.

Figure 8. Input erosion rate functions plotted against time to show representative examples. Here the maximum erosion of the oscillations is set for 1000 m/m.y., while the “pedestal” is set at 100 m/m.y. Period is set for 15 k.y. in this example. A: Full-wave rectification. Erosion increases gradually to a peak, then drops back down to the background rate before rising again. B: Square wave. Similar to A, except that rise and fall of erosion is abrupt, best modeled with a square wave. Also, the time spent at the higher erosion rate is longer. C: Exponential. Erosion rate increases exponentially with time to a peak and then abruptly drops back to the base rate. D: Sawtooth. Erosion rate increases linearly with time to a peak and then drops abruptly back to the base rate. Both of the latter two functions can also be reversed such that the increase is abrupt and the return to base conditions is more gradual. For the purposes of the modeling examples here, the conclusions do not change.

Specifically, changing the amplitude of the function defines the maximum erosion rate for the specified model scenario (units of m/m.y.). Choosing the frequency for the oscillating functions divides the maximum time into the number of cycles as defined by the input. Finally, there is the option to choose a “pedestal” for the functions, which sets the minimum erosion rate for the cycles and defines that the oscillating erosion rates are nonzero, if such is the hypothesized scenario. These three parameters guide the model to test any potential scenario for variable erosion rates. For example, drawing on the potential variations for the Oregon Coast Range site described above, the recent work on the effects of fire (Roering and Gerber, 2005) on erosion rates can be evaluated in the context of potential effects on the nuclide concentrations used to determine soil production and average erosion rates. Roering and Gerber (2005) report that postfire erosion rates exceed long-term rates by a factor of six and that fire frequency is on the order of 100–200 yr from early to late Holocene, respectively. Modeling this scenario means positing an input function for the variable erosion rate, defining the maximum possible rate under the oscillating conditions, setting the frequency of fire “events” and setting the minimum, or background, erosion rate. For this example, an exponential scenario would be chosen following a conceptual model for impact of fire on sediment flux with time (Swanson, 1981), as reported in Roering and Gerber (2005). Maximum erosion rates would be set to range from 600 m/m.y. to as high as 1800 m/m.y. to capture the range of potential background erosion rates. Minimum erosion rates would be set to range from the long-term average of ∼100 m/m.y. to 20 m/m.y., the minimum reported soil production rate from Heimsath et al. (2001b)(Fig. 3).

Motivation for the other variable erosion scenarios is also observation-based. The step function, first modeled by (Lal, 1991), then expanded by Small et al. (1997) to show the effect of exfoliation sheet size on inferred erosion rates, would, for an applicable example, be useful for the highland, southeastern Australia (Frogs Hollow) site. This function would test specifically the effect of a potential Pleistocene “stripping” event as suggested by the tor data discussed above. I do not discuss results from modeling with this scenario as do both of the above studies, as well as the Bierman and Steig (1996), flesh out the important points quite well. The step function scenario might also be modeled by using either a square wave (as in Bierman and Steig, 1996), or a rectified full-wave function with a suitably long period and different rates of maximum erosion. Both the exponential and the sawtooth functions can also be used to model the potential for variable soil cover thickness, punctuated by erosional events effectively stripping the soil. The idea here being that these stripping events would strip the soil, increase the erosion rate to the maximum soil production rate, and the site would gradually return to a local steady state as soil thickness recovered. This return of soil thickness to some steady-state value could be either exponential or linear (sawtooth), which would drive the underlying erosion rate. Reasonable maximum rates of erosion for all the functions can be set either by some governing observation—e.g., the fire example—or by positing some potentially high rate and evaluating its effects on the nuclide concentration. Similarly, natural choices for setting the frequency of events would include fire or landslide recurrence intervals, or climatic change cycles.

MODELING RESULTS

Using the input erosion rates that oscillate as a function of time, as shown for example in Figure 8, the model computes the accumulation of the cosmogenic radionuclide, 10Be, at a point. Perhaps the most relevant example of interest here, specifically for the Oregon Coast Range site, is when erosion rate can vary between 100 and 1000 m/m.y., as shown in Figure 8. The modeled 10Be accumulations in samples subject to these oscillating erosion rates are shown in Figure 9, with Figures 9A and 9B showing the response to sudden changes in erosion, potentially due to a sudden change in climate or land use. Figure 9A shows how concentration varies under the full-wave rectification scenario. Nuclide concentration at a point increases and decreases gradually to and from a peak, which corresponds to the end of the time spent under the base erosion rate of 100 m/m.y. The lowest concentrations occur under the highest erosion rate. Figure 9B shows concentration varying under the square-wave oscillation, which shows a pattern similar to the full wave, except that the drop in concentration is more abrupt and the time spent at the lowest concentration set by the highest erosion rate is longer. The increase in nuclide concentration begins with the drop from 1000 m/m.y. to 100 m/m.y., while the rapid drop in concentration begins with the sudden increase in erosion rate.

Figure 9. Modeled 10Be accumulation in point samples (either exposed rock or at the soil-rock interface) subject to the oscillating erosion rates of Figure 8. A: Full-wave rectification. Concentration increases and decreases gradually to and from a peak. B: Square wave. Similar to A, except that the drop in concentration is more abrupt and the time spent at the base concentration set by the base rate is longer. C: Exponential. The increase and decrease of concentration resembles a full-wave rectification, reflecting the exponential decrease in nuclide production with depth in the sample. D: Sawtooth. Nuclide concentration increases more rapidly and decreases exponentially.

Figure 9. Modeled 10Be accumulation in point samples (either exposed rock or at the soil-rock interface) subject to the oscillating erosion rates of Figure 8. A: Full-wave rectification. Concentration increases and decreases gradually to and from a peak. B: Square wave. Similar to A, except that the drop in concentration is more abrupt and the time spent at the base concentration set by the base rate is longer. C: Exponential. The increase and decrease of concentration resembles a full-wave rectification, reflecting the exponential decrease in nuclide production with depth in the sample. D: Sawtooth. Nuclide concentration increases more rapidly and decreases exponentially.

Figures 9C and 9D show how concentrations vary with oscillating erosion rates that are likely to be similar to changing soil thickness at a given point. The oscillating erosion conditions can also be reversed, to show an immediate increase followed by an exponential or linear decrease in erosion. In Figure 9C, the response to the exponential change in erosion rate, the increase and decrease of concentration resembles a full-wave rectification, reflecting the exponential decrease in nuclide production with depth in the sample. Specifically, the peak in nuclide concentrations is reached just after the increase in erosion rate begins and subsequently decreases smoothly to the lowest concentration corresponding to the highest erosion rate. With the immediate decrease in erosion rate, the nuclide concentration increases steadily. Figure 9D shows a similar change in concentration occurring under a sawtooth variation in erosion, though with a more rapid increase and a lower peak concentration despite the same maximum input erosion rate.

Where this modeling exercise becomes interesting is using the modeled concentrations computed under the variable input erosion rates to infer an erosion rate. The idea here is to replicate the collection of a sample with the subsequent determination of an erosion rate from the measured nuclide concentration, comparing the inferred erosion rate with the “real” erosion rate, which, in this case, is the input erosion function. This is valid specifically for a point sample. To apply this exercise to the basin-averaged problem requires coupling this model for individual points with a model keeping track of how the sediment is generated and transported out of the basin (Niemi et al., 2005). In Figure 10, the variable erosion rates shown in Figure 8 are plotted again in the same order and at the same magnitude, and are now overlain by the inferred erosion, shown by dashed red curves, rates calculated from the modeled nuclide concentrations shown in Figure 9. There is a remarkable mirroring of the input erosion rates by the inferred rates, with slight offsets evident in all the predictions. In Figures 10A–10C, there is a factor of two overestimation of even the base rate of erosion, and Figure 10B shows a 10% overestimation of the peak erosion rates under the square-wave oscillations. Figure 10D shows that under a sawtooth oscillation the base rate is never inferred from the nuclide concentrations, while the peak rate seems to be well captured. For any given time for the full-wave scenario, Figure 10A, the inferred erosion rate is off by an average of 49%, with a standard deviation of 83%. The square-wave input function, Figure 10B, is less well captured at any given time, with an average deviation of 77% in the inferred erosion rate, with a standard deviation of 154%. With the exponential function, Figure 10C, the average error is 44%, with a standard deviation of 110%. Not surprisingly, given the closer fit of the two curves in Figure 10D, the inferred erosion rate is on average 28% off from the input rate, with a standard deviation of 108%.

Figure 10. The input erosion rates modeled in Figure 8 are plotted again in the same order with the solid black curves, and are overlain by the inferred erosion rates, shown by dashed red curves, rates calculated from the modeled nuclide concentrations shown in Figure 9.

Figure 10. The input erosion rates modeled in Figure 8 are plotted again in the same order with the solid black curves, and are overlain by the inferred erosion rates, shown by dashed red curves, rates calculated from the modeled nuclide concentrations shown in Figure 9.

While this result does not bode well for point-specific sample collection, it is worth thinking about how such erosion rate variations manifest themselves in the long-term signal and could potentially be captured by sampling catchment sediments. The implication of these imposed variations in erosion rates over the long term is that the average rate in input erosion for the full wave shown in Figure 10A is 384 m/m.y., while the average rate inferred from the nuclide concentrations over the same full cycle is 406 m/m.y. This represents a surprisingly small difference of 5.6% for the inferred erosion rate. For the square-wave scenario of Figure 10B, the long-term averages are 545 and 586 m/m.y., respectively for input and inferred rates, with the inferred overestimating the input by 7.6%. The long-term average for the input exponential variation, Figure 10C, is 242 m/m.y., while the long-term average for the inferred rate is 249 m/m.y., meaning an overestimation of only 5.2%. Despite the factor of four-plus overestimation of the low erosion rates due to the sawtooth variation of Figure 10D, the difference between the long-term averages is only 2.7%, with the inferred average rate of 579 m/m.y., compared to the input average rate of 551 m/m.y. These long-term average comparisons are most meaningful when considering the potential for widely variable erosion rates across a catchment where the inferred erosion rate is seeking to capture the spatial average for the basin. These comparisons are not as meaningful for determining whether a point-specific sample is likely to have captured an average rate.

The near mirroring of even the highest erosion rates for the oscillation between 100 and 1000 m/m.y. changes with a reduction in the magnitude of the amplitude of variation, or the maximum erosion rate that the sample is subjected to. When the same oscillating erosion functions are now set to range from 10 to 100 m/m.y., potentially a more reasonable scenario for the southeastern Australia site, the erosion rates inferred from the modeled nuclide concentration do not show the same dramatic fluctuations as the input rates. An additional way that the input erosion rate scenarios can change the modeled nuclide concentrations and, therefore, the inferred erosion rates, is with the frequency of oscillation. Figure 11 shows a reduction by an order of magnitude of the range of erosion rates and also a selection of some different modeled frequencies to illustrate what can now be expressed as generalizable results. The time window plotted reflects the longer time needed to achieve steady-state concentrations under the lower input erosion rates (see Fig. 7).

Figure 11. Similar to Figure 10, except that the range of erosion rate variation is from 10 to 100 m/m.y., and the frequency of cycles is variable for each scenario. Input erosion rates are shown by solid black lines, and erosion rates inferred from the modeled nuclide concentrations are shown by the dashed red curves. A: Full-wave rectification with the same 15 k.y. period as Figure 10A. B: Square wave with a period of 10 k.y. C: Exponential variation with a period of 20 k.y. D: Sawtooth variation with a period of 2 k.y.

Figure 11. Similar to Figure 10, except that the range of erosion rate variation is from 10 to 100 m/m.y., and the frequency of cycles is variable for each scenario. Input erosion rates are shown by solid black lines, and erosion rates inferred from the modeled nuclide concentrations are shown by the dashed red curves. A: Full-wave rectification with the same 15 k.y. period as Figure 10A. B: Square wave with a period of 10 k.y. C: Exponential variation with a period of 20 k.y. D: Sawtooth variation with a period of 2 k.y.

Figure 11 plots different frequencies for the input varying erosion rates, with Figure 11A showing the same frequency for the full-wave rectification as the input for Figure 10A. Figure 11B increases the frequency for the square wave by 30% to show three full cycles in 30 k.y., while Figure 11C decreases the frequency for the exponential function to one cycle every 20 k.y., and Figure 11D increases the frequency to one cycle every 2 k.y. for the sawtooth function, all for the same range in input erosion rates. The results on the inferred erosion rates, shown also with the red dashed lines in Figure 11, are illustrative and applicable for each of the input scenarios. Most immediately striking is that the inferred erosion rates do not cycle as strongly as they do when the input erosion rate increases to a higher rate. The second most apparent aspect is that the high-frequency change of the sawtooth input has very little effect on the inferred erosion rate at any time, meaning that the steady-state nuclide concentration is changing little under the input erosion rate variations. That this is reflected in the fact that the average for both the input and the inferred rates for Figure 11D is ∼60 m/m.y. is surprising and has implications for the catchment-averaged sampling methodology, similar to the above scenario and as discussed below. The variation in rates shown by the square wave, Figure 11B, result in long-term averages that are also identical, with both the input and the inferred rates averaging 55 m/m.y. The long-term average inferred rate of 40 m/ m.y. is only 5.2% higher than the input average of 38 m/m.y. for the full-wave rectification of Figure 11A, while the long-term average for the lower-frequency exponential input function is 26 m/m.y., compared to the 29 m/m.y. inferred from the resulting nuclide concentrations, an overestimate of 11%.

As expected from visual inspection of the plots shown by Figure 11, the inferred erosion rate at any given time is, however, not as likely to capture the variable input erosion as it was for the scenarios of Figure 10. Specifically, with the full-wave variation of Figure 11A, the average error is 147% with a standard deviation of 170%. With the square-wave input, Figure 11B, the inferred erosion rate is off on average by 204%, with a standard deviation of 249%. With the exponential scenario, Figure 11C, the average error for the single cycle is 86%, with a standard deviation of 104%. Perhaps most surprisingly, the average error for the high-frequency variation of the sawtooth shown in Figure 11D is 36%, with a standard deviation of 91%. Running similar scenarios for a wide range of frequencies and amplitudes leads to several pertinent conclusions drawn from the modeling results, which are summarized below.

CONCLUSIONS

This paper reviews and builds on results from two soil-mantled landscapes that are known to push the steady-state assumption in cosmogenic nuclide-based studies of erosion and soil production rates and processes. For different reasons, the Oregon Coast Range and the southeastern Australian highlands suggest the likelihood of variable erosion rates and processes. In Oregon, stochastic landsliding and periodic forest fires can significantly impact both the processes and the rates at which any given part of the landscape is eroding. On the highlands of Australia, dramatic variations in climate during the Pleistocene led to periods of increased erosion that may have left a legacy in the cosmogenic nuclide concentrations measured today. These examples of potentially variable erosion rates are not unique. Indeed, it seems likely that most landscapes are subject to similar exceptions to the steady-state assumption. With such a view in mind, this paper also reviews the details of the CRN approach to derive a model that can numerically simulate the accumulation of CRNs under non-steady-state conditions at a specific sample location. The model presented here enables testing the effects of these potential variations in erosion rate on the nuclide concentrations of samples and, therefore, on the inferred erosion rates obtained from the samples.

Results from the modeling exercises offer some valuable conclusions. First, the greater the maximum erosion rate, the more responsive the nuclide concentration is to the erosional variations. This means that for point-specific sampling, the exposure history of the sample is especially important to constrain. For the Oregon Coast Range, there is potential for each of the individual soil production rates to be off by as much as a factor of six, if, for example, a given sampling location had experienced a recent landslide that removed several tens of centimeters of saprolite and then was rapidly filled in with upslope-derived colluvium such that there was no surficial evidence of the event. Confidence in the apparent soil production function for the region is only gained by the agreement between several samples, in this case 21 from the weathered saprolite, that suggest a similar trend. Close examination of the data shows, however, that variations in soil production rate under similar soil thicknesses can be as great as a factor of three. In a landscape that is likely to be experiencing stochastic large-magnitude events, it is critical, therefore, to analyze enough samples to ensure that any local variations in rate do not obscure the long-term average rate. The southeastern Australian highlands example offers an example of how the initial soil production function, inferred from five point samples, was statistically different from the refined function, determined by twelve samples. Determining how many samples is enough to define a function such as this is difficult to predict, but modeling the potential local variations in rate can define a range of uncertainty to be expected.

A second conclusion, also relevant to an Oregon Coast Range–like landscape, is that the higher the frequency of the erosional variation, the more likely that the inferred erosion rate is closer to the long-term average. If, for example, every thousand to two thousand years, all the trees at a site burn thoroughly and erosion rates increase by a factor of six or more, it is likely that by the time the landscape has recovered (i.e., enough time has passed that the landscape appears to be in steady state) point-specific erosion or soil production rates will capture the average. This conclusion has important implications for using catchment-averaged samples to infer rates for an entire watershed, or drainage basin. Namely, it supports the methodology and does not draw a distinction between small and large catchments. An important next step is therefore to combine point-specific modeling with a catchment-averaged model and compare model outputs with actual measurements of both point-specific and average rates.

Reducing the magnitude of the variable erosion rate led to an even better agreement between the long-term average erosion rates, irrespective of input function. Again, the higher the frequency, the lower the impact of the input variable erosion functions on the output CRN concentrations and therefore on the inferred erosion rate, supporting the use of catchment-averaged samples. Specifically, this modeling result suggests that if a catchment is experiencing highly variable rates of erosion, due to dominance of different processes across the catchment, then using an integration of sediment derived from all points is more likely to capture the long-term average than any point sample is. Under the lower erosion rate conditions, however, the impact of the variable rate on the CRN concentration is not as great, such that inferred erosion rates are often off by a factor of two or more. Results from the examples shown here to explicitly address uncertainties raised from the field examples are robust and follow similar trends for either increasing or decreasing the magnitude and frequency of the input variation. Perhaps the saving grace of this is that more slowly eroding landscapes will preserve signatures of the anomalous events such that critical field sampling can avoid them. The amending of the southeastern Australia soil production function with the addition of more samples serves as a good example of how enough CRN measurements can help insure capturing the true form of the soil production function. An obvious implication of this is determining the cost-to-benefit ratio of a given CRN sampling scheme. It appears that using catchment-averaged samples to infer erosion rates is an accurate way, in both theory and practice, to quantify erosion rates. These average erosion rates continue to show their value across a broad range of geomorphic applications and will undoubtedly be used to further couple erosion with climatic, tectonic, and anthropogenic forcings. To more fully evaluate how catchment-averaged samples might reflect variable erosion rates, the model presented here can be applied to every point in a catchment and coupled with a sediment routing model. Despite the value and apparent accuracy of average rates determined from CRNs, determining the average rate of a catchment does not untangle the way that different processes interact on the slopes and in the channel of the catchment, such that there is still value in well-planned point sampling strategies. Quantifying the sediment production and transport processes that are contributing to these average erosion rates requires higher numbers of samples and will be more dramatically affected by potential local variations in rate. Specifically, for example, quantifying landslide and bedrock failure processes continues to be elusive and will require combining a point-specific with a catchment-averaged sampling scheme.

Finally, perhaps not surprisingly, geomorphology matters. That is, constraining the sampling site such that there is little likelihood of a high-magnitude/low-frequency event affecting the nuclide concentration is critical. While this is an obvious conclusion, and one that many have articulated previously, the wider and wider application of CRN studies to quantifying surface processes warrants further emphasis of this point. This is particularly relevant to geomorphic settings that are not as well constrained as a soil-mantled, convex-up hillslope or a landslide-free first- or second-order catchment. One potential approach is to collect large numbers of samples, but that quickly becomes cost and time prohibitive. Another approach is to tackle the more complicated landscapes by examining the catchment-averaged rates, but this approach leads to empiricism rather than determinism. For example, catchments that are not uniformly soil-mantled are eroding by a range of processes. Determining the catchment-averaged rates for such diverse drainages have provided empirical support for a wide number of recent studies, but have not been able to distinguish the process-specific contributions of the sediment load. Similarly, the range of timescales over which the catchment-averaged rates are applicable has yet to be shown definitively. It is likely that here the application of multiple analytical tools—combining CRN studies with short-lived isotope studies and low-temperature thermochronometry, for example—will make the most headway. Perhaps the most difficult geomorphology to constrain is the dating of features, whether by burial-age dating in caves or thick deposits, or by profile dating on geomorphic surfaces, simply because the landscapes contributing the sediments being dated are gone or are changed significantly. Irrespective of all the uncertainties and potential pitfalls, it remains amazing how well the application of CRNs to untangling the way Earth's surface is eroding is doing.

Special thanks to Belinda Barnes for initial help with the mathematical formulation of the model and Hany Farid for invaluable help writing the Matlab® code. Kabir Heimsath, Josh Roering, Damien Kelleher, and Geoff Hunt assisted with the fieldwork across the two main study areas. Follow-up work in southeastern Australia has benefited greatly from collaboration with Ron Amundson. Ideas explored and results summarized here reflect work done in close collaboration with Bill Dietrich, John Chappell, Kuni Nishiizumi, and Bob Finkel. Thanks to Lionel Siame, Didier Bourlès, and Erik Brown, the GSA Volume Editors for their efforts to compile this volume, and to Milan Pavich and an anonymous reviewer for their valuable comments on a previous version of this manuscript. All cosmogenic results reported here were partially performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract W-7405-Eng-48. Funding for this work was from National Science Foundation grant EAR-0239655.

Ahnert

,

F.

,

1967

, The role of the equilibrium concept in the interpretation of landforms of fluvial erosion and deposition: in Macar, P., ed., L'evolution des versants: Liege, University of Liege, p.