Short communication: Increasing vertical attenuation length of cosmogenic nuclide production on steep slopes negates topographic shielding corrections for catchment erosion rates

Interpreting catchment-mean erosion rates from in situ produced cosmogenic ¹⁰Be concentrations in stream sediments requires calculating the catchment-mean ¹⁰Be surface production rate and effective mass attenuation length, both of which can vary locally due to topographic shielding and slope effects. The most common method for calculating topographic shielding accounts only for the reduction of nuclide production rates due to shielding at the surface, leading to catchment-mean corrections of up to 20 % in steep landscapes, and makes the simplifying assumption that the effective mass attenuation length for a given nuclide production mechanism is spatially uniform. Here I evaluate the validity of this assumption using a simplified catchment geometry with mean slopes ranging from 0 to 80^∘ to calculate the spatial variation in surface skyline shielding, effective mass attenuation length, and the total effective shielding factor, defined as the ratio of the shielded surface nuclide concentration to that of an unshielded horizontal surface. For flat catchments (i.e., uniform elevation of bounding ridgelines), the effect of increasing vertical attenuation length as a function of hillslope angle and skyline shielding exactly offsets the effect of decreasing surface production rate, indicating that no topographic shielding correction is needed when calculating catchment-mean vertical erosion rates. For dipping catchments (as characterized by a plane fit to the bounding ridgelines), the catchment-mean surface nuclide concentrations are also equal to that of an unshielded horizontal surface, except for cases of extremely steep range-front catchments, where the surface nuclide concentrations are counterintuitively higher than the unshielded case due to added production from oblique cosmic ray paths at depth. These results indicate that in most cases topographic shielding corrections are inappropriate for calculating catchment-mean erosion rates, and are only needed for steep catchments with nonuniform distributions of quartz and/or erosion rate. By only accounting for shielding of surface production, existing shielding approaches introduce a slope-dependent systematic error that could lead to spurious interpretations of relationships between topography and erosion rate.

1 Introduction

Measurement of in situ produced cosmogenic ¹⁰Be concentrations in stream sediments has rapidly become the primary tool for quantifying catchment-scale erosion rates over timescales of 10³–10⁵ years (Brown et al., 1995; Granger et al., 1996; von Blanckenburg, 2006; Portenga and Bierman, 2011; Codilean et al., 2018). Although requiring a number of simplifying assumptions about the steadiness of erosion and sediment transport (Bierman and Steig, 1996), erosion rates determined from ¹⁰Be concentrations in stream sediments have yielded insights into a number of key questions in tectonic geomorphology regarding the sensitivity of erosion rates to spatiotemporal patterns of climate, tectonics, and rock strength (e.g., Safran et al., 2005; Binnie et al., 2007; Ouimet et al., 2009; DiBiase et al., 2010; Bookhagen and Strecker, 2012; Miller et al., 2013; Scherler et al., 2017).

In contrast to point measurements, where a clear framework exists for converting ¹⁰Be concentrations to either a surface exposure age or steady erosion rate (e.g., Balco et al., 2008; Marrero et al., 2016), the interpretation of ¹⁰Be concentrations in stream sediment requires accounting for the spatial variation in elevation, latitude, quartz content, and erosion rate throughout a watershed (Bierman and Steig, 1996; Granger and Riebe, 2014). Additionally, topographic shielding corrections that account for the reduction of cosmic radiation flux on sloped or skyline-shielded point samples (Dunne et al., 1999) are applied to varying degrees for determining catchment-mean production rates. These shielding corrections are either applied at the pixel level (e.g., Codilean, 2006), catchment level (e.g., Binnie et al., 2006), or not at all (e.g., Portenga and Bierman, 2011). Although typically small (< 5 %), topographic corrections can be as large as 20 % for steep catchments (e.g., Norton and Vanacker, 2009). Because these corrections vary as a function of slope and relief, any systematic corrections can influence interpretations of relationships between topography and erosion rate.

