Eroding the land: Steadystate and stochastic rates and processes through a cosmogenic lens

Published:January 01, 2006
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 CRNbased methodologies are pushed to new limits. Inherent to all is an assumption of steadystate rates and processes. This paper focuses on the use of cosmogenic ^{10}Be and ^{26}Al, 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 nonsteadystate processes in shaping soilmantled landscapes. With this potential in mind, a numerical model is presented, following a review of the CRN conceptual framework, to test the effects of nonsteadystate erosion rates and processes on CRN concentrations. Results from this model focus on ^{10}Be 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 longterm erosion rate. Modeling results are applicable for pointspecific erosion rates, but endorse the catchmentaveraged 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 nonsteadystate processes. The model presented here is ready for application to catchmentaveraged processes, as well as modeling the role of muons in variable erosion rate scenarios.
INTRODUCTION
In situ–produced cosmogenic radionuclides (CRNs) are used extensively to infer erosion rates acting upon a wide variety of landforms. Original application of this methodology focused on recently glaciated bedrock where the large amount of material removed during the Last Glacial Maximum provided a wellconstrained case to determine nuclide production rates (Nishiizumi et al., 1986, 1989). The methodology is well developed (Lal, 1988, 1991) and well reviewed for geomorphic applications (Bierman, 1994; Bierman and Nichols, 2004; Cerling and Craig, 1994; Cockburn and Summerfield, 2004; Gosse and Phillips, 2001; Nishiizumi et al., 1993). Subsequent applications built upon the success of determining exposure ages of relatively unweathered rock to infer longterm, steadystate erosion rates from exposed bedrock surfaces (Bierman and Caffee, 2001; Bierman et al., 1999; Bierman and Turner, 1995; Small et al., 1997). More complicated use of CRN concentrations toward geomorphic applications involved determining river incision rates from exposure ages of strath terraces (Burbank et al., 1996; Pratt et al., 2002; Reusser et al., 2004; Weissel and Seidl, 1998), as well as from exposure ages inferred from concentrations in fluvial cobbles deposited on terraces (Anderson et al., 1996; Hancock et al., 1999; Perg et al., 2001; Repka et al., 1997). In situ concentrations of CRNs from bedrock beneath different depths of soil also led to the first determination of the rates of soil production (Heimsath et al., 1997) with subsequent and similar findings from the base of a retreating escarpment (Heimsath et al., 2000). A powerful methodology for determining average erosion rates for small watersheds was also developed, which relies on measuring CRN concentrations from stream sediments (Bierman and Steig, 1996; Brown et al., 1995b; Granger et al., 1996; Riebe et al., 2000). Each of these applications was relatively well constrained such that the exposure history of the sampled surfaces or small catchment areas could be reasonably inferred. Recent studies are pushing the limits of these constraints and are using burial age dating (Granger et al., 2001; Granger and Smith, 2000; Granger and Stock, 2004; Stock et al., 2005) as well as further applications of depthprofile dating (Perg et al., 2001; Schaller et al., 2004) to infer landscapescale uplift and erosion rates. Similarly, application of the basinaveraged methodology is being applied across larger and larger catchments (Cockburn and Summerfield, 2004; Matmon et al., 2003a; Schaller et al., 2001), in more and more complicated erosional settings (Hewawasam et al., 2003; Wobus et al., 2005), and to develop landscapescale sediment budgets (Clapp et al., 2000, 2002; Nichols et al., 2005). All in all, widespread use of CRNs enables a phenomenal amount of work on quantifying geomorphic rates and processes (Bierman and Nichols, 2004; Granger, 2002). Continued application of the CRN method of determining erosion rates and landscape evolution histories to more and more complicated geomorphic problems will depend on being able to more narrowly constrain the exposure history of the samples, a point emphasized by many of the above studies, but only tested by a few (Bierman and Steig, 1996; Lal, 1991; Small et al., 1997).
Here I briefly review the state of using CRN measurements (focusing specifically on in situ–produced ^{10}Be and ^{26}Al) 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, soilmantled 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, ^{10}Be and ^{26}Al, 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 steadystate conditions for both the overlying soil mantle and the erosion rate (equivalent to the soil production rate) of the underlying bedrock.
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 nonsteadystate 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 ^{10}Be and ^{26}Al 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 pointspecific 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 soilmantled 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 bedrocksoil 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 soilbedrock 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 soilbedrock boundary, not in the soil. The sloping soil mantle is treated like an additional buffer against cosmicray 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 ^{10}Be and ^{26}Al 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 steadystate 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 wellconstrained 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 steadystate 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, P_{S}, 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 wellmixed soils of the landscapes studied here. It is assumed that the nuclides are produced only by cosmicray 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 nonsteadystate 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, z_{0}, 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 N_{0} as a constant, Solving for N(t) thus results in an expression that only requires solving a simple integral, Under the steadystate 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 steadystate 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 z_{0} will have risen to a depth z =z_{0}– ϵT. Substituting z_{0} = Z + ϵT into equation 10 leads to The integral can now be solved analytically as follows: Equation 12 provides a closedform solution for concentration under the steadystate 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 steadystate 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 graingrain 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 pointspecific samples to determine basinwide 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 lowmagnitude/highfrequency events (e.g., processes approaching steady state) and the lowfrequency/highmagnitude 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 nonsteadystate 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 steadystate condition, the above integral does not have a closedform 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 nonsteadystate 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) = N_{0} (the initial concentration, typically equal to zero), and the depth, z(t), at time is t_{0} 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 pointspecific concentrations and inferred erosion rates. The obvious next step is to incorporate this pointspecific 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 ^{10}Be and ^{26}Al concentrations from exposed and soilmantled, weathered bedrock as well as stream sediments and bedrock. Both field sites differ in critical ways from those examined by similar studies across completely soilmantled 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 wellstudied region of the Pacific Northwest that exemplifies humidtemperate, soilmantled landscapes where stochastic as well as steadystate 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 ridgeandvalley topography (Fig. 1), and while there is some hint of longterm 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 soilmantled 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 debrisflow or landslide scars in the hollows. Creeks and rivers draining the region are mixed alluvialbedrock and are thought to be incising at roughly the same rate as the longterm 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 ^{26}Al and ^{10}Be) 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 catchmentaveraged 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 longterm 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, longterm erosion rates from cosmogenic nuclides are significantly lower than shortterm 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.
Southeastern Australian Highlands
Rolling soilmantled 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 threequarters 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 ^{26}Al and ^{10}Be) 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 soilmantled 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 soilmantled 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 ^{26}Al and ^{10}Be 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 ^{26}Al and ^{10}Be) from the fluvially sculpted and polished bedrock of the active channel bed and averaged 9 m/m.y., suggesting that the soilmantled 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.
What makes this highland site so interesting from a nonsteadystate erosional perspective are the results from a profile of bedrock samples collected from the side of an outcropping granite tor. Our interest in longterm 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 soilmantled 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 soilmantled landscape and the eroding bedrock that is outcropping. Steadystate 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.
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 nonsteadystate 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 steadystate nuclide concentration, or, in the case of an oscillating concentration, a steadystate oscillation with the imposed variable erosion rate. Conceptually, for steadystate erosion scenarios, the faster the erosion rate, the shorter the time to steadystate 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 steadystate erosion rates, functions that result in steadily decreasing steadystate 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.
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 longterm 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 longterm 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 observationbased. 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 fullwave 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 steadystate 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, ^{10}Be, 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 ^{10}Be 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 fullwave 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 squarewave 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.
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 fullwave 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 basinaveraged 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 squarewave 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 fullwave scenario, Figure 10A, the inferred erosion rate is off by an average of 49%, with a standard deviation of 83%. The squarewave 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%.
While this result does not bode well for pointspecific sample collection, it is worth thinking about how such erosion rate variations manifest themselves in the longterm 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 squarewave scenario of Figure 10B, the longterm averages are 545 and 586 m/m.y., respectively for input and inferred rates, with the inferred overestimating the input by 7.6%. The longterm average for the input exponential variation, Figure 10C, is 242 m/m.y., while the longterm average for the inferred rate is 249 m/m.y., meaning an overestimation of only 5.2%. Despite the factor of fourplus overestimation of the low erosion rates due to the sawtooth variation of Figure 10D, the difference between the longterm 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 longterm 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 pointspecific 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 steadystate concentrations under the lower input erosion rates (see Fig. 7).
Figure 11 plots different frequencies for the input varying erosion rates, with Figure 11A showing the same frequency for the fullwave 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 highfrequency change of the sawtooth input has very little effect on the inferred erosion rate at any time, meaning that the steadystate 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 catchmentaveraged sampling methodology, similar to the above scenario and as discussed below. The variation in rates shown by the square wave, Figure 11B, result in longterm averages that are also identical, with both the input and the inferred rates averaging 55 m/m.y. The longterm average inferred rate of 40 m/ m.y. is only 5.2% higher than the input average of 38 m/m.y. for the fullwave rectification of Figure 11A, while the longterm average for the lowerfrequency 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 fullwave variation of Figure 11A, the average error is 147% with a standard deviation of 170%. With the squarewave 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 highfrequency 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 soilmantled landscapes that are known to push the steadystate assumption in cosmogenic nuclidebased 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 steadystate 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 nonsteadystate 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 pointspecific 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 upslopederived 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 largemagnitude events, it is critical, therefore, to analyze enough samples to ensure that any local variations in rate do not obscure the longterm 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 longterm 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) pointspecific erosion or soil production rates will capture the average. This conclusion has important implications for using catchmentaveraged 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 pointspecific modeling with a catchmentaveraged model and compare model outputs with actual measurements of both pointspecific and average rates.
Reducing the magnitude of the variable erosion rate led to an even better agreement between the longterm 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 catchmentaveraged 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 longterm 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 costtobenefit ratio of a given CRN sampling scheme. It appears that using catchmentaveraged 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 catchmentaveraged 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 wellplanned 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 pointspecific with a catchmentaveraged sampling scheme.
Finally, perhaps not surprisingly, geomorphology matters. That is, constraining the sampling site such that there is little likelihood of a highmagnitude/lowfrequency 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 soilmantled, convexup hillslope or a landslidefree first or secondorder 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 catchmentaveraged rates, but this approach leads to empiricism rather than determinism. For example, catchments that are not uniformly soilmantled are eroding by a range of processes. Determining the catchmentaveraged rates for such diverse drainages have provided empirical support for a wide number of recent studies, but have not been able to distinguish the processspecific contributions of the sediment load. Similarly, the range of timescales over which the catchmentaveraged rates are applicable has yet to be shown definitively. It is likely that here the application of multiple analytical tools—combining CRN studies with shortlived isotope studies and lowtemperature thermochronometry, for example—will make the most headway. Perhaps the most difficult geomorphology to constrain is the dating of features, whether by burialage 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.
Contact A.M.H. for free Matlab^{®} code and instructions for model presented herein.
Contact A.M.H. for free Matlab^{®} code and instructions for model presented herein.
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. Followup 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 Bourlè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 W7405Eng48. Funding for this work was from National Science Foundation grant EAR0239655.
Figures & Tables
Contents
In SituProduced Cosmogenic Nuclides and Quantification of Geological Processes
GeoRef
 age
 Al26
 alkaline earth metals
 aluminum
 Australasia
 Australia
 Be10
 beryllium
 Coast Ranges
 cosmogenic elements
 cyclic processes
 erosion
 erosion features
 erosion rates
 exposure age
 isotopes
 landform evolution
 landscapes
 metals
 New South Wales Australia
 numerical models
 Oregon
 quantitative analysis
 quantitative geomorphology
 radioactive isotopes
 soils
 stochastic processes
 United States
 Australian Alps