The CAIRN method: Automated, reproducible calculation of catchment-averaged denudation rates from cosmogenic radionuclide concentrations

. We report a new program for calculating catchment-averaged denudation rates from cosmogenic nuclide concentrations. The method (Catchment-Averaged denudatIon Rates from cosmogenic Nuclides: CAIRN) bundles previously reported production scaling and topographic shielding algorithms. In addition, it calculates production and shielding on a pixel-by-pixel basis. We explore the effect of sampling frequency across both azimuth ( ∆ θ ) and altitude ( ∆ φ ) angles for topographic 5 shielding and show that in high relief terrain a relatively high sampling frequency is required, with a good balance achieved between accuracy and computational expense at ∆ θ = 8 ◦ and ∆ φ = 5 ◦ . CAIRN includes both internal and external uncertainty analysis, and is packaged in freely available software in order to facilitate easily reproducible denudation rate estimates. CAIRN calculates denudation rates but also automates catchment averaging of shielding and production, and thus can be 10 used to provide reproducible input parameters for the CRONUS family of online calculators.


Introduction
In-situ cosmogenic nuclides, such as 10 Be and 26 Al, are widely used to determine both exposure ages and denudation rates (e.g., Dunai, 2010;Granger et al., 2013;von Blanckenburg and Willenbring, 2014;Granger and Schaller, 2014). A denudation rate is the sum of the chemical weathering 15 rate and physical erosion rate. Since the publication of the seminal papers by Brown et al. (1995), Granger et al. (1996) and Bierman and Steig (1996), dozens of studies have used concentrations of cosmogenic nuclides in stream sediments to quantify denudation rates that are spatially averaged over eroding drainage basins. There are now more than 1000 published catchment-averaged denuda-tion rates (e.g., Portenga and Bierman, 2011;Willenbring et al., 2013a;Harel et al., 2016), with 20 many new studies published each year.
Several authors have provided standardized methods for calculating denudation rates from cosmogenic nuclide concentrations, notably the COSMOCALC package (Vermeesch, 2007) and the CRONUS-Earth online calculator (Balco et al., 2008). Here we make comparisons with the CRONUS calculator version 2.2, so we refer to it as CRONUS-2.2 for clarity. These calculators have been 25 widely adopted by the cosmogenic, quaternary science and geomorphic communities, in large part because they are easily accessible and their methods are transparent (i.e., the source files are available online). These previously published calculators are ideal for calculating denudation rates or ages from a particular site (e.g., an exposed surface or a glacial moraine). Existing calculators rely on the principle that there is an inverse relationship between denudation rate and the concentration 30 of a nuclide, because slower denudation results in more exposure to cosmic rays. In addition, these calculators make use of the fact that the concentration of a nuclide can be inverted for denudation rate if one estimates the production of the nuclide.
In the context of catchment-averaged denudation rates, nuclide production rates will vary in space, and an open-source method of calculating production and inverting nuclide concentration for de-35 nudation rate has yet to emerge. Due to the lack of an open-source tool, a wide variety of approaches to calculating catchment-averaged denudation rates are used in the literature, which makes intercomparison studies challenging (cf., Portenga and Bierman, 2011;Willenbring et al., 2013a;Harel et al., 2016).
Several factors determine the concentration of a cosmogenic nuclide in a sample. For instance, 40 elevation and latitude control the production rate of different cosmogenic nuclides (e.g., Lal, 1991;Dunai, 2000;Stone, 2000;Desilets and Zreda, 2003;Lifton et al., 2005). Production rates vary spatially, thus users of online calculators must calculate the effective production rate within a catchment using a weighted mean of nuclide production in individual pixels. The manner in which these are provided to existing calculators vary. For example, one must feed a single weighted mean produc-45 tion, after shielding corrections, to COSMOCALC. In contrast, one must calculate weighted mean shielding corrections and pass them to CRONUS-2.2, and in addition must calculate a pressure or elevation that reproduces the mean production rate before shielding.
Many authors use an averaging scheme for production wherein production is calculated in each pixel which is then passed to a calculator (e.g., Kirchner et al., 2001;Hurst et al., 2012;Munack 50 et al., 2014;Scherler et al., 2014). In addition, nuclide concentrations can be affected by partial shielding caused by snow cover, surrounding topography, and overlying layers of sediment (e.g., Balco et al., 2008). These again are spatially distributed and so authors reporting catchment-averaged denudation rates frequently report averaged shielding values. Although software packages do exist for calculating spatially averaged topographic shielding (e.g., Codilean, 2006) and snow shielding 55 (e.g., Schildgen et al., 2005), results from these models are not integrated with spatially varying production rates. Finally, in landslide dominated terrain, removal of thick layers of sediment can dilute cosmogenic nuclide concentrations in river sediment (Niemi et al., 2005;Yanites et al., 2009;West et al., 2014). This factor is often not included in denudation calculations. For these reasons, Balco et al. (2008) specifically urged development of tools dedicated to the calculation of catchment-60 averaged denudation rates from cosmogenic nuclide concentrations.
Here we present software that estimates production and shielding of the cosmogenic nuclides 10 Be and 26 Al on a pixel-by-pixel basis, and propagates uncertainty in AMS measurement and cosmogenic nuclide production. Based on these calculations the software can then calculate the expected cosmogenic nuclide concentration from a basin given a spatially homogenous denudation rate. Finally, the software uses Newton iteration to calculate the denudation rate that best reproduces the measured cosmogenic nuclide concentration. We have made this software available through an open-source platform at https://github.com/LSDtopotools/LSDTopoTools_CRNBasinwide to allow community modification and scrutiny, with the goal of enabling users to report denudation rates that can be easily reproduced by other scientists. The software distribution includes instructions for building the 70 software on a virtual machine that can function on common operating systems.

