Combined linearregression and Monte Carlo approach to modeling exposure age depth profiles
 ^{1}Department of Earth and Planetary Sciences, UC Davis, Davis, 95616, USA
 ^{2}Earth Observatory of Singapore, Nanyang Technological University, 639798, Singapore
 ^{1}Department of Earth and Planetary Sciences, UC Davis, Davis, 95616, USA
 ^{2}Earth Observatory of Singapore, Nanyang Technological University, 639798, Singapore
Correspondence: Yiran Wang (yrwwang@ucdavis.edu)
Hide author detailsCorrespondence: Yiran Wang (yrwwang@ucdavis.edu)
We introduce a set of methods for analyzing cosmogenicnuclide depth profiles that formally integrates denudation and muogenic production, while retaining the advantages of linear inversion for surfaces with inheritance and age much greater than zero. For surfaces with denudation, we present solutions for both denudation rate and total denudation depth, each with their own advantages. By combining linear inversion with Monte Carlo simulation of error propagation, our method jointly assesses uncertainty arising from measurement error and denudation constraints. Using simulated depth profiles and naturalexample depth profile data sets from the Beida River, northwest China, and Lees Ferry, Arizona, we show that our methods robustly produce accurate age and inheritance estimations for surfaces under varying circumstances. For surfaces with very low inheritance or age, it is important to apply a constrained inversion to obtain the correct result distributions. The denudationdepth approach can theoretically produce reasonably accurate age estimates even when total denudation reaches 5 times the nucleon attenuation length. The denudationrate approach, on the other hand, has the advantage of allowing direct exploration of tradeoffs between exposure age and denudation rate. Out of all the factors, lack of precise constraints for denudation rate or depth tends to be the largest contributor of age uncertainty, while negligible error results from our approximation of muogenic production using the denudationdepth approach.
 Article
(3987 KB) 
Supplement
(1739 KB)  BibTeX
 EndNote