The pixel-by-pixel skyline-shielding algorithm of Codilean (2006) results in the largest topographic shielding corrections, and has gained popularity due to its ease of implementation in the software packages TopoToolbox (Schwanghart and Scherler, 2014) and CAIRN (Mudd et al., 2016), the latter of which was used to recalculate published ¹⁰Be-derived catchment erosion rates globally as part of the OCTOPUS compilation project (Codilean et al., 2018). A key simplification of the Codilean (2006) approach is that it accounts only for the skyline shielding of surface production, and not for the change in shielding with depth, which determines the sensitivity of the effective mass attenuation length for nuclide production as a function of surface slope and skyline shielding (Dunne et al., 1999; Gosse and Phillips, 2001). Because a change in the effective mass attenuation length will directly influence the inferred erosion rate of a sample (Lal, 1991), the full depth-integrated implications of topographic shielding must be accounted for when inferring catchment erosion rates from ¹⁰Be concentrations in stream sediments.

Here I model the shielding of incoming cosmic radiation flux responsible for spallogenic production at both the surface and at depth for a simple catchment geometry to evaluate, as a function of catchment slope and relief, the total effect of topographic shielding on surface nuclide concentrations and the partitioning of shielding into surface skyline shielding and changes to the effective mass attenuation length. I then apply this framework to catchments that have a net dip (i.e., dipping plane fit to boundary ridgelines) and compare calculations of total shielding to those from typical pixel-by-pixel skyline-shielding corrections.

2 Theory

The incoming cosmic ray intensity, I(θ,d), responsible for in situ cosmogenic nuclide production by neutron spallation can be most simply described as a function of the incident ray path inclination angle above the horizon, θ, and the mass distance, d (g cm⁻²), traveled along that pathway:

$$\begin{array}{}\text{(1)}& I(\mathit{\theta },d)=I_{\mathrm{0}}{\sin }^m\mathit{\theta }{e}^{-d/\mathit{\lambda }},\end{array}$$

where I₀ is the maximum cosmic ray intensity at the surface, m is an exponent typically assumed to have a value of 2.3 (e.g., Nishiizumi et al., 1989), and λ is the mass attenuation length (g cm⁻²) for unidirectional incoming radiation (Dunne et al., 1999). The mass attenuation length for unidirectional radiation, λ, differs from the nominal mass attenuation length that describes cosmogenic nuclide production as a function of depth, Λ, due to the integration of radiation from all incident angles. Assuming m=2.3, a value of λ=1.3Λ results in a close match for horizontal unshielded surfaces with exponential production profiles typical of spallation reactions (Dunne et al., 1999; Gosse and Phillips, 2001).

For a horizontal surface sample (d=0), the unshielded total cosmic radiation flux, F₀, is described by

$$\begin{array}{}\text{(2)}& F_{\mathrm{0}}={\int }_{\mathit{\phi }=\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathit{\theta }=\mathrm{0}}^{\mathit{\pi }/\mathrm{2}}{I}_{\mathrm{0}}{\sin }^m\mathit{\theta }\cos \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }={\displaystyle \frac{\mathrm{2}\mathit{\pi }{I}_{\mathrm{0}}}{m+\mathrm{1}}},\end{array}$$

where φ is the azimuthal angle of incoming radiation, and the term cosθ accounts for the convergence of the spherical coordinate system. For point samples that are either at depth (d>0) or have an incomplete view of the sky due to topographic shielding by thick (d≫λ) objects, the total cosmic radiation flux, F, is modulated by a shielding factor, S, such that

$$\begin{array}{}\text{(3)}& S={\displaystyle \frac{F}{F_{\mathrm{0}}}}={\displaystyle \frac{m+\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{\int }_{\mathit{\phi }=\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}\left(\mathit{\phi }\right)}^{\mathit{\pi }/\mathrm{2}}{\sin }^m\mathit{\theta }{e}^{-d(\mathit{\theta },\mathit{\phi })/\mathit{\lambda }}\cos \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi },\end{array}$$

where θ₀(φ) is the inclination angle above the horizon of topographic obstructions in the direction φ and d(θ,φ) varies as a function of both ray path azimuth and inclination angle (Dunne et al., 1999; Gosse and Phillips, 2001).