Quantifying denudation rates at a single location
We derive a governing equation that tracks the concentration of a cosmogenic nuclide as it is exposed, exhumed or buried. This approach is adopted because it is the most general: specific scenarios of both steady and transient denudation and burial may therefore be derived. Our approach is 75 broadly similar to that of Parker and Perg (2005), but results are equivalent to those of more widely used derivations (e.g., Lal, 1991;Granger and Smith, 2000).
We begin by conserving the concentration of cosmogenic nuclide i through time t: where C i is the concentration of cosmogenic nuclide i (C i is typically reported in atoms g -1 ; i could 80 be 10 Be or 26 Al, for example), P i is the local production rate of cosmogenic nuclide i (in atoms g -1 yr -1 ) and λ i (yr -1 ) is the decay constant of cosmogenic nuclide i. Production can be a function of latitude, altitude (or atmospheric pressure), magnetic field strength and shielding by rock, soil, water or snow (e.g., Balco et al., 2008).
Cosmogenic nuclides can be produced by both neutrons and muons (e.g., Gosse and Phillips, 85 2001). Production by neutrons is widely modelled using a simple function in which production decays exponentially with depth (e.g., Lal, 1991). Muons, on the other hand, are modelled using a variety of schemes. The CRONUS-2.2 calculator (Balco et al., 2008) implements the scheme of Heisinger et al. (2002a, b), which requires computationally expensive integration of muon stopping over a depth profile. Field-based estimates of muon production demonstrate that Heisinger et al.
The advantage of the Heisinger et al. (2002a) scheme is that it tries to capture the physics of muon 95 passage through the near surface, and specifically the scheme models how the mean energy of muons increases as one moves to greater depths in the subsurface. This affects muon production at depth in a way that is not captured by exponential approximations. Recent work by Marrero et al. (2016) has updated the scheme of Heisinger et al. (2002a, b) to reflecting the muon production rates inferred from field studies. This method still has the disadvantage that it is computationally expensive, to the 100 extent that this computational cost is prohibitive if one is to calculate muon production in numerous pixels across a catchment.
Our approach is to approximate muon production using a sum of exponential functions (e.g., Granger and Smith, 2000;Vermeesch, 2007;Braucher et al., 2009;Schaller et al., 2009). This approach has the advantage of being computationally efficient, but it does not reflect the physics 105 of muon production and therefore does poorly at capturing muon production at depths beyond a few meters. This is unlikely to lead to large errors, however, because muon production makes up a very small percentage of the overall nuclide production at the depths where the physics-based models (Heisinger et al., 2002a, b;Marrero et al., 2016) diverge from the exponential models used in CAIRN. We specifically quantify this difference in Section 7.3, finding that the exponential approx-110 imation leads to differences between the physics-based approximation that are relatively small: for a wide range of denudation rates these differences are less than 2%.
The exponential approximation for nuclide production used in CAIRN is: where P i,SLHL is the surface production rate (atoms g -1 yr -1 ) at sea level and high latitude, F i,j is a 115 dimensionless scaling that relates the relative production of neutron spallation and muon production, S i,j is a dimensionless scaling factor that lumps the effects of production scaling and shielding of cosmic rays, d is a mass per unit area which represents the mass overlying a point under the surface (typically reported in g cm -2 ), and Λ j is the attenuation length for reaction type j (g cm -2 ). The reaction types are j = 0 for neutrons and j = 1 − 3 for muons; muons can be either slow or fast. In 120 general, production from muons relative to neutrons is greater in landscapes with a high denudation rate or at low elevation (Balco et al., 2008).
The depth d, called shielding depth, is related to depth below the surface as: where ζ (cm) is the elevation of the surface, η (cm) is the depth in the subsurface of the sample, z 125 (cm) is the elevation in a fixed reference frame and ρ (g cm -3 ) is the material density, which may be a function of depth. For a constant density, d = ρη.

Solving the governing equation
The governing equation (Eq. 1) has the general form: In our case, p(t) simply equals λ i , which is a constant in this case, and g(t) is equal to P i , which is a function of t.
Equations of this form have the solution: where const is an integration constant and which in the case of the governing equation reduces to: The term g(t) is equal to: The shielding depth, d, is a function of time: where τ is a dummy variable for time that is replaced by the limits after integration. Here t 0 is the initial time and d 0 is the initial shielding depth. In the case where denudation, denoted (g cm -2 yr -1 ), is steady in time this becomes: Here denudation is the rate of removal of mass from above the sample per unit area. If we let the concentration of the cosmogenic nuclide equal C 0 at the initial time, t 0 , and combine Eqs. (5), (7), (8), and (10), we can solve for the integration constant (const) and arrive at a solution for cosmogenic nuclide i at time t: Equation (11) is the full governing equation from which scenario-specific solutions may be derived.

Steady state solution
By convention, we consider the depth profile of cosmogenic nuclide concentration to be steady in time. This allows analytical solution of the cosmogenic nuclide concentration at any point in the 155 basin. At steady state, the particles near the surface have been removed (either through erosion or chemical weathering) at the same rate for a very long time, so we set t 0 = 0 and t = ∞. This results in a simplified form: where is the denudation rate (g cm -2 yr −1 ). If we set d = 0 (that is, we solve for material being 160 eroded from the surface, with no distributed mass loss via chemical weathering), Eq. (12) reduces to Eq. (6) from Granger and Smith (2000) for denudation only (i.e., no burial or exposure), and reduces to Eq. (8) of Lal (1991) if production is due exclusively to neutrons. If Eq. (12) is simplified to neutron only production, assumes the sample is taken from the surface (d = 0), and is solved for erosion rate, one arrives at which is equivalent to the widely used Eq. (11) from Lal (1991). However Eq. (13) requires adjustment for catchment averaged estimates of denudation rates because each point in the landscape from which sediment is derived will have its own local production and shielding factors. This is why a spatially distributed approach is required.