In situ terrestrial cosmogenic nuclide (CN) dating, especially with ^{10}Be, is a widely applied tool to estimate landform ages (e.g., Granger et al., 2013). These dates are affected by landscape processes that either remove or add CNs, lending uncertainty that may be difficult to assess without additional information. Ages of landforms constructed from sediments, such as a stream terrace, may be affected by CNs acquired by the sediments prior to deposition, termed inheritance, leading to erroneously older dates (Brocard et al., 2003; Hancock et al., 1999; Repka et al., 1997). Conversely, even a low rate of denudation of a landform after its formation will bias surfaceexposure ages younger (Lal, 1991). Under the condition of no denudation, solving a depth profile of CN concentrations via linear inversion – a technique first developed by Anderson et al. (1996) – provides a robust approach for estimating surface age and inheritance from a landform comprised of sediments. However, this method suffers from the deficit of not incorporating denudation or muogenic production process and has been succeeded by forwardmodeling approaches (e.g., Hidy et al., 2010). Despite its deficits, however, linear regression retains advantages as a robust and straightforward approach to invert for exposure age. Herein we revisit the application of linear regression for CN depth profiles with an updated approach, building upon the simplicity of the original technique, but expanding its application to surfaces with independently constrained denudation histories, and increasing the accuracy of the age and inheritance results by taking muogenic production into account.
As an inverse approach, our leastsquares linear regression directly solves for a bestfit age and inheritance, while treating the denudation rate or denudation depth as a model input, rather than an output of the model. It may be used with Monte Carlo sampling to explore the full distribution of possible ages and inheritance from the variation of input parameters. As an inverse approach, linear regression directly returns a best estimate, providing a clear way to explore how model inputs, such as erosion, affect the result. This is useful to derive an exposure age directly, or as a starting point for forward models (e.g., Hidy et al., 2010; Laloy et al., 2017; Marrero et al., 2016). In the latter case, linear inversion may help researchers tune the boundary values for the forward models to get better simulation results.
In this paper, we present a general inversion which incorporates muogenic production as a summation of exponential terms (e.g., Heisinger et al., 2002a, b), and two specific approaches for age inversion at sites with constant denudation rate. The first of these approaches applies a constant denudation rate but requires simplification by omitting muogenic production. This approach provides a useful tool to explore the tradeoffs between denudation rate and exposure age but introduces systematic errors as denudation rate increases due to the exclusion of muogenic production. The second approach introduces a solution for a constant denudation depth with a Taylorseries approximation for the muogenic production terms. With this approximation, linear regression produces robust age estimates even for surfaces with a large amount of denudation. To show the application of these techniques, we present applications to pseudorealistic depth profiles under various scenarios. We also apply our approach to two previously published example sample sites, one from the Beida River in western China from our own work (Wang et al., 2020) and one from the Lees Ferry site on the Colorado River from Hidy et al. (2010), to demonstrate model performance for realistic cases. These examples show that our revised linearregression approach is robust and can be applied to most exposure dating scenarios.
2.1 General inversion
Under conditions of constant production rate and constant denudation rate, a surface that was exposed at time t would have a concentration of a cosmogenic nuclide (N_{z}) as follows (Balco et al., 2008; Braucher et al., 2009; Lal, 1991; Lal and Arnold, 1985):
where P_{n,0}, ${P}_{{m}_{\mathrm{1}},\mathrm{0}}$, and ${P}_{{m}_{\mathrm{2}},\mathrm{0}}$ are the surface production rate induced by nucleons, negative muons, and fast muons, and Λ_{n}, ${\mathrm{\Lambda}}_{{m}_{\mathrm{1}}}$, and ${\mathrm{\Lambda}}_{{m}_{\mathrm{2}}}$ are the attenuation scale lengths (g cm^{−2}) of the nucleons and muons (negative and fast), respectively. Note that the two exponential terms for muogenic production used here are an approximation of a complex muogenic production path, because the muon energy spectrum hardens continuously with depth (Heisinger et al., 2002a, b; Marrero et al., 2016; Balco, 2017). z is the depth beneath the target surface; λ is the decay constant, and ε is a constant denudation rate, if applicable. For our purposes, we model ages using ^{10}Be, with a halflife of 1.39 Myr (Chmeleff et al., 2010; Korschinek et al., 2010), due to its wide applicability to quartzbearing sediments (Cockburn and Summerfield, 2004; Granger et al., 2013; Rixhon et al., 2017).
Based on Eq. (1), the production of cosmogenic nuclides may be simplified into two major components: the production rate at specific depth (P_{z}) and the effective exposure age of the site (T_{e}), which is the time that is required to accumulate concentration N_{z} at production rate P_{z} without erosion and radioactive decay. Therefore Eq. (1) may be rearranged into
where
The ^{10}Be concentration measured from a suite of samples (Fig. 1), C, has two components: the in situ produced concentration, N_{z}, and the inherited concentration, C_{inh},
However, ^{10}Be concentration (C) is exponential to the burial depth, based on Eqs. (2) and (3), when there is no denudation ($\mathit{\epsilon}\approx \mathrm{0}){T}_{en}={T}_{e{m}_{\mathrm{1}}}={T}_{e{m}_{\mathrm{2}}}$, and therefore Eq. (3) can be rearranged as follows:
where
This equation is an update to the linearregression approach first proposed by Anderson et al. (1996) that accounts for both nucleon and muon production, as well as radioactive decay. For the case of no denudation, CN concentration is linear to the sum of production rates via all pathways $\left(\sum _{i}{P}_{zi}\right)$, and T_{e} and C_{inh} are the slope and intercept of this linear relationship, respectively. Therefore, similar to the approach proposed by Anderson et al. (1996), we can apply leastsquares linear regression to find the slope (T_{e}) and intercept (C_{inh}) of the bestfit line to the concentration vs. production rate data of the depth profile. The exposure age, factoring in decay, may be calculated directly by rearranging Eq. (4b):
2.2 Inversion with denudation rate
For sites with constant denudation rate, ε, the effective age for each pathway (nucleons or muons) would be different, due to their different attenuation lengths. But an approximation may be made by omitting the muogenic production, on the basis that muogenic production only makes up a very small fraction of the total surface production (Braucher et al., 2003, 2011, 2013; Heisinger et al., 2002b, a; Balco et al., 2008; Balco, 2017), and Eq. (3) may be further simplified to
Using Eq. (6), a linear leastsquares regression can be applied to find the bestfit T_{en} and C_{inh}, which leads to the estimated exposure age
where
This solution illustrates the utility of separating the age model for finding T_{en} from the effect of denudation rate, contained within the parameter B. Though the error introduced to omitting muogenic production grows as the denudation rate increases, this approach (Eq. 7) provides a useful tool to examine the relationship between exposure age and denudation rate.
2.3 Inversion with denudation depth
For many practical cases, it may be more straightforward to estimate total denudation depth (D) from field evidence such as through soilprofile analysis, rather than a denudation rate (e.g., Ebert et al., 2012; Hidy et al., 2010; RuszkiczayRüdiger et al., 2016). With denudation depth, the effective age of each pathway may be rewritten as
Here we explore the application of this equation with the inclusion of muogenic production. Using a series expansion, we rewrite the effective age related to muons (negative and fast), T_{em}, into a fraction, g, of the effective age related to nucleons, T_{en}. The fraction g can be approximated solely from knowledge of the denudation depth, D (see Appendix for derivation):
Bringing g_{i} into Eq. (3), we have
where P_{ze} is the effective production rate from both nucleons and muons under the condition of a finite amount of steady denudation over the lifetime of the deposit.
Using Eq. (10), T_{en} and C_{inh} can be found by applying leastsquares linear regression with known production rates, denudation, and sample concentrations, similar to the general inversion case for no denudation described by Eq. (4).
To estimate the exposure age, we need to find the solution for
While the complicated form of Eq. (11) prohibits a direct solution, t may be found iteratively by applying the Newton's method. Using the derivative of Eq. (11),
the exposure age can then be iterated from
with initial guess, t_{0}=T_{en}.
We present a set of example applications of linear regression to CN depth profiles using both the denudationrate (Eqs. 6, 7) and the denudationdepth (Eqs. 8–13) techniques. We begin with simulated CN depth profiles to explore the impacts of scatter of sample concentration, low inheritance, denudation, and sample depth on the accuracy and precision of our approach, and to compare the performance of linear regression with the estimations using a Monte Carlo Markov chain approach. We then demonstrate the linearregression technique with two case examples. For the Beida River T2 terrace of the North Qilian Shan, China (Wang et al., 2020), we demonstrate using the denudationrate approach to explore the possible range of exposure age and the tradeoffs between age and denudation rate; we also demonstrate how to use the denudationdepth approach combined with Monte Carlo simulation to estimate the full distribution of exposure age and inheritance. For the Lees Ferry site on the Colorado River (Hidy et al., 2010), we compare results from our denudationdepth approach with the original publication.
To explore the full distribution of possible estimated age and inheritance, we consider the uncertainty propagated from five different sources: uncertainty of the ^{10}Be concentration measurements, uncertainty of sample depths, uncertainty of sediment density, uncertainty of the denudation depth or rate, and the uncertainties related to CN production and decay (i.e., the attenuation lengths, production rates, etc.). These uncertainties propagate sequentially, first from ^{10}Be concentration, sample depths, and sediment density through the leastsquares linearregression process, and second from denudation rate or depth through converting exposure age from the effective exposure age (T_{e}). Because of the limited sample sizes typical of most studies, and the variance in both concentration and depth, we apply a Monte Carlo sampling approach for the range of each input parameter. The code used here is archived in GitHub (https://github.com/YiranWangYR/10BeLeastSquares, last access: 15 August 2022; Wang, 2022). Our procedure includes the following steps.

Step 1: generate all the input parameters for one model – sampling distributions for ^{10}Be concentration, sample depth, density, denudation rate or depth, production rate, attenuation length, etc. Depending on the parameter, we represent uncertainty as either a normal or uniform distribution.

Step 2: fit these sampled input parameters with Eq. (4) or Eq. (6) or Eq. (10) to derive bestfit sample T_{en} and inheritance values.

Step 3: calculate the exposure age use T_{en} from step 2 and parameters generated from step 1 using Eq. (5) or Eq. (7) or Eqs. (11)–(13).

Repeat steps 1–3 for many times to produce a distribution of age and inheritance results.
3.1 Simulated CN depth profiles
To demonstrate how our approach will perform under different circumstances, we generate a series of simulated depth profiles to test how closely the estimated age and inheritance will reflect the true values. With these profiles, four different scenarios are tested: varying degrees of random deviation of sample concentrations, varying amount of denudation, low inheritance, and deep sample depth. For comparison, we use a forward model which operates with a Monte Carlo Markov chain algorithm (codes: https://github.com/YiranWangYR/10BeLeastSquares/blob/main/Be10_LRvsBayesian.m, last access: 15 August 2022; Wang, 2022) to estimate the exposure age and inheritance along with our linearregression approach. We choose to use our own coding following Bayesian inversion principles instead of published codes because, first, we need to ensure the input variables are the same for both approaches, and second, to be statistically significant, we need to apply both approaches several hundred times for each scenario; therefore manually doing so with published codes would be extremely time consuming.
^{a} Production rate and attenuation length are calculated (nucleon) and approximated (muon) based on a hypothetical middlelatitude, 1000 m elevation site (Stone, 2000; Balco, 2017). ^{b} The analytical uncertainty of sample concentration. ^{c} The amount that the mean sample concentration deviates from the true concentration.
For all simulated profiles, we set the true age and inheritance to be 200 kyr and 100 × 10^{3} atoms g^{−1} except for the lowinheritance scenario, where we set inheritance to 0 and 5000 atoms g^{−1}. Each simulated profile contains six samples; to mimic realistic scatter of sample concentrations, all the simulated samples are deviated from the true profile based on an assigned normal distribution of either 1 %, 2 %, 5 %, or 10 % standard deviation. Analytical uncertainties (2 % standard deviation) are assigned to each sample of the simulated profile to mimic realistic values. For each scenario, 500 simulated profiles are generated and inverted for age and inheritance using both the linearregression and forward approaches. The details of our simulation process and the parameters we used for each scenario are listed in Table 1.
3.1.1 Deviation of sample concentrations
We test three sets of simulations under different degrees of imposed deviation of sample concentrations (standard deviation of 2 %, 5 %, and 10 %) with no denudation. The 2 % deviation case matches the analytical uncertainty (2 %), while the 5 % and 10 % cases introduce excess scatter, as is typically found in CN depth profiles. For the case with 2 % of deviation, we find both linearregression and forward approaches yield results centered around the true age (200 kyr). For the linearregression approach, 95 % of the mean ages fall within 4.0 %–4.5 % error range of the true age, while for the forward approach 95 % of the mean age fall within 3.0 %–3.5 % error range of the true age. The forward approach works slightly better than linear regression for the more extreme (noisier) cases (Fig. 2 and Fig. S1 in the Supplement). For inheritance estimation, similar to the exposure age, both approaches return results centered around the true value (100 000 atoms g^{−1}), and a forward approach again performs somewhat better for noisier data. When the sample deviation increases to 5 %, then to 10 %, the estimated age and inheritance values remain centered on the correct age and inheritance while the ranges of the mean estimations for both approaches expand significantly, indicating a decrease in precision as the noise level increases (Figs. 2, S1–S3). With 5 % of deviation, 95 % of the estimated ages fall within 10 % and 8 % ranges of the true value for linearregression and forward approach, respectively. With 10 % of deviation, 95 % of the estimated age fall within 21 % and 17 % ranges of the true value for linearregression and forward approach, respectively. We also find the error range for the estimation results from the linear regression is moderately larger than a forward approach (Figs. S1–S3). Importantly, we find that while each approach provides different estimation results for the same set of samples, one model does not consistently perform better than the other (Figs. S1–S3). About 40 %–45 % of the time the linear regression returns results closer to the true age than the forward approach.
3.1.2 Low inheritance or exposure age
One important constraint for our inversion (Eqs. 6 and 10) is that both inheritance (C_{inh}) and the effective exposure age (T_{en}) should not be negative. This constraint is not important when both C_{inh} and T_{en} are much greater than zero, as none of the inverted results will fall below zero. However, when either one of these variables is small enough, inverting without a positivity constraint may lead to incorrect uncertainty distributions. We demonstrate this by using profiles generated with 5 % imposed deviation of sample concentrations and zero inheritance (Table 1). Here we compare three different scenarios: (1) not permitting negative inheritance during inversion through setting the boundary condition, C_{inh}≥0; (2) permitting negative inheritance during inversion; and (3) excluding negative inheritance outcomes after inversion. Notice that scenario 1 is a constrained inversion (nonnegative least square problem; solution in Lawson and Hanson, 1987) while scenarios 2 and 3 are not constrained. The distributions of the mean exposure age and inheritance (Fig. 3) show that, in the first scenario, the positivity constraint forces all the possible negative inheritance to be exactly zero, which leads to an age distribution skewed to the left while still centering the true age (Fig. 3a, d). In the unconstrained scenarios (2), the uncertainty distribution for inheritance incorrectly extends to negative values, though the age results were distributed normally around the true age (Fig. 3b, e). When negative inheritance results are excluded after applying the unconstrained inversion (scenario 3), the inheritance results exhibit a clipped normal distribution, while the age is not only leftskewed but also shows obvious underestimation (Fig. 3c, f), resulting from discarding the older age estimations that were related to negative inheritance. Similarly, when inverting for a surface that is sufficiently young, or the inheritance is large relative to the surface age, the unconstrained regression could also lead to an age distribution extending to negative values.
3.1.3 Denudation depth
To test the robustness of our denudationdepth approach, we tested simulated CN depth profiles with total denudation of 1, 2, 3, and 5 times the nucleonspallation attenuation length $\left(\frac{\mathit{\rho}}{{\mathrm{\Lambda}}_{n}}\right)$ and 5 % of imposed sample deviation (Table 1). The resulting distributions of the mean age (Fig. 4) are all well centered around the true age even with the largest amount of denudation tested. We also observe that when total denudation reaches 3 and 5 times the attenuation length, the ranges of the estimated mean age grow slightly wider, indicating a decrease in precision. This phenomenon can also be observed with the forward approach (Fig. S4), which suggests it is not a result of the approximation introduced through Eq. (9). Instead, we suggest this problem relates to a low concentration gradient for profiles with a large amount of denudation, making the inversion more sensitive to measurement error. We find that the distribution of the mean inheritance remains well centered around the true inheritance for all simulated profiles, even for profiles with denudation equal to 5 times the attenuation length.
3.1.4 Deep sample profiles
Samples at depths greater than ∼ 2 m are especially sensitive to muogenic production (Fig. 1). Here we test our denudationdepth approximation with depth profile samples distributed between 3 and 5 m depth to mimic the situation when near surface samples are not obtainable. Three groups of profiles, subjected to a total denudation of 0, 2, and 5 times the spallation attenuation length and with 5 % of imposed sample deviation, are tested with both linearregression and forward approach (Table 1).
We find that compared to the nearsurface profiles, results from the deep profiles show greatly reduced precision, especially for large denudation depths (Figs. 5, S5). The majority (95 % confidence) of the mean exposure ages estimated with linear regression spread between 100–300, 20–400, and 0–570 kyr, for surfaces with 0, 2, and 5 times the attenuation length of denudation, respectively (Fig. 5). This occurs with both linearregression and forward approaches, indicating the precision drop is not due to our approximation of muogenic production (Eq. 9). Similar to the examples with large denudation (Fig. 4), we suggest the low precision is a result of the very low concentration gradient at depth. Different from nearsurface profiles, a small change in the concentration gradient at depth leads to a large change of the estimated exposure age, which makes the inversion overly sensitive to random scatter of sample concentrations. As a comparison, by setting the imposed sample deviation to 1 % instead of 5 %, we find that precision increases greatly (Fig. 6). With linear regression, the majority (95 % confidence) of the mean exposure ages are distributed between 180–220, 162–242, and 131–275 kyr for 0, 2, and 5 times the attenuation length denudation, respectively (Figs. 6, S6). Thus, theoretically, both linearregression and forward approach can produce accurate estimates if the scatter of sample concentrations is low. Reducing analytical uncertainty through use of large sample masses would be crucial for the success of a deep sample profile.
3.2 Case examples
Here we present two case examples to show the application of our linearregression model to natural conditions. Wang et al. (2020) excavated a sample pit ∼ 2 m deep beneath the tread of a terrace deposit from the Beida River, within the North Qilian Shan, western China. The site is covered with 125 cm loess which was deposited continuously since 8.3 kyr. Beneath the loess cover, erosional truncation of the A and uppermost B horizon of the soil profile indicates that there may have been 20–60 cm of erosion of the terrace tread prior to the onset of loess accumulation. We collected six samples of medium to coarse sand from up to 2 m below the base of the loess. ^{10}Be concentrations for these samples, corrected for loess accumulation using the approach of Hetzel et al. (2004) are listed in Table S1.
To explore the effect of denudation rate, we use the most likely concentrations for each sample and invert for a preliminary T_{en} and inheritance values of 99.3 kyr and 8.3 × 10^{3} atoms g^{−1} respectively, using Eq. (6). Using Eq. (7), we generate a plot of exposure age vs. denudation rate based on the preliminary T_{en} value. The plot (Fig. 7) shows that if there has been no denudation, the sample site has a minimum exposure age of ∼ 102 kyr prior to loess accumulation. This exposure age increases to 160 and 240 kyr when the denudation equals to 1 and 2 times the attenuation length, respectively. When the denudation rate reaches ∼ 0.7 cm kyr^{−1}, the CN accumulation and denudation reach equilibrium and no age may be determined.
To explore the full distribution of potential age and inheritance results, we apply a Monte Carlo sampling of the sample concentrations, depths, denudation depth, etc. (Table S2). Using a normal distribution for the denudation depth with 40 cm as the mean and 10 cm as the standard deviation, we apply leastsquares linear inversion with Eq. (10). Bestfit results for ^{10}Be concentration (C_{1}) and effective production rate (P_{ze}) are shown in Fig. 8a, and the corresponding fitted depth profile curves are shown in Fig. 8b. With positivity constraints enforced, the inverted T_{en} and inheritance values are 86.3–107.2 kyr and 0–6.4 × 10^{4} atoms g^{−1}, respectively (ranges correspond to the 95 % confidence distributions for each value, Fig. 8c and d). The corresponding exposure age, calculated following Eqs. (11)–(13), is 108.6–150.6 kyr (2σ) prior to loess accumulation (Fig. 8e). The corresponding erosion rate is 0.18–0.42 cm kyr^{−1} (Fig. 8f and g). Without positivity constraints for the inheritance, the resulting exposure age and inheritance would be 108.8–155.0 kyr and −4.7 × 10^{4} to +6.4 × 10^{4} atoms g^{−1} (2σ) instead.
Hidy et al. (2010) reported a ^{10}Be depth profile from the Lees Ferry site, excavated on top of a Colorado River fill terrace. Based on the soil profile, a total erosion of 0–30 cm is estimated for the sample site. For this site Hidy et al. (2010) applied their Monte Carlo model and originally reported an exposure age and inheritance of 69.8–103 kyr, and 6.97–10.70 atoms g^{−1}, respectively (95 % confidence). After updates of their code (v1.2, Hidy et al., 2010; Mercader et al., 2012), including the incorporation of a Bayesian fitting process, their model provides a new age estimation of 76.6–96.1 kyr. See Hidy et al. (2010) for more details of the sample site, sampling and processing, and age result interpretation.
Following the original study, we apply our modeling approaches to the sand depth profile data (Table S1). In order to compare with results reported by Hidy et al. (2010), we use the same values they do for all parameters wherever possible (Table S2). Similar to the Beida River depth profile, we estimate the exposure age with the denudationdepth approach, applying Monte Carlo simulation of the uncertainty of input parameters. With a uniformly distributed 0–30 cm denudation length, the inverted bestfit lines and curves are in Fig. 9a and b. The estimated range of T_{e} and C_{inh} values are 65.2–79.9 kyr and 8.8–12.4 × 10^{4} atoms g^{−1}, respectively (95 % confidence; Fig. 9c and d). The estimated exposure age is 70.4–96.0 kyr (95 % confidence; Fig. 9e).
4.1 Case–example comparison
For the Beida River site, comparing the age determined here (108.6–151.3 kyr) to that we reported previously (107.9–164.5 kyr; Wang et al., 2020), the mean is 5 % younger, while the 95 % error range is 24 % smaller. These differences come from three different sources. First, there is a ∼ 1 % shift that arises from independently sampling depth for each measurement (the depth distributions for each sample were not sampled independently in the original paper). Second, the original paper did not take muogenic production into account. The contribution from muons leads to slightly younger age and lowers the inheritance. Third, we applied the denudationrate approach (Eq. 6) in the 2020 paper, and the corresponding denudation rate distributions are slightly different for the two approaches. Combined, the addition of muogenic contribution and using a denudationdepth instead of denudationrate approach leads to ∼ 4 % shift of the mean age and to 24 % narrowing of the error range.
For the Lees Ferry site, our age estimation of 70.4–96.0 kyr using the denudationdepth approach generally agrees with the result reported by Hidy et al. (2010), 69.8–103 kyr (or 76.6–96.1 kyr based on recalculation with the updated v1.2 code), though small discrepancies in mean age and 95 % confidence range remain between the two approaches. We suggest this arises from two different sources. First, as demonstrated in Sect. 3.1.1 and Figs. S1–S3, most of the discrepancy may be due to differences between the leastsquares and Bayesian model. Secondly, we use a twoterm approximation for muogenic process in our model instead of a 5term approximation used in the original study (Hidy et al., 2010), which may lead to a minor estimation discrepancy.
4.2 Sources of error
4.2.1 Denudation
Denudation and its uncertainty constitute a major source of error in exposure age estimation. With the same surface ^{10}Be concentration, higher denudation rate and/or larger denudation depth would result in an older effective surface age (e.g., Fig. 1). If a surface is sufficiently old, or if the denudation rate is sufficiently high, the CN buildup at the surface will reach equilibrium with nuclides removed through erosion (Lal, 1991). Figure 10a shows the relationship between denudation rate and surface age (Eq. 7). This figure suggests that once the denudation depth exceeds the mean attenuation length of nucleon spallation $\left(\frac{{\mathrm{\Lambda}}_{n}}{\mathit{\rho}}\right)$, the slope of the age versus erosion rate relationship degrades the age determination. Once the denudation depth exceeds twice the attenuation length of spallation, the age versus denudation rate curve flattens so much that it becomes effectively impossible to estimate surface age using a denudationrate approach. On the other hand, the age–denudation depth curve does not flatten as much (Fig. 10b), and therefore it is theoretically possible to use the denudation depth to determine surface age even when total denudation exceeds twice the attenuation length. Our simulated depth profile analysis also shows that the denudationdepth approach can provide accurate estimations even denudation reaches 5 times the spallation attenuation length. In practice, however, surfaces with such a large amount of denudation would subject be to large uncertainties and the denudation history may be too complex for the constant denudationrate assumption, which underlies both approaches, to be valid, casting doubt on the utility of ^{10}Be exposure dating for such cases.
Denudation affects the uncertainty of exposure age estimation in two different ways. First, the uncertainty of the final age increases with denudation rate or depth because of the nonlinear relationship between age and denudation (Fig. 10). Second, the age uncertainty will increase further through propagation of the uncertainty of the denudation. This suggests that when excavating depth profile pits in surfaces subject to denudation, it is crucial to document surface texture and analyze the soil profiles to estimate denudation depth, for a small deviation from the true denudation rate or depth would lead to a large bias in the resulting exposure age (Ebert et al., 2012; Frankel et al., 2007; Hidy et al., 2010; Mercader et al., 2012; RuszkiczayRüdiger et al., 2016).
4.2.2 Muogenic production
Muogenic production affects the accuracy of the estimated surface ages differently for the various approaches considered here. For surfaces with no denudation, muons may be fully incorporated into the inversion (Eq. 4); therefore the uncertainty only comes from the uncertainties of parameters related to muogenic production (attenuation length and production rate). When denudation is present, both denudationrate and denudationdepth approaches come with error related to muons. For the denudationrate approach, which ignores muons and only models on the relationship between ^{10}Be concentration and nucleon spallation production rates (Eq. 6), there is a slight overestimation of exposure age and inheritance. This is because the inversion process attributes a small portion of the muogenic concentration to nucleon spallation, and a larger portion is attributed to inheritance. For the denudationdepth approach, there is a negligible error related to the approximate solution for the muogenic contribution (Appendix A).
To evaluate the error related to muogenic production, we generate profiles with various exposure ages and denudation amounts and compare the estimated results with the true values (Fig. 11). We set the contribution of muogenic production to total surface production as 2 %, though in reality the contribution of muons in most places on earth is smaller (Braucher et al., 2011, 2013; Balco et al., 2008; Balco, 2017); therefore, the errors presented here can be treated as a maximum. We find the error introduced by ignoring muons with the denudationrate approach is relatively small (less than 2 % and 5 % of overestimate) when the total denudation is under 1 or 2 times the attenuation length for spallation (Fig. 11a). Above 3 times this attenuation length, the error grows drastically until no meaningful result can be found. Compared to the denudationrate approach, the denudationdepth approach reduces the error by at least 1 order of magnitude (Fig. 11b); even with very large total denudation (5 times the attenuation length), the error is smaller than 0.3 %. Like exposure age, the error in estimated inheritance is related to denudation and exposure age. As demonstrated in Fig. 11c the overestimation of inheritance from the denudationrate approach increases with the surface age but slightly decreases as the denudation increases. As a comparison, the amount of error from the denudationdepth approach (Fig. 11d) is 1 order of magnitude smaller than the denudationrate approach.
4.2.3 Other sources of error
Additional sources of error exist for CN depth profile ages that we do not consider in this study. Timedependent phenomena are not considered. Constant production rate is an important assumption needed to simplify the nuclear buildup process to apply a linearregression approach. In fact, the production rate is time dependent because the strength of Earth's magnetic field varies with time (Balco, 2017; Desilets et al., 2006; Dunai, 2001; Lifton et al., 2005; Stone, 2000). Extending our model to account for a temporally variable production rate is beyond the scope of the present study. Constant inheritance and a single value for sediment density are other key assumptions for our approach. Sediments sampled from depth profiles are assumed to be well mixed at the time of deposition and to have been deposited rapidly, such that the inherited concentration and density should be the same at every depth. This will not be true for sites with incremental deposition, and for sites where the depositional process or catchmentwide denudation rates vary with time.
We introduce a leastsquares linear inversion approach to solve cosmogenic nuclide concentration depth profiles for surface exposure age and inheritance, considering denudation rate, denudation depth, and muogenic production. This method allows propagation of error sources using Monte Carlo sampling to infer full probability distributions of age and inheritance. In addition, our model presents a straightforward way to assess the tradeoffs between exposure age and denudation rate.
Based on the inversion results of simulated profiles, we show that the leastsquares linear regression is a robust approach suitable for most exposure dating scenarios. The accuracy of linear inversion is comparable to a forward Monte Carlo modeling approach for most circumstances, except for the more deviated (noisier) sample sets. Importantly, neither inversion approach consistently outperforms the other. For surfaces with low inheritance or age, as demonstrated with the Beida River sample site and the simulated profiles, unconstrained inversion could lead to incorrect distributions for inversion results. Therefore, the linear regression we proposed here requires applying the nonnegative constraint on such depth profile data.
For surfaces with no denudation, the inversion in Eq. (4) provides an exact solution. For surfaces with denudation, the approximation of muogenic production using the denudationdepth approach (Eq. 10) introduces negligible error even for surfaces with a large amount of denudation. The denudationrate approach (Eq. 7), though less accurate, provides a useful tool to explore the exposure rate vs. denudation rate relationship. Examples of deep profiles suggest that the linear inversion approach works equally well for samples that were collected deeper than 2 m from the surface, with or without denudation. However, extra caution is needed when collecting and analyzing deep samples to minimize measurement error, as the resulting ages are much more sensitive to the scatter of concentration values relative to nearsurface profiles.
Regardless of whether employing linear regression or a forward approach, surfaces that have undergone a large amount of denudation will be subject to large uncertainties related to the denudation rate or depth. It also becomes more tenuous in such cases to assume that the denudation rate was constant throughout the history of the deposit. It is entirely possible that a sample site may have experienced episodic erosional episodes instead of constant denudation rate, which would lead to errors in the age not accounted for in the methods described here.
When ε>0, the effective exposure age T_{e} takes the following form:
This suggests values for T_{e} value are functions of both ε and t for different production pathways. Between the two variables, ε may be known, while t is unknown. Therefore, our aim is to rewrite Eq. A1 into an approximate form where t can be isolated.
We first take a natural logarithm of the T_{em} over T_{en} ratio:
The expansion of a function, $f\left(x\right)=\mathrm{ln}\left(\frac{\mathrm{1}{e}^{x}}{x}\right)$, x=B_{i}t, in Eq. (A2) may be achieved by writing a Maclaurin series with the following form:
To write the expansion, we first rewrite f(x) as
In this form, f(x) goes to zero when x goes to zero, and therefore we have
and
and
We omit higherorder derivatives from the series expansion.
Bringing Eq. (A5) into Eq. (A3), we have the expansion of f(x) as
Bringing Eq. (A6) into Eq. (A2), the natural logarithm of T_{em} over T_{en} ratio is
The contribution of the third term, t^{3}, is negligible and may be neglected. Making the substitution D = εt, the firstorder terms in Eq. (A7),
This result is independent on the exposure age, t.
With the same substitution, the secondorder terms in Eq. (A7),
In Eq. (A9), the first term of the righthand side is independent of age, t, while the second term is dependent of age. We therefore choose to omit the second term of Eq. (A9) in order to develop an ageindependent approximation. We find that this term may be omitted for two reasons. First, the absolute value of Eq. (A8) is at least 1 order of magnitude larger than Eq. (A9), and therefore omitting one term from Eq. (A9) will not lead to significant decrease in accuracy of the overall approximation. Second, for young surfaces, λt is sufficiently small that the second term of Eq. (A9) is much smaller than the first term, which means omitting it will lead to an even smaller decrease in accuracy.
Therefore, an approximate form of Eq. (A7) that is independent of t is
and the ratio between muon and nucleon effective age can be approximated as
We placed our code at GitHub (https://github.com/YiranWangYR/10BeLeastSquares, last access: 15 August 2022; https://doi.org/10.5281/zenodo.6992982, Wang, 2022).
The ^{10}Be depth profile data sets we used in this paper are published in Hidy et al. (2010, https://doi.org/10.1029/2010GC003084) and Wang et al. (2020, https://doi.org/10.1029/2019TC005901). All the data are given in the tables and in the Supplement of this paper.
The supplement related to this article is available online at: https://doi.org/10.5194/gchron45332022supplement.
YW and MEO built the linearregression models and wrote the paper together. YW developed the approximation for muogenic process and wrote the MATLAB code for both inversion approaches. YW also conducted simulations with pseudorealistic profiles.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work was supported by the US National Science Foundation (grant number EAR1524734) to Michael E. Oskin, and through Cordell Durrell Geology Field Fund to Yiran Wang. We thank the Associate Editor, Pieter Vermeesch, and the referees, Alan Hidy, Greg Balco, and an anonymous referee, for their valuable comments that helped us greatly improve the quality of our paper.
This research has been supported by the National Science Foundation (grant no. EAR1524734) and the Cordell Durrell Geology Field Fund.
This paper was edited by Pieter Vermeesch and reviewed by Alan Hidy, Greg Balco, and one anonymous referee.
Anderson, R. S., Repka, J. L., and Dick, G. S.: Explicit treatment of inheritance in dating depositional surfaces using in situ ^{10}Be and ^{26}Al, Geology, 24, 47–51, 1996.
Balco, G.: Production rate calculations for cosmicraymuonproduced ^{10}Be and ^{26}Al benchmarked against geological calibration data, Quat. Geochronol., 39, 150–173, https://doi.org/10.1016/j.quageo.2017.02.001, 2017.
Balco, G., Stone, J. O., Lifton, N. A., and Dunai, T. J.: A complete and easily accessible means of calculating surface exposure ages or erosion rates from ^{10}Be and ^{26}Al measurements, Quat. Geochronol., 3, 174–195, https://doi.org/10.1016/j.quageo.2007.12.001, 2008.
Braucher, R., Brown, E. T., Bourlès, D. L., and Colin, F.: In situ produced ^{10}Be measurements at great depths: Implications for production rates by fast muons, Earth Planet. Sci. Lett., 211, 251–258, https://doi.org/10.1016/S0012821X(03)00205X, 2003.
Braucher, R., Del Castillo, P., Siame, L., Hidy, A. J., and Bourlés, D. L.: Determination of both exposure time and denudation rate from an in situproduced ^{10}Be depth profile: A mathematical proof of uniqueness. Model sensitivity and applications to natural cases, Quat. Geochronol., 4, 56–67, https://doi.org/10.1016/j.quageo.2008.06.001, 2009.
Braucher, R., Merchel, S., Borgomano, J., and Bourlès, D. L.: Production of cosmogenic radionuclides at great depth: A multi element approach, Earth Planet. Sci. Lett., 309, 1–9, https://doi.org/10.1016/j.epsl.2011.06.036, 2011.
Braucher, R., Bourlès, D., Merchel, S., Vidal Romani, J., FernadezMosquera, D., Marti, K., Léanni, L., Chauvet, F., Arnold, M., Aumaître, G., and Keddadouche, K.: Determination of muon attenuation lengths in depth profiles from in situ produced cosmogenic nuclides, Nucl. Instrum. Meth. B, 294, 484–490, https://doi.org/10.1016/j.nimb.2012.05.023, 2013.
Brocard, G. Y., van der Beek, P. A., Bourlès, D. L., Siame, L. L., and Mugnier, J. L.: Longterm fluvial incision rates and postglacial river relaxation time in the French Western Alps from ^{10}Be dating of alluvial terraces with assessment of inheritance, soil development and wind ablation effects, Earth Planet. Sci. Lett., 209, 197–214, https://doi.org/10.1016/S0012821X(03)000311, 2003.
Chmeleff, J., von Blanckenburg, F., Kossert, K., and Jakob, D.: Determination of the ^{10}Be halflife by multicollector ICPMS and liquid scintillation counting, Nucl. Instrum. Meth. B, 268, 192–199, https://doi.org/10.1016/j.nimb.2009.09.012, 2010.
Cockburn, H. A. P. and Summerfield, M. A.: Geomorphological applications of cosmogenic isotope analysis, Prog. Phys. Geogr., 28, 1–42, https://doi.org/10.1191/0309133304pp395oa, 2004.
Desilets, D., Zreda, M., and Prabu, T.: Extended scaling factors for in situ cosmogenic nuclides: New measurements at low latitude, Earth Planet. Sci. Lett., 246, 265–276, https://doi.org/10.1016/j.epsl.2006.03.051, 2006.
Dunai, T. J.: Influence of secular variation of the geomagnetic field on production rates of in situ produced cosmogenic nuclides, Earth Planet. Sci. Lett., 193, 197–212, https://doi.org/10.1016/S0012821X(01)005039, 2001.
Ebert, K., Willenbring, J., Norton, K. P., Hall, A., and Hättestrand, C.: Meteoric ^{10}Be concentrations from saprolite and till in northern Sweden: Implications for glacial erosion and age, Quat. Geochronol., 12, 11–22, https://doi.org/10.1016/j.quageo.2012.05.005, 2012.
Frankel, K. L., Brantley, K. S., Dolan, J. F., Finkel, R. C., Klinger, R. E., Knott, J. R., Machette, M. N., Owen, L. A., Phillips, F. M., Slate, J. L., and Wernicke, B. P.: Cosmogenic ^{10}Be and 36Cl geochronology of offset alluvial fans along the northern Death Valley fault zone: Implications for transient strain in the eastern California shear zone, J. Geophys. Res.Sol. Ea., 112, 1–18, https://doi.org/10.1029/2006JB004350, 2007.
Granger, D. E., Lifton, N. A., and Willenbring, J. K.: A cosmic trip: 25 years of cosmogenic nuclides in geology, Bull. Geol. Soc. Am., 125, 1379–1402, https://doi.org/10.1130/B30774.1, 2013.
Hancock, G. S., Anderson, R. S., Chadwick, O. A., and Finkel, R. C.: Dating fluvial terraces with ^{10}Be and ^{26}Al profiles: Application to the Wind River, Wyoming, Geomorphology, 27, 41–60, https://doi.org/10.1016/S0169555X(98)000890, 1999.
Heisinger, B., Lal, D., Jull, A. J. T., Kubik, P., IvyOchs, S., Knie, K., and Nolte, E.: Production of selected cosmogenic radionuclides by muons: 2. Capture of negative muons, Earth Planet. Sci. Lett., 200, 357–369, https://doi.org/10.1016/S0012821X(02)006416, 2002a.
Heisinger, B., Lal, D., Jull, A. J. T., Kubik, P., IvyOchs, S., Neumaier, S., Knie, K., Lazarev, V., and Nolte, E.: Production of selected cosmogenic radionuclides by muons 1. Fast muons, Earth Planet. Sci. Lett., 200, 345–355, https://doi.org/10.1016/S0012821X(02)006404, 2002b.
Hetzel, R., Tao, M., Stokes, S., Niedermann, S., IvyOchs, S., Gao, B., Strecker, M. R., and Kubik, P. W.: Late Pleistocene/Holocene slip rate of the Zhangye thrust (Qilian Shan, China) and implications for the active growth of the northeastern Tibetan Plateau, Tectonics, 23, TC6006, https://doi.org/10.1029/2004TC001653, 2004.
Hidy, A. J., Gosse, J. C., Pederson, J. L., Mattern, J. P., and Finkel, R. C.: A geologically constrained Monte Carlo approach to modeling exposure ages from profiles of cosmogenic nuclides: An example from Lees Ferry, Arizona, Geochem. Geophy. Geosy., 11, Q0AA10, https://doi.org/10.1029/2010GC003084, 2010.
Korschinek, G., Bergmaier, A., Faestermann, T., Gerstmann, U. C., Knie, K., Rugel, G., Wallner, A., Dillmann, I., Dollinger, G., von Gostomski, C. L., Kossert, K., Maiti, M., Poutivtsev, M., and Remmert, A.: A new value for the halflife of ^{10}Be by HeavyIon Elastic Recoil Detection and liquid scintillation counting, Nucl. Instrum. Meth. B, 268, 187–191, https://doi.org/10.1016/j.nimb.2009.09.020, 2010.
Lal, D.: Cosmic ray labeling of erosion surfaces: in situ nuclide production rates and erosion models, Earth Planet. Sci. Lett., 104, 424–439, https://doi.org/10.1016/0012821X(91)90220C, 1991.
Lal, D. and Arnold, J. R.: Tracing quartz through the environment, Proc. Indian Acad. Sci. – Earth Planet. Sci., 94, 1–5, https://doi.org/10.1007/BF02863403, 1985.
Laloy, E., Beerten, K., Vanacker, V., Christl, M., Rogiers, B., and Wouters, L.: Bayesian inversion of a CRN depth profile to infer Quaternary erosion of the northwestern Campine Plateau (NE Belgium), Earth Surf. Dynam., 5, 331–345, https://doi.org/10.5194/esurf53312017, 2017.
Lawson, C. L. and Hanson, R. J.: Solving least squares problems, Society for Industrial and Applied Mathematics, 160–165, https://doi.org/10.1137/1.9781611971217, 1987.
Lifton, N. A., Bieber, J. W., Clem, J. M., Duldig, M. L., Evenson, P., Humble, J. E., and Pyle, R.: Addressing solar modulation and longterm uncertainties in scaling secondary cosmic rays for in situ cosmogenic nuclide applications, Earth Planet. Sci. Lett., 239, 140–161, https://doi.org/10.1016/j.epsl.2005.07.001, 2005.
Marrero, S. M., Phillips, F. M., Borchers, B., Lifton, N., Aumer, R., and Balco, G.: Cosmogenic nuclide systematics and the CRONUScalc program, Quat. Geochronol., 31, 160–187, https://doi.org/10.1016/j.quageo.2015.09.005, 2016.
Martin, L. C. P., Blard, P. H., Balco, G., Lavé, J., Delunel, R., Lifton, N., and Laurent, V.: The CREp program and the ICED production rate calibration database: A fully parameterizable and updated online tool to compute cosmicray exposure ages, Quat. Geochronol., 38, 25–49, https://doi.org/10.1016/j.quageo.2016.11.006, 2017.
Mercader, J., Gosse, J. C., Bennett, T., Hidy, A. J., and Rood, D. H.: Cosmogenic nuclide age constraints on Middle Stone Age lithics from Niassa, Mozambique, Quaternary Sci. Rev., 47, 116–130, https://doi.org/10.1016/j.quascirev.2012.05.018, 2012.
Repka, J. L., Anderson, R. S., and Finkel, R. C.: Cosmogenic dating of fluvial terraces, Fremont River, Utah, Earth Planet. Sci. Lett., 152, 59–73, https://doi.org/10.1016/s0012821x(97)001490, 1997.
Rixhon, G., Briant, R. M., Cordier, S., Duval, M., Jones, A., and Scholz, D.: Revealing the pace of river landscape evolution during the Quaternary: recent developments in numerical dating methods, Quaternary Sci. Rev., 166, 91–113, https://doi.org/10.1016/j.quascirev.2016.08.016, 2017.
RuszkiczayRüdiger, Z., Braucher, R., Novothny, Á., Csillag, G., Fodor, L., Molnár, G., and Madarász, B.: Tectonic and climatic control on terrace formation: Coupling in situ produced ^{10}Be depth profiles and luminescence approach, Danube River, Hungary, Central Europe, Quaternary Sci. Rev., 131, 127–147, https://doi.org/10.1016/j.quascirev.2015.10.041, 2016.
Stone, J. O.: Air pressure and cosmogenic isotope production, J. Geophys. Res.Sol. Ea., 105, 23753–23759, https://doi.org/10.1029/2000JB900181, 2000.
Wang, Y.: YiranWangYR/10BeLeastSquares: Combined LinearRegression and Monte Carlo Approach to Modeling Exposure Age Depth ProfilesMatlab Code (v2.1), Zenodo [code], https://doi.org/10.5281/zenodo.6992982, 2022.
Wang, Y., Oskin, M. E., Zhang, H., Li, Y., Hu, X., and Lei, J.: Deducing CrustalScale ReverseFault Geometry and Slip Distribution From Folded River Terraces, Qilian Shan, China, Tectonics, 39, e2019TC005901, https://doi.org/10.1029/2019TC005901, 2020.