Equation (3) has two implications for interpreting exposure ages or erosion rates from cosmogenic nuclide concentrations of samples partially shielded by skyline topography (θ₀(φ)>0). First, skyline shielding will reduce the surface production rate of cosmogenic nuclides by a factor of S₀:

$$\begin{array}{}\text{(4)}& S_{\mathrm{0}}={\displaystyle \frac{m+\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{\int }_{\mathit{\phi }=\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}\left(\mathit{\phi }\right)}^{\mathit{\pi }/\mathrm{2}}{\sin }^m\mathit{\theta }\cos \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }.\end{array}$$

Second, due to shielding of low-intensity cosmic radiation below incident angles of θ₀(φ), the effective mass attenuation length, Λ_eff, will increase relative to the nominal mass attenuation length for describing cosmogenic nuclide production as a function of depth, Λ (Dunne et al., 1999; Gosse and Phillips, 2001). For calculating surface exposure ages, only the reduction in surface production rate due to skyline shielding needs to be taken into account, and Eq. (4) is easily calculated for single points in the landscape (e.g., Balco et al., 2008). However, for determining erosion rates both the surface shielding and changing effective attenuation length must be accounted for, which requires solving Eq. (3) numerically as a function of vertical depth below the surface, as described in Sect. 3 below.

The importance of accounting for both changes in surface production rate, P, and changes in the effective mass attenuation length, Λ_eff, is illustrated by the analytical solution for nuclide concentration, C, measured on a steadily eroding surface for a stable nuclide with an exponential decrease in production rate with depth:

$$\begin{array}{}\text{(5)}& C=P{\mathrm{\Lambda }}_{\mathrm{eff}}/E,\end{array}$$

where E is erosion rate (g cm⁻² yr⁻¹; Lal, 1991). From Eq. (5) it is clear that increasing Λ_eff counters the effect of decreasing P in determining the surface nuclide concentration (or alternatively for inferring erosion rate).

3 Topographic shielding model for a simplified catchment geometry

For stream sediment samples that require calculating cosmogenic nuclide production rates across an entire catchment, solving Eq. (3) as a function of depth is presently too computationally intensive to be practical. Consequently, numerical implementations of topographic shielding calculations at the catchment scale make the simplifying assumption that Λ_eff=Λ, and thus S=S₀ (Codilean, 2006; Schwanghart and Scherler, 2014; Mudd et al., 2016), accounting only for the effect of decreasing surface production rate, P. Here I use a simplified catchment geometry to solve Eq. (3) and directly calculate the impact of topographic shielding and surface slope on interpretations of catchment erosion rates from cosmogenic nuclide concentrations in stream sediments. For simplicity, I assume that cosmogenic nuclides are only produced by neutron spallation (i.e., Λ=160 g cm⁻²) and that the erosion rate is high enough that radioactive decay is negligible (i.e., E > 0.01 g cm⁻² yr⁻¹ for ¹⁰Be).

Throughout the analysis below, both the effective mass attenuation length, Λ_eff, and erosion rate, E, are defined in the vertical, rather than slope-normal direction. The vertical (with respect to the geoid) reference frame was chosen for three reasons. First, most studies report erosion rate as a vertical lowering rate and assume primarily vertical exhumation pathways. Second, treatment of slope-normal processes introduces a grid-scale dependence of erosion and shielding calculations that varies with topographic roughness (Norton and Vanacker, 2009). Third, for the case of uniform erosion rate, the resulting shielding calculations do not depend on the choice of reference frame, as long as the orientation of Λ_eff and E are defined similarly.

Figure 1 Model catchment setup, showing (a) map view, (b) cross section along azimuthal angle φ (note that |γ|=α for φ=0), and (c) close-up of hillslope cross section.