Snow and self shielding
Equation (12) is restrictive in that it only considers material removed from a specific depth, i.e. removed for a single value of d. In reality samples may come from a zone of finite thickness. This finite thickness can contribute some shielding to the sample, i.e. the bottom of a sample is shielded by the mass of the sample that overlies. This shielding is called self shielding and is generally im-175 plemented by assuming that self shielding can simply be approximated by a reduction in neutron production (e.g., Vermeesch, 2007;Balco et al., 2008). Snow can also reduce production of cosmogenic nuclides (e.g., Gosse and Phillips, 2001). Typically these two forms of shielding (snow and self) are incorporated in denudation rate calculators as a scaling coefficient calculated before solving the governing equations (e.g., Vermeesch, 2007;Balco et al., 2008), i.e. snow and self shielding are 180 incorporated into the S i,j term.
Our strategy is slightly different: we calculate snow and self shielding by integrating the cosmogenic nuclide concentration over a finite depth in eroded material. For example, if there is no snow, the concentration of cosmogenic nuclides at a given location is obtained by depth-averaging the steady concentrations from zero depth (the surface) to the thickness of eroded material. If snow is 185 present, the concentration is determined by depth-averaging from the mean snow depth (d s ) to the thickness of the removed material (d t ). Both d s and d t are shielding thicknesses, therefore they are in units of g cm -2 and thus differences in material density are taken into account. The depth-averaged concentration is then: 190 In most applications, the thickness of the removed material will be 0, i.e. the particles from which nuclide concentrations are measured in detrital sediment are derived from a thin layer removed from the surface of the catchment. However, the solution described by Eq. (14) allows some flexibility so that future users can explore different erosion scenarios, for example removal of sediment through mass wasting. We discuss this in Sect. 6, but for the current contribution we focus on steady-state 195 scenarios.

Topographic shielding
In addition to snow and self shielding, locations in hilly or mountainous areas can also receive a reduced flux of cosmic rays because these have been shielded by surrounding topography (Dunne et al., 1999). We adopt the method of Codilean (2006), in which both the effect of dipping sample As ∆θ and ∆φ values decrease, the accuracy with which the shielding is calculated is expected to increase, as we are modelling shielding at finer resolutions. However, this benefit is attenuated by increasing computational cost when these values tend towards (1 • , 1 • ). Codilean (2006) compared the accuracy of different ∆θ and ∆φ by comparing them to a minimum step size of (5 • , 5 • ). Here we 210 exploit the efficiency of our software and the considerable increase in computing power since 2006 to explore smaller step sizes. We make the assumption that a step size of (1 • , 1 • ), corresponding to 32 400 iterations of the shielding algorithm, is an accurate representation of the true shielding factor to the extent that any further refinement in the measurements would not yield a significant change in the results of the cosmogenic nuclide calculations.

215
In order to determine the optimal balance between measurement accuracy and computational efficiency, the full range of (∆θ, ∆φ) pairs were used to derive shielding values for each cell of a worst-case scenario: a high-relief section of the Himalaya (650 km 2 with a 7000 m range in elevation). Table 1 presents the maximum absolute residual value (the error of the pixel with the greatest error) for topographic shielding of the corresponding step sizes when compared to the shielding 220 derived for (1 • , 1 • ). Using values below Codilean (2006)'s suggested threshold of (5 • , 5 • ) gives increasingly small returns for a larger computational burden. We suggest that a (∆θ, ∆φ) pair of (8 • , 5 • ), requiring 810 iterations, is an optimal value for any high relief landscape, yielding a maximum absolute error in our test site of 0.018. On lower relief landscapes the (∆θ, ∆φ) values could be increased to achieve the same level of accuracy. We note that these data are determined using a 90 m 225 resolution DEM, and errors will be higher for finer resolution DEMs (Norton and Vanacker, 2009).
Our topographic shielding calculations rely on two approximations that can lead to some uncertainty. First, the method of Codilean (2006) assumes the horizon attenuates all cosmic rays, and secondly the production of cosmogenic nuclides obeys a power law relationship between the cosine of the zenith angle. Argento et al. (2015) have shown these assumptions to be inaccurate. In addition, 230 the Codilean (2006) method does not include changes to the flux penetration distance on the gradient of the topographic surface (e.g., Dunne et al., 1999;Balco, 2014). Thus our method, while precise, reflects a simplified model of the true physics of topographic shielding.

