Phosphorus (P) is an essential element for all life on Earth through its role in forming ATP and as a structural component of DNA (Nelson et al., 2008). Consequently, the P cycle in terrestrial and marine environments has been studied extensively (Filippelli, 2002;Paytan and McLaughlin, 2007;Turner and Condron, 2013;Walker and Syers, 1976). Ecological research has shown that P fertility of terrestrial ecosystems is strongly linked to the weathering trajectory of soils with time: on geomorphically stable landforms, increasingly chemically altered soils lead to a declining pool of plant-available P, which can cause a decline of primary production and biomass, and strongly influence species and functional diversity (Crews et al., 1995;Eger et al., 2013b;Peltzer et al., 2010;Zemunik et al., 2015). The depletion of plantavailable P, however, is not simply a result of P weathering loss but also due to intensive biochemical transformations and recycling (Frossard et al., 2000).
Our current understanding of long-term P transformations is largely based on soil chronosequence studies; a study concept that takes advantage of a set of landforms that formed at different but known times in the past that have been minimally rejuvenated by erosion or deposition. In this framework, all other soil forming factors since cessation of erosion or deposition are assumed to have been similar between sites, allowing for isolation of the influence of time on soil development. Synthesising multiple soil chronosequences in New Zealand, Walker and Syers (1976) established the seminal soil P development concept: with increasing time, bio-available P declines as a result of leaching and the transformation of primary, rock-derived apatite P into less directly bio-available P forms such as organic P and P adsorbed to or occluded into secondary oxides. Whereas apatite P can be made directly bio-available as PO 4
3-through mineral dissolution in an acidic soil environment, the physically occluded P fraction, in particular, comprises P forms that are highly stabilized (Smeck, 1985) and hence not readily accessible by biota as a result of physical protection in mineral structures (primary or secondary silicate minerals, oxides, oxyhydroxides), organic matter and soil micro-aggregates (Blake et al., 2003;Guo and Yost, 1998). The Walker and Syers paradigm of P development has been found to be generally valid for a range of soils in different climatic and lithologic settings (Crews et al., 1995;Eger et al., 2011;Selmants and Hart, 2010;Turner and Laliberté, 2015). However, the nominally non-eroding setting of a chronosequence is a special case, as most of Earth's surface undergoes either net erosion (Larsen et al., 2014b) or deposition. Hillslopes are predominantly erosional landforms, where gravity and physical disturbances facilitated by water or bioturbation drive the downslope movement of soil, which is then delivered to fluvial systems or deposited on convergent sections of slopes or at slope-valley transitions. As mass is physically and chemically lost from a soil profile on an eroding hillslope, soil cover is maintained over time by the counterbalancing process of soil production (Gilbert, 1877;Heimsath et al., 1997), the conversion of parent material to soil. Soil production is regarded as a natural rejuvenator of soil nutrients by the replacement of weathered, nutrient-poor material with unweathered substrate (Amundson et al., 2015;Porder and Hilley, 2011;Porder et al., 2007b;Vitousek et al., 2003). The 'fertilisation' through soil production on slopes could be especially significant for soil P because in most terrestrial settings P is supplied to the biogeochemical cycle by weathering of the P-bearing mineral apatite and hence is delivered to the base of the soil by the parent material, unless there are external sources of P, such as atmospheric input. Dust has a major impact on soil P budgets in sufficiently P-depleted soils and/or where dust deposition rates are high (e.g., Chadwick et al., 1999;Eger et al., 2013a). Atmospheric input may even play an important role in P cycling at younger stages of ecosystem development in some locations (Arvin et al., 2017;Boyle et al., 2013).
The role of hillslope topography and soil erosion processes need to be considered when evaluating soil P pools and fractionation as it will affect the time soil material is residing on the slope before removal by chemical or physical processes (Agbenin and Tiessen, 1994;Amundson et al., 2015;Porder and Hilley, 2011;Porder et al., 2007b;Vitousek et al., 2003). For example, in Hawaii lower proportions of occluded P but more organic P were found on a hillslope in comparison to the geomorphically stable shield surface, indicating rejuvenation via slope dynamics (erosion and deposition) (Vitousek et al., 2003). However, no clear trends of P fractionation existed across the hillslope itself, from the shoulder (younger soils) to the toeslope (older soils). P fractionation data from ridge-slope-valley transects in Puerto Rico demonstrated the dominant control on the spatial distribution of more labile P forms was topography; labile P was lowest on the ridge and generally increased downslope toward the valley (Mage and Porder, 2013). In contrast, parent material was the main control on occluded and total P, with the highest values in the valleys, and apatite P (< 5% of total P in all soils) was unrelated to either topography or parent material (Mage and Porder, 2013). Selected soil P fractions (total P, apatite P, labile P and occluded P at 0-20 cm depth) on ridgetops in Puerto Rico were not significantly controlled by erosion rates or soil residence time, however, erosion rates and residence times varied little between sites (McClintock et al., 2015). Data from slope transects in Brazil showed that young upper slope soils (Entisols) have higher apatite P and lower labile P concentrations than Inceptisols in mid and lower slope positions (Agbenin and Tiessen, 1994). Differences in relative soil residence times induced by erosion were deemed the likely reason for the behaviour of apatite P. With only the study from Brazil adhering to the Pdevelopment concept derived from chronosequences, the relationship between P fractions and the relative soil age on slopes is less clear.
The divergence in P fractionation on eroding slopes relative to what is predicted from chronosequence studies highlights the need to reconcile the apparently different behaviour of P observed in different topographic settings. We suggest that comparing these findings in the context of soil P evolution as proposed by Walker and Syers (1976) is the most promising approach. Amundson et al. (2015) proposed a unifying concept in which temporal shifts from N to P nutrient limitation in terrestrial ecosystems are related to the continuum of residence times of minerals within the soil. The concept of Amundson et al. (2015) builds on new appreciation of tectonic uplift as a driver of erosion and thus P supply in the otherwise P-depleted tropical soils (Porder et al., 2007b). Uplift is typically associated with tectonic plate margins and a major control of erosion rates that are inversely related to soil residence times. Soil residence time in these studies is defined as the length of time that is required for soil material to be removed by erosion and replaced by soil production, and during which soil particles experience physical and biogeochemical conditions at the top of the weathering profile (Almond et al., 2007;Dere et al., 2013;McClintock et al., 2015). Compared to chronosequences developed in flat landforms, Amundson et al. (2015) suggested that residence times for most hillslope soils in temperate climates give rise to neither N nor P limitation. In other words, soils on eroding hillslopes are not too young to have N limitation or too old to be depleted in mineral P.
Whether eroding hillslope soils indeed occupy an optimal residence time window with respect to P limitation remains to be tested. There are few data that directly link individual P fractions to absolute soil residence times (McClintock et al., 2015). Additionally, previous studies of soil P on eroding hillslopes are largely limited to tropical landscapes (Abekoe and Tiessen, 1998;Agbenin and Tiessen, 1994;Araújo et al., 2004;Mage and Porder, 2013;McClintock et al., 2015;Porder and Hilley, 2011;Porder et al., 2007b;Vitousek et al., 2003). In these actively eroding tropical systems, deep chemical alteration of bedrock causes soils to be depleted in apatite P, which provides the first indication that the optimal window hypothesis may not be applicable globally. However, the applicability of these studies from tropical landscapes to extra-tropical regions may also be limited. In contrast to temperate climate regions, in the tropics, deep and more completely weathered profiles prevail, mineralisation rates of organic matter are higher, low-reactivity clays and pedogenic oxide/hydroxides increasingly dominate the residual soils, and the legacy of glacial/periglacial conditions during the Pleistocene is largely absent.
Here we present new P fractionation data quantitatively linked to hillslope soil residence times across two gradients of erosion rates in temperate ecosystems and compare them to published results regarding patterns and rates of P transformation. We initially hypothesised, based on the proposal by Amundson et al. (2015), that higher soil production and erosion rates and hence shorter residence times result in high total soil P concentrations and high proportions of primary mineral P as expected for immature soils, whereas lower erosion rates and longer residence times result in low total soil P due to the intensive weathering of older soil particles, and a high proportion of secondary P forms as expected in more mature soils. However, our data do not support this hypothesis and instead, somewhat distinct from the conceptual framework laid out in Walker and Syers (1976), highlight the significance of weathering below the base of the soil in temperate climates, biological uptake of P and potential dust accretion.
We first require a consistent framework for the measure of time for our soils. As we will show, soil residence time and soil age provide consistent temporal references to which soil P dynamics from geomorphically active and stable landscapes can be compared. We conceptualize that the mass balance of a hillslope soil (Fig. 1) is largely determined by the difference between the mass losses via physical and chemical erosion and mass input via soil production; our conceptualization assumes aeolian inputs are negligible. In this simplified view, a hillslope soil is defined as a part of a weathering profile that has been not only chemically weathered but also physically disturbed and mixed. In contrast, the lower part of the weathering profile is considered chemically weathered but physically undisturbed (i.e., saprolite). We assume soil mass removal via physical and chemical erosion is balanced by soil production, such that a steady-state is reached (Heimsath et al., 1997). Since our focus is on transformations and losses of P in hillslope soils (i.e., soil as the residual of the weathering process), we are concerned with the ages of the particles with respect to their initial incorporation within the soil, as exposure to weathering and leaching increases as a function of particle age. Hillslope soil particles have a distribution of age that is unknowable in all but the simplest case of steady state soil thickness, together with either complete mixing or plug flow (Mudd and Yoo, 2010) and the absence of chemical weathering. Hence, pragmatically, we seek only a metric to rank soils according their exposure to weathering. Under aforementioned simplifying conditions, mean particle age of the soil (ψ), mean soil residence time (χ, average age of particles leaving the soil) and soil turnover time (φ, the length of time that it takes for a soil particle to be completely depleted by the outgoing flux) are equal (Almond et al., 2007;Mudd and Yoo, 2010;Yoo and Mudd, 2008). We adopt the mean particle age, and, assuming perfect mixing and steady state in the absence of chemical weathering, we estimate it by the soil turnover time. Soil turnover time is calculated as the ratio of the mass of the soil and the outgoing mass flux from that soil (Mudd and Yoo, 2010). Assuming steady state, the outgoing mass flux (erosion) equals the rate of conversion of parent material to soil (i.e., soil production rate) as determined by cosmogenic nuclide measurements at each of our sites (see below) corrected for chemical mass loss.
where ψ is mean particle age, φ residence time and χ is turnover time (T).
[Zr] represents the mass concentration of the immobile element zirconium (MM -1 ), ρ is bulk density (ML -3 ), h is soil thickness (L), D is soil production/erosion rate (LT -1 ), and subscripts s and r indicate soil and parent material, respectively. The term [Zr] s /[Zr] r converts the soil erosion rate, which includes a chemical weathering component, into a physical erosion rate (e.g., Riebe et al., 2003). Since soil thickness and soil production rate units are given in length, the inclusion of ρ s /ρ r accounts for the dilation between parent material and soil. Here [Zr] s / [Zr] r is typically larger than 1 because of Zr enrichment in soils as a result of leaching of other more soluble elements. In contrast, the bulk density ratio between soils and the parent material is typically < 1, contributing to cancelling the effect of Zr enrichment in soils. Thus we further simplify our soil particle age metric to h/D, similar to other studies (e.g., Amundson et al., 2015;Porder et al., 2007b). In the literature, mean soil particle age, mean soil residence time, or soil turnover time have been used interchangeably or authors simply referred to soil residence time without strict definitions based on reservoir theory (e.g., Almond et al., 2007;Amundson et al., 2015;Green et al., 2006;Porder and Hilley, 2011). We follow this convention and use the term soil residence time instead of the mean particle age or turnover time.
While an approximate steady state is a useful concept for investigating eroding soils, soils developing on geomorphically stable landforms are only minimally affected by physical erosion. Still, a mean age of soil mineral particles can be defined (Yoo and Mudd, 2008). A soil consists of mineral particles that have a range of time lengths (i.e., ages) since their physical incorporation into the soil from the underlying parent materials. The maximum age of mineral grains cannot be older than the age of the soil, however. Soil age is defined as time length since cessation of erosion or deposition. Additionally, in noneroding chronosequences, soils become increasingly thicker with time as chemically more inert soil material residually accumulates, which slows the downward propagation of soil development into the parent material (Lebedeva et al., 2010), and hence the rate of incorporation of nutrient-replenishing parent material (Yoo and Mudd, 2008). Thus, the number of mineral grains introduced to the soil from the underlying parent material exponentially decreases over time. As a consequence, it is expected that the mean age of the mineral particles is less than soil age, but the distribution of individual mineral grains' ages is skewed toward the early phase of soil formation. Thus, the mean age of mineral grains in a soil is proportional to the soil age (Yoo and Mudd, 2008).
The P data come from two temperate locations that differ substantially in rainfall and soil production/erosion rates. Soil thicknesses, soil production/erosion rates, and calculated soil residence times (Eq. ( 1)) for each soil are reported in Tables 1 and 2. Soil thicknesses and soil production/erosion rates were reported previously and the methods and discussion concerning these data, and the range of parameters are described by Larsen et al. (2014a) for the Western Southern Alps (WSA) sites, and Hurst et al. (2012) and Yoo et al. (2011) for the Feather River (FR) sites. Soil residence times at FR sites were also reported in Wang et al. (2018).
The first study area is located in the western Southern Alps (WSA) of New Zealand at the collisional boundary of the Australian and Pacific Plates (Fig. 2), resulting in up to 10 mm y -1 tectonic uplift (Little et al., 2005;Tippett and Kamp, 1993). Soil parent material is schist derived from a greywacke protolith. The Southern Alps form an orographic barrier against the prevailing westerly airstream resulting in a mean annual precipitation of 10,391 mm (1979-2015, maximum Dec 1099 mm, minimum July 643 mm), with a mean annual temperature of 5.5 °C (NIWA, 2016;Tonkin and Basher, 2001) at ~900 m asl. The natural vegetation cover is a podocarp-hardwood forest and subalpine, dense scrub/low tree communities (Wardle, 1977). Topography is heavily dissected by a dense drainage network of steep, V-shaped valleys including waterfalls, gorges, and narrow ridge lines (Whitehouse, 1988). Landslides are frequent as a result of earthquakes and high rainfall (Hilton et al., 2008;Hovius et al., 1997;Korup et al., 2004), but return intervals are long enough to allow the formation of thin soil and regolith cover at any point on the landscape between failures (Larsen et al., 2014a;Whitehouse, 1988).
The soil production/erosion rates in the WSA (Table 1) are amongst the highest in the world (Larsen et al., 2014a) and the soils are weakly developed Entisols or Inceptisols (Soil Survey Staff, 2014). All individual soil sampling sites were located on the main ridges or in local, meter-scale convexities on smaller divides emanating from the main ridges to avoid effects from landsliding. As such the site selection aimed at fulfilling the steady state assumption required by the in-situ cosmogenic nuclide method to yield reliable soil production rates at each site. We do not necessarily expect these sites to be representative of the average soil thickness in each of the WSA catchments. Local slope at the soil sites ranged between 24°to 50°. The second study area is in the Feather River catchment (FR) in the Northern Sierra Nevada of California, USA (Fig. 2). The FR site is within the lower reaches of the Middle Fork Feather River, where mean annual precipitation is 1750 mm and the mean annual temperature 12.5 °C (PRISM Climate Group, www.prism.oregonstate.edu). The bedrock at the study site is granodiorite, but the adjacent area features a complex intrusion of granitoid plutons into metamorphic and ophilitic rocks (Saucedo and Wagner, 1992). Erosion rates vary with topography, with lower erosion rates of 20-40 mm ky -1 for a relatively flat relict upland surface and high erosion rates of 200-250 mm ky -1 on the steep slopes draining to the deeply incised canyon of the Feather River (Hurst et al., 2012;Riebe et al., 2000;Wakabayashi and Sawyer, 2001).
The FR study sites are located within the Bald Rock tributary basin that descends from a relict surface (850 m asl) to the Middle Fork Feather River (310 m asl). Spatially detailed rainfall data are lacking in the region. However, the region's precipitation map (Western Regional Climate Center, https://wrcc.dri.edu/Climate/maps.php) suggests that the elevation difference within the tributary basin causes only ~15% of variation in the annual mean precipitation. Relatively constant climate within the basin is also reflected in homogenous presence of mixed conifer forest (Milodowski et al., 2014). The overall slope gradients of the FR sites within the tributary basin increase from approximately 15°t o 31°toward the Middle Fork Feather River. A knick-point, which has been initiated by the incision of the Middle Fork Feather River, has been migrating upward through the tributary basin (Attal et al., 2015). Our sites comprise three eroding hillslope transects: POMD is near a low relief plateau and located above the knick-point, BRC is below the knick-point, and FTA is between the knick-point and the plateau. According to Hurst et al. (2012), catchment scale erosion rates adequately represent the spatial variability of erosion rates within the tributary basins and vary from 35.7 mm ky -1 at POMD to 250 mm ky -1 at BRC, with intermediate rates at the FTA sites. Since the soil thicknesses of the FTA soils do not differ significantly from those of POMD and BRC, they will have soil residence times that are between those of POMD and BRC (Table 2). Consistent with the range of residence times, all soils are Inceptisols. The BRC soils with highest erosion rates have substantially more coarse grain sizes and are more heterogeneous in their thicknesses as compared to POMD and FTA (Wang et al., 2018). Unlike POMD and FTA, which have continuous soil cover, BRC is also characterised by patchy bedrock outcrops (Milodowski et al., 2015). Though a generally negative relationship between soil thickness and erosion rate has been observed at an adjacent ridge line (Gabet et al., 2015), within the Bald Rock Basin soil thickness is relatively insensitive to erosion rate (Yoo et al., 2011).
With respect to our conceptual framework of hillslope soils (Fig. 1), we define soil as the sum of pedogenic A and B horizons (Soil Survey Staff, 2014). Our field observations clearly indicated the effects of physical disturbance by tree roots and tree throws in mixing these horizons, qualifying the sum of A and B horizons as the mobile soil. The zone of chemical weathering between the soil and fresh bedrock, which is termed saprolite (Fig. 1), is characterised in our study areas by wellpreserved rock fabric indicative of minimal physical disturbance. At WSA, the thin (< 0.5 m) saprolite zone can be at times better described as R horizon that comprises in-situ (i.e., physically connected to bedrock below), minimally weathered bedrock, and mm-sized cracks containing material from the overlying B horizon. At the WSA sites we took bulk-samples (combined A and B horizons) of each soil profile (one sample per site, total 22 samples from 22 soil sites). We took great care to obtain bulk samples that represented the true proportions of each soil horizon in the soils (i.e., no preferential sampling of one horizon) by cutting back the profile face with a spade over the entire depth of the soil and collecting the cut-back material. At the FR site each hillslope (POMD, FTA, BRC) was sampled at the summit, shoulder, and backslope for soil and saprolite material (convex to straight slopes). Each soil pit was excavated to the depth of 20-30 cm below the soil-saprolite boundary and soil samples were collected by horizons and depth intervals. Because little differences in soil geochemistry and morphology were observed as a function of topographic locations within each hillslope group (Yoo et al., 2011), our detailed P fractionation measurements were limited to two soil profiles from each hillslope.
We primarily present P data of soil samples and saprolite (FR only) for the following fractions: total P (P total ), the primary, apatite-derived P fraction (P apatite ), organic P as the organically bound P (P org ), the nonoccluded, iron and aluminium oxide-bound P (P Fe/Al ), and the occluded/recalcitrant/residual P fraction (P occ ). The extraction procedures differed between the WSA and FR sites due to the dates when the analyses were conducted (WSA in 2013, FR in 2016). For the WSA sites, P total was extracted by NaOH fusion in nickel crucibles (Blakemore et al., 1987;Smith and Bain, 1982), and the extracts analysed following Murphy and Riley (1962). P org was extracted following the ignition method of Saunders and Williams (1955). The modified Hedley sequential fractionation with 0.1 M NaOH and 1 M HCl after Tiessen and Moir (1993) yielded inorganic P Fe/Al and P apatite , respectively. P org , P Fe/ Al , and P apatite extracts were quantified also following Murphy and Riley (1962). The difference between total P and the sum of P org , P apatite , and P Fe/Al is regarded as the occluded/recalcitrant/residual P (P occ ). We note that our P occ fraction does not discriminate between inorganic and organic P occ .
The FR samples underwent a more detailed fractionation than the WSA samples following the scheme by Condron et al. (1996). This scheme involves a sequential extraction of 6 consecutive steps on the same soil sample: 1) extraction of labile inorganic P with 1 M ammonium chloride (Pi NH4Cl ); 2) inorganic and organic P (Pi bic and Po bic ) with 0.5 M sodium bicarbonate (NaHCO 3 at pH 8.5); 3) inorganic and organic P (Pi OH _I and Po OH _I) with 0.1 M NaOH; 4) P apatite (Pi HCl ) with 1 M HCl; 5) a second extraction with 0.1 M NaOH (Pi OH _II and Po OH _II); and a final digestion with concentrated H 2 SO 4 and 30% H 2 O 2 to yield the residual P (Olsen and Sommers, 1982). The inorganic P concentration in acid extracts was quantified following Murphy and Riley (1962). Inorganic P in alkaline extracts followed Dick and Tabatabai (1977), whereas the organic P was obtained by the difference between the inorganic P and total P concentrations after digestion with ammonium persulfate and H 2 SO 4 in an autoclave.
To allow for comparability, the FR P fractions were combined to be equivalent to the WSA fractions: organic P (P org ) is the sum of Po bic , Po OH_I and Po OH_II ; apatite P (P apatite ) is equivalent to Pi HCl ; P Fe/Al is Pi OH_I ; and occluded P (P occ ) equals the sum of Pi OH_II and P residual (here we account for the fact that the simpler P fractionation of the WSA samples does not include a second NaOH extraction that was performed on the Feather River samples). P data are only reported for the mineral horizons (Table 1 for WSA and Table 2 for FR). At both sites we measured pH of the bulk samples (WSA) and selected depth increments (FR) at a soil/water mass ratio of 1:2.5. To compare results between the two sites and previous work, we focus on the ratios of P apatite , P occ , P org , and P Fe/Al to P total rather than absolute P concentrations, as this approach allows comparison of sites with varying concentrations of P in the parent material (Hahm et al., 2014;Mage and Porder, 2013;Porder and Ramachandran, 2013). Regression analysis and derivation of regression model parameters was conducted using R (R Core Team, 2017).
Soils are very acidic (Table 1) with pH values as low as 3.2, similar to other published data from the region (e.g., Almond and Tonkin, 1999;Stevens, 1968;Tonkin and Basher, 2001). Since these are composite values representing the substrate over the entire depth of each soil, the values are likely to be lower for the topsoils and higher for the subsoils alone. Secondary P is the predominant form of P in the WSA samples (83-97%), whereas apatite P remains between 3% to 17% of P total over the entire range of soil residence times (Fig. 3, Table 1). The P apatite /P total ratio is weakly inversely correlated with soil residence time (χ) (P apatite /P total = 0.1144 -0.0002χ , R 2 = 0.18, p = 0.045). In contrast, neither P org nor P mainly associated with pedogenic oxides (P occ and P Fe/Al ) are statistically significantly correlated with residence time.
Soils at FR are slightly acidic throughout. At the FR sites (Fig. 3), depth-weighted contributions of P fractions to P total of each soil show P apatite is always < 3% of P total , whereas P org and P occ are clearly dominant. The high proportions of P org and P occ are present across the range of soil residence times (Fig. 3). Only P org decreases slightly with increasing soil residence time (χ) following a statistically significant power-law model (P org /P total = 0.76χ -0.113 R 2 = 0.89, p < 0.005). Since P fractions are measured in consecutive depth intervals and soils are deeper at the FR sites than the WSA sites, a more detailed picture of the P fractions across soil depths was obtained (Fig. 4). Despite the rudimentary morphological development of the Inceptisols, there are major depth gradients in P chemistry. Nearly all P fractions and P total have highest concentrations in the topsoil (Fig. 4, Table 2). High topsoil concentrations are most strongly expressed in the more bio-available secondary forms of P (NH 4 Cl, NaHCO 3 , and first NaOH extractions) but also P apatite . The decline for most P fractions with depth continues beyond the soil and reaches deep (> 200 cm) into the saprolite (Cr) with the exception of P apatite concentrations, which increase again in concentration at depths > 150 cm.
Soil residence times vary by up to four orders of magnitude but residence times exert little control on the concentrations of total P; P apatite /P total ratios remain low (< 18% at WSA, < 3% at FR) for all soils (Tables 1 and 2). The statistically significant model linking soil residence time to P apatite /P total , albeit statistically weak, indicates that the paradigm of P apatite loss with increasing soil development time cannot be rejected at least for WSA. However, at WSA the proportion of P apatite to P total at the shortest soil residence time is < 18% in contrast to the classic Walker & Syers paradigm that predicts the dominance of P apatite over secondary P forms in such young soils with rudimentary profile morphology (Entisols, Inceptisols).
In Fig. 5, we compared our own ratios of P apatite /P total as a function of soil residence time to the published ratios from soil chronosequence studies. It appears that soils at both of our study sites, despite not being morphologically developed beyond Inceptisols, have already reached the late stage of soil P development with a very low and largely invariant proportion of P apatite and very high secondary P forms typical for older chronosequence sites (Fig. 5). For instance, contrasting the WSA sites against the Spodosols developed on the nearby Franz Josef chronosequence shows that the 10% average proportion of P apatite in the WSA soils is at best similar to the top 30 cm (to stay within our range of WSA soil depths) of > 1000-to 5000-year-old soils of the Franz Josef chronosequence (Stevens, 1968). Additionally, the P occ /P total ratios of the hillslope soils at WSA are so high that they are only replicated at the 120,000 y-old, retrogressive stage of the Franz Josef chronosequence (data from Stevens, 1968).
Contrasting the FR sites to the Merced River chronosequence (Harden, 1987), developed on granite-derived alluvium and located just west of the Sierra Nevada, reveals that, like the WSA, the P chemistry of hillslope soils is comparable to that of old soils. We note here that the Merced River sites developed in a dryer climate than the FR sites (mean annual precipitation: 300 mm). Although P fractionation data are not available for the Merced River chronosequence, the site's apatite concentrations can serve as a proxy for the depletion of primary mineral P (P apatite ) (Harden, 1987). Apatite concentrations in Merced River soils decrease 10-fold within the first 40 ky of soil formation with little change thereafter (> 40 ky to 600 ky). The initially rapid decline of apatite observed at the Merced River chronosequence is similar to the trend in P apatite depletion at Franz Josef and other chronosequences (Fig. 5). Comparing the FR sites to Merced River chronosequence, the low and invariant contributions of P apatite to P total at the FR site signals that FR soils have already reached that stage of severe apatite depletion only observed in Merced River soils older than 40 ky that exhibit much greater morphological maturity and chemical differentiation (e.g. layers of illuvial clay-enrichment in the soil) than the FR hillslope soils.
The only other published P fractionation data in Fig. 5 from eroding hillslopes with temporal data are from Puerto Rico (McClintock et al., 2015). McClintock et al. (2015) reported soil residence time for the top 20 cm of soils and every sample contained < 5% of P apatite (their HCl- The FTA soil residence times are used for illustration purposes are derived from the average of the erosion rates of 35.7 mm ky -1 (POMD) and 250 mm ky -1 (BRC). Soil residence time has little influence on P fractions and secondary P fractions clearly dominate over P apatite at all times. extractable P). P occ (residual-P), P org (NaHCO 3 -P o + NaOH-P o ) and P Fe/ Al (NaOH-P i ) contribute on average 55%, 30%, and 13% to P total , respectively. The residence times of their soils are comparable to those we studied (Fig. 5), and there is also similarity to our average P inventory that shows the contributions to P total of P occ > P org = P Fe/Al for WSA and P occ > P org > P Fe/Al for FR (Fig. 3). These similar patterns in soil P persist despite the large climatic difference between WSA, FR and Puerto Rico (Fig. 6), indicating that climate is not a driver of such patterns observed in P fractions.
Our results differ from the patterns reported for soil chronosequence studies (Fig. 5) but are consistent with soil P fractionation studies conducted on eroding hillslopes underlain by crystalline bedrock (Fig. 7). These published data are compiled from studies that explicitly describe soil sampling locations on hillslopes and from in-situ soil parent materials (or local regolith). Most sites in Fig. 7 are upslope locations where soil production from the underlying bedrock maintains soil cover and colluvial deposition is limited. All soils from eroding hillslopes in the published literature, despite their presumably short soil residence times due to erosion, have very low percentages of P apatite (Fig. 7). Eroding soils also contain a high proportion of P occ relative to P total (Abekoe and Tiessen, 1998;Araújo et al., 2004;Homyak et al., 2014;Mage and Porder, 2013;McClintock et al., 2015;Vitousek et al., 2003), consistent with results from the WSA and FR sites.
In Fig. 7, the soils from semi-arid northern Brazil (Agbenin and Tiessen, 1994) are an exception to the low contribution of P apatite to P total in eroding soils. The P apatite contribution to P total in A horizons of these upslope soils reaches 60 ± 18% in shallow Entisols, but that quickly decreases to 17 ± 18% lower on the slope where thick depositional Inceptisols are found. However, this northern Brazilian hillslope (Agbenin and Tiessen, 1994) is underlain by apatite-rich syenite that is unusual for the region (Araújo et al., 2004). In summary, excluding the study site underlain by apatite-rich syenite, no previous work from eroding hillslopes we examined documents P apatite contributing > 30% of P total .
Despite the extreme rainfall rates, it seems unlikely that the strong depletion of P apatite at WSA is simply a reflection of the high rainfall in accelerating the weathering and transformation of rock/soil P. Fig. 7 shows the range of precipitation from the published hillslope P studies including those from this study. It is clear that the exceedingly low contribution of apatite to total P on eroding hillslopes is not limited to regions of high rainfall but rather is a norm across a wide range of precipitation rates. Although they did not measure phosphorus, Dixon et al. ( 2009) linked erosion rates and chemical weathering state of soils and saprolite. They found that the chemical weathering state of the saprolite determined the chemical weathering of the soil: when the saprolite was highly weathered, additional weathering in the soil was low, whereas when the saprolite was only weakly chemically altered, the contribution of soil weathering to the overall chemical weathering of the weathering column would be high. Given that Dixon et al. (2009) found a strong relationship between erosion rate and weathering rate of the saprolite (but not the soil), we expected that the P inventory of the saprolite at FR would respond to erosion rates. However, this (i.e., higher erosion rates/lower residence times = less strongly depleted apatite P in the saprolite) does not seem to be the case at FR (see saprolite samples in Table 2).
One potential contribution to the P depth profiles we see (e.g., the increase in total P from saprolite to soil observed at FR; Fig. 4) is dust deposition. Substantial dust deposition is highly unlikely for the WSA sites, as studies have shown that even in favourable conditions of local dust mobilisation (e.g. close to an unvegetated braided river in the coastal plain) any effect of dust on soil P is limited to areas close to the dust source (< 2 km) (Eger et al., 2013a). There is no local dust-producing source in the vicinity of our WSA sites and long-range deposition from Australia (Holocene dust deposition rate 0.6 g m -2 y -1 ; Marx et al., 2009) will have little impact at such high erosion rates (lowest rate of all WSA soils 307 g m -2 y -1 ; Larsen et al., 2014a). For the FR sites, as indicated by the peaks of most P fractions and P total in the topsoil, deposition of dust may be more significant even in eroding (and thus rejuvenating) soils. Aciego et al. (2017) extrapolated a threemonth dust trapping record from the driest months in the Sierra Nevada to annual deposition rates of 3 to 36 g m -2 . Hence, although this extrapolation might be an overestimation due to limiting the measurements to the dry season, we acknowledge the likely accretion of Pbearing dust in the FR study area. However, dust deposition has little effect on our interpretation. With increasing soil depth, closer to the parent material source, and decreasing potential impact of atmospheric deposition, the P apatite remains low and secondary P forms remain clearly dominant regardless of soil residence time (Fig. 4). Alternative
Soil Residence Time or Age (Years) 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 Papatite/Ptotal (%) 0 20 40 60 80 100 WSA (This Study) Feather River (This Study)
Frantz Joseph (Stevens, 1968) Hawaii (Crews et al., 1995) Western Australia (Turner and Laliberté 2015) New Mexico (Lajtha and Schlesinger 1988) Haast River, NZ (Eger and others 2011) Mendocino California (Izquierdo et al., 2013) Cooloola, Australia (Chen and others 2015) Arizona Desert (Selmants and Hart 2010) Puerto Rico (McClintock et al., 2015) Fig. 5. Relative contributions of P apatite to P total across soil residence time and soil age gradients (filled solid points denote studies on hillslope soils, whereas hollow points are chronosequence studies).
Fig. 6. Climate data from the sites discussed in this study (Figs. 8 and 9). References as follows: 1, 2 Franz Josef, NZ (Richardson et al., 2004) 3 Haast River, NZ (Eger et al., 2011); 4 Mendocino, California (Izquierdo et al., 2013); 5 Arizona desert (Selmants and Hart, 2010); 6 Hawaii (Crews et al., 1995); 7 New Mexico (Lajtha and Schlesinger, 1988); 8, 9 Western Australia (Turner and Laliberté, 2015); 10 Northern Brazil (Agbenin and Tiessen, 1994); 11 Hawaii (Vitousek et al., 2003);12 Puerto Rico (McClintock et al., 2015); 13 Puerto Rico (Mage and Porder, 2013); 14 Cooloola, AUS (Chen et al., 2015); 15, 16 Northern Brazil (Araújo et al., 2004); 17 Ghana (Abekoe and Tiessen, 1998); 18 Sierra Nevada, California (Homyak et al., 2014).
or complementary explanations for the surface peak in P concentrations include P uplift by plants (Jobbágy and Jackson, 2004) or bioturbation within the soil (e.g., frequently observed tree throw in the study area).
All of the soils at our field sites are either Entisols or Inceptisols. Nevertheless, they are highly depleted in P apatite . Our data are largely consistent with other soil P studies conducted for eroding hillslopes (Fig. 8a). These observations clearly deviate from the general relationship of P fractions and soil orders postulated first by Smeck (1985), which is a pedological extension of the Walker and Syers (1976) paradigm, such that the progressive change of soil orders is aligned with the predictable changes in soil P. Soils are assumed to develop in a sequence from Entisols to Inceptisols to Alfisols to Ultisols (or Spodosols) to Oxisols (Smeck, 1985). The concept of correlation between soil P fractions and soil orders was later confirmed through global data compilations (Cross and Schlesinger, 1995;Lajtha and Schlesinger, 1988;Yang and Post, 2011). However, most of the data sets used to build these relationships between P fractions and soil order are from geomorphically stable landforms. In contrast, data from eroding hillslopes, regardless of soil order, show low P apatite contributions to P total (Fig. 8a). Consequently, neither soil residence time nor soil order is able to predict the low P apatite contributions on eroding, soil-mantled hillslopes.
4.5. Why do soil order and soil residence time fail to explain the contribution of P apatite ?
Whereas soil residence time on eroding hillslopes explains the dominance of soil orders typical of young geomorphic surfaces, it fails to account for the low contribution of P apatite to P total .
Soil order is determined largely by field observations of soil morphology including soil colour coatings, texture, structure, and horizons. The vertical depth distribution of these morphological properties is particularly diagnostic for several soil orders. For example, vertical distribution of soil texture and B horizon development are critical for determining a series of soils from Inceptisols to Ultisols (Soil Survey Staff, 2014). Continual mixing and/or consequent rejuvenation of a soil by erosion and soil production, for instance, would physically prevent the development of such vertical properties, similar to the effects of bioturbation (Johnson and Watson-Stegner, 1987). Because of this, soil orders characteristic for young geomorphic surfaces can develop from strongly weathered parent material as long as soil residence time is short and thus prevent significant vertical horizonation within the soils. It is also notable that the only soil order that is associated with mature soil development in Fig. 8a is the Oxisol. This is because in soil taxonomy, Oxisols, unlike Ultisols, do not require strong vertical stratification in clay contents and the classification is largely dependent on heavily-weathered soil minerals. The observation that soil order can be decoupled from weathering state of the parent material is not limited to eroding hillslopes. At the Cooloola sand dune soil chronosequence in Australia (Chen et al., 2015) (Fig. 8b), unlike at most other soil chronosequences, the young Entisols exhibits low levels of apatite P, simply because highly weathered sand deposits constitute the site's soil parent material.
The insensitivity of soil order to pre-weathering in parent material is consistent with our sites where Entisols and Inceptisols have formed from already chemically weathered saprolite. This is evident from the FR data (Fig. 4), where the soils have not formed from fresh bedrock but from saprolite overlying unweathered granodiorite. The saprolite weathering is evident by the dominance of secondary P forms (Fig. 4) and enrichment in biogeochemically conservative elements such as Zr (Yoo et al., 2011). We did not reach the depth to fresh bedrock despite hand augering to depths of 2 to 9 m below the soil-saprolite boundary.
Therefore at least for eroding, soil mantled hillslopes, available data suggest that soil P dynamics neither proceed in tandem with the general developmental sequence of soil orders as proposed by Smeck (1985) nor make soil residence time a good predictor of soil P dynamics.
We believe the reason for the discrepancy between chronosequences and hillslopes is derived from a fundamental difference between 'erosional' soils and 'depositional' soils. In contrast to hillslope soils, most chronosequences are originally developed in relatively unweathered parent material of water-, and glacier-transported origin. These transport mechanisms usually comprise comminution and particle size sorting. Deposition of lighter and more weathered mineral particles (clays, oxides) in the lowlands becomes less likely since these particles offer less resistance to physical transport (Dellinger et al., 2014;Kautz and Martin, 2007). Instead, the less weathered particles of mostly larger size fractions (silt, sand, > 2 mm) preferentially accumulate and ultimately form the parent material of lowland chronosequences (e.g., see parent material of chronosequence soils from NZ and California: Eger et al., 2011;Harden, 1987;Ross et al., 1977;Stevens, 1968;Wells and Goff, 2007). Fig. 7. Contribution of P apatite to P total versus mean annual precipitation for eroding hillslope soils from this study and from the literature. In plotting published data, we did not attempt to average the reported values or combine results from different depths for calculating soil profile integrated values. Instead all of the reported values are included in this figure. For the study conducted in the high Sierra Nevada (Homyak et al., 2014), we note that their reported P values are averaged over several soil profiles that include soils on hillslopes and adjacent depositional settings.
Chronosequences that are formed from volcanic rocks, like in Hawaii (Crews et al., 1995;Vitousek, 2004), behave similar to chronosequences developed from sedimentary lithologies: lava flows in Hawaii create new, minimally eroding geomorphic surfaces from initially unweathered, P apatite -rich parent material, conceptually similar to chronosequences on sedimentary deposits that involve particle sizedifferentiating transport (Fig. 5). To our knowledge, the Cooloola coastal dune sequence is the only published soil chronosequence with P fractionation data derived from a pre-weathered allochthonous parent material. Not unlike our residence time gradients, it shows low and invariant P apatite values across the entire sequence (Chen et al., 2015). Thus, both concepts, soil residence time and soil age in their narrow definitions do not consider any pre-weathering of parent material. However, soil age is often able to structure the evolution of soils on chronosequences because the parent material at the start of soil formation is usually minimally weathered.
It has been proposed (Porder et al., 2007a;Uhlig and von Blanckenburg, 2016) that a 'vertically oriented' version of the Walker and Syers (1976) model of P evolution applies to the changes of P fractions with depth (Fig. 9A). Such model recognises the inverse relationship between soil depth and mineral age following the incorporation of minerals into the active weathering zone of saprolite and soil.
However, the FR data and the review of existing studies allow this model to be modified at multiple fronts. In contrast to the expectation from the vertically oriented Walker and Syers model, total P does not gradually increase with increasing soil depth (Fig. 9B). Our data from FR (Fig. 4) indicate that total P decreases as bedrock chemically weathers to saprolite but that total P is greater in soils than in saprolite, albeit soil P dominated by secondary P forms. We attribute higher soil P concentrations to two processes: 1) dust deposition, and 2) biological nutrient redistribution (nutrient uplift) (Jobbágy and Jackson, 2004), whereby roots propagate into the saprolite and take up bio-available P from the saprolite zone. Plant-bound P is then returned to the soil as organic P, and partly transformed into other secondary (inorganic) P forms. Enrichment of P in surface soils associated with biological nutrient uptake and/or atmospheric deposition has been commonly observed (Chadwick and Asner, 2016;Merritts et al., 1992;Yoo et al., 2015) and does not seem to be necessarily limited to a particular P fraction or geomorphic setting (slopes vs. chronosequences) (Agbenin and Tiessen, 1994;Homyak et al., 2014;Lajtha and Schlesinger, 1988;Fig. 8. Soil P apatite /P total plotted against soil orders for a) eroding hillslopes, and b) soil chronosequences. Mage and Porder, 2013;Stevens, 1968;Turner and Laliberté, 2015). Additionally, organic matter is concentrated in the top of the weathering profile which together with Fe/Al oxides and secondary silicate clays partly protects P from leaching through the formation of P org , P occ and P Fe/Al . With increasing depth, the P-depleted saprolite zone beneath the enriched soil will eventually transition into more unweathered parent material with higher P concentrations and an increase of P apatite .
Therefore, it is not erosion directly that rejuvenates P. It is instead plant uptake of P at depth and dust deposition that rejuvenates soil P. Indirectly, erosion is required to maintain an ongoing supply of P to the root exploration zone of plants.
We characterised P fractionation in soils from eroding hillslopes across two soil residence time gradients and compared these new results against published soil P data from hillslopes and soil chronosequences on non-eroding landforms. We tested the Walker and Syers paradigm of soil P development as derived from soil chronosequences against hillslope soils through the conceptual link between soil residence times and soil ages. A naive application of this P model to eroding hillslopes predicts dominance of P apatite over secondary P in soils with very short residence times. However, we find the majority of soil P we and others have measured exists in the form of secondary P (83-97% in our data) regardless of soil residence time. Furthermore, soil residence time also does not explain the distribution of the secondary P forms. We conclude that the fundamental difference between chronosequence and hillslope soil derives from the weathering occurring in the bedrock (formation of saprolite) before it becomes part of the mobile soil. During initial stages of chronosequence development P apatite almost always dominates and P-depleted saprolite is normally not present. In contrast, on hillslopes weathered bedrock or saprolite appears to be common, combined with soils of short residence times and immature soil development. The legacy of pre-soil weathering of the underlying saprolite effectively counteracts the fertilising potential of the tectonic upliftsoil erosionsoil production feedback. Our data also indicate that plants may play an important role in redistributing P by uplift from the saprolite zone to the soil. Together with external dust, this redistribution increases soil P concentrations relative to the saprolite. Our work suggests that limits on ecosystem development through a decline of bio-available soil P forms may be more relevant to eroding hillslope soils than previously thought. Fig. 9. A) Vertically oriented Walker and Syers model and B) a new P dynamics model for weathering profiles based on the results from the Feather River sites. The WSA data show a similar pattern but lack the same detailed depth resolution due to the different sampling protocols. Although, we did not measure P total for fresh bedrock at the Feather River, globally compiled P contents of granodiorite bedrock types show mean value of 810 mg kg -1 with 25% value of 480 and 75% value of 1004 mg kg -1 (Porder and Ramachandran, 2013).