Figure 2 Dipping catchment shielding geometry, illustrated using example from the San Gabriel Mountains, California, USA. Image is centered on 34.20^∘ N, 117.61^∘ W. Colored lines indicate planes fit through bounding ridgelines dipping at angle β. $\stackrel{\mathrm{‾}}{S_{\mathrm{0}}}$ indicates mean surface skyline-shielding parameter calculated using algorithm of Codilean (2006), and ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$ indicates the mean total effective shielding factor calculated from the simplified catchment geometry.

3.1 Simplified catchment geometry and model setup

Catchment geometry is simplified as an infinitely long v-shaped valley with width 2L_h and uniform hillslope angle α (Fig. 1). Because the ridgelines have uniform elevation, there is no net dip to the catchment; the effect of valley inclination will be assessed in Sect. 3.3. At a horizontal distance from the ridgeline x and vertical depth below the surface z, the shielding factor, S(x,z), is defined as

$$\begin{array}{ll}{\displaystyle }S\left(x,z\right)& {\displaystyle }={\displaystyle \frac{m+\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{\int }_{\mathit{\phi }=\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}(x,L_{\mathrm{h}},z,\mathit{\phi },\mathit{\alpha })}^{\mathit{\pi }/\mathrm{2}}{\sin }^m\mathit{\theta }{e}^{-d(z,\mathit{\rho },\mathit{\theta },\mathit{\gamma }\left(\mathit{\alpha },\mathit{\phi }\right))/\mathit{\lambda }}\\ \text{(6)}& {\displaystyle }& {\displaystyle }\cos \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi },\end{array}$$

where ρ is rock density and γ is the apparent dip of the hillslope in the azimuthal direction φ (Fig. 1b). The inclination angle integration limit, θ₀, is a function of topographic skyline-shielding inclination, and can be determined geometrically (Fig. 1) as