Production scaling
Production of cosmogenic nuclides varies as a function of both elevation (defined via atmospheric 235 pressure) and latitude and these variations are accounted for by using one of several possible scaling schemes. The classic scaling model of Lal (1991), later modified by Stone (2000), is the simplest and is referred to herein as Lal/Stone. Later scaling models (Dunai, 2000(Dunai, , 2001Desilets and Zreda, 2003;Lifton et al., 2005Lifton et al., , 2014 have incorporated other parameters such as time-dependent geomag-netic field variations, solar modulation, and nuclide-specific information, resulting in a total of seven possible scaling models in the most recent CRONUS calculator (Marrero et al., 2016).
These scaling schemes vary in complexity and therefore computational expense. Time-dependent scaling schemes are far more computationally expensive than the time-independent scheme of Lal/Stone, which does not consider variations in geomagnetic field strength. Recent calibration results (Borchers et al., 2016;Phillips et al., 2016a), including a low-latitude, high-altitude site in Peru (Kelly et al.,245 2015; Phillips et al., 2016b) suggest that the time-independent Lal/Stone scheme performs similarly to the physics-based schemes presented in Lifton et al. (2014) and fits the data better than several other scaling schemes (Dunai, 2000;Desilets and Zreda, 2003;Lifton et al., 2005). For these reasons, we scale production rates using the Lal/Stone scheme. This may lead to some uncertainty because production rates are scaled by the intensity of the Earth's geomagnetic field (e.g., Dunai, 2010), and 250 this intensity has been relatively high over the last 20 kyrs (Valet et al., 2005;Lifton et al., 2014), meaning that this approximation could lead to some uncertainty in samples with slow denudation rates. For example, a rock removal rate of 0.03 mm/yr would remove 60 cm in 20 kyrs, and most production of nuclides occurs in the top 60 cm of rock (Lal, 1991). However, in cases with faster denudation rates, the uncertainty introduced by assuming time-invariant production rates is likely to 255 be much smaller than other sources of uncertainty.
The Lal/Stone scaling scheme requires air pressure, whereas most published studies include only elevation information. We follow the approach of Balco et al. (2008) and convert latitude and elevation data to pressure using the NCEP2 climate reanalysis data (Compo et al., 2011). In certain areas, the ERA-40 reanalysis (Uppala et al., 2005) has been shown to provide more accurate results and 260 due to CAIRN's open source design new models can be readily incorporated into the software. Here we retain the NCEP2 reanalysis to better compare our results with CRONUS-2.2. We note that if users deploy CAIRN as a spatial averaging front end to online calculators, they should be vigilant to use the same air pressure conversion method in both CAIRN and the online calculator.

265
To calculate the concentration of a cosmogenic nuclide, the scaling factors for each production pathway (S i,j ) must be computed. Both topographic shielding and production rate scaling are subsumed within the scaling terms (S i,j ), whereas snow and self shielding are computed separately (see Sect. 2.3). These scaling terms are not computed for each production pathway, but rather are lumped into a single value. We therefore need to compute the values of the individual scaling factors, S i,j .

270
To do this, we follow the method of Vermeesch (2007) and calculate scaling factors using an effective attenuation depth. This is necessary because, when considering multiple production pathways, the scaling terms for individual production mechanisms may vary depending on elevation, shielding, sample thickness, or denudation rates. For example, muogenic pathways will contribute relatively more to production when there is more shielding since muogenic reactions penetrate deeper than 275 spallation.
To determine the scaling terms for the individual production mechanisms (S i,j ), we first compute the total scaling at a location (S tot ), which we define as the product of the production rate scaling (S p ) and the topographic shielding (S t ), that is S tot = S t S p . Production scaling (S p ) is estimated using the Lal/Stone scaling scheme and S t is calculated using our topographic shielding algorithms. 280 We then derive the scaling factors for the individual production mechanisms, S i , j, by employing a virtual attenuation length, Λ v , in units of g cm -2 , following the method of Vermeesch (2007): We must therefore calculate Λ v based on S tot . The individual production mechanisms must be set such that: In Eq. (16), S tot and F i,j are known, whereas S i,j are functions of Λ v . We thus iterate upon Λ v , calculating S i,j using Eq. (15) using Newton's method until Eq. (16) converges on a solution for Λ v .
Once the virtual attenuation length is solved, the S i,j terms are then used in Eq. (14).
3 Denudation rates across a catchment 290 So far we have described the calculations that predict the concentration of a cosmogenic nuclide at one specific location in a basin. All existing cosmogenic nuclide calculators contain some form of these calculations. A wide variety of approaches to scale calculations of cosmogenic nuclide concentrations within a single location to the concentration across entire catchments have been used in the literature. Some authors have averaged production rates on a pixel-by-pixel basis but have 295 not considered topographic shielding (e.g., Belmont et al., 2007;DiBiase et al., 2010;Portenga and Bierman, 2011). Others have calculated an average scaling by integrating the product of topographic shielding and production on a pixel-by-pixel basis (e.g., Ouimet et al., 2009;Hurst et al., 2012;Lupker et al., 2012;Scherler et al., 2014). Another strategy is to calculate both averaged topographic shielding and production scaling values for a basin (e.g., Abbühl et al., 2010). All of these approaches 300 involve some degree of spatial averaging of production, shielding, or a combination of the two before catchment-averaged denudation rates can be estimated.
The approach we take in CAIRN differs in that shielding and production rates are not averaged: these are calculated locally at each pixel. For a given denudation rate, , the concentration of cosmogenic nuclides from each pixel is calculated, then the catchment-averaged concentration is the 305 average of the concentrations from all pixels. This concentration requires no weighting because the denudation rate is considered to be spatially homogenous. The denudation rate for the basin is then iterated upon with Newton's method until the predicted concentration of cosmogenic nuclides emerging from the catchment matches the measured concentration (see Algorithm 1).
We should note here that the version of CAIRN reported in this contribution calculates the de-310 nudation rate across an entire catchment required to produce the observed concentration of the target cosmogenic nuclide. That is, CAIRN assumes denudation rates and target mineral concentrations are the same everywhere in the catchment. Users can explore the effect of instantaneously removing mass by modifying d t Eq. 14, and d t can be spatially heterogeneous. However, even if users choose spatially heterogeneous d t , CAIRN will still calculate the spatially homogenous background 315 denudation rate in light of dilution by mass wasting or stripping of material from the landscape.
Future adaptations of the code could account for nested basins, as this sampling strategy is common in many studies of basin averaged erosion rates, or changes in the concentration of target minerals as employed by, for example Safran et al. (2006) or Carretier et al. (2015). Our software is open source so other groups can make adjustments to CAIRN to suit their needs. These potential future 320 developments, however, are beyond the scope of this contribution.

Uncertainty propagation
We calculate uncertainty from both internal (nuclide concentration uncertainties from accelerator mass spectrometry (AMS) measurements) and external (shielding and production rate) sources using Gaussian propagation of uncertainty following Balco et al. (2008). We do note that some authors 325 have used a Monte Carlo approach in determining cosmogenic nuclide-derived denudation rates because parameter uncertainties can have non-gaussian distributions (e.g., West et al., 2015). CAIRN, at present, does not implement a Monte Carlo uncertainty approach but rather follows conventional Gaussian propagation of uncertainty.

Gaussian propagation of uncertainty 330
Uncertainties are calculated in terms of the denudation rate, , in units of g cm -2 yr -1 , so that no assumption about material density is necessary. The standard deviation of the denudation rate, s , is calculated with where s x is the standard deviation of x, s y is the standard deviation of y, and so on. The variables x 335 and y can represent any uncertain parameter, such as the measurement uncertainty or the production rate of the nuclide. All uncertainties (e.g., nuclide concentration) are assumed to be at the one sigma level unless otherwise stated. The derivatives in Eq. (17) are calculated using the nominal value plus the associated uncertainty and then recalculating the denudation rate in the original, pixel-by-pixel fashion.

340
Three uncertainties are included in the calculation: i) the uncertainty in cosmogenic nuclide concentration, ii) the uncertainty in the production rate at sea level, high latitude (P i,SLHL ), and iii) uncertainty in muon production. Uncertainty in cosmogenic nuclide concentration is reported by authors alongside concentrations. For the cosmogenic nuclide concentration uncertainty, the concentration is used directly to determine the denudation rate uncertainty. For all other parameters, the 345 uncertainty values help to predict a new concentration in each pixel, which is then used to determine denudation rate uncertainty. It is important to note here that we do not calculate uncertainties inherent in the basin-averaging approach which assumes spatial homogeneity in source material and denudation rates, and denudation that is steady in time; we address these uncertainties in Sect. 6.
The uncertainty on the production rate (P i,SLHL ) is based on that used in the CRONUS-2.2 cal-350 culator (Balco et al., 2008): in CRONUS-2.2 the uncertainty is 0.39 atoms cm -2 yr -1 for 10 Be based on a production rate of 4.49 atoms cm -2 yr -1 . This means the uncertainty in CRONUS-2.2 is 8.7% of P i,SLHL for 10 Be. We use this uncertainty for both 10 Be and 26 Al based on our the production rates reported in Table (2). Although the recent CRONUS-Earth calibration (Borchers et al., 2016) has produced new production rates for both 10 Be and 26 Al, the production rate uncertainties remain 355 in the same range as those used here (Phillips et al., 2016a).

