The sensitivity of cosmogenic radionuclide analysis to soil bulk density: Implications for soil formation rates

Improving our knowledge of soil formation is critical so that we can better understand the first‐order controls on soil thickness and more effectively inform land‐management decisions. Cosmogenic radionuclide analysis has allowed soil scientists to more accurately constrain the rates at which soils form from bedrock. In such analysis, the concentration of an isotope, such as Beryllium‐10, is measured from a sample of bedrock. Because this concentration is partly governed by the lowering of the bedrock‐soil interface, a cosmogenic depth‐profile model can be fitted to infer the bedrock and surface lowering rates compatible with the measured concentrations. Given that the bedrock‐soil interface is shielded by soil, the cosmic rays responsible for the in‐situ production of the radionuclide are attenuated, with attenuation rates dependent on the density profile of this soil. Many studies have assumed that soil bulk density is either equal to that of the bedrock or constant with depth. The failure to acknowledge the variations in soil bulk density means that cosmogenically derived soil formation rates previously published may be under‐ or overestimates. Here, we deploy a new model called “CoSOILcal” to a global compilation of cosmogenic analyses of soil formation and, by making use of estimated bulk density profiles, recalculate rates of soil formation to assess the sensitivity to this important parameter. We found that where a soil mantle >0.25 m overlies the soil‐bedrock interface, accounting for the soil bulk density profile brings about a significantly slower rate of soil formation than that previously published. Moreover, the impact of using bulk density profiles on cosmogenically derived soil formation rates increases as soil thickens. These findings call into question the accuracy of our existing soil formation knowledge and we suggest that future cosmogenic radionuclide analysis must consider the bulk density profile of the overlying soil.

Beryllium-10, is measured from a sample of bedrock. Because this concentration is partly governed by the lowering of the bedrock-soil interface, a cosmogenic depth-profile model can be fitted to infer the bedrock and surface lowering rates compatible with the measured concentrations. Given that the bedrock-soil interface is shielded by soil, the cosmic rays responsible for the insitu production of the radionuclide are attenuated, with attenuation rates dependent on the density profile of this soil. Many studies have assumed that soil bulk density is either equal to that of the bedrock or constant with depth. The failure to acknowledge the variations in soil bulk density means that cosmogenically derived soil formation rates previously published may be under-or overestimates. Here, we deploy a new model called "CoSOILcal" to a global compilation of cosmogenic analyses of soil formation and, by making use of estimated bulk density profiles, recalculate rates of soil formation to assess the sensitivity to this important parameter. We found that where a soil mantle >0.25 m overlies the soil-bedrock interface, accounting for the soil bulk density profile brings about a significantly slower rate of soil formation than that previously published. Moreover, the impact of using bulk density profiles on cosmogenically derived soil formation rates increases as soil thickens. These findings call into question the accuracy of our existing soil formation knowledge and we suggest that future cosmogenic radionuclide analysis must consider the bulk density profile of the overlying soil.

Highlights
• The effect of heterogeneities in soil bulk density on cosmogenically derived soil formation rates is unknown.
• Soil formation rates are recalculated using a new model to analyse the effect of density variations.
• Accounting for density in soils >0.25 m thickness brings about significantly slower soil formation rates.
• Measuring soil bulk density is essential when cosmogenically deriving soil formation rates.