$$\begin{array}{}\text{(7)}& \tan {\mathit{\theta }}_{\mathrm{0}}=\left\{\begin{array}{l}{\displaystyle \frac{\left(x\tan \mathit{\alpha }+z\right)\cos \mathit{\phi }}{\mathrm{2}{L}_{\mathrm{h}}-x}},\mathrm{0}\le \mathit{\phi }<{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\\ -\tan \mathit{\alpha }\cos \mathit{\phi }-{\displaystyle \frac{z}{x}}\cos \mathit{\phi },{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\le \mathit{\phi }\le \mathit{\pi }\end{array}\right..\end{array}$$

The apparent dip, γ, can be derived from the model geometry in Fig. 1 as

$$\begin{array}{}\text{(8)}& \tan \mathit{\gamma }=-\tan \mathit{\alpha }\cos \mathit{\phi },\end{array}$$

and the mass distance traveled through rock by a given incident ray as

$$\begin{array}{}\text{(9)}& d={\displaystyle \frac{\mathit{\rho }z\cos \mathit{\gamma }}{\sin \left(\mathit{\theta }-\mathit{\gamma }\right)}}.\end{array}$$

Equation (5) was numerically solved for a series of hillslopes over a grid of $(x/L_{\mathrm{h}}=[\mathrm{0},\mathrm{1}];\mathit{\rho }z/\mathrm{\Lambda }=[\mathrm{0},\mathrm{40}\left]\right)$ with horizontal spacing $\mathrm{d}x=L_{\mathrm{h}}/\mathrm{500}$ and vertical spacing $\mathrm{d}z=\mathrm{\Lambda }/\mathrm{500}\mathit{\rho }$. To characterize mean slope controls on the shielding factor, S(x,z), the above calculation was applied to nine hillslopes with mean slope, α, ranging from 0 to 80^∘ in 10^∘ increments. Because $L_{\mathrm{h}}\gg \mathrm{\Lambda }/\mathit{\rho }$ for most natural landscapes, the resulting distribution of shielding factors is independent of hillslope scale.

Figure 3 Total shielding factor, S(x,z), as a function of normalized vertical depth and distance from ridgeline for varying hillslope angle, α.

3.2 Calculation of shielding parameters from model results

After applying Eq. (6) to a hillslope, it is straightforward to calculate the surface skyline-shielding component, $S_{\mathrm{0}}\left(x\right)=S(x,\mathrm{0})$. This skyline-shielding component should match the topographic shielding factor determined from the algorithm of Codilean (2006); for comparison this parameter was calculated at each pixel in the model catchment using TopoToolbox (Schwanghart and Scherler, 2014). Two additional parameters were calculated at each slope position using Eq. (5): the effective vertical mass attenuation length, Λ_eff(x), and the total effective shielding factor, C_eff(x).

Although spallogenic production of cosmogenic nuclides following Eq. (1) is well-described by an exponential decrease with depth for horizontal unshielded surfaces, this is not true in general for shielded samples (Dunne et al., 1999). The effective vertical mass attenuation length, Λ_eff(x), is approximated by the vertical depth below the surface at which the shielding factor is 5 % of the surface shielding (i.e., 3 e-folding lengths) such that

$$\begin{array}{}\text{(10)}& S\left(x,{\displaystyle \frac{\mathrm{3}{\mathrm{\Lambda }}_{\mathrm{eff}}\left(x\right)}{\mathit{\rho }}}\right)=\mathrm{0.05}S(x,\mathrm{0}).\end{array}$$

If nuclide production as a function of depth deviates from an exponential decline, it is inaccurate to use the analytical relationship between surface sample concentration, C(x) (atoms g⁻¹), and steady-state vertical erosion rate, E (g cm⁻² yr⁻¹), typically applied to eroding samples

$$\begin{array}{}\text{(11)}& C\left(x\right)={\displaystyle \frac{S\left(x\right)P_{\mathrm{0}}\left(x\right){\mathrm{\Lambda }}_{\mathrm{eff}}\left(x\right)}{E}},\end{array}$$

where P₀(x) is the unshielded surface production rate, corrected for latitude and air pressure (Lal, 1991). Equation (11) derives from integrating the path history of a particle being exhumed vertically at a steady rate E and emerging at the surface with an accumulated nuclide concentration C(x):

$$\begin{array}{}\text{(12)}& C\left(x\right)=P_{\mathrm{0}}\left(x\right){\int }_{t_{\mathrm{0}}}^{t_{\mathrm{surface}}}S(x,z(t\left)\right)\mathrm{d}t,\end{array}$$

which can be parameterized in terms of vertical depth below the surface, z, according to

$$\begin{array}{}\text{(13)}& C\left(x\right)={\displaystyle \frac{P_{\mathrm{0}}\left(x\right)}{E/\mathit{\rho }}}{\int }_{\mathrm{0}}^{z_{\mathrm{0}}}S(x,z)\mathrm{d}z,\end{array}$$

where the depth of a rock parcel below the surface z₀ at time t₀ is deep enough such that there is no cosmogenic nuclide production ($z_{\mathrm{0}}=\mathrm{40}\mathrm{\Lambda }/\mathit{\rho }$ for the calculations below) and $t_{\mathrm{surface}}=t_{\mathrm{0}}+\mathit{\rho }{z}_{\mathrm{0}}/E$ is the time it takes for a rock parcel to travel from depth z₀ to the surface (assuming a vertical exhumation pathway). Because there is no analytical solution for Eq. (13), the integral needs to be solved numerically. A total effective shielding factor, C_eff(x), acts as a correction factor to interpret local erosion rate from a sample concentration, defined by

$$\begin{array}{}\text{(14)}& C_{\mathrm{eff}}\left(x\right)={\displaystyle \frac{C_{\mathrm{shielded}}\left(x\right)}{C_{\mathrm{unshielded}}\left(x\right)}}={\displaystyle \frac{{\sum }_{z=\mathrm{0}}^{z=z_{\mathrm{0}}}S(x,z)}{{\sum }_{z=\mathrm{0}}^{z=z_{\mathrm{0}}}S{}^{\prime }(x,z)}},\end{array}$$

where ${\sum }_{z=\mathrm{0}}^{z=z_{\mathrm{0}}}S{}^{\prime }(x,z)$ is the integrated shielding depth profile for the case α=0 (i.e., no slope or skyline shielding), and C_eff(x) does not depend on spatial variations in latitude or air pressure corrections. Finally, a mean effective shielding factor, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$, is defined for the whole hillslope as

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}={\displaystyle \frac{\mathrm{1}}{L_{\mathrm{h}}}}{\sum }_{x=\mathrm{0}}^{x=L_h}{C}_{\mathrm{eff}}\left(x\right),\end{array}$$

which is equivalent to the catchment-mean shielding factor for the simplified valley geometry shown in Fig. 1.

Figure 4 Plots of (a) surface skyline-shielding factor (b) normalized effective vertical attenuation length, and (c) total effective shielding factor as a function of normalized distance from ridgeline for model runs with α=0–80^∘. Dashed line in (a) indicates topographic shielding calculation using algorithm of Codilean (2006) applied to a digital elevation model of the case α=80^∘.

Figure 5 Plot of normalized production rate relative to horizontal unshielded surface as a function of normalized vertical depth for a 60^∘ slope with no additional skyline shielding. Near the surface, production rates are decreased due to slope shielding of incoming cosmic radiation; however, production rates at depth increase relative to the unshielded case due to additional radiation along shorter oblique pathways (Fig. 1c).

3.3 Approximation for dipping catchments

Although the above framework accounts for variations in catchment relief and hillslope angle, α, in all cases there is no net dip to the entire catchment (i.e., ridgeline elevations are uniform), which is not the case for natural watersheds. To simplify the geometry of a dipping catchment, I use a similar approach as Binnie et al. (2006) to model the catchment as a plane fit through the bounding ridgelines with dip β. I focus on two end-member cases, using examples from the San Gabriel Mountains, California, USA, for illustration (Fig. 2). First, for an “interior” catchment that is tributary to a larger valley within a mountain range, the catchment will have a net shielding similar to the geometry of the hillslope in Fig. 1. Consequently, the shielding geometry can be approximated by Eqs. (6)–(9) with α=β. For the case of an “exterior” catchment that has a net dip β but no opposing skyline shielding, Eq. (7) becomes

$$\begin{array}{}\text{(16)}& \tan {\mathit{\theta }}_{\mathrm{0}}=\left\{\begin{array}{l}\mathrm{0},\mathrm{0}\le \mathit{\phi }<{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},\\ -\tan \mathit{\alpha }\cos \mathit{\theta }-{\displaystyle \frac{z}{x}}\cos \mathit{\phi },{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\le \mathit{\phi }\le \mathit{\pi }.\end{array}\right.\end{array}$$

For both examples, I compared the catchment-mean shielding factor, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$, to the mean surface skyline-shielding factor, $\stackrel{\mathrm{‾}}{S_{\mathrm{0}}}$, as calculated using the commonly applied topographic shielding algorithm of Codilean (2006) in TopoToolbox (Schwanghart and Scherler, 2014).

Figure 6 Plots showing mean total shielding factor, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$, for (a) simple horizontal catchment case (Fig. 1) for varying mean hillslope angle, α, which is equivalent to the “interior” dipping catchment case as a function of catchment dip, β (Fig. 2), and (b), the mean total shielding factor, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$, for the “exterior” dipping catchment model as a function of catchment dip, β (Fig. 2). Red points in (a) indicate the relationship between the mean surface skyline-shielding factor, ${\stackrel{\mathrm{‾}}{S}}_{\mathrm{0}}$, as a function of mean hillslope angle for compilation of catchment ¹⁰Be data in the Himalaya and eastern Tibet as reported by Scherler et al. (2017). Red and yellow squares indicate the mean surface skyline factor, ${\stackrel{\mathrm{‾}}{S}}_{\mathrm{0}}$, calculated for example catchments from the San Gabriel Mountains (Fig. 2). Arrows indicate the difference between mean surface skyline-shielding factor and mean total shielding factor, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$.

4 Model results

For the catchment geometry shown in Fig. 1, the local shielding factor, S(x,z), decreases with increasing depth, z; distance downslope, x; and increasing slope, α (Fig. 3). The surface skyline-shielding factor, S₀(x), decreases with distance downslope, x, and increasing hillslope angle, α, with the greatest shielding occurring in the valley bottoms of steep catchments (Fig. 4a). For the case α=80^∘, comparison of S₀(x) with the topographic shielding algorithm of Codilean (2006) shows that the two are equivalent (Fig. 4a).

The normalized effective attenuation length, Λ_eff∕Λ, decreases as a function of distance downslope and increases with increasing hillslope angle (Fig. 4b). Although cosmogenic nuclide production is concentrated at depths of $\mathit{\rho }z/\mathrm{\Lambda }=[\mathrm{0},\mathrm{3}]$ for low slopes, production rates at depth for very steep slopes can be greater than those of flat landscapes, despite lower surface production rates (Fig. 5). This effect emerges in part due to the increased effective attenuation length for collimated radiation in skyline-shielded samples (up to a factor of 1.3 – Dunne et al., 1999; Gosse and Phillips, 2001), but mainly because on steep slopes a point at vertical depth z below the surface is receiving incident radiation from oblique pathways that can be shorter than those overhead (Fig. 1c). Consequently, there is an additional radiation flux that increases the effective vertical mass attenuation length, Λ_eff, an effect that is most pronounced near ridgelines ($x/L_{\mathrm{h}}<\sim \mathrm{0.4})$ where skyline shielding is minimized (Figs. 3, 4b).

The combined effect of the decrease in surface production (Fig. 4a) and the increase in effective attenuation length (Fig. 4b) leads to a pattern whereby the total effective shielding factor, C_eff(x), is greater than 1 along the upper portion of hillslopes and less than 1 along the lower portion of hillslopes near the valley bottom (Fig. 4c). Although there may be considerable variation in shielding depending on slope position for steep slopes ($\mathit{\alpha }>\mathrm{60}{}^{\circ })$, the mean effective shielding parameter, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$, is unity for all cases (Fig. 6a).

For the case of dipping catchments (Fig. 2), the sensitivity of the mean effective shielding parameter to catchment dip, β, depends on whether catchments are “interior” (i.e., shielded by an opposing catchment) or “exterior” (i.e., no external skyline shielding). For “interior” catchments, the shielding calculations are identical to the analysis above, and thus ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$ is again unity for all cases (Fig. 6a). For “exterior” catchments, the increase in effective attenuation length at steep slopes due to shorter oblique radiation pathways (Fig. 1c) is larger than the decrease in surface production due to skyline shielding, and ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$ is greater than 1 (Fig. 6b). However, for all but the most extreme catchment dips (β≤40^∘), ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$ is effectively 1 (within 1 %).

For the two example catchments in the San Gabriel Mountains (Fig. 2), the mean total effective shielding factor, ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{eff}}$, is 1.00, despite steep catchment dips (β=17 and 32^∘) and high mean surface skyline shielding, $\stackrel{\mathrm{‾}}{S_{\mathrm{0}}}$ ($\stackrel{\mathrm{‾}}{S_{\mathrm{0}}}=\mathrm{0.87}$ and 0.84 as calculated by the algorithm of Codilean, 2006; Fig. 6a).

5 Implications for interpreting catchment erosion rates from ¹⁰Be concentrations in stream sediment

The above results indicate that no correction factor for topographic shielding is needed to infer catchment-mean erosion rates from ¹⁰Be concentrations in stream sediments for most cases, as long as the assumptions of spatially uniform quartz content and steady uniform erosion rate are valid. Only in the extreme case of an “exterior” catchment with mean dip β > 40^∘ will such corrections be necessary. Although the approach of only calculating the surface skyline-shielding component of the total effective shielding factor is appropriate for calculating surface exposure ages, neglecting the slope and shielding controls on the effective mass attenuation length leads to a systematic underprediction of the actual erosion rate. The magnitude of this underprediction increases with increasing catchment mean slope, as highlighted by a compilation of catchment erosion rates from steep catchments in the Himalaya and eastern Tibetan Plateau (red data points, Fig. 5a).

For catchments with spatially variable quartz content or erosion rate, a spatially distributed total effective shielding factor, C_eff, must be calculated at each pixel. Although calculating the surface skyline-shielding component is straightforward (Codilean, 2006), solving Eq. (3) at depth for arbitrary catchment geometries is presently too computationally intensive to be practical. However, while not entirely transferable to arbitrarily rough topography (e.g., Norton and Vanacker, 2009), Fig. 4c suggests that for slopes less than 40^∘, the total effective shielding factor does not vary significantly across the hillslope. For steep catchments with spatially variable quartz content or erosion rate, direct calculation of shielding at depth is likely needed to calculate the spatially distributed total effective shielding parameter. In particular, shielding calculations in landscapes dominated by cliff retreat are poorly suited for treatment in a vertical reference frame (e.g., Ward and Anderson, 2011).

The modeling approach above assumes a simplified angular distribution of cosmic radiation flux (Eq. 1) and only accounts for cosmogenic nuclide production via spallation. In actuality, the cosmic radiation flux does not go to zero at the horizon, and becomes increasingly collimated (higher m) with increasing atmospheric depth (Argento et al., 2015). Thus, the sensitivity of the effective mass attenuation length to shielding will increase with increasing elevation. However, the magnitude of changes in the effective mass attenuation length due to shielding-induced collimation is at most 30 % (Dunne et al., 1999), compared to the potential factor of 3 or more increase due to shorter oblique radiation pathways on very steep slopes (Figs. 1c; 4b). For hillslope gradients commonly observed in cosmogenic nuclide studies of steep landscapes (30–40^∘), the increase in effective mass attenuation length due to shielding-induced collimation and slope effects are 2 %–5 % and 6 %–15 %, respectively (Dunne et al., 1999; Fig. 4b). The dependence of Λ on atmospheric depth, which is typically not accounted for in catchment erosion studies, is minor (< 10 % for extreme case of catchment with 4 km of relief Marrero et al., 2016) compared to the above slope effect for most landscapes. Treatment of cosmogenic nuclide production by muons is less constrained than spallogenic production, but the angular distribution of production by muons is likely similar to that for spallation reactions and also sensitive to latitude and atmospheric depth (Heisinger et al., 2002a, b).

Overall, the effect of topographic shielding corrections on interpreting catchment erosion rates is small compared to typical assumptions inherent to detrital cosmogenic nuclide methods. In particular, the assumption of steady lowering is likely to be increasingly inappropriate for rapidly eroding landscapes characterized by a significant contribution of muonogenic production or slowly eroding landscapes where ¹⁰Be concentrations integrate over glacial–interglacial climate cycles. Steep landscapes characterized by stochastic mass wasting present additional complications (Niemi et al., 2005; Yanites et al., 2009), requiring the nontrivial calculation of spatially distributed shielding parameters for an arbitrary catchment geometry. Nonetheless, in all cases accounting only for surface skyline shielding (e.g., Codilean, 2006) without including its concurrent influence on the effective attenuation length yields incorrect results.

6 Conclusions

The simplified model presented here for catchment-scale topographic shielding of incoming cosmic radiation highlights the two competing effects of slope and skyline shielding. As catchment relief increases, surface production rates are reduced due to increased skyline shielding. However, for shielded samples radiation is increasingly collimated, and for sloped surfaces oblique radiation pathways increase nuclide production at depth. Both of these effects lead to deeper effective vertical mass attenuation lengths, which offset the reduction in surface production when inferring erosion rates from cosmogenic nuclide concentrations. At the catchment scale, the mean total effective shielding factor is 1 for a large range of catchment geometries, suggesting that topographic shielding corrections for catchment samples are generally not needed, and that applying commonly used topographic shielding algorithms leads to underestimation of true erosion rates by up to 20 %. Although these corrections are typically small compared to other methodological uncertainties, they vary systematically with slope and relief. Consequently, misapplication of shielding correction factors could influence interpretations of relationships between topography and erosion rate.