Uncertainty from snow shielding
Uncertainties from nuclide concentration, muon production, and production rates are calculated in-365 ternally by our software. Uncertainties from snow and self shielding rely on user-supplied information and therefore must be estimated separately.
Snow shielding can be supplied as a constant effective snow thickness (in g cm -2 ) or spatially distributed information in the form of a raster. Most snow shielding calculations reported in the literature are based on an effective attenuation estimated by the thickness of snow (e.g., Balco et al.,370 2008), but recent field-based measurements indicate that snow may attenuate fluxes of cosmic rays to a greater extent than assumed in simple mass-based snow shielding calculations (Zweck et al., 2013;Delunel et al., 2014). However these uncertainties are small compared to the extreme uncertainties of the thickness, extent and duration of snow over millennial timescales, which are unlikely to ever be well constrained. If no snow shielding values are provided, the software assumes that there is no 375 snow cover.
To calculate uncertainties, users must supply two scenarios for these shielding factors. For example, the user could provide two snow thickness rasters representing variation in snow thickness with 1σ uncertainty (how an author might calculate this could fill another paper and is beyond the scope of our study). The denudation rates of these two scenarios would then be calculated, and the square 380 of the difference in these two denudation rates would then be inserted into Eq. (17). In this way users can calculate shielding uncertainties manually.

Summary of CAIRN parameters for denudation calculations
To summarize, CAIRN predicts cosmogenic nuclide production from neutrons and muons using a four exponential approximation of data from Braucher et al. (2009). These production rates are 385 scaled using Lal/Stone time-independent scaling. Production is calculated at every pixel, with atmospheric pressure calculated via interpolation from the NCEP2 reanalysis data (Compo et al., 2011).
Topographic shielding is calculated using the method of Codilean (2006), and scaled production rates are multiplied by topographic, snow, and self shielding at each pixel. Decay rates, attenuation lengths, and parameters for production are reported in Table 2. Denudation rates are reported in g 390 cm −2 yr −1 because in these units no assumptions about density, which is spatially heterogeneous, are required. In addition, users must report the AMS standard when supplying nuclide concentrations to CAIRN and the concentrations are then normalized following the same scheme as Balco et al. (2008). The CAIRN software prints these parameters to a file so that if they change in the future based on new calibration datasets, users will be able to both view and report these updated 395 values.

Spatial averaging for ingestion by other denudation rate calculators
In addition to producing denudation rates, CAIRN also provides spatially-averaged production rates Self shielding used for spatial averaging is calculated for each pixel k with: where S self,k is the self shielding correction for the k th pixel, d t,k is the shielding thickness for the 415 k th pixel (in g cm -2 ). Equation (18) is used in both COSMOCALC and CRONUS. In the CRONUS calculators, snow shielding is lumped with topographic shielding, therefore the CRONUS calculators presume the user will determine the product of snow and topographic shielding at a site with a method of their choice. COSMOCALC includes a snow shielding calculator which assumes that the equivalent depth of snow (in g cm -2 ) attenuates neutron production following the formula: where S snow,k is the snow shielding correction of the k th pixel and d s,k is the time-averaged depth of snow water equivalent in g cm -2 . We adopt this approximation when performing spatial averaging. Recent work suggests snow may attenuate spallation to a greater degree than predicted by Eq.
(19) (Delunel et al., 2014), and Zweck et al. (2013) suggest that the attenuation length for snow is 425 reduced compared to rock (they report an attenuation length of 109 g cm −2 for snow). However, the uncertainty in historic snow thickness vastly outweighs uncertainties from the snow shielding equation. Although there have been methods suggested to model the evolution of snow thickness through time (e.g., Beniston et al., 2003), the averaging time for eroded particles that accumulate cosmogenic nuclides is on the order of thousands to tens of thousands of years (e.g., Lal, 1991), and 430 reconstructing snow thickness over this timescale is highly uncertain. Users wishing to approximate the Zweck et al. (2013) attenuation lengths can feed CAIRN snow rasters with a thicker apparent snow layer. Overall, we therefore recommend that users include a large range of snow thickness in their uncertainty analysis, guided by historical observations of snow depth.

Spatial averaging for COSMOCALC
In COSMOCALC's erosion calculator (which calculates denudation), the required inputs are a combined shielding and scaling term, the cosmogenic nuclide concentration and the uncertainty in the cosmogenic nuclide concentration. That is, scaling and shielding are combined in a single, spatially averaged term. We calculate the scaling factor S CCtot , which is a lumped shielding and scaling term, where terms are calculated on a pixel-by-pixel basis. Snow shielding is calculated from Eq. (19)