K E Y W O R D S
Beryllium-10, bulk density, cosmogenic radionuclide analysis, CoSOILcal, soil formation 1 | INTRODUCTION How, where and why do soils form? These questions represent some of the oldest scholarly enquiries within soil science (Dokuchaev, 1879). Being able to identify the processes of, and the factors that influence, soil formation can help to inform our understanding of the soil system: its processes and functions, and the delivery of multiple ecosystem services. Given the diverse range of environments in which soils form, the subdiscipline of pedogenesis has a wide focus. An inexhaustive list of the types of enquiries undertaken by soil formation scholars includes the study of the accumulation and transformation of parent material (Hurni, 1983;Jenny, 1941;Minasny & McBratney, 1999;Simonson, 1959), the horizonisation of soil profiles (Bockheim & Gennadiyev, 2000;McAuliffe, 1994; and the factors that influence the evolution of soil properties (Richter, Oh, Fimmen, & Jackson, 2007;Schaetzl & Thompson, 2005;Vereecken et al., 2016). For this paper, we define "soil formation" here and hereafter as the process by which bedrock material converts into soil (Targulian & Krasilnikov, 2007;Egli, Dahms, & Norton, 2014).
One of the most important questions asked by soil scientists is: how fast does soil form? (Stockmann, Minasny, & McBratney, 2014). Knowledge of the first-order balance between rates of soil formation and erosion is integral if we are to ensure the long-term sustainability of global soil resources (Montgomery, 2007). Although measuring soil erosion is a long-established practice within soil science (Poesen, 2017;Quinton, Govers, Van Oost, & Bardgett, 2010), quantifying the rates at which soils form from bedrock has received less widespread attention (Duan, Shi, Li, Rong, & Fen, 2017;Schertz, 1983). Only within the past 20 to 30 years have technological advancements and interdisciplinary liaisons allowed soil scientists to more precisely constrain the rates at which soils form from bedrock (Heimsath, Dietrich, Nishiizumi, & Finkal, 1997;Román-Sánchez, Laguna, et al., 2019;Román-Sánchez, Reimann, Wallinga, & Vanwalleghem, 2019). By conducting analyses across a range of climatological and lithological contexts, it has also become possible to assess the extent to which the state factors of soil formationclimate, organisms, relief, parent material and time (Jenny, 1941)influence these soil formation rates (Stockmann et al., 2014).
The development of cosmogenic radionuclide analysis has demonstrated that soil thickness exerts a significant internal control on these state factors, and by extension, soil formation rates (Larsen et al., 2014). Many authors have observed that soil thickening leads to an exponential decline in soil formation rates . Thicker soil more effectively insulates the parent material against temperature and precipitation variations that drive weathering processes (Heimsath, Fink, & Hancock, 2009;Minasny & McBratney, 1999).
However, soil formation rates are not solely determined by this relationship with soil thickness. Moreover, Yu, Faybishenko, Hunt, and Ghanbarian (2017) show that rates of bedrock weathering are instead constrained by the transmission of water and solutes down the soil profile. A major determinant in the dynamics of this process is the bulk density of the soil, with greater bulk densities limiting the volume of water and solutes and slowing their infiltration to the bedrock (Gabet, Eldelman, & Langner, 2006).
Despite the fact that soil bulk density influences soil formation (Neely, DiBiase, Corbett, Bierman, & Caffee, 2019;Price & Velbel, 2003), precise density profiles are usually not integrated into the cosmogenic nuclide's production models. The cosmogenic nuclide concentrations in bedrock samples under the soil are fundamentally dependent upon two factors. One of these is the duration that the bedrock has been exposed to cosmic rays, as the cosmic bombardment of the minerals in the uppermost metres of bedrock produces these nuclides. Therefore, longer exposure times give rise to greater cosmogenic nuclide concentrations. The second factor is the evolution of the effective depth (lithostatic pressure) to the bedrock with time, which can be numerically related to the rate at which the bedrock weathers into soil (Lal, 1991;Stockmann et al., 2014). Given that the bedrock weathering rate is the desired dependent variable, concentrations of the radionuclide -N in Equation (1)can be measured using accelerated mass spectrometry (AMS) and interpolated to solve for bedrock weathering rates (ε): where P i is the annual production rate of the radionuclide by spallation; fast muons and stopping muons (sp, μ f and μ − ) at a surface with slope ϴ; z is the sample depth; ρ is the mean density of material overlying the sample; λ is the decay constant of the radionuclide; and Λ i are the mean attenuation lengths of the cosmic radiations (Lal, 1991). Cosmic rays are attenuated when they pass through the soil to reach the underlying bedrock. Accounting for the factors that drive this attenuation is critical so that accurate bedrock weathering rates can be determined. Two terms in Equation (1) directly address this attenuation: the depth of the sample and the density of the overburden material (in this case, the soil) (Balco, Stone, Lifton, & Dunai, 2008). Many studies (Heimsath et al., 1997;Owen et al., 2011;Riggins, Anderson, Anderson, & Tye, 2011) have assumed that the density of the soil is either equal to the bedrock density or is constant with depth (but see Larsen et al., 2014). This fails to acknowledge the heterogeneities of the soil profile and, in particular, the spatial variation in bulk density (Evans et al., 2019). As a result, all previous cosmogenic radionuclide analyses estimating soil formation rates for bedrock overlain by soil that have not measured and/or accounted for the variation in the soil density may have yielded data which are under-or overestimates.
Here, we amass an inventory of cosmogenically derived soil formation rates previously reported for bedrock underlain by soil, where spatial changes in soil bulk density have hitherto not been employed. Employing the CoSOILcal programme, the first of its kind that considers overburden density (Rodés & Evans, 2019), we aim to assess the sensitivity of these soil formation rates when bulk density data are accounted for.
From each of the shortlisted studies, raw data were extracted, including sample latitude, longitude and elevation, 10 Be concentration, the concentration uncertainty and the soil formation rate. The density assumed by the authors (which was generally either bedrock density or average soil density) was also recorded. Although the production rate of 10 Be is influenced by topographic obstructions, this can be addressed by calculating a shielding factor. This represents the ratio of the 10 Be production rate at the obstructed site to that at a site where the surface is flat and the horizon is clear (Balco et al., 2008). These shielding factors were recorded from each study. Some studies did not report data for all of the above criteria. As a result, the inventory was truncated to only analyse entries with a complete dataset. This resulted in the removal of 101 entries, permitting 163 for analysis.
The resulting inventory of soil formation rates used in this analysis was collated from 12 studies, representing 10 unique locations across Australia, the USA, Chile and the UK, five different climates (according to the Köppen classification system) and all three major rock types (n = 163; see Supplementary Information for the full dataset). The median depth for the inventory was 0.35 m, with 72% between 0 and 0.5 m, 25% between 0.5 and 1.0 m, and 12% > 1.0 m.
The variations in the bulk density of the soil above the bedrock were not provided in any of the inventoried studies or their accompanying supplementary information files. Therefore, in the absence of these bulk density data, fine earth bulk densities were estimated for five depths (0, 50, 100, 150 and 200 cm) down the soil profile at each site using the International Soil Reference Information Centre (ISRIC) Global Soil Information System "SoilGrids" (250 m resolution; June 2016 update; see Hengl et al., 2017). We acknowledge the fact that the ISRIC 250 m raster was not intended for this type of sitespecific analysis; a better approach would be to measure the bulk densities down the soil profiles at each site studied in the inventory. However, in the absence of bulk densities measured at the site scale, we use the ISRIC data here solely as a means by which to demonstrate the sensitivity of soil formation rates to bulk density.
For this sensitivity analysis, we employed the CoSOILcal model (Rodés & Evans, 2019) first applied in Evans et al. (2019). The main objective of the CoSOILcal model is to calculate a "best fit" bedrock lowering rate and its associated uncertainty at a site with known latitude, longitude, elevation and shielding using measured concentrations of in-situ cosmogenic radionuclides (and their uncertainties) at or below the bedrock-soil interface of known depth, taking into account the overlying soil bulk density (Equation 2).
A sequence of modelled bulk densities are generated and logarithmically distributed between 1 cm and 100 m, including those which are measured "in the field" (z,x). Densities shallower and deeper than those measured are extrapolated using the shallowest and deepest measurement, respectively. The remaining densities to be calculated are those that lie in between each measured density; these are linearly interpolated from the nearest neighbours. Second, a sequence of erosion rates (ε) are generated and logarithmically distributed between 1 cm y −1 and 100 m My −1 .
Next, the surface production rates (P) of the radionuclide are calculated, based on the inputted latitude, longitude and elevation data, as well as the apparent attenuation lengths of fast (Λ μf ) and stopping muons (Λ μ -) under the soil surface. The model then uses these surface production rate and attenuation data, as well as the generated bulk density profile (z,x) and the sequence of erosion rates (ε), to calculate the concentrations of the cosmogenic isotope at several depths (zs) down the soil profile for a given landscape age. The landscape age here refers to the time when the production and accumulation of the radionuclide began. In many cases this refers to the last major erosion event (e.g., the orogenic uplift of the bedrock in the study area or the erosion associated with glaciation). The landscape age is inputted by the user. Time (t) is discretized in a sequence of 100 values logarithmically distributed between 100 years and the landscape age. For each of these 100 timesteps, the model calculates an effective depth (x) by an interpolation of zs + ε Á t down the profile. The concentration of the cosmogenic isotope that accumulates during each timestep is then calculated using the following equation: where N is the cosmogenic isotope concentration accumulated during Δt timestep at the effective depth x; ε is the bedrock weathering rate; P is the annual production rate of the radionuclide by spallation, fast muons and stopping muons (sp, μ f and μ − ) at a surface with slope ϴ; λ is the decay constant of the radionuclide; Λ is the mean attenuation lengths of the cosmic radiations; and ρ is the density of overburden material for the time frame t -Δt to t. All of the modelled concentrations for the 100 time steps are then summed: where T is the landscape age. Because isotope concentrations are measured using AMS, bedrock weathering rates (ε) can be found by the simple interpolation of N. "Best fit" bedrock weathering rates can be computed by first calculating the deviation (s) of these modelled isotope concentrations (C) from those that were measured using AMS: where C is the modelled isotope concentration; M is the measured isotope concentration; and σ M is the uncertainty of the measured concentrations. Chi-squared values are then computed as the sum of the squared deviations. Where chi-squared values are smaller than the minimum chi-squared value plus the number of samples, the modelled bedrock weathering rates are considered to fit the data within a one-sigma confidence level. Further details about the model can be found in Rodés and Evans (2019). This model was run for each entry within the soil formation rate inventory (n = 163). Soil formation rates from both the original inventory and those recalculated using the CoSOILcal model were not normally distributed (the Anderson Darling test statistics were − 5.1 and − 12.2, respectively; p > .05 for both tests), so the Mann-Whitney U-test (a non-parametric statistical test for difference) was run. All statistical analyses were completed at 95% significance on a standard Excel workbook.
Employing the CoSOILcal model, the median rate of soil formation for the total inventory increased by 16% to 0.034 mm y −1 (range, 0.001-0.47 mm y −1 ) (Figure 1). However, this was not found to be statistically significant (p > 0.05). Moreover, 68% of the inventory reported instances where soil formation rates decreased after the adoption of the CoSOILcal model. Within this subset, the mean reduction in soil formation rates was 0.02 mm y −1 (range, 0-0.21 mm y −1 ). In just under a third of instances, soil formation rates increased after applying CoSOILcal, with the mean increase for this subset being 0.03 mm y −1 (range, 0.001-0.22 mm y −1 ). With respect to the effect of climate, the fastest rates were still associated with Warm Summer Mediterranean (Csb) contexts and, although the median had decreased by 0.02 mm y −1 in comparison to the original data, this was not found to be statistically significant (p > .05). On the contrary, the CoSOILcal output showed that the slowest rates were not associated with Arid Cold (Bwk) conditions, as was the case previously, but with Humid Subtropical (Cfa) climates (n = 26) (Figure 2; top panel). The median soil formation rate for the Humid Subtropical subset was 0.02 mm y −1 , representing a decrease of 0.006 mm y −1 in comparison to the original data, which was not statistically significant (p > .05). For Arid Cold climates, the CoSOILcal output reported an increase of 0.04 mm y −1 in the median soil formation rate, which was statistically significant (p < .05).
The CoSOILcal output had less of an impact on the lithological influence on soil formation (Figure 2; bottom panel). The fastest rates were again associated with metamorphic lithologies. For this subset, the median soil formation rate had increased by 0.01 mm y −1 but this was not found to be significant (p > .05). Likewise, the slowest rates were associated with igneous lithologies. For this subset, the CoSOILcal output reported a small decrease of 0.002 mm y −1 in the median soil formation rate, which again was not found to be significant (p > .05).

| DISCUSSION
The lack of a significant difference between soil formation rates prior to the application of CoSOILcal and those after the model was run may be explained by soil thickness. Moreover, the basis of the CoSOILcal model is that more of the heterogeneity in the density of the soil overlying the bedrock is accounted for. It is, therefore, reasonable to hypothesize that there is a threshold soil thickness under which the employment of these density data does not bring about a significant difference to soil formation rates. This was tested using the global inventory collated for this paper (Figure 2). It was found that the CoSOILcal model brought about a statistically significant difference in soil formation rates for soils >0.25 m F I G U R E 1 Soil formation rates from the global inventory previously published (x axis) and those calculated using the CoSOILcal model (y axis). The diagonal line represents y = x. The inset shows a zoomed projection of rates between 0 and 0.1 mm y −1 . The full dataset can be found in the Supplementary Information [Color figure can be viewed at wileyonlinelibrary.com] (p < .05). For those soils within the inventory >0.25 m (n = 95), the median soil formation rate derived using the CoSOILcal model was 1.2 times slower than the median of the rates previously published. Furthermore, as soil thickness increased, the difference between the CoSOILcal median and that of the original inventory also increased. For instance, for soils >0.4 m (n = 64) and > 0.5 m (n = 46) in thickness, the medians calculated using CoSOILcal were 1.3 and 1.7 times slower than those previously published, respectively.
Given that over half of the soil formation rates in this study's inventory are attributed to soils that are >0.25 m (n = 95), of which 48% of these are >0.5 m (n = 46), it calls into question the accuracy of these data. Moreover, it suggests that for these deeper soils, cosmogenically derived soil formation rates may be slower than we have previously estimated. This may have wider implications for some of the land-management decisions that have been based around these rates. For example, soil formation rates have been previously used to derive rates of soil loss tolerance. Although there are multiple methods of calculating soil loss tolerance values, one of those that has been popularly used is based on the premise that rates of soil erosion should be curtailed to those of soil formation (Di Stefano & Ferro, 2016). A review of the soil formation rates used hitherto to calculate soil loss tolerance values is beyond the scope of this study. However, our findings suggest that if the soil formation rates used stem from studies where the soil depth exceeded 0.25 m, the soil loss tolerance values may have been overestimated.
To address the potential inaccuracies of the existing soil formation rate data, the CoSOILcal model can be applied post hoc so that the measurement of radionuclide concentrations does not have to be repeated. However, this would require bulk density values being measured from the positions where sampling for cosmogenic radionuclide analysis took place. To ensure that these F I G U R E 2 Difference in the soil formation rate when the CoSOILcal model was employed in comparison to the original dataset (n = 163) for different climatic regions (top panel) and lithologies (bottom panel). The full dataset can be found in the Supplementary Information recalculated rates are as accurate as possible, we argue that the approach taken in this paper (that is, the use of the ISRIC 250 m raster for soil bulk density) should not be used to make these corrections. The use of the ISRIC data here was solely to demonstrate the sensitivity of soil formation rates to bulk density. Nevertheless, following the preliminary analyses presented in this paper, a global effort is now required to revisit our network of soil formation studies and recalculate rates using soil bulk density data measured at the site scale.
Furthermore, we argue that future attempts to derive soil formation rates using cosmogenic radionuclide analysis should encompass the measurement of bulk density down the soil profile and the use of these data in the CoSOILcal model when calculating bedrock lowering rates. Accurately quantifying soil formation rates is essential given that these are often used to guide policy decisions on soil conservation and erosion mitigation (Montgomery, 2007;Verheijen, Jones, Rickson, & Smith, 2009).
Although we have focused this paper on the implications of bulk density for soil formation rates, it is also important to acknowledge that these findings are equally impactful beyond soil science. Cosmogenic radionuclide analysis has been used to derive rates of bedrock weathering for a range of geomorphological discourses and Earth System models (Cockburn & Summerfield, 2004). These include calculating the long-term rates of landscape evolution (Heimsath et al., 1997) and quantifying the mobilization of bedrock-derived petrogenic carbon (Hemingway et al., 2018). The bulk density profiles of the unconsolidated material overlying the bedrock in these studies will have a similar effect on attenuating cosmic rays and should be accounted for. As a result, we argue that both our findings and the CoSOILcal model represent a significant contribution to multiple communities across Earth Sciences.

| CONCLUSION
We have demonstrated that applying higher resolution estimates of soil bulk density when using cosmogenic radionuclide analysis to calculate rates of soil formation is an important consideration. Applying the CoSOILcal model to a global inventory of previously published analyses, we found that soil formation rates, modelled from measured concentrations of 10 Be in the bedrock, were significantly different when high-resolution bulk density data for soils >0.25 m in thickness were applied. Furthermore, the impact of soil bulk density on cosmogenically derived rates increases with soil thickness. These findings highlight potentially important implications for the use of cosmogenic radionuclide analysis both within and beyond soil science. Not only does our work suggest that our existing soil formation rate inventory might require revisiting, but it also implies that future work that uses cosmogenic radionuclide analysis to derive soil formation rates should account for the bulk densities down the soil profile and employ the CoSOILcal model to calculate the rates of bedrock lowering. This is especially important given that soil formation rates have been, and continue to be, employed when making land-management policies and decisions.