Spatial averaging for the CRONUS calculators
The CRONUS calculators (CRONUS-2.2 and CRONUScalc) require a lumped shielding value and information about either the elevation or pressure of the sample. Spatial averaging of the lumped 450 shielding value, S CRshield , is calculated with: Note that we fold the self shielding into the lumped shielding term so that when transferring data to the CRONUS calculator the sample thickness should be set to 0.
The CRONUS calculators then calculate production using either an elevation or pressure. Produc-455 tion rates are nonlinear with either elevation or pressure, so we must compute an effective pressure that reproduces the mean production rate in the catchment. This is because the arithmetic average of either elevations or pressures within the catchment, when converted to production rate, will not result in the average production rate due to this nonlinearity. CAIRN calculates an effective pressure that reproduces the effective production rate over the catchment. The average production rate is 460 calculated with: We then use the Newton iteration on the Lal/Stone scaling scheme to find the pressure which reproduces the basin average production rate (S ef f p 6 Uncertainties introduced by spatial and temporal variability CAIRN provides uncertainty estimates based on uncertainties in the measurement of nuclide concentrations, and uncertainties in production rates. It does, however, make an assumption of steady 470 erosion, and also makes assumptions likely to be violated almost everywhere on Earth due to the long timescales of geomorphic adjustment, which are on the order of tens of thousands to millions of years (e.g., Fernandes and Dietrich, 1997;Roering et al., 2001;Whipple, 2001;Mudd and Furbish, 2007;Braun et al., 2015) versus climate oscillations that are tens to hundreds of thousands of years (e.g., Lisiecki and Raymo, 2005). In addition, spatial heterogeneity in lithology and target 475 mineral concentrations can lead to additional uncertainty to denudation rate estimates (e.g., Safran et al., 2006;Carretier et al., 2015). Mass wasting can also perturb the concentration of cosmogenic nuclides (e.g., Niemi et al., 2005;Yanites et al., 2009), leading to further uncertainties. Finally, as noted in Sect. 4.2, if snow shielding is to be taken into account, one must estimate the shielding provided by snow over millennial timescales, which, to put it mildly, are difficult to constrain.

480
For the problem of spatially heterogeneous lithology, careful geologic mapping, such as that done by a handful of recent authors (e.g., Safran et al., 2006;West et al., 2014;McPhillips et al., 2014;Carretier et al., 2015), can alleviate some of the uncertainty, but such mapping is logistically challenging. For landsliding, mass removal can be measured in the field, modelled (e.g., Niemi et al., 2005;Yanites et al., 2009), or approximated using mapped landslide inventories (e.g., Hovius et al., 485 1997;Korup, 2005). These may be combined with data on landslide area-volume relationships (e.g., Guzzetti et al., 2009). The main difficulty here is that it takes some time for the cosmogenic nuclide concentration to readjust after mass removal (e.g., Schaller and Ehlers, 2006;Muzikar, 2009;Mudd, 2016) and thus one must make some estimate of not only the spatial distribution of landslides but their evolution through time (Yanites et al., 2009). Simulating nuclide concentrations in settings 490 where denudation rates vary in space and time is possible (Mudd, 2016), but computationally intensive and one must have some confidence that one can accurately reconstruct the temporal evolution of denudation rates. Although recent progress has been made in deriving time series of denudation rates from current topography (e.g., Whittaker et al., 2008;Pritchard et al., 2009;Hurst et al., 2013;Goren et al., 2014;Fox et al., 2014;Croissant and Braun, 2014;Rudge et al., 2015), these methods 495 still suffer from the fact that we lack devices for time travel and struggle to test such reconstructions.
Ultimately, uncertainties in the spatial distribution of denucation and source material, and temporal uncertainties in denudation rates, mean that the uncertainties reported by CAIRN are the minimum uncertainties: they do not take into account landscape transience, lithology, or variation in snow shielding. The fact that catchment-averaged denudation rates carry additional uncertainties is 500 well known, and Dunai (2010)  Our spatially-averaged production scaling and shielding estimates are approximations of spatial averaging reported by other authors. We compare our data to both published denudation rate estimates, 510 and to estimates of denudation rates generated by the CRONUS calculator given the spatial averaging described in Section 5.3. In our comparisons we use seven published cosmogenic datasets (Table 3). These datasets were chosen to span a wide range of locations (i.e., differing latitudes and elevations) and denudation rates. The parameters used by CAIRN for these comparisons are reported in Table 2.

515
It will perhaps aid the reader if we explain how denudation rate estimates may vary between methods. Firstly, production rates are nonlinearly related to elevation, and thus spatial averaging of the product of production scaling and shielding is not the same as the product of the spatial averages of production scaling and shielding. In addition, previous studies and other calculators have chosen different parameters for cosmogenic nuclide production and shielding. For example, 520 past publications have used a wide variety of methods for estimating topographic shielding (e.g., see Table 3). Choices of spallation and muon production rates also affect the final denudation rate.
Consider a measured nuclide concentration that one uses to infer a denudation rate. If one assumes a high production rate (via either muons or spallation), it means that for a given denudation rate the predicted nuclide concentration is higher. Thus, for a given nuclide concentration, the inferred 525 denudation rate is higher if the assumed production rate is higher (see dashed lines in Figure 1). If the inferred shielding is higher, then for a given denudation rate the production is lower, and the inferred denudation for a given concentration will be lower.

Spatial averaging of production and shielding vs pixel-by-pixel calculations
First, we compare results of two methods using the exponential approximation of muon produc-530 tion (Eq. 12), used in both COSMOCALC and the CAIRN calculator. The difference in calculating denudation rates by iterating upon cosmogenic nuclide concentration from all pixels in a basin (the CAIRN method) and calculating it by using a spatial average of the production of scaling and production terms (Eq. 20) is virtually zero if snow and self shielding are spatially homogenous (Figure 2a).
Thus we find that combining all scaling and shielding terms in a single lumped term is adequate for 535 calculating denudation rates if computational power is limited.
Separating production rate scaling from shielding leads to slightly larger uncertainty (Figure 2b), but in terms of the total uncertainty this averaging also leads to small uncertainties (on the order of 1-2% compared to 10-20% from other sources of uncertainty). We suspect that many users will want to compare rates determined by our software with the popular CRONUS calculators (Balco et al., 2008;540 Marrero et al., 2016). The CRONUS calculators internally scale production rates while shielding is supplied by the user. Consequently, the uncertainties plotted in Figure 2b approximate uncertainties arising from the spatial averaging process that users must pass to the CRONUS calculators. Some users may wish to calculate denudation rates using time-dependent scaling schemes, which is not possible in CAIRN, but CAIRN can be used as a front end to the CRONUS calculators via its spatial 545 averaging capabilities with the confidence that this will only introduce relatively small errors.

Comparison with existing denudation rate estimates
Denudation rates reported in the literature from catchment-averaged cosmogenic nuclide concentrations are calculated using a wide variety of methods. The term erosion rate is often substituted for denudation rate although few studies attempt to account for chemical weathering (cf., Kirchner et al.,550 2001; Riebe et al., 2001). Studies differ in their strategies for production rate scaling, topographic, snow, and self shielding, and the manner in which spatial averaging is performed. In many cases there is insufficient detail reported that might enable other groups to reproduce reported denudation rates.
A primary motivation behind CAIRN is to provide an open-source means of computing denudation rates that may then be reproduced by other groups. We have incorporated reported snow shielding 555 from previous studies by inverting Eq. 19 for an annual average snow thickness and then distributing this thickness over the entire DEM. We acknowledge this is a poor representation of snow thickness but snow shielding rasters are rarely available and in most cases there is little reported snow shielding.
The diversity in methods for calculating denudation rates reported in the literature means that it is 560 difficult to compare denudation rates when they come from different studies. This problem has been highlighted by previous data intercomparison studies (Portenga and Bierman, 2011;Willenbring et al., 2013a;Harel et al., 2016). High latitude production rates under Lal/Stone scaling of 10 Be have changed in the last 10 years due to an ever increasing number of calibration sites (e.g., Phillips et al., 2016a) and changing AMS standards (Nishiizumi, 2004). In some cases, muons are not con-565 sidered, whereas other studies use a variety of different muon production schemes (e.g., Table 3).
Topographic shielding is occasionally not considered (particularly in older studies). In some cases the horizon elevation is recorded from a limited number of directions (e.g., COSMOCALC includes a calculator using 8 directions), and in other instances the computational method of Codilean (2006) is used. Studies also cite Dunne et al. (1999) for shielding but this paper lists several methods for 570 calculating shielding: the equations therein depend on the number and geometry of shielding objects and this information is seldom reported. Even when the more robust method of Codilean (2006) is used, the spacing of azimuth and angle of elevation is often not reported.
Studies typically report erosion or denudation rates in dimensions of length per time, but this requires an assumption about density, which can vary spatially and is sometimes not reported. Most 575 studies use a rock equivalent denudation rate (as opposed to a regolith or soil denudation rate) and thus densities assumed are typically rock densities (see Table 3). Because denudation rates are traditionally reported in dimensions of length per time, we do not suggest future authors cease reporting denudation in these dimensions, but we do recommend also reporting denudation rates in dimensions of mass per area per time (e.g., g cm −2 yr −1 ) because these units allow simpler comparison between 580 sites as they require no assumptions about spatially heterogeneous density.
Of our seven example datasets (Table 3) do not give enough information to reproduce their shielding calculations, but we note that authors that employ the equations of Dunne et al. (1999) use a limited number of horizon measurements to calculate shielding. For example in COSMOCALC (Vermeesch, 2007), users are expected to input horizon values at 45 • intervals. Our calculations suggest that this can lead to lower maximum shielding differences between this method and the CAIRN method (Table 1). An example of the 590 potential underestimates of topographic shielding is shown in Figure 4. is lower than the reported denudation rate. Reasons for this vary since the method used to calculate denudation rates vary in each example study, but differences are likely to be due to the higher pro-595 duction rates used in previous studies (Table 3) and slightly greater topographic shielding in CAIRN (see Figure 1).
One component of CAIRN that requires caution is that the snapping of cosmogenic samples to channels is automated: if errors in the DEM place the main channel in the wrong location, or GPS coordinates of the sampling location contain large errors (common in older datasets), there is a chance the basin selected by CAIRN will not be the same as the sampled basin. This can result in large errors as production rates vary significantly with elevation. We have provided a tool in the github repository that allows users to check the basins that are associated with cosmogenic nuclide samples. If these do not match the expected basins, then users will need to manually change the latitude and longitude of the samples until they are located near the correct channel. 605 We wish to emphasize that the relative denudation rates do not change significantly between CAIRN and reported values (as evidenced by a clustering about the 1:1 line in Figure 5). In addition previous studies contain elements modulating denudation rates that are not contained within the current version of CAIRN. For example, Kirchner et al. (2001) reports true physical erosion rather than denudation and Safran et al. (2006) modified their denudation rates based on the quartz 610 content of the source areas.

Comparison with the CRONUS calculators
The results from CAIRN are compared to results from both CRONUS calculators. When comparing output from CAIRN with output from the online CRONUS-2.2 calculator, far larger uncertainties (up to to 40% of the denudation rate) occur. These differences are not controlled by denudation 615 rate ( Figure 6a) but are instead mainly a function of the production rate (Figure 6b). In the previous section, we found that differences due to spatial averaging and separation of shielding from production scaling are small. The large difference is primarily due to the difference in spallation production rates and the over-production of muons in CRONUS version 2.2, as described by Balco et al. (2013). According to Balco et al. (2013), future versions of this CRONUS calculator will be 620 updated to have significantly reduced muogenic production consistent with recent studies (Braucher et al., 2003(Braucher et al., , 2011(Braucher et al., , 2013Phillips et al., 2016a). If production rates in CRONUS are changed to reflect the production rates from Braucher et al. (2011), we find that differences are quite small (Figure 7).
We see from this figure that in locations with high production rates just under half of these differences between CAIRN and CRONUS-2.2 are from the different spallation rates, whereas in locations 625 with low production rates, most of the differences are due to the higher muon production present in CRONUS-2.2.
The other CRONUS calculator, CRONUScalc, incorporates new spallation production rates and muon production is calculated using production rates based on a deep core from Antarctica (Marrero et al., 2016;Phillips et al., 2016a). In order to examine the underlying source of discrepan-630 cies between the three calculators, we plot the total and muon production rates for the CAIRN, CRONUS-2.2, and CRONUScalc calculators in Figure (8). The production rates for CRONUS-2.2 are calculated directly from the MATLAB scripts available online. The CRONUScalc production rates are approximated as a three exponential analytical function with parameters shown in Table 4.
Although total production rates appear relatively similar, CRONUScalc and CAIRN predict signifi-635 cantly smaller muon contributions that CRONUS-2.2. The result is that for the same denudation rate, the CRONUS-2.2 calculator produces significantly more (in some cases 40% more) atoms than using CAIRN or CRONUScalc (Figure 9) leading to a large discrepancy in calculated denudation rates between CRONUS-2.2 and the other two calculators (CAIRN and CRONUScalc), which both incorporate more recent muon production rates. The CAIRN outputs of topographic shielding, as well as 640 the spatial averaging of both production scaling and shielding, are independent of these calculators and will still provide spatial averaging for use with future calculator versions, even as production rates and mechanisms are updated.
We have used the spatially averaged shielding and scaling outputs from CAIRN to determine differences between CAIRN and CRONUScalc. We find that there is a 2.5% to 5% difference between 645 the denudation rates predicted by CAIRN and those predicted by CRONUScalc ( Figure 10). Currently CRONUScalc is not able to calculate very high denudation rates (for rates greater than~0.06 g cm −2 yr −1 the current version of CRONUScalc crashes; it was designed for exposure ages and becomes computationally unstable at high erosion rates) so we cannot compare CAIRN to CRONUScalc for all of the example datasets. The differences in Figure Table 4). It is important to note that the CAIRN implementation of muons from Marrero et al. (2016) assumes that Λ = 160 g cm −2 for spallation, whereas in CRONUScalc this attenuation length can vary as a function of latitude and pressure. We compare the denudation rates from CAIRN using the production parameters in Table 4 ( CAIRN −CRC ) with the default production scheme of Braucher et al. (2009) 660 in Figure (11). The differences here are smaller (mostly less than 2%) suggesting that much of the difference seen in Figure (10) is due to spatial averaging.

Conclusions
We present an automated, open-source method for calculating catchment-averaged denudation rates based on the concentrations of in-situ cosmogenic nuclides collected in stream sediment. Our catchment-665 averaged denudation rate method (CAIRN) predicts cosmogenic nuclide concentrations based on pixel-by-pixel scaling and shielding. These concentrations are then averaged to predict the catchmentaveraged concentration. Newton iteration is then used to find the denudation rate for which the predicted concentration matches the measured concentration and to derive associated uncertainties. In addition, CAIRN provides spatially averaged shielding and scaling values that can be used by other 670 popular calculators (which do not provide spatial averaging, e.g., CRONUS and COSMOCALC).
The CAIRN method is provided as open-source software so that reported denudation rates can be easily reproduced.
The CAIRN method is intended to streamline the computation and reporting of catchment-averaged denudation rates, but it has limitations that may be the subject of future developments. At the mo-675 ment CAIRN assumes steady erosion; there is no facility for incorporating transient erosion rates which might affect nuclide concentrations in transient landscapes (e.g., Willenbring et al., 2013b;Mudd, 2016). In addition, the method does not include a facility for nesting basins in which the denudation rate in a large basin incorporates the denudation rates from smaller basins that it contains.
The calculator cannot account for differing source areas of material, so at the moment it is not capa-680 ble of using different particle size fractions to identify denudation hot spots (e.g., Riebe et al., 2015;Carretier et al., 2016). Despite these limitations, the CAIRN method addresses the need to provide transparent, reproducible estimates of denudation rates.
Our open source framework allows other users to update the algorithms (e.g., a nesting function could be built on top of the current CAIRN architecture) and different atmospheric reanalysis data 685 or new muon scaling schemes can be added as needed in the future. Thus we hope it will provide a platform for more nuanced estimates of denudation rates from cosmogenic nuclides in the future.

Software and data availability
The software is available at the LSDTopoTools Github website (https://github.com/LSDtopotools/).

Denudation rate
Greater shielding Greater production Figure 1. A schematic drawing of the predicted concentration of a nuclide as a function of denudation rate. If production rates are assumed to be higher, the predicted concentration will be higher for a given denudation rate.
If shielding is greater, the predicted concentration is lower for a predicted denudation rate. Thus assumptions about production and shielding will affect the inferred denudation rate given a sample with fixed concentration, shown with the dashed lines.

a.
b. Figure 2. Differences between the denudation rate calculated by CAIRN ( CAIRN ) and the denudation rate using the production factor (SCCtot) (which includes production scaling and shielding) passed to COSMOCALC ( CC ) (a.), and differences between the denudation rate calculated by CAIRN ( CAIRN ) and the denudation rate using separate spatial averages for shielding and production scaling that are then averaged ( CC−CRON U S ) as a function of production factor (b.). In this case the production factor is calculated by multiplying the separately averaged shielding (S CRShield ) and scaling (S ef f p ) factors. This approach emulates the data requirements for CRONUS-2.2, which calculates production scaling and accepts a single shielding factor (for snow and topography combined). Although the shielding and scaling emulate data requirements for CRONUS-2.2, the denudation rate is calculated using the exponential production method of CAIRN and COSMOCALC.    Figure 8. Production rates of 10 Be as a function of depth for muons only (a.) and total production (b.). These production rates are calculated using the Lal/Stone scaling at 70 degrees North and with a pressure of 1007 hPa (near sea level). Note the logarithmic depth scale: eroding particles spend a large amount of their exposure history below 100 g cm −2 and so increased muon production at these depths, despite being a small fraction of the total production, plays a significant role in determining the total nuclide concentration (see Figure 9).  (Balco et al., 2008) and CRONUScalc (Marrero et al., 2016). These concentrations are calculated for a hypothetical site at 70 degrees North and near sea level (1007 hPa). Note that although the default production scheme in CAIRN is the Braucher et al. (2009) scheme, the production from CRONUScalc (Marrero et al., 2016) can also be used (see Table 4). Figure 10. Differences between the denudation rate calculated by CAIRN ( CAIRN ) and the denudation rate calculated with CRONUScalc ( CRCalc ) as a function of CAIRN denudation rate for selected studies. Figure 11. Differences between the denudation rate calculated by CAIRN using the parameters in Table 4 to approximate CRONUScalc production ( CAIRN −CRCalc ) and the denudation rate calculated with CRONUScalc ( CRCalc ) as a function of CAIRN denudation rate for selected studies.